]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
New functions output_indices and set_polynomial_ordering.
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 9 Mar 2004 14:39:31 +0000 (14:39 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 9 Mar 2004 14:39:31 +0000 (14:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@8698 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial_space.h
deal.II/base/source/polynomial_space.cc

index b837533f08dc4bbcc76bc112f6505e0191d09fdb..c3b375d5fcaba6313800558d56cb867b19cae89c 100644 (file)
  * Given a vector of <i>n</i> one-dimensional polynomials
  * <i>P<sub>0</sub></i> to <i>P<sub>n</sub></i>, where
  * <i>P<sub>i</sub></i> has degree <i>i</i>, this class generates all
- * polynomials of the form <i> P<sub>ijk</sub>(x,y,z) =
+ * multi-dimensional polynomials of the form <i>
+ * P<sub>ijk</sub>(x,y,z) =
  * P<sub>i</sub>(x)P<sub>j</sub>(y)P<sub>k</sub>(z)</i>, where the sum
  * of <i>i</i>, <i>j</i> and <i>k</i> is less than or equal <i>n</i>.
  *
- * @author Guido Kanschat, 2002, Wolfgang Bangerth, 2003
+ * The @ref{output_indices} function prints the ordering of the
+ * polynomials, i.e. for each multi-dimensional polynomial in the
+ * polynomial space it gives the indices i,j,k of the one-dimensional
+ * polynomials in x,y and z direction. The ordering of the
+ * multi-dimensional polynomials can be changed by using the
+ * @p{set_polynomial_ordering} function.
+ *
+ * @author Guido Kanschat, 2002, Wolfgang Bangerth, 2003, Ralf Hartmann 2004
  */
 template <int dim>
 class PolynomialSpace
@@ -58,6 +66,20 @@ class PolynomialSpace
     template <class Pol>
     PolynomialSpace (const std::vector<Pol> &pols);
 
+                                    /**
+                                     * Prints the list of the indices
+                                     * to <tt>out</tt>.
+                                     */
+    void output_indices(std::ostream &out) const;
+
+                                    /**
+                                     * Sets the ordering of the
+                                     * polynomials. Requires
+                                     * <tt>index_map.size()==n()</tt>. Stores
+                                     * a copy of <tt>index_map</tt>.
+                                     */
+    void set_polynomial_ordering(const vector<unsigned int> &index_map);
+    
                                     /**
                                      * Computes the value and the
                                      * first and second derivatives
@@ -152,7 +174,20 @@ class PolynomialSpace
                    int, int, int,
                    << "Dimension " << arg1 << " not equal to " << arg2 << " nor to " << arg3);
 
-           
+  protected:
+    
+                                    /**
+                                     * Compute numbers in x, y and z
+                                     * direction. Given an index
+                                     * <tt>n</tt> in the d-dimensional
+                                     * polynomial space, compute the
+                                     * indices i,j,k such that
+                                     * <i>p<sub>n</sub>(x,y,z) =
+                                     * p<sub>i</sub>(x)p<sub>j</sub>(y)p<sub>k</sub>(z)</i>.
+                                     */
+    void compute_index (const unsigned int n,
+                        unsigned int      (&index)[dim]) const;
+
   private:
                                     /**
                                      * Copy of the vector <tt>pols</tt> of
@@ -169,16 +204,16 @@ class PolynomialSpace
     const unsigned int n_pols;
 
                                     /**
-                                     * Compute numbers in x, y and z
-                                     * direction. Given an index
-                                     * <tt>n</tt> in the d-dimensional
-                                     * polynomial space, compute the
-                                     * indices i,j,k such that
-                                     * <i>p<sub>n</sub>(x,y,z) =
-                                     * p<sub>i</sub>(x)p<sub>j</sub>(y)p<sub>k</sub>(z)</i>.
+                                     * Index map for reordering the
+                                     * polynomials.
                                      */
-    void compute_index (const unsigned int n,
-                        unsigned int      (&index)[dim]) const;
+    std::vector<unsigned int> index_map;
+
+                                    /**
+                                     * Index map for reordering the
+                                     * polynomials.
+                                     */
+    std::vector<unsigned int> reverse_index_map;
     
                                     /**
                                      * Static function used in the
@@ -208,12 +243,23 @@ void PolynomialSpace<3>::compute_index(const unsigned int n,
 
 template <int dim>
 template <class Pol>
-PolynomialSpace<dim>::
-PolynomialSpace (const std::vector<Pol> &pols)
+PolynomialSpace<dim>::PolynomialSpace (const std::vector<Pol> &pols)
                :
                polynomials (pols.begin(), pols.end()),
-               n_pols (compute_n_pols(polynomials.size()))
-{}
+               n_pols (compute_n_pols(polynomials.size())),
+               index_map(n_pols),
+               reverse_index_map(n_pols)
+{
+                                  // per default set this index map
+                                  // to identity. This map can be
+                                  // changed by the user through the
+                                  // set_polynomial_ordering function
+  for (unsigned int i=0; i<n_pols; ++i)
+    {
+      index_map[i]=i;
+      reverse_index_map[i]=i;
+    }
+}
 
 
 template<int dim>
index aefd6e5ff34b83e39f4e852ea0b15464b2119cb8..eaa6ef02ffdfea8f5c97b58ba8988441ab153610 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2000, 2001, 2002, 2003 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -17,6 +17,7 @@
 #include <base/table.h>
 
 
+
 template <int dim>
 unsigned int
 PolynomialSpace<dim>::compute_n_pols (const unsigned int n)
@@ -34,9 +35,12 @@ PolynomialSpace<dim>::compute_n_pols (const unsigned int n)
 template <>
 void
 PolynomialSpace<1>::
-compute_index (const unsigned int n,
+compute_index (const unsigned int i,
                unsigned int      (&index)[1]) const
 {
+  Assert(i<index_map.size(),
+        ExcIndexRange(i,0,index_map.size()));
+  const unsigned int n=index_map[i];
   index[0] = n;
 }
 
@@ -45,9 +49,12 @@ compute_index (const unsigned int n,
 template <>
 void
 PolynomialSpace<2>::
-compute_index (const unsigned int n,
+compute_index (const unsigned int i,
                unsigned int      (&index)[2]) const
 {
+  Assert(i<index_map.size(),
+        ExcIndexRange(i,0,index_map.size()));
+  const unsigned int n=index_map[i];
                                    // there should be a better way to
                                    // write this function (not
                                    // linear in n_1d), someone
@@ -70,9 +77,12 @@ compute_index (const unsigned int n,
 template <>
 void
 PolynomialSpace<3>::
-compute_index (const unsigned int n,
+compute_index (const unsigned int i,
                unsigned int      (&index)[3]) const
 {
+  Assert(i<index_map.size(),
+        ExcIndexRange(i,0,index_map.size()));
+  const unsigned int n=index_map[i];
                                    // there should be a better way to
                                    // write this function (not
                                    // quadratic in n_1d), someone
@@ -96,6 +106,37 @@ compute_index (const unsigned int n,
 }
 
 
+template <int dim>
+void
+PolynomialSpace<dim>::output_indices(std::ostream &out) const
+{
+  unsigned int ix[dim];
+  for (unsigned int i=0; i<n_pols; ++i)
+    {
+      compute_index(i,ix);
+      out << i << "\t";
+      for (unsigned int d=0; d<dim; ++d)
+       out << ix[d] << " ";
+      out << endl;
+    }
+}
+
+
+
+template <int dim>
+void
+PolynomialSpace<dim>::set_polynomial_ordering(
+  const vector<unsigned int> &imap)
+{
+  Assert(imap.size()==index_map.size(),
+        ExcDimensionMismatch(imap.size(), index_map.size()));
+
+  index_map=imap;
+  for (unsigned int i=0; i<index_map.size(); ++i)
+    reverse_index_map[index_map[i]]=i;
+}
+
+
 
 template <int dim>
 double
@@ -104,7 +145,6 @@ PolynomialSpace<dim>::compute_value (const unsigned int i,
 {
   unsigned int ix[dim];
   compute_index(i,ix);
-
                                    // take the product of the
                                    // polynomials in the various space
                                    // directions
@@ -235,9 +275,10 @@ PolynomialSpace<dim>::compute (const Point<dim>            &p,
       for (unsigned int iz=0;iz<((dim>2) ? n_1d : 1);++iz)
        for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy)
          for (unsigned int ix=0; ix<n_1d-iy-iz; ++ix)
-           values[k++] = v[0][ix][0]
-                         * ((dim>1) ? v[1][iy][0] : 1.)
-                         * ((dim>2) ? v[2][iz][0] : 1.);
+           values[reverse_index_map[k++]] =
+             v[0][ix][0]
+             * ((dim>1) ? v[1][iy][0] : 1.)
+             * ((dim>2) ? v[2][iz][0] : 1.);
     }
   
   if (update_grads)
@@ -248,11 +289,11 @@ PolynomialSpace<dim>::compute (const Point<dim>            &p,
        for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy)
          for (unsigned int ix=0; ix<n_1d-iy-iz; ++ix)
            {
+             const unsigned int k2=reverse_index_map[k++];
              for (unsigned int d=0;d<dim;++d)
-               grads[k][d] = v[0][ix][(d==0) ? 1 : 0]
-                              * ((dim>1) ? v[1][iy][(d==1) ? 1 : 0] : 1.)
-                              * ((dim>2) ? v[2][iz][(d==2) ? 1 : 0] : 1.);
-             ++k;
+               grads[k2][d] = v[0][ix][(d==0) ? 1 : 0]
+                  * ((dim>1) ? v[1][iy][(d==1) ? 1 : 0] : 1.)
+                  * ((dim>2) ? v[2][iz][(d==2) ? 1 : 0] : 1.);
            }
     }
 
@@ -264,6 +305,7 @@ PolynomialSpace<dim>::compute (const Point<dim>            &p,
        for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy)
          for (unsigned int ix=0; ix<n_1d-iy-iz; ++ix)
            {
+             const unsigned int k2=reverse_index_map[k++];
              for (unsigned int d1=0; d1<dim; ++d1)
                for (unsigned int d2=0; d2<dim; ++d2)
                  {
@@ -277,11 +319,11 @@ PolynomialSpace<dim>::compute (const Point<dim>            &p,
                    const unsigned int
                      j2 = ((d1==2) ? 1 : 0) + ((d2==2) ? 1 : 0);
                    
-                   grad_grads[k][d1][d2] = v[0][ix][j0]
-                                           * ((dim>1) ? v[1][iy][j1] : 1.)
-                                           * ((dim>2) ? v[2][iz][j2] : 1.);
+                   grad_grads[k2][d1][d2] =
+                     v[0][ix][j0]
+                     * ((dim>1) ? v[1][iy][j1] : 1.)
+                     * ((dim>2) ? v[2][iz][j2] : 1.);
                  }
-             ++k;
            }
     }
 }

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.