From 39184813db91ac0e62ecaa9f2f852b2d93c222bd Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 21 Apr 2003 16:11:54 +0000
Subject: [PATCH] Add anisotropic polynomials.

git-svn-id: https://svn.dealii.org/trunk@7416 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../include/base/tensor_product_polynomials.h | 232 +++++++++++++++
 .../base/source/tensor_product_polynomials.cc | 268 ++++++++++++++++++
 2 files changed, 500 insertions(+)

diff --git a/deal.II/base/include/base/tensor_product_polynomials.h b/deal.II/base/include/base/tensor_product_polynomials.h
index 7dbcd9d70a..066d66520e 100644
--- a/deal.II/base/include/base/tensor_product_polynomials.h
+++ b/deal.II/base/include/base/tensor_product_polynomials.h
@@ -222,6 +222,223 @@ class TensorProductPolynomials
 
 
 
+
+/**
+ * Anisotropic tensor product of given polynomials.
+ *
+ * Given one-dimensional polynomials @{Px1}, @{Px2}, ... in
+ * x-direction, @{Py1}, @{Py2}, ... in y-direction, and so on, this
+ * class generates polynomials of the form @p{ Qijk(x,y,z) =
+ * Pxi(x)Pyj(y)Pzk(z)}. If the base polynomials are mutually
+ * orthogonal on the interval $[-1,1]$ or $[0,d], then the tensor
+ * product polynomials are orthogonal on $[-1,1]^d$ or $[0,1]^d$,
+ * respectively.
+ *
+ * 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{Px1(x)Py1(y)}, @p{Px2(x)Py1(y)},
+ * @p{Px3(x)Py1(y)}, ..., @p{Px1(x)Py2(y)}, @p{Px2(x)Py2(y)},
+ * @p{Px3(x)Py2(y)}, ..., and likewise in 3d.
+ * 
+ * @author Wolfgang Bangerth 2003
+ */
+template <int dim>
+class AnisotropicPolynomials
+{
+  public:
+				     /**
+				      * Constructor. @p{pols} is a
+				      * table of one-dimensional
+				      * polynomials. The number of
+				      * rows in this table should be
+				      * equal to the space dimension,
+				      * with the elements of each row
+				      * giving the polynomials that
+				      * shall be used in this
+				      * particular coordinate
+				      * direction. These polynomials
+				      * may vary between coordinates,
+				      * as well as their number.
+				      */
+    AnisotropicPolynomials (const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols);
+
+				     /**
+				      * Computes the value and the
+				      * first and second derivatives
+				      * of each tensor product
+				      * polynomial at @p{unit_point}.
+				      *
+				      * 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,
+                  std::vector<Tensor<1,dim> > &grads,
+                  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.
+				      */
+    Tensor<2,dim> compute_grad_grad (const unsigned int i,
+                                     const Point<dim> &p) const;
+
+				     /**
+				      * Returns the number of tensor
+				      * product polynomials. It is the
+				      * product of the number of
+				      * polynomials in each coordinate
+				      * direction.
+				      */
+    unsigned int n() const;
+
+				     /**
+				      * Exception.
+				      */
+    DeclException3 (ExcDimensionMismatch2,
+		    int, int, int,
+		    << "Dimension " << arg1 << " not equal to " << arg2 << " nor to " << arg3);
+                                     /**
+                                      * Exception
+                                      */
+    DeclException1 (ExcInvalidDim,
+                    int,
+                    << "The number of rows in this table must be equal to the "
+                    << "space dimension, but is " << arg1);
+	    
+  private:
+				     /**
+				      * Copy of the vector @p{pols} of
+				      * polynomials given to the
+				      * constructor.
+				      */
+    const std::vector<std::vector<Polynomials::Polynomial<double> > > polynomials;
+
+				     /**
+				      * Number of tensor product
+				      * polynomials. This is
+				      * @p{Nx*Ny*Nz}, or with terms
+				      * dropped if the number of space
+				      * dimensions is less than 3.
+				      */
+    const unsigned int n_tensor_pols;
+
+                                     /**
+                                      * 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;
+    
+				     /**
+				      * Given the input to the
+				      * constructor, compute
+				      * @p{n_tensor_pols}.
+				      */
+    static
+    unsigned int
+    get_n_tensor_pols (const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols);
+};
+
+
+
+
 /* -------------- declaration of explicit specializations --- */
 
 template <>
@@ -265,4 +482,19 @@ TensorProductPolynomials(const std::vector<Pol> &pols)
 
 
 
+template <>
+void
+AnisotropicPolynomials<1>::compute_index(const unsigned int n,
+                                         unsigned int      (&index)[1]) const;
+template <>
+void
+AnisotropicPolynomials<2>::compute_index(const unsigned int n,
+                                         unsigned int      (&index)[2]) const;
+template <>
+void
+AnisotropicPolynomials<3>::compute_index(const unsigned int n,
+                                         unsigned int      (&index)[3]) const;
+
+
+
 #endif
diff --git a/deal.II/base/source/tensor_product_polynomials.cc b/deal.II/base/source/tensor_product_polynomials.cc
index d60ecdea40..83e6dbc8c6 100644
--- a/deal.II/base/source/tensor_product_polynomials.cc
+++ b/deal.II/base/source/tensor_product_polynomials.cc
@@ -18,6 +18,8 @@
 
 
 
