--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2009 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.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__polynomials_adini_h
+#define __deal2__polynomials_adini_h
+
+#include <base/point.h>
+#include <base/tensor.h>
+#include <base/table.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * The cubic polynomial space for the Adini element
+ *
+ * This space consists of the cubic space <i>P<sub>3</sub></i>
+ * augmented by the functions <i>xy<sup>3</sup></i> and
+ * <i>x<sup>3</sup>y</i>.
+ *
+ * The basis of the space is chosen to match the node functionals of
+ * the Adini element.
+ *
+ * @todo This polynomial space is implemented in 2D only.
+ *
+ * @author Bärbel Janssen, 2007
+ */
+
+class AdiniPoly
+{
+ public:
+ /**
+ * Constructor for
+ * the polynomials of
+ * the described space
+ */
+ AdiniPoly ();
+ /**
+ * Computes the value and the
+ * first and second derivatives
+ * of each polynomial at
+ * <tt>unit_point</tt>.
+ *
+ * The size of the vectors must
+ * either be equal 0 or equal
+ * n(). In the first case,
+ * the function will not compute
+ * these values, i.e. you
+ * indicate what you want to have
+ * computed by resizing those
+ * vectors which you want filled.
+ *
+ * If you need values or
+ * derivatives of all polynomials
+ * then use this function, rather
+ * than using any of the
+ * compute_value(),
+ * compute_grad() or
+ * compute_grad_grad()
+ * functions, see below, in a
+ * loop over all polynomials.
+ */
+
+ void compute (const Point<2> &unit_point,
+ std::vector<double> &values,
+ std::vector<Tensor<1,2> > &grads,
+ std::vector< Tensor<2,2> > &grad_grads) const;
+
+ /**
+ * Computes the value of the
+ * <tt>i</tt>th polynomial at
+ * <tt>unit_point</tt>.
+ *
+ * Consider using compute() instead.
+ */
+
+ double compute_value (const unsigned int i,
+ const Point<2> &p) const;
+
+ /**
+ * Computes the gradient of the
+ * <tt>i</tt>th polynomial at
+ * <tt>unit_point</tt>.
+ *
+ * Consider using compute() instead.
+ */
+
+ Tensor<1,2> compute_grad (const unsigned int i,
+ const Point<2> &p) const;
+ /**
+ * Computes the second derivative
+ * (grad_grad) of the <tt>i</tt>th
+ * polynomial at
+ * <tt>unit_point</tt>.
+ *
+ * Consider using compute() instead.
+ */
+
+ Tensor<2,2> compute_grad_grad (const unsigned int i, const Point<2> &p) const;
+ Tensor<2,2> compute_grad_grad_2 (const unsigned int i, const Point<2> &p) const;
+
+ private:
+ /**
+ * Store the coefficients of the
+ * polynominals in the order
+ * \f$1,x,y,x^2,y^2,xy,x^3,y^3,xy^2,x^2y,x^3y,xy^3\f$
+ */
+ Table<2, double> coef;
+
+ /**
+ * Store the coefficients of the x-derivative
+ * of the polynominals in the order
+ * \f$1,x,y,x^2,y^2,xy,x^3,y^3,xy^2,x^2y,x^3y,xy^3\f$
+ */
+
+ Table<2, double> dx;
+ /**
+ * Store the coefficients of the y-derivative
+ * of the polynominals in the order
+ * \f$1,x,y,x^2,y^2,xy,x^3,y^3,xy^2,x^2y,x^3y,xy^3\f$
+ */
+
+ Table<2, double> dy;
+ /**
+ * Store the coefficients of the second x-derivative
+ * of the polynominals in the order
+ * \f$1,x,y,x^2,y^2,xy,x^3,y^3,xy^2,x^2y,x^3y,xy^3\f$
+ */
+ Table<2, double> dxx;
+ /**
+ * Store the coefficients of the second y-derivative
+ * of the polynominals in the order
+ * \f$1,x,y,x^2,y^2,xy,x^3,y^3,xy^2,x^2y,x^3y,xy^3\f$
+ */
+ Table<2, double> dyy;
+ /**
+ * Store the coefficients of the second mixed derivative
+ * of the polynominals in the order
+ * \f$1,x,y,x^2,y^2,xy,x^3,y^3,xy^2,x^2y,x^3y,xy^3\f$
+ */
+ Table<2, double> dxy;
+
+};
+
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 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.
+//
+//---------------------------------------------------------------------------
+
+#include <base/polynomials_adini.h>
+
+#define ENTER_COEFFICIENTS(koefs,z,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11) \
+ koefs(0, z)= a0; koefs(1, z)=a1; koefs(2, z)=a2; koefs(3, z)=a3; koefs( 4, z)=a4 ; koefs( 5, z)=a5 ; \
+ koefs(6, z)= a6; koefs(7, z)=a7; koefs(8, z)=a8; koefs(9, z)=a9; koefs(10, z)=a10; koefs(11, z)=a11;
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+AdiniPoly::AdiniPoly() :
+ coef(12,12),
+ dx(12,12),
+ dy(12,12),
+ dxx(12,12),
+ dyy(12,12),
+ dxy(12,12)
+{
+ // 1 x y xx yy xy 3x 3y xyy xxy 3xy x3y
+ // 0 1 2 3 4 5 6 7 8 9 10 11
+ ENTER_COEFFICIENTS(coef, 0, 1, 0, 0,-3,-3,-1, 2, 2, 3, 3,-2,-2);
+ ENTER_COEFFICIENTS(coef, 1, 0, 1, 0,-2, 0,-1, 1, 0, 0, 2,-1, 0);
+ ENTER_COEFFICIENTS(coef, 2, 0, 0, 1, 0,-2,-1, 0, 1, 2, 0, 0,-1);
+ ENTER_COEFFICIENTS(coef, 3, 0, 0, 0, 3, 0, 1,-2, 0,-3,-3, 2, 2);
+ ENTER_COEFFICIENTS(coef, 4, 0, 0, 0,-1, 0, 0, 1, 0, 0, 1,-1, 0);
+ ENTER_COEFFICIENTS(coef, 5, 0, 0, 0, 0, 0, 1, 0, 0,-2, 0, 0, 1);
+ ENTER_COEFFICIENTS(coef, 6, 0, 0, 0, 0, 3, 1, 0,-2,-3,-3, 2, 2);
+ ENTER_COEFFICIENTS(coef, 7, 0, 0, 0, 0, 0, 1, 0, 0, 0,-2, 1, 0);
+ ENTER_COEFFICIENTS(coef, 8, 0, 0, 0, 0,-1, 0, 0, 1, 1, 0, 0,-1);
+ ENTER_COEFFICIENTS(coef, 9, 0, 0, 0, 0, 0,-1, 0, 0, 3, 3,-2,-2);
+ ENTER_COEFFICIENTS(coef, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0);
+ ENTER_COEFFICIENTS(coef, 11, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 1);
+
+ ENTER_COEFFICIENTS( dx, 0, 0,-6,-1, 6, 3, 6, 0,-2, 0,-6, 0, 0);
+ ENTER_COEFFICIENTS( dx, 1, 1,-4,-1, 3, 0, 4, 0, 0, 0,-3, 0, 0);
+ ENTER_COEFFICIENTS( dx, 2, 0, 0,-1, 0, 2, 0, 0,-1, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dx, 3, 0, 6, 1,-6,-3,-6, 0, 2, 0, 6, 0, 0);
+ ENTER_COEFFICIENTS( dx, 4, 0,-2, 0, 3, 0, 2, 0, 0, 0,-3, 0, 0);
+ ENTER_COEFFICIENTS( dx, 5, 0, 0, 1, 0,-2, 0, 0, 1, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dx, 6, 0, 0, 1, 0,-3,-6, 0, 2, 0, 6, 0, 0);
+ ENTER_COEFFICIENTS( dx, 7, 0, 0, 1, 0, 0,-4, 0, 0, 0, 3, 0, 0);
+ ENTER_COEFFICIENTS( dx, 8, 0, 0, 0, 0, 1, 0, 0,-1, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dx, 9, 0, 0,-1, 0, 3, 6, 0,-2, 0,-6, 0, 0);
+ ENTER_COEFFICIENTS( dx, 10, 0, 0, 0, 0, 0,-2, 0, 0, 0, 3, 0, 0);
+ ENTER_COEFFICIENTS( dx, 11, 0, 0, 0, 0,-1, 0, 0, 1, 0, 0, 0, 0);
+
+ ENTER_COEFFICIENTS( dy, 0, 0,-1,-6, 3, 6, 6,-2, 0,-6, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 1, 0,-1, 0, 2, 0, 0,-1, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 2, 1,-1,-4, 0, 3, 4, 0, 0,-3, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 3, 0, 1, 0,-3, 0,-6, 2, 0, 6, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 4, 0, 0, 0, 1, 0, 0,-1, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 5, 0, 1, 0, 0, 0,-4, 0, 0, 3, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 6, 0, 1, 6,-3,-6,-6, 2, 0, 6, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 7, 0, 1, 0,-2, 0, 0, 1, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 8, 0, 0,-2, 0, 3, 2, 0, 0,-3, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 9, 0,-1, 0, 3, 0, 6,-2, 0,-6, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 10, 0, 0, 0,-1, 0, 0, 1, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dy, 11, 0, 0, 0, 0, 0,-2, 0, 0, 3, 0, 0, 0);
+
+ // 0 1 2 3 4 5 6 7 8 9 10 11
+ ENTER_COEFFICIENTS( dxx, 0,-6, 12, 6, 0, 0,-12, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 1,-4, 6, 4, 0, 0, -6, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 3, 6,-12,-6, 0, 0, 12, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 4,-2, 6, 2, 0, 0, -6, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 6, 0, 0,-6, 0, 0, 12, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 7, 0, 0,-4, 0, 0, 6, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 9, 0, 0, 6, 0, 0,-12, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 10, 0, 0,-2, 0, 0, 6, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxx, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+
+ ENTER_COEFFICIENTS( dyy, 0,-6, 6, 12, 0, 0,-12, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 2,-4, 4, 6, 0, 0, -6, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 3, 0,-6, 0, 0, 0, 12, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 5, 0,-4, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 6, 6,-6,-12, 0, 0, 12, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 8,-2, 2, 6, 0, 0, -6, 0,-0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 9, 0, 6, 0, 0, 0,-12, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dyy, 11, 0,-2, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0);
+
+ ENTER_COEFFICIENTS( dxy, 0,-1, 6, 6,-6,-6, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 1,-1, 4, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 2,-1, 0, 4, 0,-3, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 3, 1,-6,-6, 6, 6, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 4, 0, 2, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 5, 1, 0,-4, 0, 3, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 6, 1,-6,-6, 6, 6, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 7, 1,-4, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 8, 0, 0, 2, 0,-3, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 9,-1, 6, 6,-6,-6, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 10, 0,-2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0);
+ ENTER_COEFFICIENTS( dxy, 11, 0, 0,-2, 0, 3, 0, 0, 0, 0, 0, 0, 0);
+}
+
+void
+AdiniPoly::compute (const Point<2> &unit_point,
+ std::vector<double> &values,
+ std::vector<Tensor<1,2> > &grads,
+ std::vector< Tensor<2,2> > &grad_grads) const
+{
+ if(values.size()>0)//nur dann was tun
+ {
+ for(unsigned int i=0; i<values.size(); ++i)
+ {
+ values[i] = compute_value(i, unit_point);
+ }
+ }
+
+ if(grads.size()>0)//nur dann was tun
+ {
+ for(unsigned int i=0; i<grads.size(); ++i)
+ {
+ grads[i] = compute_grad(i, unit_point);
+ }
+ }
+
+ if(grad_grads.size()>0)//nur dann was tun
+ {
+ for(unsigned int i=0; i<grad_grads.size(); ++i)
+ {
+ grad_grads[i] = compute_grad_grad(i, unit_point);
+ }
+ }
+ return;
+}
+
+double AdiniPoly::compute_value (const unsigned int i, const Point<2> &p) const
+{
+ const double x = p(0);
+ const double y = p(1);
+ return coef( 0,i) +
+ coef( 1,i)*x +
+ coef( 2,i)*y +
+ coef( 3,i)*x*x +
+ coef( 4,i)*y*y +
+ coef( 5,i)*x*y +
+ coef( 6,i)*x*x*x +
+ coef( 7,i)*y*y*y +
+ coef( 8,i)*x*y*y +
+ coef( 9,i)*x*x*y +
+ coef(10,i)*x*x*x*y +
+ coef(11,i)*x*y*y*y;
+}
+
+Tensor<1,2> AdiniPoly::compute_grad (const unsigned int i, const Point<2> &p) const
+{
+ const double x = p(0);
+ const double y = p(1);
+ Tensor<1,2> tensor;
+ tensor[0]= dx( 0,i) +
+ dx( 1,i)*x +
+ dx( 2,i)*y +
+ dx( 3,i)*x*x +
+ dx( 4,i)*y*y +
+ dx( 5,i)*x*y +
+ dx( 6,i)*x*x*x +
+ dx( 7,i)*y*y*y +
+ dx( 8,i)*x*y*y +
+ dx( 9,i)*x*x*y +
+ dx(10,i)*x*x*x*y +
+ dx(11,i)*x*y*y*y;
+
+ tensor[1]= dy( 0,i) +
+ dy( 1,i)*x +
+ dy( 2,i)*y +
+ dy( 3,i)*x*x +
+ dy( 4,i)*y*y +
+ dy( 5,i)*x*y +
+ dy( 6,i)*x*x*x +
+ dy( 7,i)*y*y*y +
+ dy( 8,i)*x*y*y +
+ dy( 9,i)*x*x*y +
+ dy(10,i)*x*x*x*y +
+ dy(11,i)*x*y*y*y;
+ return tensor;
+}
+
+Tensor<2,2> AdiniPoly::compute_grad_grad (const unsigned int i, const Point<2> &p) const
+{
+ const double x = p(0);
+ const double y = p(1);
+ Tensor<2,2> tensor;
+ tensor[0][0]= dxx( 0,i) +
+ dxx( 1,i)*x +
+ dxx( 2,i)*y +
+ dxx( 3,i)*x*x +
+ dxx( 4,i)*y*y +
+ dxx( 5,i)*x*y +
+ dxx( 6,i)*x*x*x +
+ dxx( 7,i)*y*y*y +
+ dxx( 8,i)*x*y*y +
+ dxx( 9,i)*x*x*y +
+ dxx(10,i)*x*x*x*y +
+ dxx(11,i)*x*y*y*y;
+ tensor[0][1]= dxy( 0,i) +
+ dxy( 1,i)*x +
+ dxy( 2,i)*y +
+ dxy( 3,i)*x*x +
+ dxy( 4,i)*y*y +
+ dxy( 5,i)*x*y +
+ dxy( 6,i)*x*x*x +
+ dxy( 7,i)*y*y*y +
+ dxy( 8,i)*x*y*y +
+ dxy( 9,i)*x*x*y +
+ dxy(10,i)*x*x*x*y +
+ dxy(11,i)*x*y*y*y;
+ tensor[1][0]= tensor[0][1];
+ tensor[1][1]= dyy( 0,i) +
+ dyy( 1,i)*x +
+ dyy( 2,i)*y +
+ dyy( 3,i)*x*x +
+ dyy( 4,i)*y*y +
+ dyy( 5,i)*x*y +
+ dyy( 6,i)*x*x*x +
+ dyy( 7,i)*y*y*y +
+ dyy( 8,i)*x*y*y +
+ dyy( 9,i)*x*x*y +
+ dyy(10,i)*x*x*x*y +
+ dyy(11,i)*x*y*y*y;
+ return tensor;
+}
+
+
+DEAL_II_NAMESPACE_CLOSE