]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Clean up some real old and unintelligible cruft.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 Apr 2003 16:09:21 +0000 (16:09 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 Apr 2003 16:09:21 +0000 (16:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@7415 0785d39b-7218-0410-832d-ea1e28bc413d

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

index ad06c30f7c73e733252f3d5706356501bddb6a81..0151d48daef9d3802633d5458e533980b907fa4a 100644 (file)
@@ -33,7 +33,7 @@
  * of the form @p{ Pijk(x,y,z) = Pi(x)Pj(y)Pk(z)}, where the sum of
  * @p{i}, @p{j} and @p{k} is less than or equal @p{n}.
  *
- * @author Guido Kanschat, 2002
+ * @author Guido Kanschat, 2002, Wolfgang Bangerth, 2003
  */
 template <int dim>
 class PolynomialSpace
@@ -55,7 +55,7 @@ class PolynomialSpace
                                      * @p{Polynomial<double>}.
                                      */
     template <class Pol>
-    PolynomialSpace(const std::vector<Pol> &pols);
+    PolynomialSpace (const std::vector<Pol> &pols);
 
                                     /**
                                      * Computes the value and the
@@ -65,9 +65,12 @@ class PolynomialSpace
                                      *
                                      * The size of the vectors must
                                      * either be equal @p{0} or equal
-                                     * @p{n()}.  In the
-                                     * first case, the function will
-                                     * not compute these values.
+                                     * @p{n()}.  In the first case,
+                                     * the function will not compute
+                                     * these values, i.e. you
+                                     * indicate what you want to have
+                                     * computed by resizing those
+                                     * vectors which you want filled.
                                      *
                                      * If you need values or
                                      * derivatives of all polynomials
@@ -79,8 +82,8 @@ class PolynomialSpace
                                      * functions, see below, in a
                                      * loop over all polynomials.
                                      */
-    void compute (const Point<dim>                     &unit_point,
-                 std::vector<double>                  &values,
+    void compute (const Point<dim>            &unit_point,
+                 std::vector<double>         &values,
                  std::vector<Tensor<1,dim> > &grads,
                  std::vector<Tensor<2,dim> > &grad_grads) const;
     
@@ -112,8 +115,8 @@ class PolynomialSpace
                                      *
                                      * Consider using @p{compute} instead.
                                      */
-    Tensor<2,dim> compute_grad_grad(const unsigned int i,
-                                   const Point<dim> &p) const;
+    Tensor<2,dim> compute_grad_grad (const unsigned int i,
+                                     const Point<dim> &p) const;
 
                                     /**
                                      * Return the number of
@@ -127,7 +130,7 @@ class PolynomialSpace
                                      * 2d, and @p{N(N+1)(N+2)/6 in
                                      * 3d.
                                      */
-    unsigned int n() const;
+    unsigned int n () const;
 
                                     /**
                                      * Exception.
@@ -161,10 +164,8 @@ class PolynomialSpace
                                      * @p{p_n(x,y,z) =
                                      * p_i(x)p_j(y)p_k(z)}.
                                      */
-    void compute_index(unsigned int n,
-                      unsigned int& nx,
-                      unsigned int& ny,
-                      unsigned int& nz) const;
+    void compute_index (const unsigned int n,
+                        unsigned int      (&index)[dim]) const;
     
                                     /**
                                      * Static function used in the
@@ -175,6 +176,21 @@ class PolynomialSpace
 };
 
 
+/* -------------- declaration of explicit specializations --- */
+
+template <>
+void PolynomialSpace<1>::compute_index(const unsigned int n,
+                                       unsigned int      (&index)[1]) const;
+template <>
+void PolynomialSpace<2>::compute_index(const unsigned int n,
+                                       unsigned int      (&index)[2]) const;
+template <>
+void PolynomialSpace<3>::compute_index(const unsigned int n,
+                                       unsigned int      (&index)[3]) const;
+
+
+
+/* -------------- inline and template functions ------------- */
 
 template <int dim>
 template <class Pol>
index 92d05e3caaf152f45b6099fa15a405530d9df22c..7dbcd9d70ab05d9189ce75d36bf313cab5f64069 100644 (file)
@@ -19,7 +19,6 @@
 #include <base/tensor.h>
 #include <base/point.h>
 #include <base/polynomial.h>
-#include <base/smartpointer.h>
 
 #include <vector>
 
  * $[0,d], then the tensor product polynomials are orthogonal on
  * $[-1,1]^d$ or $[0,1]^d$, respectively.
  *
- * @author Ralf Hartmann, 2000, documentation Guido Kanschat
+ * Indexing is as following: the order of dim-dimensional polynomials
+ * is x-coordinates running fastest, then y-coordinate, etc. The first
+ * few polynomials are thus @p{P1(x)P1(y)}, @p{P2(x)P1(y)},
+ * @p{P3(x)P1(y)}, ..., @p{P1(x)P2(y)}, @p{P2(x)P2(y)},
+ * @p{P3(x)P2(y)}, ..., and likewise in 3d.
+ * 
+ * @author Ralf Hartmann, Guido Kanschat, 2000, Wolfgang Bangerth 2003
  */
 template <int dim>
 class TensorProductPolynomials
@@ -42,13 +47,15 @@ class TensorProductPolynomials
   public:
                                     /**
                                      * Constructor. @p{pols} is a
-                                     * vector of pointers to
-                                     * one-dimensional polynomials
-                                     * and will be copied into the
-                                     * member variable @p{polynomials}.
+                                     * vector of objects that should
+                                     * be derived or otherwise
+                                     * convertible to one-dimensional
+                                     * polynomial objects and will be
+                                     * copied into the member
+                                     * variable @p{polynomials}.
                                      */
     template <class Pol>
-    TensorProductPolynomials(const std::vector<Pol> &pols);
+    TensorProductPolynomials (const std::vector<Pol> &pols);
 
                                     /**
                                      * Computes the value and the
@@ -74,10 +81,10 @@ class TensorProductPolynomials
                                      * loop over all tensor product
                                      * polynomials.
                                      */
-    void compute(const Point<dim>                     &unit_point,
-                std::vector<double>                  &values,
-                std::vector<Tensor<1,dim> > &grads,
-                std::vector<Tensor<2,dim> > &grad_grads) const;
+    void compute (const Point<dim>            &unit_point,
+                  std::vector<double>         &values,
+                  std::vector<Tensor<1,dim> > &grads,
+                  std::vector<Tensor<2,dim> > &grad_grads) const;
     
                                     /**
                                      * Computes the value of the
@@ -158,8 +165,8 @@ class TensorProductPolynomials
                                      * and in a much more efficient
                                      * way.
                                      */
-    Tensor<2,dim> compute_grad_grad(const unsigned int i,
-                                   const Point<dim> &p) const;
+    Tensor<2,dim> compute_grad_grad (const unsigned int i,
+                                     const Point<dim> &p) const;
 
                                     /**
                                      * Returns the number of tensor
@@ -191,42 +198,70 @@ class TensorProductPolynomials
                                      */
     const unsigned int n_tensor_pols;
 
-                                    /**
-                                     * @p{n_pols_to[n]=polynomials.size()^n}
-                                     * Filled by the constructor.
-                                     *
-                                     * For internal use only. 
-                                     */
-    std::vector<unsigned int> n_pols_to;
+                                     /**
+                                      * Each tensor product polynomial
+                                      * @รพ{i} is a product of
+                                      * one-dimensional polynomials in
+                                      * each space direction. Compute
+                                      * the indices of these
+                                      * one-dimensional polynomials
+                                      * for each space direction,
+                                      * given the index @p{i}.
+                                      */
+    void compute_index (const unsigned int i,
+                        unsigned int       (&indices)[dim]) const;
     
                                     /**
                                      * Computes @p{x} to the power of
-                                     * @p{y} for unsigned int @p{x}
-                                     * and @p{y}. It is a private
-                                     * function as it is only used in
-                                     * this class.
+                                     * @p{dim} for unsigned int @p{x}.
+                                     * Used in the constructor.
                                      */
-    static unsigned int power(const unsigned int x, const unsigned int y);
+    static
+    unsigned int x_to_the_dim (const unsigned int x);
 };
 
 
 
+/* -------------- declaration of explicit specializations --- */
+
+template <>
+void
+TensorProductPolynomials<1>::compute_index(const unsigned int n,
+                                           unsigned int      (&index)[1]) const;
+template <>
+void
+TensorProductPolynomials<2>::compute_index(const unsigned int n,
+                                           unsigned int      (&index)[2]) const;
+template <>
+void
+TensorProductPolynomials<3>::compute_index(const unsigned int n,
+                                           unsigned int      (&index)[3]) const;
+
+
+/* ---------------- template and inline functions ---------- */
+
+template <int dim>
+inline
+unsigned int
+TensorProductPolynomials<dim>::
+x_to_the_dim (const unsigned int x)
+{
+  unsigned int y = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    y *= x;
+  return y;
+}
+
+
+
 template <int dim>
 template <class Pol>
 TensorProductPolynomials<dim>::
 TensorProductPolynomials(const std::vector<Pol> &pols)
                :
                polynomials (pols.begin(), pols.end()),
-               n_tensor_pols(power(pols.size(), dim)),
-               n_pols_to(dim+1)
-{
-  const unsigned int n_pols=polynomials.size();
-
-  n_pols_to[0]=1;
-  for (unsigned int i=0; i<dim; ++i)
-    n_pols_to[i+1]=n_pols_to[i]*n_pols;
-  Assert(n_pols_to[dim]==n_tensor_pols, ExcInternalError());
-}
+               n_tensor_pols(x_to_the_dim(pols.size()))
+{}
 
 
 
index 85b084e084b4adad1bb67298bd7b3337d14296d7..1debf26e89aafbceeb58242f3f28bdad35d3e52e 100644 (file)
@@ -22,7 +22,7 @@ unsigned int
 PolynomialSpace<dim>::compute_n_pols (const unsigned int n)
 {
   unsigned int n_pols = n;
-  for (unsigned int i=1;i<dim;++i)
+  for (unsigned int i=1; i<dim; ++i)
     {
       n_pols *= (n+i);
       n_pols /= (i+1);
@@ -31,65 +31,110 @@ PolynomialSpace<dim>::compute_n_pols (const unsigned int n)
 }
 
 
-template <int dim>
+template <>
+void
+PolynomialSpace<1>::
+compute_index(const unsigned int n,
+              unsigned int      (&index)[1]) const
+{
+  index[0] = n;
+}
+
+
+
+template <>
 void
-PolynomialSpace<dim>::compute_index(unsigned int n,
-                                   unsigned int& nx,
-                                   unsigned int& ny,
-                                   unsigned int& nz) const
+PolynomialSpace<2>::
+compute_index (const unsigned int n,
+               unsigned int      (&index)[2]) const
 {
+                                   // there should be a better way to
+                                   // write this function (not
+                                   // linear in n_1d), someone
+                                   // should think about this...
   const unsigned int n_1d=polynomials.size();
   unsigned int k=0;
-  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)
-       if (k++ == n)
-         {
-           nz = iz;
-           ny = iy;
-           nx = ix;
-           return;
-         }
+  for (unsigned int iy=0; iy<n_1d; ++iy)
+    if (n < k+n_1d-iy)
+      {
+        index[0] = n-k;
+        index[1] = iy;
+        return;
+      }  
+    else
+      k+=n_1d-iy;
 }
 
 
+
+template <>
+void
+PolynomialSpace<3>::
+compute_index (const unsigned int n,
+               unsigned int      (&index)[3]) const
+{
+                                   // there should be a better way to
+                                   // write this function (not
+                                   // quadratic in n_1d), someone
+                                   // should think about this...
+                                   //
+                                   // (ah, and yes: the original
+                                   // algorithm was even cubic!)
+  const unsigned int n_1d=polynomials.size();
+  unsigned int k=0;
+  for (unsigned int iz=0; iz<n_1d; ++iz)
+    for (unsigned int iy=0; iy<n_1d-iz; ++iy)
+      if (n < k+n_1d-iy-iz)
+        {
+          index[0] = n-k;
+          index[1] = iy;
+          index[2] = iz;
+          return;
+        }
+      else
+        k += n_1d-iy-iz;
+}
+
+
+
 template <int dim>
 double
-PolynomialSpace<dim>::compute_value(const unsigned int i,
-                                   const Point<dim> & p) const
+PolynomialSpace<dim>::compute_value (const unsigned int i,
+                                     const Point<dim>  &p) const
 {
-  unsigned int ix = 0;
-  unsigned int iy = 0;
-  unsigned int iz = 0;
-  compute_index(i,ix,iy,iz);
-  
-  double result = polynomials[ix].value(p(0));
-  if (dim>1)
-    result *= polynomials[iy].value(p(1));
-  if (dim>2)
-    result *= polynomials[iz].value(p(2));
+  unsigned int ix[dim];
+  compute_index(i,ix);
+
+                                   // take the product of the
+                                   // polynomials in the various space
+                                   // directions
+  double result = 1.;
+  for (unsigned int d=0; d<dim; ++d)
+    result *= polynomials[ix[d]].value(p(d));
   return result;
 }
 
-  
+
+
 template <int dim>
 Tensor<1,dim>
-PolynomialSpace<dim>::compute_grad(const unsigned int i,
-                                  const Point<dim> &p) const
+PolynomialSpace<dim>::compute_grad (const unsigned int i,
+                                    const Point<dim>  &p) const
 {
-  unsigned int ix[3];
-  compute_index(i,ix[0],ix[1],ix[2]);
+  unsigned int ix[dim];
+  compute_index(i,ix);
   
   Tensor<1,dim> result;
-  for (unsigned int d=0;d<dim;++d)
+  for (unsigned int d=0; d<dim; ++d)
     result[d] = 1.;
-  
+
+                                   // get value and first derivative
   std::vector<double> v(2);
-  for (unsigned int d=0;d<dim;++d)
+  for (unsigned int d=0; d<dim; ++d)
     {
       polynomials[ix[d]].value(p(d), v);
       result[d] *= v[1];
-      for (unsigned int d1=0;d1<dim;++d1)
+      for (unsigned int d1=0; d1<dim; ++d1)
        if (d1 != d)
          result[d1] *= v[0];
     }
@@ -99,29 +144,31 @@ PolynomialSpace<dim>::compute_grad(const unsigned int i,
 
 template <int dim>
 Tensor<2,dim>
-PolynomialSpace<dim>::compute_grad_grad(const unsigned int i,
-                                       const Point<dim> &p) const
+PolynomialSpace<dim>::compute_grad_grad (const unsigned int i,
+                                         const Point<dim>  &p) const
 {
-  unsigned int ix[3];
-  compute_index(i,ix[0],ix[1],ix[2]);
+  unsigned int ix[dim];
+  compute_index(i,ix);
   
   Tensor<2,dim> result;
-  for (unsigned int d=0;d<dim;++d)
-    for (unsigned int d1=0;d1<dim;++d1)
+  for (unsigned int d=0; d<dim; ++d)
+    for (unsigned int d1=0; d1<dim; ++d1)
       result[d][d1] = 1.;
   
+                                   // get value, first and second
+                                   // derivatives
   std::vector<double> v(3);
-  for (unsigned int d=0;d<dim;++d)
+  for (unsigned int d=0; d<dim; ++d)
     {
       polynomials[ix[d]].value(p(d), v);
       result[d][d] *= v[2];
-      for (unsigned int d1=0;d1<dim;++d1)
+      for (unsigned int d1=0; d1<dim; ++d1)
        {
          if (d1 != d)
            {
              result[d][d1] *= v[1];
              result[d1][d] *= v[1];
-             for (unsigned int d2=0;d2<dim;++d2)
+             for (unsigned int d2=0; d2<dim; ++d2)
                if (d2 != d)
                  result[d1][d2] *= v[0];
            }
@@ -134,11 +181,11 @@ PolynomialSpace<dim>::compute_grad_grad(const unsigned int i,
 
 
 template <int dim>
-void PolynomialSpace<dim>::compute(
-  const Point<dim>                     &p,
-  std::vector<double>                  &values,
-  std::vector<Tensor<1,dim> > &grads,
-  std::vector<Tensor<2,dim> > &grad_grads) const
+void
+PolynomialSpace<dim>::compute (const Point<dim>            &p,
+                               std::vector<double>         &values,
+                               std::vector<Tensor<1,dim> > &grads,
+                               std::vector<Tensor<2,dim> > &grad_grads) const
 {
   const unsigned int n_1d=polynomials.size();
   
index 6c9a62c1d3854a4eb47cbd3f0843d2af083ccb0e..d60ecdea401c0d58f5d39159c47ec8b1bc604221 100644 (file)
 
 
 
+template <>
+void
+TensorProductPolynomials<1>::
+compute_index (const unsigned int i,
+               unsigned int       (&indices)[1]) const
+{
+  Assert (i<polynomials.size(), ExcInternalError());
+  indices[0] = i;
+}
 
-template <int dim>
-unsigned int TensorProductPolynomials<dim>::power(const unsigned int x,
-                                                 const unsigned int y)
+
+
+template <>
+void
+TensorProductPolynomials<2>::
+compute_index (const unsigned int i,
+               unsigned int       (&indices)[2]) const
 {
-  unsigned int value=1;
-  for (unsigned int i=0; i<y; ++i)
-    value*=x;
-  return value;
+  const unsigned int n_pols = polynomials.size();
+  Assert (i<n_pols*n_pols, ExcInternalError());
+
+  indices[0] = i % n_pols;
+  indices[1] = i / n_pols;
+}
+
+
+
+template <>
+void
+TensorProductPolynomials<3>::
+compute_index (const unsigned int i,
+               unsigned int       (&indices)[3]) const
+{
+  const unsigned int n_pols = polynomials.size();
+  Assert (i<n_pols*n_pols*n_pols, ExcInternalError());
+
+  indices[0] = i % n_pols;
+  indices[1] = (i/n_pols) % n_pols;
+  indices[2] = i / (n_pols*n_pols);
 }
 
 
 
 template <int dim>
 double
-TensorProductPolynomials<dim>::compute_value(const unsigned int i,
-                                            const Point<dim> &p) const
+TensorProductPolynomials<dim>::compute_value (const unsigned int i,
+                                              const Point<dim> &p) const
 {
-  const unsigned int n_pols=polynomials.size();
+  unsigned int indices[dim];
+  compute_index (i, indices);
   
   double value=1.;
   for (unsigned int d=0; d<dim; ++d)
-    value *= polynomials[(i/n_pols_to[d])%n_pols].value(p(d));
+    value *= polynomials[indices[d]].value(p(d));
   
   return value;
 }
@@ -48,10 +79,11 @@ TensorProductPolynomials<dim>::compute_value(const unsigned int i,
   
 template <int dim>
 Tensor<1,dim>
-TensorProductPolynomials<dim>::compute_grad(const unsigned int i,
-                                           const Point<dim> &p) const
+TensorProductPolynomials<dim>::compute_grad (const unsigned int i,
+                                             const Point<dim> &p) const
 {
-  const unsigned int n_pols=polynomials.size();
+  unsigned int indices[dim];
+  compute_index (i, indices);
 
                                    // compute values and
                                    // uni-directional derivatives at
@@ -59,15 +91,15 @@ TensorProductPolynomials<dim>::compute_grad(const unsigned int i,
                                    // co-ordinate direction
   std::vector<std::vector<double> > v(dim, std::vector<double> (2));
   for (unsigned int d=0; d<dim; ++d)
-    polynomials[(i/n_pols_to[d])%n_pols].value(p(d), v[d]);
+    polynomials[indices[d]].value(p(d), v[d]);
   
   Tensor<1,dim> grad;
   for (unsigned int d=0; d<dim; ++d)
-    grad[d]=1.;
-  
-  for (unsigned int d=0; d<dim; ++d)
-    for (unsigned int x=0; x<dim; ++x)
-      grad[d]*=v[x][d==x];
+    {
+      grad[d] = 1.;
+      for (unsigned int x=0; x<dim; ++x)
+        grad[d] *= v[x][d==x];
+    }
   
   return grad;
 }
@@ -75,35 +107,34 @@ TensorProductPolynomials<dim>::compute_grad(const unsigned int i,
 
 template <int dim>
 Tensor<2,dim>
-TensorProductPolynomials<dim>::compute_grad_grad(const unsigned int i,
-                                                const Point<dim> &p) const
+TensorProductPolynomials<dim>::compute_grad_grad (const unsigned int i,
+                                                  const Point<dim> &p) const
 {
-  const unsigned int n_pols=polynomials.size();
-    
+  unsigned int indices[dim];
+  compute_index (i, indices);
+
   std::vector<std::vector<double> > v(dim, std::vector<double> (3));
   for (unsigned int d=0; d<dim; ++d)
-    polynomials[(i/n_pols_to[d])%n_pols].value(p(d), v[d]);
+    polynomials[indices[d]].value(p(d), v[d]);
   
   Tensor<2,dim> grad_grad;
-
   for (unsigned int d1=0; d1<dim; ++d1)
     for (unsigned int d2=0; d2<dim; ++d2)
-      grad_grad[d1][d2]=1.;
-  
-  for (unsigned int x=0; x<dim; ++x)
-    for (unsigned int d1=0; d1<dim; ++d1)
-      for (unsigned int d2=0; d2<dim; ++d2)
-       {
-         unsigned int derivative=0;
-         if (d1==x || d2==x)
-           {
-             if (d1==d2)
-               derivative=2;
-             else
-               derivative=1;
-           } 
-         grad_grad[d1][d2]*=v[x][derivative];
-       }
+      {
+        grad_grad[d1][d2] = 1.;
+        for (unsigned int x=0; x<dim; ++x)
+          {
+            unsigned int derivative=0;
+            if (d1==x || d2==x)
+              {
+                if (d1==d2)
+                  derivative=2;
+                else
+                  derivative=1;
+              } 
+            grad_grad[d1][d2] *= v[x][derivative];
+          }
+      }
 
   return grad_grad;
 }
@@ -112,96 +143,96 @@ TensorProductPolynomials<dim>::compute_grad_grad(const unsigned int i,
 
 
 template <int dim>
-void TensorProductPolynomials<dim>::compute(
-  const Point<dim>                     &p,
-  std::vector<double>                  &values,
-  std::vector<Tensor<1,dim> > &grads,
-  std::vector<Tensor<2,dim> > &grad_grads) const
+void
+TensorProductPolynomials<dim>::
+compute (const Point<dim>            &p,
+         std::vector<double>         &values,
+         std::vector<Tensor<1,dim> > &grads,
+         std::vector<Tensor<2,dim> > &grad_grads) const
 {
-  const unsigned int n_pols=polynomials.size();
-  
-  Assert(values.size()==n_tensor_pols || values.size()==0,
-        ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
-  Assert(grads.size()==n_tensor_pols|| grads.size()==0,
-        ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
-  Assert(grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
-        ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
-
-  unsigned int v_size=0;
-  bool update_values=false, update_grads=false, update_grad_grads=false;
-  if (values.size()==n_tensor_pols)
-    {
-      update_values=true;
-      v_size=1;
-    }
-  if (grads.size()==n_tensor_pols)
-    {
-      update_grads=true;
-      v_size=2;
-    }
-  if (grad_grads.size()==n_tensor_pols)
-    {
-      update_grad_grads=true;
-      v_size=3;
-    }
+  Assert (values.size()==n_tensor_pols || values.size()==0,
+          ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
+  Assert (grads.size()==n_tensor_pols|| grads.size()==0,
+          ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
+  Assert (grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
+          ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
+
+  const bool update_values     = (values.size() == n_tensor_pols),
+             update_grads      = (grads.size()==n_tensor_pols),
+             update_grad_grads = (grad_grads.size()==n_tensor_pols);
+
+                                   // check how many
+                                   // values/derivatives we have to
+                                   // compute
+  unsigned int n_values_and_derivatives = 0;
+  if (update_values)
+    n_values_and_derivatives = 1;
+  if (update_grads)
+    n_values_and_derivatives = 2;
+  if (update_grad_grads)
+    n_values_and_derivatives = 3;
+
 
-  Table<2,std::vector<double> > v(dim, n_pols);
-  for (unsigned int d=0; d<v.size()[0]; ++d)
-    for (unsigned int i=0; i<v.size()[1]; ++i)
+                                   // compute the values (and
+                                   // derivatives, if necessary) of
+                                   // all polynomials at this
+                                   // evaluation point
+  Table<2,std::vector<double> > v(dim, polynomials.size());
+  for (unsigned int d=0; d<dim; ++d)
+    for (unsigned int i=0; i<polynomials.size(); ++i)
       {
-        v(d,i).resize (v_size, 0.);
+        v(d,i).resize (n_values_and_derivatives, 0.);
         polynomials[i].value(p(d), v(d,i));
       };
   
-  if (update_values)
+  for (unsigned int i=0; i<n_tensor_pols; ++i)
     {
-      for (unsigned int i=0; i<n_tensor_pols; ++i)
-       values[i]=1;
+                                       // first get the
+                                       // one-dimensional indices of
+                                       // this particular tensor
+                                       // product polynomial
+      unsigned int indices[dim];
+      compute_index (i, indices);
       
-      for (unsigned int x=0; x<dim; ++x)
-       for (unsigned int i=0; i<n_tensor_pols; ++i)
-         values[i]*=v[x][(i/n_pols_to[x])%n_pols][0];
-    }
+      if (update_values)
+        {
+          values[i] = 1;
+          for (unsigned int x=0; x<dim; ++x)
+            values[i] *= v(x,indices[x])[0];
+        }
   
-  if (update_grads)
-    {
-      for (unsigned int i=0; i<n_tensor_pols; ++i)
-       for (unsigned int d=0; d<dim; ++d)
-         grads[i][d]=1.;
-
-      for (unsigned int x=0; x<dim; ++x)
-       for (unsigned int i=0; i<n_tensor_pols; ++i)
-         for (unsigned int d=0; d<dim; ++d)
-           grads[i][d]*=v[x][(i/n_pols_to[x])%n_pols][d==x];
-    }
-
-  if (update_grad_grads)
-    {
-      for (unsigned int i=0; i<n_tensor_pols; ++i)
-       for (unsigned int d1=0; d1<dim; ++d1)
-         for (unsigned int d2=0; d2<dim; ++d2)
-           grad_grads[i][d1][d2]=1.;
-
-      for (unsigned int x=0; x<dim; ++x)
-       for (unsigned int i=0; i<n_tensor_pols; ++i)
-         for (unsigned int d1=0; d1<dim; ++d1)
-           for (unsigned int d2=0; d2<dim; ++d2)
-             {
-               unsigned int derivative=0;
-               if (d1==x || d2==x)
-                 {
-                   if (d1==d2)
-                     derivative=2;
-                   else
-                     derivative=1;
-                 } 
-               grad_grads[i][d1][d2]*=
-                 v[x][(i/n_pols_to[x])%n_pols][derivative];
-             }
+      if (update_grads)
+        for (unsigned int d=0; d<dim; ++d)
+          {
+            grads[i][d] = 1.;            
+            for (unsigned int x=0; x<dim; ++x)
+              grads[i][d] *= v(x,indices[x])[d==x];
+          }
+
+      if (update_grad_grads)
+        for (unsigned int d1=0; d1<dim; ++d1)
+          for (unsigned int d2=0; d2<dim; ++d2)
+            {
+              grad_grads[i][d1][d2] = 1.;
+              for (unsigned int x=0; x<dim; ++x)
+                {
+                  unsigned int derivative=0;
+                  if (d1==x || d2==x)
+                    {
+                      if (d1==d2)
+                        derivative=2;
+                      else
+                        derivative=1;
+                    } 
+                  grad_grads[i][d1][d2]
+                    *= v(x,indices[x])[derivative];
+                }
+            }
     }
 }
 
 
+
 template<int dim>
 unsigned int
 TensorProductPolynomials<dim>::n() const

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.