From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 7 May 2001 13:20:03 +0000 (+0000)
Subject: Implementation of TensorProductPolynomials::compute_value, compute_grad and compute_g... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9a5aac8c3420066f44105c238126e779e60497ec;p=dealii-svn.git

Implementation of TensorProductPolynomials::compute_value, compute_grad and compute_grad_grad function.


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 3f091d8807..69b7fdf377 100644
--- 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 67443cb4dc..fcec54f1b4 100644
--- 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));