+/* ------------------- TensorProductPolynomials -------------- */
+
 template <>
 void
 TensorProductPolynomials<1>::
@@ -240,7 +242,273 @@ TensorProductPolynomials<dim>::n() const
   return n_tensor_pols;
 }
 
+
+
+
+/* ------------------- AnisotropicPolynomials -------------- */
+
+
+template <int dim>
+AnisotropicPolynomials<dim>::
+AnisotropicPolynomials(const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols)
+		:
+		polynomials (pols),
+		n_tensor_pols(get_n_tensor_pols(pols))
+{
+  Assert (pols.size() == dim, ExcInvalidDim(pols.size()));
+  for (unsigned int d=0; d<dim; ++d)
+    Assert (pols[d].size() > 0,
+            ExcMessage ("The number of polynomials must be larger than zero "
+                        "for all coordinate directions."));
+}
+
+
+
+
+template <>
+void
+AnisotropicPolynomials<1>::
+compute_index (const unsigned int i,
+               unsigned int       (&indices)[1]) const
+{
+  Assert (i<polynomials[0].size(), ExcInternalError());
+  indices[0] = i;
+}
+
+
+
+template <>
+void
+AnisotropicPolynomials<2>::
+compute_index (const unsigned int i,
+               unsigned int       (&indices)[2]) const
+{
+  Assert (i < polynomials[0].size()*polynomials[1].size(),
+          ExcInternalError());
+
+  indices[0] = i % polynomials[0].size();
+  indices[1] = i / polynomials[0].size();
+}
+
+
+
+template <>
+void
+AnisotropicPolynomials<3>::
+compute_index (const unsigned int i,
+               unsigned int       (&indices)[3]) const
+{
+  Assert (i < polynomials[0].size()*polynomials[1].size()*polynomials[2].size(),
+          ExcInternalError());
+
+  indices[0] = i % polynomials[0].size();
+  indices[1] = (i/polynomials[0].size()) % polynomials[1].size();
+  indices[2] = i / (polynomials[0].size()*polynomials[1].size());
+}
+
+
+
+template <int dim>
+double
+AnisotropicPolynomials<dim>::compute_value (const unsigned int i,
+                                              const Point<dim> &p) const
+{
+  unsigned int indices[dim];
+  compute_index (i, indices);
+  
+  double value=1.;
+  for (unsigned int d=0; d<dim; ++d)
+    value *= polynomials[d][indices[d]].value(p(d));
+  
+  return value;
+}
+
+  
+template <int dim>
+Tensor<1,dim>
+AnisotropicPolynomials<dim>::compute_grad (const unsigned int i,
+                                             const Point<dim> &p) const
+{
+  unsigned int indices[dim];
+  compute_index (i, indices);
+
+                                   // compute values and
+                                   // uni-directional derivatives at
+                                   // the given point in each
+                                   // co-ordinate direction
+  std::vector<std::vector<double> > v(dim, std::vector<double> (2));
+  for (unsigned int d=0; d<dim; ++d)
+    polynomials[d][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 x=0; x<dim; ++x)
+        grad[d] *= v[x][d==x];
+    }
+  
+  return grad;
+}
+
+
+template <int dim>
+Tensor<2,dim>
+AnisotropicPolynomials<dim>::compute_grad_grad (const unsigned int i,
+                                                  const Point<dim> &p) const
+{
+  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[d][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)
+          {
+            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
+AnisotropicPolynomials<dim>::
+compute (const Point<dim>            &p,
+         std::vector<double>         &values,
+         std::vector<Tensor<1,dim> > &grads,
+         std::vector<Tensor<2,dim> > &grad_grads) const
+{
+  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;
+
+
+                                   // compute the values (and
+                                   // derivatives, if necessary) of
+                                   // all polynomials at this
+                                   // evaluation point
+  std::vector<std::vector<std::vector<double> > > v(dim);
+  for (unsigned int d=0; d<dim; ++d)
+    {
+      v[d].resize (polynomials[d].size());
+      for (unsigned int i=0; i<polynomials[d].size(); ++i)
+        {
+          v[d][i].resize (n_values_and_derivatives, 0.);
+          polynomials[d][i].value(p(d), v[d][i]);
+        };
+    }
   
+  for (unsigned int i=0; i<n_tensor_pols; ++i)
+    {
+                                       // first get the
+                                       // one-dimensional indices of
+                                       // this particular tensor
+                                       // product polynomial
+      unsigned int indices[dim];
+      compute_index (i, indices);
+      
+      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 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 ? 1 : 0];
+          }
+
+      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
+AnisotropicPolynomials<dim>::n() const
+{
+  return n_tensor_pols;
+}
+
+
+template <int dim>
+unsigned int
+AnisotropicPolynomials<dim>::
+get_n_tensor_pols (const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols)
+{
+  unsigned int y = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    y *= pols[d].size();
+  return y;
+}
+
+
+
+/* ------------------- explicit instantiations -------------- */
 template class TensorProductPolynomials<1>;
 template class TensorProductPolynomials<2>;
 template class TensorProductPolynomials<3>;
+
+template class AnisotropicPolynomials<1>;
+template class AnisotropicPolynomials<2>;
+template class AnisotropicPolynomials<3>;
-- 
2.39.5