From fed6625325bb78704737a9cd4de32ea7d778a9d4 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Mon, 11 May 2009 10:33:02 +0000 Subject: [PATCH] add first implementation of Adini polynomials git-svn-id: https://svn.dealii.org/trunk@18832 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomials_adini.h | 158 ++++++++++++ deal.II/base/source/polynomials_adini.cc | 243 ++++++++++++++++++ 2 files changed, 401 insertions(+) create mode 100644 deal.II/base/include/base/polynomials_adini.h create mode 100644 deal.II/base/source/polynomials_adini.cc diff --git a/deal.II/base/include/base/polynomials_adini.h b/deal.II/base/include/base/polynomials_adini.h new file mode 100644 index 0000000000..ffa00fd920 --- /dev/null +++ b/deal.II/base/include/base/polynomials_adini.h @@ -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 +#include +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * The cubic polynomial space for the Adini element + * + * This space consists of the cubic space P3 + * augmented by the functions xy3 and + * x3y. + * + * 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 + * unit_point. + * + * 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 &values, + std::vector > &grads, + std::vector< Tensor<2,2> > &grad_grads) const; + + /** + * Computes the value of the + * ith polynomial at + * unit_point. + * + * Consider using compute() instead. + */ + + double compute_value (const unsigned int i, + const Point<2> &p) const; + + /** + * Computes the gradient of the + * ith polynomial at + * unit_point. + * + * 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 ith + * polynomial at + * unit_point. + * + * 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 index 0000000000..3685dfda0b --- /dev/null +++ b/deal.II/base/source/polynomials_adini.cc @@ -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 + +#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 &values, + std::vector > &grads, + std::vector< Tensor<2,2> > &grad_grads) const +{ + if(values.size()>0)//nur dann was tun + { + for(unsigned int i=0; i0)//nur dann was tun + { + for(unsigned int i=0; i0)//nur dann was tun + { + for(unsigned int i=0; i &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 -- 2.39.5