dune-localfunctions  2.2.0
pyramidp2localbasis.hh
Go to the documentation of this file.
00001 // -*- tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=2 sw=2 sts=2:
00003 #ifndef DUNE_PYRAMID_P2_LOCALBASIS_HH
00004 #define DUNE_PYRAMID_P2_LOCALBASIS_HH
00005 
00006 #include <dune/common/fmatrix.hh>
00007 
00008 #include <dune/localfunctions/common/localbasis.hh>
00009 
00010 namespace Dune
00011 {
00026   template<class D, class R>
00027   class PyramidP2LocalBasis
00028   {
00029   public:
00031     typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
00032                              Dune::FieldMatrix<R,1,3> > Traits;
00033     
00035     unsigned int size () const
00036     {
00037       return 14;
00038     }
00039     
00041     inline void evaluateFunction (const typename Traits::DomainType& in,
00042                                   std::vector<typename Traits::RangeType>& out) const
00043     { 
00044       out.resize(14);
00045       
00046       // transform to reference element with base [-1,1]^2
00047       const R x = 2.0*in[0] + in[2] - 1.0;
00048       const R y = 2.0*in[1] + in[2] - 1.0;
00049       const R z = in[2];
00050       
00051       if (x > y)
00052         {
00053           // vertices
00054           out[0] = 0.25*(x + z)*(x + z - 1)*(y - z - 1)*(y - z);
00055           out[1] = -0.25*(x + z)*(y - z)*((x + z + 1)*(-y + z + 1) - 4*z) - z*(x - y);
00056           out[2] = 0.25*(x + z)*(y - z)*(y - z + 1)*(x + z - 1);
00057           out[3] = 0.25*(y - z)*(x + z)*(y - z + 1)*(x + z + 1);
00058           out[4] = z*(2*z - 1);
00059           
00060           // lower edges
00061           out[5] = -0.5*(y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
00062           out[6] = -0.5*(y - z + 1)*(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1));
00063           out[7] = -0.5*(x + z - 1)*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1));
00064           out[8] = -0.5*(y - z + 1)*(x + z - 1)*(x + 1)*y;
00065           
00066           // upper edges
00067           out[9] = z*(x + z - 1)*(y - z - 1);
00068           out[10] = -z*((x + z + 1)*(y - z - 1) + 4*z);
00069           out[11] = -z*(y - z + 1)*(x + z - 1);
00070           out[12] = z*(y - z + 1)*(x + z + 1);
00071           
00072           // base face
00073           out[13] = (y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
00074         }
00075       else
00076         {
00077           // vertices
00078           out[0] = 0.25*(y + z)*(y + z - 1)*(x - z - 1)*(x - z);
00079           out[1] = -0.25*(x - z)*(y + z)*(x - z + 1)*(-y - z + 1);
00080           out[2] = 0.25*(x - z)*(y + z)*((x - z - 1)*(y + z + 1) + 4*z) + z*(x - y);
00081           out[3] = 0.25*(y + z)*(x - z)*(x - z + 1)*(y + z + 1);
00082           out[4] = z*(2*z - 1);
00083           
00084           // lower edges
00085           out[5] = -0.5*(y + z - 1)*(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1));
00086           out[6] = -0.5*(x - z + 1)*(y + z - 1)*(y + 1)*x;
00087           out[7] = -0.5*(x - z + 1)*(y + z - 1)*(x - 1)*y;
00088           out[8] = -0.5*(x - z + 1)*(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1));
00089           
00090           // upper edges
00091           out[9] = z*(y + z - 1)*(x - z - 1);
00092           out[10] = -z*(x - z + 1)*(y + z - 1);
00093           out[11] = -z*((y + z + 1)*(x - z - 1) + 4*z);
00094           out[12] = z*(x - z + 1)*(y + z + 1);
00095           
00096           // base face
00097           out[13] = (x - z + 1)*(y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1));
00098         }
00099     }
00100     
00102     inline void 
00103     evaluateJacobian (const typename Traits::DomainType& in,         // position
00104                       std::vector<typename Traits::JacobianType>& out) const      // return value
00105     {  
00106       out.resize(14);
00107       
00108       // transform to reference element with base [-1,1]^2
00109       const R x = 2.0*in[0] + in[2] - 1.0;
00110       const R y = 2.0*in[1] + in[2] - 1.0;
00111       const R z = in[2];
00112       
00113       // transformation of the gradient leads to a multiplication
00114       // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
00115       if (x > y)
00116         {
00117           // vertices
00118           out[0][0][0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
00119           out[0][0][1] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
00120           out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
00121             + 0.25*((2*x + 2*z - 1)*(y - z - 1)*(y - z)
00122                     + (x + z)*(x + z - 1)*(-2*y + 2*z + 1));
00123           
00124           out[1][0][0] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
00125                                    + (x + z)*(y - z)*(-y + z + 1)) - z);
00126           out[1][0][1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
00127                                    + (x + z)*(y - z)*(-(x + z + 1))) + z);
00128           out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
00129             - 0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
00130                     - (x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
00131                     + (x + z)*(y - z)*(x - y + 2*z - 2))
00132             - (x - y);
00133           
00134           out[2][0][0] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
00135           out[2][0][1] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
00136           out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
00137             + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z - 1)
00138                     + (x + z)*(y - z)*(y - x - 2*z + 2));
00139           
00140           out[3][0][0] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
00141           out[3][0][1] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
00142           out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
00143             + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z + 1)
00144                     + (y - z)*(x + z)*(y - x - 2*z));
00145           
00146           out[4][0][0] = 0;
00147           out[4][0][1] = 0;
00148           out[4][0][2] = 4*z - 1;
00149           
00150           // lower edges
00151           out[5][0][0] = -((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
00152                            + (y - z + 1)*(x + z - 1)*((y - 1) + z));
00153           out[5][0][1] = -((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
00154                            + (y - z + 1)*(x + z - 1)*((x + 1) - z));
00155           out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
00156             - 0.5*((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
00157                    + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
00158           
00159           out[6][0][0] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
00160           out[6][0][1] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1)
00161                            + (y - z + 1)*((x + z + 1)*x + 2*z));
00162           out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
00163             - 0.5*(-(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1))
00164                    + (y - z + 1)*(((y - 1)*x - 1) + 2*y + 1));
00165           
00166           out[7][0][0] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
00167                            + (x + z - 1)*((y - z - 1)*y + 2*z));
00168           out[7][0][1] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
00169           out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
00170             - 0.5*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
00171                    + (x + z - 1)*((-(x + 1)*y - 1) + 2*x + 1));
00172           
00173           out[8][0][0] = -(y - z + 1)*(2*x + z)*y;
00174           out[8][0][1] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
00175           out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
00176             - 0.5*(-x + y - 2*z + 2)*(x + 1)*y;
00177           
00178           // upper edges
00179           out[9][0][0] = 2*z*(y - z - 1);
00180           out[9][0][1] = 2*z*(x + z - 1);
00181           out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
00182             + (x + z - 1)*(y - z - 1) + z*(-x + y - 2*z);
00183           
00184           out[10][0][0] = -2*z*(y - z - 1);
00185           out[10][0][1] = -2*z*(x + z + 1);
00186           out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
00187             - ((x + z + 1)*(y - z - 1) + 4*z)
00188             - z*(-x + y - 2*z + 2);
00189           
00190           out[11][0][0] = -2*z*(y - z + 1);
00191           out[11][0][1] = -2*z*(x + z - 1);
00192           out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
00193             - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
00194           
00195           out[12][0][0] = 2*z*(y - z + 1);
00196           out[12][0][1] = 2*z*(x + z + 1);
00197           out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
00198             + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
00199           
00200           // base face
00201           out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
00202                              + (y - z + 1)*(x + z - 1)*(y - 1 + z));
00203           out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
00204                              + (y - z + 1)*(x + z - 1)*(x + 1 - z));
00205           out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
00206             + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
00207                + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
00208         }
00209       else
00210         {
00211           // vertices
00212           out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
00213           out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
00214           out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
00215             + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
00216                     + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
00217           
00218           out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
00219           out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
00220           out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
00221             - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
00222                     + (x - z)*(y + z)*(-x + y + 2*z - 2));
00223           
00224           out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
00225                               + (x - z)*(y + z)*(y + z + 1) + 4*z);
00226           out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
00227                               + (x - z)*(y + z)*(x - z - 1) - 4*z);
00228           out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
00229             + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
00230                     + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
00231           
00232           out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
00233           out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
00234           out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
00235             + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
00236                     + (y + z)*(x - z)*(x - y - 2*z));
00237           
00238           out[4][0][0] = 0;
00239           out[4][0][1] = 0;
00240           out[4][0][2] = 4*z - 1;
00241           
00242           // lower edges
00243           out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
00244           out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
00245                            + (y + z - 1)*((x - z - 1)*x + 2*z));
00246           out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
00247             - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
00248                    + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
00249           
00250           out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
00251           out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
00252           out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
00253             - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
00254           
00255           out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
00256           out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
00257           out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
00258             - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
00259           
00260           out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
00261                            + (x - z + 1)*((y + z + 1)*y + 2*z));
00262           out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
00263           out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
00264             - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
00265                    + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
00266           
00267           // upper edges
00268           out[9][0][0] = 2*z*(y + z - 1);
00269           out[9][0][1] = 2*z*(x - z - 1);
00270           out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
00271               + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
00272           
00273           out[10][0][0] = -2*z*(y + z - 1);
00274           out[10][0][1] = -2*z*(x - z + 1);
00275           out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
00276             - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
00277           
00278           out[11][0][0] = -2*z*(y + z + 1);
00279           out[11][0][1] = -2*z*(x - z - 1);
00280           out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
00281             - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
00282           
00283           out[12][0][0] = 2*z*(y + z + 1);
00284           out[12][0][1] = 2*z*(x - z + 1);
00285           out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
00286             + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
00287           
00288           // base face
00289           out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
00290                              + (x - z + 1)*(y + z - 1)*(y + 1 - z));
00291           out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
00292                              + (x - z + 1)*(y + z - 1)*(x - 1 + z));
00293           out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
00294             + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
00295             + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
00296         }
00297     }
00298     
00300     unsigned int order () const
00301     {
00302       return 2;
00303     }
00304   };
00305 }
00306 #endif