From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 21 Apr 2003 16:09:21 +0000 (+0000)
Subject: Clean up some real old and unintelligible cruft.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c060d0364c910fb54f4ebfba64c8ca3cfbde9b9a;p=dealii-svn.git

Clean up some real old and unintelligible cruft.


git-svn-id: https://svn.dealii.org/trunk@7415 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/base/include/base/polynomial_space.h b/deal.II/base/include/base/polynomial_space.h
index ad06c30f7c..0151d48dae 100644
--- a/deal.II/base/include/base/polynomial_space.h
+++ b/deal.II/base/include/base/polynomial_space.h
@@ -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>
diff --git a/deal.II/base/include/base/tensor_product_polynomials.h b/deal.II/base/include/base/tensor_product_polynomials.h
index 92d05e3caa..7dbcd9d70a 100644
--- a/deal.II/base/include/base/tensor_product_polynomials.h
+++ b/deal.II/base/include/base/tensor_product_polynomials.h
@@ -19,7 +19,6 @@
 #include <base/tensor.h>
 #include <base/point.h>
 #include <base/polynomial.h>
-#include <base/smartpointer.h>
 
 #include <vector>
 
@@ -34,7 +33,13 @@
  * $[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()))
+{}
 
 
 
diff --git a/deal.II/base/source/polynomial_space.cc b/deal.II/base/source/polynomial_space.cc
index 85b084e084..1debf26e89 100644
--- a/deal.II/base/source/polynomial_space.cc
+++ b/deal.II/base/source/polynomial_space.cc
@@ -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();
   
diff --git a/deal.II/base/source/tensor_product_polynomials.cc b/deal.II/base/source/tensor_product_polynomials.cc
index 6c9a62c1d3..d60ecdea40 100644
--- a/deal.II/base/source/tensor_product_polynomials.cc
+++ b/deal.II/base/source/tensor_product_polynomials.cc
@@ -18,29 +18,60 @@
 
 
 
+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