]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove lexicographic_to_hierarchic_numbering. Implement hierarchic_to_fe_q_hierarchic...
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Fri, 5 Aug 2005 08:47:47 +0000 (08:47 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Fri, 5 Aug 2005 08:47:47 +0000 (08:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@11227 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 040bb30fe5296f03843c6a34896541b25838e2fe..62dea1042d8b5f19b8fedf765c3a8f900ed93089 100644 (file)
@@ -228,7 +228,7 @@ template <int dim> 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 <int dim>
 class FE_Q_Hierarchical : public FE_Poly<TensorProductPolynomials<dim>,dim>
@@ -315,75 +315,63 @@ class FE_Q_Hierarchical : public FE_Poly<TensorProductPolynomials<dim>,dim>
     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:
+                                     * 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<unsigned int>
-    lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
-                                          const unsigned int            degree);
+    std::vector<unsigned int> hierarchic_to_fe_q_hierarchical_numbering (
+      const FiniteElementData<dim> &fe);
 
                                     /**
                                      * This is an analogon to the
@@ -392,7 +380,7 @@ class FE_Q_Hierarchical : public FE_Poly<TensorProductPolynomials<dim>,dim>
                                      */
     static
     std::vector<unsigned int>
-    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<unsigned int>
-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
index bb6c9145e6d18ab604132004b19f39f0202deea9..6bc2343da6f52bc5a90f48db1cff70466d240a69 100644 (file)
@@ -46,10 +46,10 @@ FE_Q_Hierarchical<dim>::FE_Q_Hierarchical (const unsigned int degree)
                    get_dpo_vector(degree),1, degree).dofs_per_cell, false),
                  std::vector<std::vector<bool> >(FiniteElementData<dim>(
                    get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector<bool>(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<FullMatrix<double> > &do
   const std::vector<unsigned int> &renumber=
     this->poly_space.get_numbering();
   
-
   for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
     {
       this->prolongation[c].reinit (this->dofs_per_cell, this->dofs_per_cell);
@@ -575,251 +574,186 @@ FE_Q_Hierarchical<dim>::get_dpo_vector(const unsigned int deg)
 template <int dim>
 std::vector<unsigned int>
 FE_Q_Hierarchical<dim>::
-lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
-                                      const unsigned int            degree)
+hierarchic_to_fe_q_hierarchical_numbering (const FiniteElementData<dim> &fe)
 {
-  std::vector<unsigned int> renumber (fe_data.dofs_per_cell, 0);
-  
-  const unsigned int n = degree+1;
+  Assert (fe.n_components() == 1, ExcInternalError());
+  std::vector<unsigned int> 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<fe.dofs_per_cell; ++i)
+         h2l[i] = i;
 
-  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, 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<fe.dofs_per_line; ++i)
+         h2l[next_index++] = 2+i;
+                                        // right line
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = (2+i)*n+1;
+                                        // top line
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = n+2+i;
+                                        // left line
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = (2+i)*n;
+                                        // inside quad
+       Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
+               ExcInternalError());
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+           h2l[next_index++] = (2+i)*n+2+j;
+
+       Assert (next_index == fe.dofs_per_cell, ExcInternalError());
+
+       break;
+      }
+
+      case 3:
+      {
+       unsigned int next_index = 0;
+       const unsigned int n2=n*n;
+                                        // first the eight vertices
+                                        // front face, counterclock wise
+       h2l[next_index++] = 0;
+       h2l[next_index++] = 1;
+       h2l[next_index++] = n2+1;
+       h2l[next_index++] = n2;
+                                        // back face, counterclock wise
+       h2l[next_index++] = n;
+       h2l[next_index++] = n+1;
+       h2l[next_index++] = n2+n+1;
+       h2l[next_index++] = n2+n;
+
+                                        // now the lines
+                                        // front face, counterclock wise
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = 2+i;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = (2+i)*n2+1;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = n2+2+i;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = (2+i)*n2;
+                                        // back face, counterclock wise
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = n+2+i;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = (2+i)*n2+n+1;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = n2+n+2+i;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = (2+i)*n2+n;
+                                        // lines in y-direction,
+                                        // counterclock wise
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = (2+i)*n;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = (2+i)*n+1;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = n2+(2+i)*n+1;
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         h2l[next_index++] = n2+(2+i)*n;
+
+                                        // inside quads
+       Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
+               ExcInternalError());
+                                        // front face
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+           h2l[next_index++] = (2+i)*n2+2+j;
+                                        // back face
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+           h2l[next_index++] = (2+i)*n2+n+2+j;
+                                        // bottom face
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+           h2l[next_index++] = (2+i)*n+2+j;
+                                        // right face
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+           h2l[next_index++] = (2+i)*n2+(2+j)*n+1;
+                                        // top face
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+           h2l[next_index++] = n2+(2+i)*n+2+j;
+                                        // left face
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+           h2l[next_index++] = (2+i)*n2+(2+j)*n;
+
+                                        // inside hex
+       Assert (fe.dofs_per_hex == fe.dofs_per_quad*fe.dofs_per_line,
+               ExcInternalError());
+       for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+         for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+           for (unsigned int k=0; k<fe.dofs_per_line; ++k)
+             h2l[next_index++] = (2+i)*n2+(2+j)*n+2+k;
+
+       Assert (next_index == fe.dofs_per_cell, ExcInternalError());
        
-             default:
-                   Assert(false, ExcNotImplemented());
-           }
+       break;
+      }       
 
-         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 (unsigned int i=0; i<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 = 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<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)
-           {
-                                              // 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<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 = 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 <int dim>
 std::vector<unsigned int>
 FE_Q_Hierarchical<dim>::
-face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
+face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int degree)
 {
   FiniteElementData<dim-1> fe_data(FE_Q_Hierarchical<dim-1>::get_dpo_vector(degree),1);
-  return FE_Q_Hierarchical<dim-1>::lexicographic_to_hierarchic_numbering (fe_data,
-                                                                         degree); 
+  return invert_numbering(FE_Q_Hierarchical<dim-1>::
+                         hierarchic_to_fe_q_hierarchical_numbering (fe_data));
 }
 
 
@@ -827,7 +761,7 @@ face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
 
 template <>
 std::vector<unsigned int>
-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<unsigned int> ();
 }

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.