Implementation of TensorProductPolynomials::compute_value, compute_grad and compute_g...

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

diff --git a/deal.II/base/include/base/tensor_product_polynomials.h b/deal.II/base/include/base/tensor_product_polynomials.h
index 3f091d88077cb371e209d70043e423d6dacc8ce0..69b7fdf3773146c26469f46a910f49d1fed03a02 100644 (file)
--- a/deal.II/base/include/base/tensor_product_polynomials.h
+++ b/deal.II/base/include/base/tensor_product_polynomials.h
@@ -47,20 +47,115 @@ class TensorProductPolynomials
     TensorProductPolynomials(const std::vector<Pol> &pols);
 
                                     /**
-                                     * Calculates the polynomials
-                                     * and their derivatives at
-                                     * @p{unit_point}.
+                                     * Computes the value and the
+                                     * first and second derivatives
+                                     * of each tensor product
+                                     * polynomial at @p{unit_point}.
                                      *
-                                     * The vectors must either have
-                                     * length @p{0} or number of
-                                     * polynomials. In the first
-                                     * case, the function will not
-                                     * compute these values.
+                                     * The size of the vectors must
+                                     * either be equal @p{0} or equal
+                                     * @p{n_tensor_pols}.  In the
+                                     * first case, the function will
+                                     * not compute these values.
+                                     *
+                                     * If you need values or
+                                     * derivatives of all tensor
+                                     * product polynomials then use
+                                     * this function, rather than
+                                     * using any of the
+                                     * @p{compute_value},
+                                     * @p{compute_grad} or
+                                     * @p{compute_grad_grad}
+                                     * functions, see below, in a
+                                     * loop over all tensor product
+                                     * polynomials.
+                                     */
+    void compute(const Point<dim>                     &unit_point,
+                std::vector<double>                  &values,
+                typename std::vector<Tensor<1,dim> > &grads,
+                typename std::vector<Tensor<2,dim> > &grad_grads) const;
+    
+                                    /**
+                                     * Computes the value of the
+                                     * @p{i}th tensor product
+                                     * polynomial at
+                                     * @p{unit_point}. Here @p{i} is
+                                     * given in tensor product
+                                     * numbering.
+                                     *
+                                     * Note, that using this function
+                                     * within a loop over all tensor
+                                     * product polynomials is not
+                                     * efficient, because then each
+                                     * point value of the underlying
+                                     * (one-dimensional) polynomials
+                                     * is (unnecessarily) computed
+                                     * several times.  Instead use
+                                     * the @p{compute} function, see
+                                     * above, with
+                                     * @p{values.size()==n_tensor_pols}
+                                     * to get the point values of all
+                                     * tensor polynomials all at once
+                                     * and in a much more efficient
+                                     * way.
+                                     */
+    double compute_value (const unsigned int i,
+                         const Point<dim> &p) const;
+
+                                    /**
+                                     * Computes the grad of the
+                                     * @p{i}th tensor product
+                                     * polynomial at
+                                     * @p{unit_point}. Here @p{i} is
+                                     * given in tensor product
+                                     * numbering.
+                                     *
+                                     * Note, that using this function
+                                     * within a loop over all tensor
+                                     * product polynomials is not
+                                     * efficient, because then each
+                                     * derivative value of the
+                                     * underlying (one-dimensional)
+                                     * polynomials is (unnecessarily)
+                                     * computed several times.
+                                     * Instead use the @p{compute}
+                                     * function, see above, with
+                                     * @p{grads.size()==n_tensor_pols}
+                                     * to get the point value of all
+                                     * tensor polynomials all at once
+                                     * and in a much more efficient
+                                     * way.
+                                     */
+    Tensor<1,dim> compute_grad (const unsigned int i,
+                               const Point<dim> &p) const;
+
+                                    /**
+                                     * Computes the second
+                                     * derivative (grad_grad) of the
+                                     * @p{i}th tensor product
+                                     * polynomial at
+                                     * @p{unit_point}. Here @p{i} is
+                                     * given in tensor product
+                                     * numbering.
+                                     *
+                                     * Note, that using this function
+                                     * within a loop over all tensor
+                                     * product polynomials is not
+                                     * efficient, because then each
+                                     * derivative value of the
+                                     * underlying (one-dimensional)
+                                     * polynomials is (unnecessarily)
+                                     * computed several times.
+                                     * Instead use the @p{compute}
+                                     * function, see above, with
+                                     * @p{grad_grads.size()==n_tensor_pols}
+                                     * to get the point value of all
+                                     * tensor polynomials all at once
+                                     * and in a much more efficient
+                                     * way.
                                      */
-    void compute (const Point<dim>                     &unit_point,
-                 std::vector<double>                  &values,
-                 typename std::vector<Tensor<1,dim> > &grads,
-                 typename std::vector<Tensor<2,dim> > &grad_grads) const;
+    Tensor<2,dim> compute_grad_grad(const unsigned int i,
+                                   const Point<dim> &p) const;
 
                                     /**
                                      * Returns the number of tensor
@@ -92,21 +187,35 @@ 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;
+    
+    
     static unsigned int power(const unsigned int x, const unsigned int 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_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());
+}
 
 
 

diff --git a/deal.II/base/source/tensor_product_polynomials.cc b/deal.II/base/source/tensor_product_polynomials.cc
index 67443cb4dcb956ec71d022fc5062cbe14ff7bbd0..fcec54f1b4bc74f7fa548942ff13ef8da6d12bc1 100644 (file)
--- a/deal.II/base/source/tensor_product_polynomials.cc
+++ b/deal.II/base/source/tensor_product_polynomials.cc
@@ -17,6 +17,7 @@
 
 
 
+
 template <int dim>
 unsigned int TensorProductPolynomials<dim>::power(const unsigned int x,
                                                  const unsigned int y)
@@ -28,6 +29,84 @@ unsigned int TensorProductPolynomials<dim>::power(const unsigned int x,
 }
 
 
+
+template <int dim>
+double
+TensorProductPolynomials<dim>::compute_value(const unsigned int i,
+                                            const Point<dim> &p) const
+{
+  const unsigned int n_pols=polynomials.size();
+  
+  double value=1.;
+  for (unsigned int d=0; d<dim; ++d)
+    value *= polynomials[(i/n_pols_to[d])%n_pols].value(p(d));
+  
+  return value;
+}
+
+  
+template <int dim>
+Tensor<1,dim>
+TensorProductPolynomials<dim>::compute_grad(const unsigned int i,
+                                           const Point<dim> &p) const
+{
+  const unsigned int n_pols=polynomials.size();
+  
+  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]);
+  
+  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];
+  
+  return grad;
+}
+
+
+template <int dim>
+Tensor<2,dim>
+TensorProductPolynomials<dim>::compute_grad_grad(const unsigned int i,
+                                                const Point<dim> &p) const
+{
+  const unsigned int n_pols=polynomials.size();
+    
+  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]);
+  
+  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];
+       }
+
+  return grad_grad;
+}
+
+
+
+
 template <int dim>
 void TensorProductPolynomials<dim>::compute(
   const Point<dim>                     &p,
@@ -35,12 +114,7 @@ void TensorProductPolynomials<dim>::compute(
   typename std::vector<Tensor<1,dim> > &grads,
   typename std::vector<Tensor<2,dim> > &grad_grads) const
 {
-  unsigned int n_pols=polynomials.size();
-  std::vector<unsigned int> n_pols_to(dim+1);
-  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());
+  const unsigned int n_pols=polynomials.size();
   
   Assert(values.size()==n_tensor_pols || values.size()==0,
         ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));