dune-localfunctions
2.2.0
|
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