]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add first implementation of Adini polynomials
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 11 May 2009 10:33:02 +0000 (10:33 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 11 May 2009 10:33:02 +0000 (10:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@18832 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomials_adini.h [new file with mode: 0644]
deal.II/base/source/polynomials_adini.cc [new file with mode: 0644]

diff --git a/deal.II/base/include/base/polynomials_adini.h b/deal.II/base/include/base/polynomials_adini.h
new file mode 100644 (file)
index 0000000..ffa00fd
--- /dev/null
@@ -0,0 +1,158 @@
+//---------------------------------------------------------------------------
+//    $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
diff --git a/deal.II/base/source/polynomials_adini.cc b/deal.II/base/source/polynomials_adini.cc
new file mode 100644 (file)
index 0000000..3685dfd
--- /dev/null
@@ -0,0 +1,243 @@
+//---------------------------------------------------------------------------
+//    $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

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.