--- /dev/null
+//---------------------- 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
* product polynomials. For $n$
* 1d polynomials this is $n^dim$.
*/
- unsigned int n_tensor_product_polynomials() const;
+ unsigned int n() const;
/**
* Exception.
// 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.
--- /dev/null
+//---------------------- 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>;
template<int dim>
unsigned int
-TensorProductPolynomials<dim>::n_tensor_product_polynomials() const
+TensorProductPolynomials<dim>::n() const
{
return n_tensor_pols;
}
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());
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);
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;
}
+
+DEAL::Legendre
DEAL::P0 * P0 =2.00
DEAL::P1 * P0 =0.00
DEAL::P1 * P1 =0.667
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
#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()
{
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);
}
-
-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::