]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new class PolynomialSpace
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 11 Feb 2002 12:58:09 +0000 (12:58 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 11 Feb 2002 12:58:09 +0000 (12:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@5500 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial_space.h [new file with mode: 0644]
deal.II/base/include/base/tensor_product_polynomials.h
deal.II/base/source/legendre.cc
deal.II/base/source/polynomial_space.cc [new file with mode: 0644]
deal.II/base/source/tensor_product_polynomials.cc
deal.II/deal.II/source/fe/mapping_q.cc
tests/base/polynomial1d.cc
tests/base/polynomial1d.checked
tests/base/polynomial_test.cc
tests/base/polynomial_test.checked

diff --git a/deal.II/base/include/base/polynomial_space.h b/deal.II/base/include/base/polynomial_space.h
new file mode 100644 (file)
index 0000000..6fd8423
--- /dev/null
@@ -0,0 +1,169 @@
+//----------------------  polynomials.h  -------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------  polynomials.h  -------------
+#ifndef __deal2__polynomial_space_h
+#define __deal2__polynomial_space_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+#include <base/tensor.h>
+#include <base/point.h>
+#include <base/polynomial.h>
+#include <base/smartpointer.h>
+
+#include <vector>
+
+
+/**
+ * Polynomial space of degree at most n in higher dimensions.
+ *
+ * Given a vector of @{n} one-dimensional polynomials @{P0} to @{Pn},
+ * where @{Pi} has degree @p{i}, this class generates all polynomials
+ * 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 below/equal @p{n}.
+ *
+ * @author Guido Kanschat, 2002
+ */
+template <int dim>
+class PolynomialSpace
+{
+  public:
+                                    /**
+                                     * Constructor. @p{pols} is a
+                                     * vector of pointers to
+                                     * one-dimensional polynomials
+                                     * and will be copied into the
+                                     * member variable @p{polynomials}.
+                                     */
+    template <class Pol>
+    PolynomialSpace(const typename std::vector<Pol> &pols);
+
+                                    /**
+                                     * Computes the value and the
+                                     * first and second derivatives
+                                     * of each polynomial at
+                                     * @p{unit_point}.
+                                     *
+                                     * 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.
+                                     *
+                                     * If you need values or
+                                     * derivatives of all 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 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 polynomial at
+                                     * @p{unit_point}.
+                                     *
+                                     * Consider using @p{compute} instead.
+                                     */
+    double compute_value (const unsigned int i,
+                         const Point<dim> &p) const;
+
+                                    /**
+                                     * Computes the gradient of the
+                                     * @p{i}th polynomial at
+                                     * @p{unit_point}.
+                                     *
+                                     * Consider using @p{compute} instead.
+                                     */
+    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
+                                     * polynomial at
+                                     * @p{unit_point}.
+                                     *
+                                     * Consider using @p{compute} instead.
+                                     */
+    Tensor<2,dim> compute_grad_grad(const unsigned int i,
+                                   const Point<dim> &p) const;
+
+                                    /**
+                                     * Returns the number of tensor
+                                     * product polynomials. For $n$
+                                     * 1d polynomials this is $n^dim$.
+                                     */
+    unsigned int n() const;
+
+                                    /**
+                                     * Exception.
+                                     */
+    DeclException3 (ExcDimensionMismatch2,
+                   int, int, int,
+                   << "Dimension " << arg1 << " not equal to " << arg2 << " nor to " << arg3);
+
+           
+  private:
+                                    /**
+                                     * Copy of the vector @p{pols} of
+                                     * polynomials given to the
+                                     * constructor.
+                                     */
+    std::vector<Polynomial<double> > polynomials;
+
+                                    /**
+                                     * Number of tensor product
+                                     * polynomials. For $n$ 1d
+                                     * polynomials this is $n^dim$.
+                                     */
+    unsigned int n_pols;
+    
+                                    /**
+                                     * 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.
+                                     */
+    static unsigned int power(const unsigned int x, const unsigned int y);
+};
+
+
+
+template <int dim>
+template <class Pol>
+PolynomialSpace<dim>::PolynomialSpace(
+  const typename std::vector<Pol> &pols):
+               polynomials (pols.begin(), pols.end())
+{
+  const unsigned int n=polynomials.size();
+
+  n_pols = n;
+  for (unsigned int i=1;i<dim;++i)
+    {
+      n_pols *= (n+i);
+      n_pols /= (i+1);
+    }
+}
+
+
+
+#endif
index f897c5dd8f6afc8298a270075221ced3fdf09d41..24a5786992554d4292ff82839a1151dbb1d3ba44 100644 (file)
@@ -163,7 +163,7 @@ class TensorProductPolynomials
                                      * product polynomials. For $n$
                                      * 1d polynomials this is $n^dim$.
                                      */
-    unsigned int n_tensor_product_polynomials() const;
+    unsigned int n() const;
 
                                     /**
                                      * Exception.
index d6810aa59d50ff48cdbe4a4b360a8d7a8aed22cb..230dcc9a5cc38ead6f2ab4092117f0df224414b0 100644 (file)
@@ -23,7 +23,8 @@
 // a memory checked such as "purify". Maybe, this should be handled somehow
 // to avoid this confusion in future.
 
-
+//TODO:[GK] These polynomials are orthogonal on [-1,1], but the
+//integral over p^2 is not 1. Find out if this is ok.
 
 // Reserve space for polynomials up to degree 19. Should be sufficient
 // for the start.
diff --git a/deal.II/base/source/polynomial_space.cc b/deal.II/base/source/polynomial_space.cc
new file mode 100644 (file)
index 0000000..f42e51a
--- /dev/null
@@ -0,0 +1,186 @@
+//----------------------  polynomial_space.cc  ------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------  polynomial_space.cc  ------------
+
+
+#include <base/exceptions.h>
+#include <base/polynomial_space.h>
+
+
+
+template <int dim>
+unsigned int PolynomialSpace<dim>::power(const unsigned int x,
+                                    const unsigned int y)
+{
+  unsigned int value=1;
+  for (unsigned int i=0; i<y; ++i)
+    value*=x;
+  return value;
+}
+
+
+
+template <int dim>
+double
+PolynomialSpace<dim>::compute_value(const unsigned int i,
+                               const Point<dim> &p) const
+{
+  Assert(false, ExcNotImplemented());
+  return 0.;
+}
+
+  
+template <int dim>
+Tensor<1,dim>
+PolynomialSpace<dim>::compute_grad(const unsigned int i,
+                              const Point<dim> &p) const
+{
+  Assert(false, ExcNotImplemented());
+  return Tensor<1,dim>();
+}
+
+
+template <int dim>
+Tensor<2,dim>
+PolynomialSpace<dim>::compute_grad_grad(const unsigned int i,
+                                   const Point<dim> &p) const
+{
+  Assert(false, ExcNotImplemented());
+  return Tensor<2,dim>();
+}
+
+
+
+
+template <int dim>
+void PolynomialSpace<dim>::compute(
+  const Point<dim>                     &p,
+  std::vector<double>                  &values,
+  typename std::vector<Tensor<1,dim> > &grads,
+  typename std::vector<Tensor<2,dim> > &grad_grads) const
+{
+  const unsigned int n_1d=polynomials.size();
+  
+  Assert(values.size()==n_pols || values.size()==0,
+        ExcDimensionMismatch2(values.size(), n_pols, 0));
+  Assert(grads.size()==n_pols|| grads.size()==0,
+        ExcDimensionMismatch2(grads.size(), n_pols, 0));
+  Assert(grad_grads.size()==n_pols|| grad_grads.size()==0,
+        ExcDimensionMismatch2(grad_grads.size(), n_pols, 0));
+
+  unsigned int v_size=0;
+  bool update_values=false, update_grads=false, update_grad_grads=false;
+  if (values.size()==n_pols)
+    {
+      update_values=true;
+      v_size=1;
+    }
+  if (grads.size()==n_pols)
+    {
+      update_grads=true;
+      v_size=2;
+    }
+  if (grad_grads.size()==n_pols)
+    {
+      update_grad_grads=true;
+      v_size=3;
+    }
+
+                                  // Store data in a single
+                                  // vector. Access is by
+                                  // v[d][n][o]
+                                  //  d: coordinate direction
+                                  //  n: number of 1d polynomial
+                                  //  o: order of derivative
+  std::vector<std::vector<std::vector<double> > >
+    v(dim,
+      std::vector<std::vector<double> > (n_1d,
+                                        std::vector<double> (v_size, 0.)));
+
+  for (unsigned int d=0; d<dim; ++d)
+    {
+      std::vector<std::vector<double> >& v_d=v[d];
+      Assert(v_d.size()==n_1d, ExcInternalError());
+      for (unsigned int i=0; i<n_1d; ++i)
+       polynomials[i].value(p(d), v_d[i]);
+    }
+
+  if (update_values)
+    {
+      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)
+           values[k++] = v[0][ix][0]
+                         * ((dim>1) ? v[1][iy][0] : 1.)
+                         * ((dim>2) ? v[2][iz][0] : 1.);
+    }
+  
+  if (update_grads)
+    {
+      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)
+           {
+             for (unsigned int d=0;d<dim;++d)
+               grads[k][d] = v[0][ix][(d==0) ? 1 : 0]
+    * ((dim>1) ? v[1][iy][(d==1) ? 1 : 0] : 1.)
+    * ((dim>2) ? v[2][iz][(d==2) ? 1 : 0] : 1.);
+             ++k;
+           }
+    }
+
+  if (update_grad_grads)
+    {
+      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)
+           {
+             for (unsigned int d1=0; d1<dim; ++d1)
+               for (unsigned int d2=0; d2<dim; ++d2)
+                 {
+                                                    // Derivative
+                                                    // order for each
+                                                    // direction
+                   const unsigned int
+                     j0 = ((d1==0) ? 1 : 0) + ((d2==0) ? 1 : 0);
+                   const unsigned int
+                     j1 = ((d1==1) ? 1 : 0) + ((d2==1) ? 1 : 0);
+                   const unsigned int
+                     j2 = ((d1==2) ? 1 : 0) + ((d2==2) ? 1 : 0);
+                   
+                   grad_grads[k][d1][d2] = v[0][ix][j0]
+                                           * ((dim>1) ? v[1][iy][j1] : 1.)
+                                           * ((dim>2) ? v[2][iz][j2] : 1.);
+                 }
+             ++k;
+           }
+    }
+}
+
+
+template<int dim>
+unsigned int
+PolynomialSpace<dim>::n() const
+{
+  return n_pols;
+}
+
+  
+template class PolynomialSpace<1>;
+template class PolynomialSpace<2>;
+template class PolynomialSpace<3>;
index 56c71b7a8d11691345e331c35e788adaecd4ea92..e6522ac3fdce6b438acee35f2e6f0d5d9b0ea631 100644 (file)
@@ -202,7 +202,7 @@ void TensorProductPolynomials<dim>::compute(
 
 template<int dim>
 unsigned int
-TensorProductPolynomials<dim>::n_tensor_product_polynomials() const
+TensorProductPolynomials<dim>::n() const
 {
   return n_tensor_pols;
 }
index d062606494bbe5e360690a6600d5b8884a09190d..bdad90bbfe9b2c269358554a365d5fa003b2c655 100644 (file)
@@ -120,7 +120,7 @@ MappingQ<dim>::MappingQ (const unsigned int p):
     v.push_back(LagrangeEquidistant(degree,i));
 
   tensor_pols = new TensorProductPolynomials<dim> (v);
-  Assert (n_shape_functions==tensor_pols->n_tensor_product_polynomials(),
+  Assert (n_shape_functions==tensor_pols->n(),
          ExcInternalError());
   Assert(n_inner+n_outer==n_shape_functions, ExcInternalError());
   
index b40f43618a3407d4f6e3666aa0a3f94a26bead3c..1e9636d5f398948c6d809152e98199ad3299de2a 100644 (file)
@@ -45,6 +45,9 @@ int main ()
   deallog.depth_console(0);
 
   std::vector<Polynomial<double> > p (15);
+
+  deallog << "Legendre" << endl;
+  
   for (unsigned int i=0;i<p.size();++i)
     p[i] = Legendre<double>(i);
   
@@ -52,4 +55,19 @@ int main ()
     for (unsigned int j=0;j<=i;++j)
       deallog << 'P' << i << " * P" << j
              << " =" << scalar_product(p[i], p[j]) << std::endl;
+
+
+  deallog << "LagrangeEquidistant" << endl;
+  
+  p.resize(6);
+  for (unsigned int i=0;i<p.size();++i)
+    p[i] = LagrangeEquidistant(p.size(), i);
+
+                                  // We add 1.0001 bacuse of bugs in
+                                  // the ostream classes
+  for (unsigned int i=0;i<p.size();++i)
+    for (unsigned int j=0;j<p.size();++j)
+      deallog << 'P' << i << "(x" << j
+             << ") =" << p[i].value((double) j/p.size())+1.0001 << std::endl;
 }
+
index 9704110fcf326cdf27dd291b1920ef7279b1947a..31c5c15030532a749bf1ec48bc5b0647458ae46b 100644 (file)
@@ -1,4 +1,5 @@
 
+DEAL::Legendre
 DEAL::P0 * P0 =2.00
 DEAL::P1 * P0 =0.00
 DEAL::P1 * P1 =0.667
@@ -119,3 +120,40 @@ DEAL::P14 * P11 =0.00
 DEAL::P14 * P12 =0.00
 DEAL::P14 * P13 =0.00
 DEAL::P14 * P14 =0.0690
+DEAL::LagrangeEquidistant
+DEAL::P0(x0) =2.00
+DEAL::P0(x1) =1.00
+DEAL::P0(x2) =1.00
+DEAL::P0(x3) =1.00
+DEAL::P0(x4) =1.00
+DEAL::P0(x5) =1.00
+DEAL::P1(x0) =1.00
+DEAL::P1(x1) =2.00
+DEAL::P1(x2) =1.00
+DEAL::P1(x3) =1.00
+DEAL::P1(x4) =1.00
+DEAL::P1(x5) =1.00
+DEAL::P2(x0) =1.00
+DEAL::P2(x1) =1.00
+DEAL::P2(x2) =2.00
+DEAL::P2(x3) =1.00
+DEAL::P2(x4) =1.00
+DEAL::P2(x5) =1.00
+DEAL::P3(x0) =1.00
+DEAL::P3(x1) =1.00
+DEAL::P3(x2) =1.00
+DEAL::P3(x3) =2.00
+DEAL::P3(x4) =1.00
+DEAL::P3(x5) =1.00
+DEAL::P4(x0) =1.00
+DEAL::P4(x1) =1.00
+DEAL::P4(x2) =1.00
+DEAL::P4(x3) =1.00
+DEAL::P4(x4) =2.00
+DEAL::P4(x5) =1.00
+DEAL::P5(x0) =1.00
+DEAL::P5(x1) =1.00
+DEAL::P5(x2) =1.00
+DEAL::P5(x3) =1.00
+DEAL::P5(x4) =1.00
+DEAL::P5(x5) =2.00
index 841a5aea7b390f7698a8a18ba416f5a8bad5e584..ff68d9e637d058d191f4667fb61e5f40a8391c0f 100644 (file)
 
 #include <base/logstream.h>
 #include <base/tensor_product_polynomials.h>
+#include <base/polynomial_space.h>
 
+//using std;
 
 extern "C"
 void abort()
 {}
 
 
+template<int dim, class POLY>
+void check_poly(const Point<dim>& x,
+               const POLY& p)
+{
+  const unsigned int n = p.n();
+  vector<double> values (n);
+  vector<Tensor<1,dim> > gradients(n);
+  vector<Tensor<2,dim> > second(n);
+  
+  p.compute (x, values, gradients, second);
+  
+  for (unsigned int k=0;k<n;++k)
+    {
+      values[k] *= pow(10, dim);
+      gradients[k] *= pow(10, dim);
+      
+      deallog << 'P' << k << "\t= " << values[k]
+             << "\tgradient\t";
+      for (unsigned int d=0;d<dim;++d)
+       deallog << gradients[k][d] << '\t';
+      deallog << "\t2nd\t";
+      for (unsigned int d1=0;d1<dim;++d1)
+       for (unsigned int d2=0;d2<dim;++d2)
+         deallog << second[k][d1][d2] << '\t';
+      deallog << endl;
+    }
+  deallog << endl;
+}
 
-bool equals_delta_ij(double value, unsigned int i, unsigned int j)
+
+template <int dim>
+void
+check_tensor (const vector<Polynomial<double> >& v,
+             const Point<dim>& x)
 {
-  double eps=1e-14;
-  if ((i==j && std::fabs(value-1)<eps) || (i!=j && std::fabs(value)<eps))
-    return true;
-  else
-    return false;
+  deallog.push("Tensor");
+  TensorProductPolynomials<dim> p(v);
+  check_poly (x, p);
+  deallog.pop();
 }
 
-void Q3_4th_shape_function_values_and_grads_dim2(
-  const Point<2> &point, double &v_exact,
-  Tensor<1,2> &grad_exact, Tensor<2,2> &grad2);
 
+template <int dim>
+void
+check_poly (const vector<Polynomial<double> >& v,
+           const Point<dim>& x)
+{
+  deallog.push("Polyno");
+  PolynomialSpace<dim> p(v);
+  check_poly (x, p);
+  deallog.pop();
+}
 
 
+void
+check_dimensions (const vector<Polynomial<double> >& p)
+{
+  deallog.push("1d");
+  check_tensor(p, Point<1>(.5));
+  check_poly(p, Point<1>(.5));
+  deallog.pop();
+  deallog.push("2d");
+  check_tensor(p, Point<2>(.5, .2));
+  check_poly(p, Point<2>(.5, .2));
+  deallog.pop();
+  deallog.push("3d");
+  check_tensor(p, Point<3>(.5, .2, .3));
+  check_poly(p, Point<3>(.5, .2, .3));
+  deallog.pop();
+}
 
 int main()
 {
@@ -47,185 +103,15 @@ int main()
   logfile.precision(4);
   deallog.attach(logfile);
   deallog.depth_console(0);
-  
-  std::vector<double> values(1);
-  deallog << "LagrangeEquidistant polynoms:" << std::endl;
-  for (unsigned int order=1; order<=4; ++order)
-    {
-      deallog << "Polynomial p of order " << order << std::endl;
-      for (unsigned int s_point=0; s_point<=order; ++s_point)
-       {
-         LagrangeEquidistant polynom(order, s_point);
-
-                                          // support points in vertices
-         for (unsigned int i=0; i<=order; ++i)
-           {
-             double x=static_cast<double>(i)/order;
-             polynom.value(x, values);
-             deallog << " p_" << s_point << "(" << x << ")";
-//           deallog << "=" << values[0];
-             if (equals_delta_ij(values[0], s_point, i))
-               deallog << "   ok";
-             else
-               deallog << "   false";
-             deallog << std::endl;
-
-                                              // now also check
-                                              // whether the other
-                                              // @p{value} function
-                                              // returns the same
-                                              // result
-             if (polynom.value(x) != values[0])
-               {
-                 deallog << "The two `value' functions return different results!"
-                         << std::endl;
-                 abort ();
-               };
-           }
-       }
-    }
-
-  deallog << std::endl << "Test derivatives computed by the Horner scheme:" << std::endl;
-  LagrangeEquidistant pol(4, 2);
-  std::vector<double> v_horner(6);
-  for (unsigned int i=0; i<=10; ++i)
-    {
-      double xi=i*0.1;
-      deallog << "x=" << xi << ",    all derivatives: ";
-      std::vector<double> v_exact(6);
-      
-      v_exact[0]=64.0*xi*xi*xi*xi-128.0*xi*xi*xi+76.0*xi*xi-12.0*xi;
-      v_exact[1]=256.0*xi*xi*xi-384.0*xi*xi+152.0*xi-12.0;
-      v_exact[2]=768.0*xi*xi-768.0*xi+152.0;
-      v_exact[3]=1536*xi-768;
-      v_exact[4]=1536;
-      v_exact[5]=0;
-
-      pol.value(xi, v_horner);
-
-      bool ok=true;
-      for (unsigned int i=0; i<v_exact.size(); ++i)
-       {
-//       deallog << "v_horner[i]=" << v_horner[i]
-//            << "   v_exact[i]=" << v_exact[i] << std::endl;
-         if (std::fabs(v_horner[i]-v_exact[i])>1e-12)
-           ok=false;
-       }
-
-      if (ok)
-       deallog << "ok";
-      else
-       deallog << "false";
-
-      deallog << std::endl;
-    }
-
-  if (true)
-    {
-      deallog << std::endl << "Derivatives of a polynomial of degree 0 (a constant function)." << std::endl;      
-      std::vector<double> a_const(1,1.);
-      const Polynomial<double> pol_const(a_const);
-      std::vector<double> exact_values(5,0.);
-      exact_values[0]=1.;
-      std::vector<double> computed_values(5);
-
-      pol_const.value(0.24, computed_values);
-      bool ok=true;
-      for (unsigned int i=0; i<exact_values.size(); ++i)
-       {
-         if (std::fabs(computed_values[i]-exact_values[i])>1e-15)
-           ok=false;
-       }
-
-      if (ok)
-       deallog << "ok";
-      else
-       deallog << "false";
-
-      deallog << std::endl;
-    }
-  
-  
-
-
-  deallog << std::endl << "Test of TensorProductPolynomials:" << std::endl;
-  deallog << "2D Example:" << std::endl;
-  unsigned int p=3,
-   n_tensor_pols=(p+1)*(p+1);
-  std::vector<Polynomial<double> > pols;
-  
-  for (unsigned int i=0; i<=p; ++i)
-    pols.push_back(LagrangeEquidistant(p, i));
-  
-  TensorProductPolynomials<2> tp_pol(pols);
 
-  double v_exact;
-  Tensor<1,2> grad_exact;
-  Tensor<2,2> grad_grad_exact;
-  
-  Point<2> point(0.35,0.62);
-                                  // 4th shape function of Q3<2> is
-                                  // equivalent to its 1st shape
-                                  // function in lexicographical
-                                  // order.
-  Q3_4th_shape_function_values_and_grads_dim2(point, v_exact, grad_exact, grad_grad_exact);
-  
-  unsigned int i=1;
-  double v=tp_pol.compute_value(i, point);
-  Tensor<1,2> grad=tp_pol.compute_grad(i, point);
-  Tensor<2,2> grad_grad=tp_pol.compute_grad_grad(i, point);
-
-  std::vector<double> vs(n_tensor_pols);
-  std::vector<Tensor<1,2> > grads(n_tensor_pols);
-  std::vector<Tensor<2,2> > grad_grads(n_tensor_pols);
-  tp_pol.compute(point, vs, grads, grad_grads);
-
-
-  deallog << "v=" << v << std::endl;
-  deallog << "vs[" << i << "]=" << vs[i] << std::endl;
-  deallog << "v_exact=" << v_exact << std::endl;
-  deallog << "grad=" << grad << std::endl;
-  deallog << "grads[" << i << "]=" << grads[i] << std::endl;
-  deallog << "grad_exact=" << grad_exact << std::endl;
-  for (unsigned int j=0; j<grad_grads[i].dimension; ++j)
-    for (unsigned int k=0; k<grad_grads[i].dimension; ++k)
-      {
-       deallog << "grad_grad[" << j << "][" << k << "]="
-               << grad_grad[j][k] << std::endl;
-       deallog << "grad_grads[" << i<< "][" << j << "][" << k << "]="
-               << grad_grads[i][j][k] << std::endl;
-       deallog << "grad_grad_exact[" << j << "][" << k << "]="
-               << grad_grad_exact[j][k] << std::endl;
-      }
-}
+  vector<Polynomial<double> > p(3);
+  for (unsigned int i=0;i<p.size();++i)
+    p[i] = LagrangeEquidistant(p.size(), i);
 
+  check_dimensions(p);
 
+  for (unsigned int i=0;i<p.size();++i)
+    p[i] = Legendre<double>(i);
 
-void Q3_4th_shape_function_values_and_grads_dim2(
-  const Point<2> &point, double &v_exact, Tensor<1,2> &grad_exact, Tensor<2,2> &grad2)
-{
-                                  // the following functions
-                                  // are taken from fe_lib.cubic.cc
-  const double xi=point(0),
-             eta=point(1);
-  
-  v_exact=9.0*xi-45.0/2.0*xi*xi+27.0/2.0*xi*xi*xi+(
-    -99.0/2.0*xi+495.0/4.0*xi*xi-297.0/4.0*xi*xi*xi)*eta+(81.0*xi-405.0/2.0*xi*xi+243.0/2.0*xi*xi*xi)*
-                eta*eta+(-81.0/2.0*xi+405.0/4.0*xi*xi-243.0/4.0*xi*xi*xi)*eta*eta*eta;
-  
-  grad_exact[0]=9.0-45.0*xi+81.0/2.0*xi*xi+(-99.0/2.0+495.0/2.0*xi-891.0/4.0*xi*xi)*eta+
-               (81.0-405.0*xi+729.0/2.0*xi*xi)*eta*eta+(-81.0/2.0+405.0/2.0*xi-729.0/4.0*xi*xi)*eta*eta*eta;
-  grad_exact[1]=-99.0/2.0*xi+495.0/4.0*xi*xi-297.0/4.0*xi*xi*xi+2.0*(
-    81.0*xi-405.0/2.0*xi*xi+243.0/2.0*xi*xi*xi)*eta+3.0*(
-      -81.0/2.0*xi+405.0/4.0*xi*xi-243.0/4.0*xi*xi*xi)*eta*eta;
-  
-  grad2[0][0] = -45.0+81.0*xi+(495.0/2.0-891.0/2.0*xi)*eta+(-405.0+729.0*xi)*eta*eta+
-               (405.0/2.0-729.0/2.0*xi)*eta*eta*eta;
-  grad2[0][1] = -99.0/2.0+495.0/2.0*xi-891.0/4.0*xi*xi+2.0*(81.0-405.0*xi+729.0/2.0*xi*xi)*eta+
-               3.0*(-81.0/2.0+405.0/2.0*xi-729.0/4.0*xi*xi)*eta*eta;
-  grad2[1][0] = -99.0/2.0+495.0/2.0*xi-891.0/4.0*xi*xi+2.0*(81.0-405.0*xi+729.0/2.0*xi*xi)*eta+
-               3.0*(-81.0/2.0+405.0/2.0*xi-729.0/4.0*xi*xi)*eta*eta;
-  grad2[1][1] = 162.0*xi-405.0*xi*xi+243.0*xi*xi*xi+6.0*(-81.0/2.0*xi+405.0/4.0*xi*xi-
-                                                        243.0/4.0*xi*xi*xi)*eta;
+  check_dimensions(p);
 }
-
index e8f11416e41d3cdfe5266125b5f8c942f50fabb9..414c0120a69dfe02b57a8db1c1eecac62938fd50 100644 (file)
 
-DEAL::LagrangeEquidistant polynoms:
-DEAL::Polynomial p of order 1
-DEAL:: p_0(0.000)   ok
-DEAL:: p_0(1.000)   ok
-DEAL:: p_1(0.000)   ok
-DEAL:: p_1(1.000)   ok
-DEAL::Polynomial p of order 2
-DEAL:: p_0(0.000)   ok
-DEAL:: p_0(0.5000)   ok
-DEAL:: p_0(1.000)   ok
-DEAL:: p_1(0.000)   ok
-DEAL:: p_1(0.5000)   ok
-DEAL:: p_1(1.000)   ok
-DEAL:: p_2(0.000)   ok
-DEAL:: p_2(0.5000)   ok
-DEAL:: p_2(1.000)   ok
-DEAL::Polynomial p of order 3
-DEAL:: p_0(0.000)   ok
-DEAL:: p_0(0.3333)   ok
-DEAL:: p_0(0.6667)   ok
-DEAL:: p_0(1.000)   ok
-DEAL:: p_1(0.000)   ok
-DEAL:: p_1(0.3333)   ok
-DEAL:: p_1(0.6667)   ok
-DEAL:: p_1(1.000)   ok
-DEAL:: p_2(0.000)   ok
-DEAL:: p_2(0.3333)   ok
-DEAL:: p_2(0.6667)   ok
-DEAL:: p_2(1.000)   ok
-DEAL:: p_3(0.000)   ok
-DEAL:: p_3(0.3333)   ok
-DEAL:: p_3(0.6667)   ok
-DEAL:: p_3(1.000)   ok
-DEAL::Polynomial p of order 4
-DEAL:: p_0(0.000)   ok
-DEAL:: p_0(0.2500)   ok
-DEAL:: p_0(0.5000)   ok
-DEAL:: p_0(0.7500)   ok
-DEAL:: p_0(1.000)   ok
-DEAL:: p_1(0.000)   ok
-DEAL:: p_1(0.2500)   ok
-DEAL:: p_1(0.5000)   ok
-DEAL:: p_1(0.7500)   ok
-DEAL:: p_1(1.000)   ok
-DEAL:: p_2(0.000)   ok
-DEAL:: p_2(0.2500)   ok
-DEAL:: p_2(0.5000)   ok
-DEAL:: p_2(0.7500)   ok
-DEAL:: p_2(1.000)   ok
-DEAL:: p_3(0.000)   ok
-DEAL:: p_3(0.2500)   ok
-DEAL:: p_3(0.5000)   ok
-DEAL:: p_3(0.7500)   ok
-DEAL:: p_3(1.000)   ok
-DEAL:: p_4(0.000)   ok
-DEAL:: p_4(0.2500)   ok
-DEAL:: p_4(0.5000)   ok
-DEAL:: p_4(0.7500)   ok
-DEAL:: p_4(1.000)   ok
-DEAL::
-DEAL::Test derivatives computed by the Horner scheme:
-DEAL::x=0.000,    all derivatives: ok
-DEAL::x=0.1000,    all derivatives: ok
-DEAL::x=0.2000,    all derivatives: ok
-DEAL::x=0.3000,    all derivatives: ok
-DEAL::x=0.4000,    all derivatives: ok
-DEAL::x=0.5000,    all derivatives: ok
-DEAL::x=0.6000,    all derivatives: ok
-DEAL::x=0.7000,    all derivatives: ok
-DEAL::x=0.8000,    all derivatives: ok
-DEAL::x=0.9000,    all derivatives: ok
-DEAL::x=1.000,    all derivatives: ok
-DEAL::
-DEAL::Derivatives of a polynomial of degree 0 (a constant function).
-DEAL::ok
-DEAL::
-DEAL::Test of TensorProductPolynomials:
-DEAL::2D Example:
-DEAL::v=-0.02225
-DEAL::vs[1]=-0.02225
-DEAL::v_exact=-0.02225
-DEAL::grad=0.04092 0.4577
-DEAL::grads[1]=0.04092 0.4577
-DEAL::grad_exact=0.04092 0.4577
-DEAL::grad_grad[0][0]=0.3809
-DEAL::grad_grads[1][0][0]=0.3809
-DEAL::grad_grad_exact[0][0]=0.3809
-DEAL::grad_grad[0][1]=-0.8418
-DEAL::grad_grads[1][0][1]=-0.8418
-DEAL::grad_grad_exact[0][1]=-0.8418
-DEAL::grad_grad[1][0]=-0.8418
-DEAL::grad_grads[1][1][0]=-0.8418
-DEAL::grad_grad_exact[1][0]=-0.8418
-DEAL::grad_grad[1][1]=1.225
-DEAL::grad_grads[1][1][1]=1.225
-DEAL::grad_grad_exact[1][1]=1.225
+DEAL:1d:Tensor::P0     = -0.6250       gradient        1.250           2nd     4.500   
+DEAL:1d:Tensor::P1     = 5.625 gradient        -33.75          2nd     -4.500  
+DEAL:1d:Tensor::P2     = 5.625 gradient        33.75           2nd     -4.500  
+DEAL:1d:Tensor::
+DEAL:1d:Polyno::P0     = -0.6250       gradient        1.250           2nd     4.500   
+DEAL:1d:Polyno::P1     = 5.625 gradient        -33.75          2nd     -4.500  
+DEAL:1d:Polyno::P2     = 5.625 gradient        33.75           2nd     -4.500  
+DEAL:1d:Polyno::
+DEAL:2d:Tensor::P0     = -1.400        gradient        2.800   15.25           2nd     1.008   -0.3050 -0.3050 -0.7875 
+DEAL:2d:Tensor::P1     = 12.60 gradient        -75.60  -137.2          2nd     -1.008  8.235   8.235   7.087   
+DEAL:2d:Tensor::P2     = 12.60 gradient        75.60   -137.2          2nd     -1.008  -8.235  -8.235  7.087   
+DEAL:2d:Tensor::P3     = -6.300        gradient        12.60   -10.12          2nd     4.536   0.2025  0.2025  1.800   
+DEAL:2d:Tensor::P4     = 56.70 gradient        -340.2  91.12           2nd     -4.536  -5.467  -5.467  -16.20  
+DEAL:2d:Tensor::P5     = 56.70 gradient        340.2   91.12           2nd     -4.536  5.467   5.467   -16.20  
+DEAL:2d:Tensor::P6     = 1.800 gradient        -3.600  -6.750          2nd     -1.296  0.1350  0.1350  -1.238  
+DEAL:2d:Tensor::P7     = -16.20        gradient        97.20   60.75           2nd     1.296   -3.645  -3.645  11.14   
+DEAL:2d:Tensor::P8     = -16.20        gradient        -97.20  60.75           2nd     1.296   3.645   3.645   11.14   
+DEAL:2d:Tensor::
+DEAL:2d:Polyno::P0     = -1.400        gradient        2.800   15.25           2nd     1.008   -0.3050 -0.3050 -0.7875 
+DEAL:2d:Polyno::P1     = 12.60 gradient        -75.60  -137.2          2nd     -1.008  8.235   8.235   7.087   
+DEAL:2d:Polyno::P2     = 12.60 gradient        75.60   -137.2          2nd     -1.008  -8.235  -8.235  7.087   
+DEAL:2d:Polyno::P3     = -6.300        gradient        12.60   -10.12          2nd     4.536   0.2025  0.2025  1.800   
+DEAL:2d:Polyno::P4     = 56.70 gradient        -340.2  91.12           2nd     -4.536  -5.467  -5.467  -16.20  
+DEAL:2d:Polyno::P5     = 1.800 gradient        -3.600  -6.750          2nd     -1.296  0.1350  0.1350  -1.238  
+DEAL:2d:Polyno::
+DEAL:3d:Tensor::P0     = -0.5390       gradient        1.078   5.871   18.41           2nd     0.03881 -0.01174        -0.03682        -0.01174        -0.03032        -0.2005 -0.03682        -0.2005 -0.1386 
+DEAL:3d:Tensor::P1     = 4.851 gradient        -29.11  -52.84  -165.7          2nd     -0.03881        0.3170  0.9941  0.3170  0.2729  1.805   0.9941  1.805   1.247   
+DEAL:3d:Tensor::P2     = 4.851 gradient        29.11   -52.84  -165.7          2nd     -0.03881        -0.3170 -0.9941 -0.3170 0.2729  1.805   -0.9941 1.805   1.247   
+DEAL:3d:Tensor::P3     = -2.426        gradient        4.851   -3.898  82.84           2nd     0.1746  0.007796        -0.1657 0.007796        0.06930 0.1331  -0.1657 0.1331  -0.6237 
+DEAL:3d:Tensor::P4     = 21.83 gradient        -131.0  35.08   -745.6          2nd     -0.1746 -0.2105 4.474   -0.2105 -0.6237 -1.198  4.474   -1.198  5.613   
+DEAL:3d:Tensor::P5     = 21.83 gradient        131.0   35.08   -745.6          2nd     -0.1746 0.2105  -4.474  0.2105  -0.6237 -1.198  -4.474  -1.198  5.613   
+DEAL:3d:Tensor::P6     = 0.6930        gradient        -1.386  -2.599  -23.67          2nd     -0.04990        0.005198        0.04734 0.005198        -0.04764        0.08876 0.04734 0.08876 0.1782  
+DEAL:3d:Tensor::P7     = -6.237        gradient        37.42   23.39   213.0           2nd     0.04990 -0.1403 -1.278  -0.1403 0.4288  -0.7989 -1.278  -0.7989 -1.604  
+DEAL:3d:Tensor::P8     = -6.237        gradient        -37.42  23.39   213.0           2nd     0.04990 0.1403  1.278   0.1403  0.4288  -0.7989 1.278   -0.7989 -1.604  
+DEAL:3d:Tensor::P9     = -14.55        gradient        29.11   158.5   11.97           2nd     1.048   -0.3170 -0.02394        -0.3170 -0.8186 -0.1304 -0.02394        -0.1304 0.2898  
+DEAL:3d:Tensor::P10    = 131.0 gradient        -785.9  -1427.  -107.7          2nd     -1.048  8.560   0.6464  8.560   7.367   1.173   0.6464  1.173   -2.608  
+DEAL:3d:Tensor::P11    = 131.0 gradient        785.9   -1427.  -107.7          2nd     -1.048  -8.560  -0.6464 -8.560  7.367   1.173   -0.6464 1.173   -2.608  
+DEAL:3d:Tensor::P12    = -65.49        gradient        131.0   -105.2  53.86           2nd     4.715   0.2105  -0.1077 0.2105  1.871   0.08657 -0.1077 0.08657 1.304   
+DEAL:3d:Tensor::P13    = 589.4 gradient        -3536.  947.2   -484.8          2nd     -4.715  -5.683  2.909   -5.683  -16.84  -0.7791 2.909   -0.7791 -11.74  
+DEAL:3d:Tensor::P14    = 589.4 gradient        3536.   947.2   -484.8          2nd     -4.715  5.683   -2.909  5.683   -16.84  -0.7791 -2.909  -0.7791 -11.74  
+DEAL:3d:Tensor::P15    = 18.71 gradient        -37.42  -70.17  -15.39          2nd     -1.347  0.1403  0.03078 0.1403  -1.286  0.05771 0.03078 0.05771 -0.3726 
+DEAL:3d:Tensor::P16    = -168.4        gradient        1010.   631.5   138.5           2nd     1.347   -3.789  -0.8311 -3.789  11.58   -0.5194 -0.8311 -0.5194 3.353   
+DEAL:3d:Tensor::P17    = -168.4        gradient        -1010.  631.5   138.5           2nd     1.347   3.789   0.8311  3.789   11.58   -0.5194 0.8311  -0.5194 3.353   
+DEAL:3d:Tensor::P18    = 1.323 gradient        -2.646  -14.41  -37.17          2nd     -0.09526        0.02882 0.07434 0.02882 0.07442 0.4049  0.07434 0.4049  -0.1638 
+DEAL:3d:Tensor::P19    = -11.91        gradient        71.44   129.7   334.5           2nd     0.09526 -0.7782 -2.007  -0.7782 -0.6698 -3.644  -2.007  -3.644  1.474   
+DEAL:3d:Tensor::P20    = -11.91        gradient        -71.44  129.7   334.5           2nd     0.09526 0.7782  2.007   0.7782  -0.6698 -3.644  2.007   -3.644  1.474   
+DEAL:3d:Tensor::P21    = 5.954 gradient        -11.91  9.568   -167.3          2nd     -0.4287 -0.01914        0.3345  -0.01914        -0.1701 -0.2688 0.3345  -0.2688 -0.7371 
+DEAL:3d:Tensor::P22    = -53.58        gradient        321.5   -86.11  1505.           2nd     0.4287  0.5167  -9.032  0.5167  1.531   2.419   -9.032  2.419   6.634   
+DEAL:3d:Tensor::P23    = -53.58        gradient        -321.5  -86.11  1505.           2nd     0.4287  -0.5167 9.032   -0.5167 1.531   2.419   9.032   2.419   6.634   
+DEAL:3d:Tensor::P24    = -1.701        gradient        3.402   6.379   47.79           2nd     0.1225  -0.01276        -0.09558        -0.01276        0.1169  -0.1792 -0.09558        -0.1792 0.2106  
+DEAL:3d:Tensor::P25    = 15.31 gradient        -91.85  -57.41  -430.1          2nd     -0.1225 0.3445  2.581   0.3445  -1.052  1.613   2.581   1.613   -1.895  
+DEAL:3d:Tensor::P26    = 15.31 gradient        91.85   -57.41  -430.1          2nd     -0.1225 -0.3445 -2.581  -0.3445 -1.052  1.613   -2.581  1.613   -1.895  
+DEAL:3d:Tensor::
+DEAL:3d:Polyno::P0     = -0.5390       gradient        1.078   5.871   18.41           2nd     0.03881 -0.01174        -0.03682        -0.01174        -0.03032        -0.2005 -0.03682        -0.2005 -0.1386 
+DEAL:3d:Polyno::P1     = 4.851 gradient        -29.11  -52.84  -165.7          2nd     -0.03881        0.3170  0.9941  0.3170  0.2729  1.805   0.9941  1.805   1.247   
+DEAL:3d:Polyno::P2     = 4.851 gradient        29.11   -52.84  -165.7          2nd     -0.03881        -0.3170 -0.9941 -0.3170 0.2729  1.805   -0.9941 1.805   1.247   
+DEAL:3d:Polyno::P3     = -2.426        gradient        4.851   -3.898  82.84           2nd     0.1746  0.007796        -0.1657 0.007796        0.06930 0.1331  -0.1657 0.1331  -0.6237 
+DEAL:3d:Polyno::P4     = 21.83 gradient        -131.0  35.08   -745.6          2nd     -0.1746 -0.2105 4.474   -0.2105 -0.6237 -1.198  4.474   -1.198  5.613   
+DEAL:3d:Polyno::P5     = 0.6930        gradient        -1.386  -2.599  -23.67          2nd     -0.04990        0.005198        0.04734 0.005198        -0.04764        0.08876 0.04734 0.08876 0.1782  
+DEAL:3d:Polyno::P6     = -14.55        gradient        29.11   158.5   11.97           2nd     1.048   -0.3170 -0.02394        -0.3170 -0.8186 -0.1304 -0.02394        -0.1304 0.2898  
+DEAL:3d:Polyno::P7     = 131.0 gradient        -785.9  -1427.  -107.7          2nd     -1.048  8.560   0.6464  8.560   7.367   1.173   0.6464  1.173   -2.608  
+DEAL:3d:Polyno::P8     = -65.49        gradient        131.0   -105.2  53.86           2nd     4.715   0.2105  -0.1077 0.2105  1.871   0.08657 -0.1077 0.08657 1.304   
+DEAL:3d:Polyno::P9     = 1.323 gradient        -2.646  -14.41  -37.17          2nd     -0.09526        0.02882 0.07434 0.02882 0.07442 0.4049  0.07434 0.4049  -0.1638 
+DEAL:3d:Polyno::
+DEAL:1d:Tensor::P0     = 10.00 gradient        0.000           2nd     0.000   
+DEAL:1d:Tensor::P1     = 5.000 gradient        10.00           2nd     0.000   
+DEAL:1d:Tensor::P2     = -1.250        gradient        15.00           2nd     3.000   
+DEAL:1d:Tensor::
+DEAL:1d:Polyno::P0     = 10.00 gradient        0.000           2nd     0.000   
+DEAL:1d:Polyno::P1     = 5.000 gradient        10.00           2nd     0.000   
+DEAL:1d:Polyno::P2     = -1.250        gradient        15.00           2nd     3.000   
+DEAL:1d:Polyno::
+DEAL:2d:Tensor::P0     = 100.0 gradient        0.000   0.000           2nd     0.000   0.000   0.000   0.000   
+DEAL:2d:Tensor::P1     = 50.00 gradient        100.0   0.000           2nd     0.000   0.000   0.000   0.000   
+DEAL:2d:Tensor::P2     = -12.50        gradient        150.0   0.000           2nd     3.000   0.000   0.000   0.000   
+DEAL:2d:Tensor::P3     = 20.00 gradient        0.000   100.0           2nd     0.000   0.000   0.000   0.000   
+DEAL:2d:Tensor::P4     = 10.00 gradient        20.00   50.00           2nd     0.000   1.000   1.000   0.000   
+DEAL:2d:Tensor::P5     = -2.500        gradient        30.00   -12.50          2nd     0.6000  1.500   1.500   0.000   
+DEAL:2d:Tensor::P6     = -44.00        gradient        0.000   60.00           2nd     0.000   0.000   0.000   3.000   
+DEAL:2d:Tensor::P7     = -22.00        gradient        -44.00  30.00           2nd     0.000   0.6000  0.6000  1.500   
+DEAL:2d:Tensor::P8     = 5.500 gradient        -66.00  -7.500          2nd     -1.320  0.9000  0.9000  -0.3750 
+DEAL:2d:Tensor::
+DEAL:2d:Polyno::P0     = 100.0 gradient        0.000   0.000           2nd     0.000   0.000   0.000   0.000   
+DEAL:2d:Polyno::P1     = 50.00 gradient        100.0   0.000           2nd     0.000   0.000   0.000   0.000   
+DEAL:2d:Polyno::P2     = -12.50        gradient        150.0   0.000           2nd     3.000   0.000   0.000   0.000   
+DEAL:2d:Polyno::P3     = 20.00 gradient        0.000   100.0           2nd     0.000   0.000   0.000   0.000   
+DEAL:2d:Polyno::P4     = 10.00 gradient        20.00   50.00           2nd     0.000   1.000   1.000   0.000   
+DEAL:2d:Polyno::P5     = -44.00        gradient        0.000   60.00           2nd     0.000   0.000   0.000   3.000   
+DEAL:2d:Polyno::
+DEAL:3d:Tensor::P0     = 1000. gradient        0.000   0.000   0.000           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P1     = 500.0 gradient        1000.   0.000   0.000           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P2     = -125.0        gradient        1500.   0.000   0.000           2nd     3.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P3     = 200.0 gradient        0.000   1000.   0.000           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P4     = 100.0 gradient        200.0   500.0   0.000           2nd     0.000   1.000   0.000   1.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P5     = -25.00        gradient        300.0   -125.0  0.000           2nd     0.6000  1.500   0.000   1.500   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P6     = -440.0        gradient        0.000   600.0   0.000           2nd     0.000   0.000   0.000   0.000   3.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P7     = -220.0        gradient        -440.0  300.0   0.000           2nd     0.000   0.6000  0.000   0.6000  1.500   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P8     = 55.00 gradient        -660.0  -75.00  0.000           2nd     -1.320  0.9000  0.000   0.9000  -0.3750 0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P9     = 300.0 gradient        0.000   0.000   1000.           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Tensor::P10    = 150.0 gradient        300.0   0.000   500.0           2nd     0.000   0.000   1.000   0.000   0.000   0.000   1.000   0.000   0.000   
+DEAL:3d:Tensor::P11    = -37.50        gradient        450.0   0.000   -125.0          2nd     0.9000  0.000   1.500   0.000   0.000   0.000   1.500   0.000   0.000   
+DEAL:3d:Tensor::P12    = 60.00 gradient        0.000   300.0   200.0           2nd     0.000   0.000   0.000   0.000   0.000   1.000   0.000   1.000   0.000   
+DEAL:3d:Tensor::P13    = 30.00 gradient        60.00   150.0   100.0           2nd     0.000   0.3000  0.2000  0.3000  0.000   0.5000  0.2000  0.5000  0.000   
+DEAL:3d:Tensor::P14    = -7.500        gradient        90.00   -37.50  -25.00          2nd     0.1800  0.4500  0.3000  0.4500  0.000   -0.1250 0.3000  -0.1250 0.000   
+DEAL:3d:Tensor::P15    = -132.0        gradient        0.000   180.0   -440.0          2nd     0.000   0.000   0.000   0.000   0.9000  0.6000  0.000   0.6000  0.000   
+DEAL:3d:Tensor::P16    = -66.00        gradient        -132.0  90.00   -220.0          2nd     0.000   0.1800  -0.4400 0.1800  0.4500  0.3000  -0.4400 0.3000  0.000   
+DEAL:3d:Tensor::P17    = 16.50 gradient        -198.0  -22.50  55.00           2nd     -0.3960 0.2700  -0.6600 0.2700  -0.1125 -0.07500        -0.6600 -0.07500        0.000   
+DEAL:3d:Tensor::P18    = -365.0        gradient        0.000   0.000   900.0           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   3.000   
+DEAL:3d:Tensor::P19    = -182.5        gradient        -365.0  0.000   450.0           2nd     0.000   0.000   0.9000  0.000   0.000   0.000   0.9000  0.000   1.500   
+DEAL:3d:Tensor::P20    = 45.62 gradient        -547.5  0.000   -112.5          2nd     -1.095  0.000   1.350   0.000   0.000   0.000   1.350   0.000   -0.3750 
+DEAL:3d:Tensor::P21    = -73.00        gradient        0.000   -365.0  180.0           2nd     0.000   0.000   0.000   0.000   0.000   0.9000  0.000   0.9000  0.6000  
+DEAL:3d:Tensor::P22    = -36.50        gradient        -73.00  -182.5  90.00           2nd     0.000   -0.3650 0.1800  -0.3650 0.000   0.4500  0.1800  0.4500  0.3000  
+DEAL:3d:Tensor::P23    = 9.125 gradient        -109.5  45.62   -22.50          2nd     -0.2190 -0.5475 0.2700  -0.5475 0.000   -0.1125 0.2700  -0.1125 -0.07500        
+DEAL:3d:Tensor::P24    = 160.6 gradient        0.000   -219.0  -396.0          2nd     0.000   0.000   0.000   0.000   -1.095  0.5400  0.000   0.5400  -1.320  
+DEAL:3d:Tensor::P25    = 80.30 gradient        160.6   -109.5  -198.0          2nd     0.000   -0.2190 -0.3960 -0.2190 -0.5475 0.2700  -0.3960 0.2700  -0.6600 
+DEAL:3d:Tensor::P26    = -20.07        gradient        240.9   27.38   49.50           2nd     0.4818  -0.3285 -0.5940 -0.3285 0.1369  -0.06750        -0.5940 -0.06750        0.1650  
+DEAL:3d:Tensor::
+DEAL:3d:Polyno::P0     = 1000. gradient        0.000   0.000   0.000           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Polyno::P1     = 500.0 gradient        1000.   0.000   0.000           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Polyno::P2     = -125.0        gradient        1500.   0.000   0.000           2nd     3.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Polyno::P3     = 200.0 gradient        0.000   1000.   0.000           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Polyno::P4     = 100.0 gradient        200.0   500.0   0.000           2nd     0.000   1.000   0.000   1.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Polyno::P5     = -440.0        gradient        0.000   600.0   0.000           2nd     0.000   0.000   0.000   0.000   3.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Polyno::P6     = 300.0 gradient        0.000   0.000   1000.           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   
+DEAL:3d:Polyno::P7     = 150.0 gradient        300.0   0.000   500.0           2nd     0.000   0.000   1.000   0.000   0.000   0.000   1.000   0.000   0.000   
+DEAL:3d:Polyno::P8     = 60.00 gradient        0.000   300.0   200.0           2nd     0.000   0.000   0.000   0.000   0.000   1.000   0.000   1.000   0.000   
+DEAL:3d:Polyno::P9     = -365.0        gradient        0.000   0.000   900.0           2nd     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   3.000   
+DEAL:3d:Polyno::

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.