1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 2 // vi: set et ts=4 sw=2 sts=2: 3 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS4_CUBE2D_LOCALBASIS_HH 4 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS4_CUBE2D_LOCALBASIS_HH 5 6 #include <bitset> 7 #include <numeric> 8 #include <vector> 9 10 #include <dune/common/fmatrix.hh> 11 12 #include "../../common/localbasis.hh" 13 14 namespace Dune 15 { 16 /** 17 * \ingroup LocalBasisImplementation 18 * \brief Second order Raviart-Thomas shape functions on the reference quadrilateral. 19 * 20 * \tparam D Type to represent the field in the domain. 21 * \tparam R Type to represent the field in the range. 22 * 23 * \nosubgrouping 24 */ 25 template<class D, class R> 26 class RT4Cube2DLocalBasis 27 { 28 29 public: 30 typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,2,Dune::FieldVector<R,2>, 31 Dune::FieldMatrix<R,2,2> > Traits; 32 33 /** 34 * \brief Make set number s, where 0 <= s < 16 35 * 36 * \param s Edge orientation indicator 37 */ RT4Cube2DLocalBasis(std::bitset<4> s=0)38 RT4Cube2DLocalBasis (std::bitset<4> s = 0) 39 { 40 sign0 = (s[0]) ? -1.0 : 1.0; 41 sign1 = (s[1]) ? -1.0 : 1.0; 42 sign2 = (s[2]) ? -1.0 : 1.0; 43 sign3 = (s[3]) ? -1.0 : 1.0; 44 } 45 46 //! \brief number of shape functions size() const47 unsigned int size () const 48 { 49 return 60; 50 } 51 52 /** 53 * \brief Evaluate all shape functions 54 * 55 * \param in Position 56 * \param out return value 57 */ evaluateFunction(const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const58 inline void evaluateFunction (const typename Traits::DomainType& in, 59 std::vector<typename Traits::RangeType>& out) const 60 { 61 out.resize(60); 62 63 auto const& x = in[0], y = in[1]; 64 65 const auto l1_x = 2*x - 1; 66 const auto l2_x = x*(6*x - 6) + 1; 67 const auto l3_x = x*(x*(20*x - 30) + 12) - 1; 68 const auto l4_x = x*(x*(x*(70*x - 140) + 90) - 20) + 1; 69 const auto l5_x = x*(x*(x*(x*(252*x - 630) + 560) - 210) + 30) - 1; 70 const auto l1_y = 2*y - 1; 71 const auto l2_y = y*(6*y - 6) + 1; 72 const auto l3_y = y*(y*(20*y - 30) + 12) - 1; 73 const auto l4_y = y*(y*(y*(70*y - 140) + 90) - 20) + 1; 74 const auto l5_y = y*(y*(y*(y*(252*y - 630) + 560) - 210) + 30) - 1; 75 76 out[0][0]=sign0*(0.5*(-l4_x)+0.5*l5_x); 77 out[0][1]=0.0; 78 out[1][0]=-(1.5)*l4_x*l1_y+1.5*l5_x*l1_y; 79 out[1][1]=0.0; 80 out[2][0]=sign0*(-(2.5)*l4_x*l2_y+2.5*l5_x*l2_y); 81 out[2][1]=0.0; 82 out[3][0]=-(3.5)*l4_x*l3_y+3.5*l5_x*l3_y; 83 out[3][1]=0.0; 84 out[4][0]=sign0*(-(4.5)*l4_x*l4_y+4.5*l5_x*l4_y); 85 out[4][1]=0.0; 86 87 out[5][0]=sign1*(0.5*l4_x+0.5*l5_x); 88 out[5][1]=0.0; 89 out[6][0]=-(1.5)*l4_x*l1_y-1.5*l5_x*l1_y; 90 out[6][1]=0.0; 91 out[7][0]=sign1*(2.5*l4_x*l2_y+2.5*l5_x*l2_y); 92 out[7][1]=0.0; 93 out[8][0]=-(3.5)*l4_x*l3_y-3.5*l5_x*l3_y; 94 out[8][1]=0.0; 95 out[9][0]=sign1*(4.5*l4_x*l4_y+4.5*l5_x*l4_y); 96 out[9][1]=0.0; 97 98 out[10][0]=0.0; 99 out[10][1]=sign2*(0.5*(-l4_y)+0.5*l5_y); 100 out[11][0]=0.0; 101 out[11][1]=1.5*l1_x*l4_y-1.5*l1_x*l5_y; 102 out[12][0]=0.0; 103 out[12][1]=sign2*(-(2.5)*l2_x*l4_y+2.5*l2_x*l5_y); 104 out[13][0]=0.0; 105 out[13][1]=3.5*l3_x*l4_y-3.5*l3_x*l5_y; 106 out[14][0]=0.0; 107 out[14][1]=sign2*(-(4.5)*l4_x*l4_y+4.5*l4_x*l5_y); 108 109 out[15][0]=0.0; 110 out[15][1]=sign3*(0.5*l4_y+0.5*l5_y); 111 out[16][0]=0.0; 112 out[16][1]=1.5*l1_x*l4_y+1.5*l1_x*l5_y; 113 out[17][0]=0.0; 114 out[17][1]=sign3*(2.5*l2_x*l4_y+2.5*l2_x*l5_y); 115 out[18][0]=0.0; 116 out[18][1]=3.5*l3_x*l4_y+3.5*l3_x*l5_y; 117 out[19][0]=0.0; 118 out[19][1]=sign3*(4.5*l4_x*l4_y+4.5*l4_x*l5_y); 119 120 out[20][0]=1.0-l4_x; 121 out[20][1]=0.0; 122 out[21][0]=3.0*l1_y-3.0*l4_x*l1_y; 123 out[21][1]=0.0; 124 out[22][0]=5.0*l2_y-5.0*l4_x*l2_y; 125 out[22][1]=0.0; 126 out[23][0]=7.0*l3_y-7.0*l4_x*l3_y; 127 out[23][1]=0.0; 128 out[24][0]=9.0*l4_y-9.0*l4_x*l4_y; 129 out[24][1]=0.0; 130 out[25][0]=3.0*l1_x-3.0*l5_x; 131 out[25][1]=0.0; 132 out[26][0]=9.0*l1_x*l1_y-9.0*l5_x*l1_y; 133 out[26][1]=0.0; 134 out[27][0]=15.0*l1_x*l2_y-15.0*l5_x*l2_y; 135 out[27][1]=0.0; 136 out[28][0]=21.0*l1_x*l3_y-21.0*l5_x*l3_y; 137 out[28][1]=0.0; 138 out[29][0]=27.0*l1_x*l4_y-27.0*l5_x*l4_y; 139 out[29][1]=0.0; 140 out[30][0]=5.0*l2_x-5.0*l4_x; 141 out[30][1]=0.0; 142 out[31][0]=15.0*l2_x*l1_y-15.0*l4_x*l1_y; 143 out[31][1]=0.0; 144 out[32][0]=25.0*l2_x*l2_y-25.0*l4_x*l2_y; 145 out[32][1]=0.0; 146 out[33][0]=35.0*l2_x*l3_y-35.0*l4_x*l3_y; 147 out[33][1]=0.0; 148 out[34][0]=45.0*l2_x*l4_y-45.0*l4_x*l4_y; 149 out[34][1]=0.0; 150 out[35][0]=7.0*l3_x-7.0*l5_x; 151 out[35][1]=0.0; 152 out[36][0]=21.0*l3_x*l1_y-21.0*l5_x*l1_y; 153 out[36][1]=0.0; 154 out[37][0]=35.0*l3_x*l2_y-35.0*l5_x*l2_y; 155 out[37][1]=0.0; 156 out[38][0]=49.0*l3_x*l3_y-49.0*l5_x*l3_y; 157 out[38][1]=0.0; 158 out[39][0]=63.0*l3_x*l4_y-63.0*l5_x*l4_y; 159 out[39][1]=0.0; 160 out[40][0]=0.0; 161 out[40][1]=1.0-l4_y; 162 out[41][0]=0.0; 163 out[41][1]=3.0*l1_y-3.0*l5_y; 164 out[42][0]=0.0; 165 out[42][1]=5.0*l2_y-5.0*l4_y; 166 out[43][0]=0.0; 167 out[43][1]=7.0*l3_y-7.0*l5_y; 168 out[44][0]=0.0; 169 out[44][1]=3.0*l1_x-3.0*l1_x*l4_y; 170 out[45][0]=0.0; 171 out[45][1]=9.0*l1_x*l1_y-9.0*l1_x*l5_y; 172 out[46][0]=0.0; 173 out[46][1]=15.0*l1_x*l2_y-15.0*l1_x*l4_y; 174 out[47][0]=0.0; 175 out[47][1]=21.0*l1_x*l3_y-21.0*l1_x*l5_y; 176 out[48][0]=0.0; 177 out[48][1]=5.0*l2_x-5.0*l2_x*l4_y; 178 out[49][0]=0.0; 179 out[49][1]=15.0*l2_x*l1_y-15.0*l2_x*l5_y; 180 out[50][0]=0.0; 181 out[50][1]=25.0*l2_x*l2_y-25.0*l2_x*l4_y; 182 out[51][0]=0.0; 183 out[51][1]=35.0*l2_x*l3_y-35.0*l2_x*l5_y; 184 out[52][0]=0.0; 185 out[52][1]=7.0*l3_x-7.0*l3_x*l4_y; 186 out[53][0]=0.0; 187 out[53][1]=21.0*l3_x*l1_y-21.0*l3_x*l5_y; 188 out[54][0]=0.0; 189 out[54][1]=35.0*l3_x*l2_y-35.0*l3_x*l4_y; 190 out[55][0]=0.0; 191 out[55][1]=49.0*l3_x*l3_y-49.0*l3_x*l5_y; 192 out[56][0]=0.0; 193 out[56][1]=9.0*l4_x-9.0*l4_x*l4_y; 194 out[57][0]=0.0; 195 out[57][1]=27.0*l4_x*l1_y-27.0*l4_x*l5_y; 196 out[58][0]=0.0; 197 out[58][1]=45.0*l4_x*l2_y-45.0*l4_x*l4_y; 198 out[59][0]=0.0; 199 out[59][1]=63.0*l4_x*l3_y-63.0*l4_x*l5_y; 200 } 201 202 /** 203 * \brief Evaluate Jacobian of all shape functions 204 * 205 * \param in Position 206 * \param out return value 207 */ evaluateJacobian(const typename Traits::DomainType & in,std::vector<typename Traits::JacobianType> & out) const208 inline void evaluateJacobian (const typename Traits::DomainType& in, 209 std::vector<typename Traits::JacobianType>& out) const 210 { 211 out.resize(60); 212 auto const& x = in[0], y = in[1]; 213 214 const auto l1_x = 2*x - 1; 215 const auto l2_x = x*(6*x - 6) + 1; 216 const auto l3_x = x*(x*(20*x - 30) + 12) - 1; 217 const auto l4_x = x*(x*(x*(70*x - 140) + 90) - 20) + 1; 218 const auto l5_x = x*(x*(x*(x*(252*x - 630) + 560) - 210) + 30) - 1; 219 const auto l1_y = 2*y - 1; 220 const auto l2_y = y*(6*y - 6) + 1; 221 const auto l3_y = y*(y*(20*y - 30) + 12) - 1; 222 const auto l4_y = y*(y*(y*(70*y - 140) + 90) - 20) + 1; 223 const auto l5_y = y*(y*(y*(y*(252*y - 630) + 560) - 210) + 30) - 1; 224 225 const auto dxl1_x = 2.0; 226 const auto dxl2_x = 12*x - 6; 227 const auto dxl3_x = x*(60*x - 60) + 12; 228 const auto dxl4_x = x*(x*(280*x - 420) + 180) - 20; 229 const auto dxl5_x = x*(x*(x*(1260*x - 2520) + 1680) - 420) + 30; 230 const auto dyl1_y = 2.0; 231 const auto dyl2_y = 12*y - 6; 232 const auto dyl3_y = y*(60*y - 60) + 12; 233 const auto dyl4_y = y*(y*(280*y - 420) + 180) - 20; 234 const auto dyl5_y = y*(y*(y*(1260*y - 2520) + 1680) - 420) + 30; 235 236 // x-component 237 out[0][0][0]=sign0*(0.5*(-dxl4_x)+0.5*dxl5_x); 238 out[0][1][0]=0.0; 239 out[1][0][0]=-(1.5)*dxl4_x*l1_y+1.5*dxl5_x*l1_y; 240 out[1][1][0]=0.0; 241 out[2][0][0]=sign0*(-(2.5)*dxl4_x*l2_y+2.5*dxl5_x*l2_y); 242 out[2][1][0]=0.0; 243 out[3][0][0]=-(3.5)*dxl4_x*l3_y+3.5*dxl5_x*l3_y; 244 out[3][1][0]=0.0; 245 out[4][0][0]=sign0*(-(4.5)*dxl4_x*l4_y+4.5*dxl5_x*l4_y); 246 out[4][1][0]=0.0; 247 248 out[5][0][0]=sign1*(0.5*dxl4_x+0.5*dxl5_x); 249 out[5][1][0]=0.0; 250 out[6][0][0]=-(1.5)*dxl4_x*l1_y-1.5*dxl5_x*l1_y; 251 out[6][1][0]=0.0; 252 out[7][0][0]=sign1*(2.5*dxl4_x*l2_y+2.5*dxl5_x*l2_y); 253 out[7][1][0]=0.0; 254 out[8][0][0]=-(3.5)*dxl4_x*l3_y-3.5*dxl5_x*l3_y; 255 out[8][1][0]=0.0; 256 out[9][0][0]=sign1*(4.5*dxl4_x*l4_y+4.5*dxl5_x*l4_y); 257 out[9][1][0]=0.0; 258 259 out[10][0][0]=0.0; 260 out[10][1][0]=0.0; 261 out[11][0][0]=0.0; 262 out[11][1][0]=1.5*dxl1_x*l4_y-1.5*dxl1_x*l5_y; 263 out[12][0][0]=0.0; 264 out[12][1][0]=sign2*(-(2.5)*dxl2_x*l4_y+2.5*dxl2_x*l5_y); 265 out[13][0][0]=0.0; 266 out[13][1][0]=3.5*dxl3_x*l4_y-3.5*dxl3_x*l5_y; 267 out[14][0][0]=0.0; 268 out[14][1][0]=sign2*(-(4.5)*dxl4_x*l4_y+4.5*dxl4_x*l5_y); 269 270 out[15][0][0]=0.0; 271 out[15][1][0]=0.0; 272 out[16][0][0]=0.0; 273 out[16][1][0]=1.5*dxl1_x*l4_y+1.5*dxl1_x*l5_y; 274 out[17][0][0]=0.0; 275 out[17][1][0]=sign3*(2.5*dxl2_x*l4_y+2.5*dxl2_x*l5_y); 276 out[18][0][0]=0.0; 277 out[18][1][0]=3.5*dxl3_x*l4_y+3.5*dxl3_x*l5_y; 278 out[19][0][0]=0.0; 279 out[19][1][0]=sign3*(4.5*dxl4_x*l4_y+4.5*dxl4_x*l5_y); 280 281 out[20][0][0]=-dxl4_x; 282 out[20][1][0]=0.0; 283 out[21][0][0]=-3.0*dxl4_x*l1_y; 284 out[21][1][0]=0.0; 285 out[22][0][0]=-5.0*dxl4_x*l2_y; 286 out[22][1][0]=0.0; 287 out[23][0][0]=-7.0*dxl4_x*l3_y; 288 out[23][1][0]=0.0; 289 out[24][0][0]=-9.0*dxl4_x*l4_y; 290 out[24][1][0]=0.0; 291 out[25][0][0]=3.0*dxl1_x-3.0*dxl5_x; 292 out[25][1][0]=0.0; 293 out[26][0][0]=9.0*dxl1_x*l1_y-9.0*dxl5_x*l1_y; 294 out[26][1][0]=0.0; 295 out[27][0][0]=15.0*dxl1_x*l2_y-15.0*dxl5_x*l2_y; 296 out[27][1][0]=0.0; 297 out[28][0][0]=21.0*dxl1_x*l3_y-21.0*dxl5_x*l3_y; 298 out[28][1][0]=0.0; 299 out[29][0][0]=27.0*dxl1_x*l4_y-27.0*dxl5_x*l4_y; 300 out[29][1][0]=0.0; 301 out[30][0][0]=5.0*dxl2_x-5.0*dxl4_x; 302 out[30][1][0]=0.0; 303 out[31][0][0]=15.0*dxl2_x*l1_y-15.0*dxl4_x*l1_y; 304 out[31][1][0]=0.0; 305 out[32][0][0]=25.0*dxl2_x*l2_y-25.0*dxl4_x*l2_y; 306 out[32][1][0]=0.0; 307 out[33][0][0]=35.0*dxl2_x*l3_y-35.0*dxl4_x*l3_y; 308 out[33][1][0]=0.0; 309 out[34][0][0]=45.0*dxl2_x*l4_y-45.0*dxl4_x*l4_y; 310 out[34][1][0]=0.0; 311 out[35][0][0]=7.0*dxl3_x-7.0*dxl5_x; 312 out[35][1][0]=0.0; 313 out[36][0][0]=21.0*dxl3_x*l1_y-21.0*dxl5_x*l1_y; 314 out[36][1][0]=0.0; 315 out[37][0][0]=35.0*dxl3_x*l2_y-35.0*dxl5_x*l2_y; 316 out[37][1][0]=0.0; 317 out[38][0][0]=49.0*dxl3_x*l3_y-49.0*dxl5_x*l3_y; 318 out[38][1][0]=0.0; 319 out[39][0][0]=63.0*dxl3_x*l4_y-63.0*dxl5_x*l4_y; 320 out[39][1][0]=0.0; 321 out[40][0][0]=0.0; 322 out[40][1][0]=0.0; 323 out[41][0][0]=0.0; 324 out[41][1][0]=0.0; 325 out[42][0][0]=0.0; 326 out[42][1][0]=0.0; 327 out[43][0][0]=0.0; 328 out[43][1][0]=0.0; 329 out[44][0][0]=0.0; 330 out[44][1][0]=3.0*dxl1_x-3.0*dxl1_x*l4_y; 331 out[45][0][0]=0.0; 332 out[45][1][0]=9.0*dxl1_x*l1_y-9.0*dxl1_x*l5_y; 333 out[46][0][0]=0.0; 334 out[46][1][0]=15.0*dxl1_x*l2_y-15.0*dxl1_x*l4_y; 335 out[47][0][0]=0.0; 336 out[47][1][0]=21.0*dxl1_x*l3_y-21.0*dxl1_x*l5_y; 337 out[48][0][0]=0.0; 338 out[48][1][0]=5.0*dxl2_x-5.0*dxl2_x*l4_y; 339 out[49][0][0]=0.0; 340 out[49][1][0]=15.0*dxl2_x*l1_y-15.0*dxl2_x*l5_y; 341 out[50][0][0]=0.0; 342 out[50][1][0]=25.0*dxl2_x*l2_y-25.0*dxl2_x*l4_y; 343 out[51][0][0]=0.0; 344 out[51][1][0]=35.0*dxl2_x*l3_y-35.0*dxl2_x*l5_y; 345 out[52][0][0]=0.0; 346 out[52][1][0]=7.0*dxl3_x-7.0*dxl3_x*l4_y; 347 out[53][0][0]=0.0; 348 out[53][1][0]=21.0*dxl3_x*l1_y-21.0*dxl3_x*l5_y; 349 out[54][0][0]=0.0; 350 out[54][1][0]=35.0*dxl3_x*l2_y-35.0*dxl3_x*l4_y; 351 out[55][0][0]=0.0; 352 out[55][1][0]=49.0*dxl3_x*l3_y-49.0*dxl3_x*l5_y; 353 out[56][0][0]=0.0; 354 out[56][1][0]=9.0*dxl4_x-9.0*dxl4_x*l4_y; 355 out[57][0][0]=0.0; 356 out[57][1][0]=27.0*dxl4_x*l1_y-27.0*dxl4_x*l5_y; 357 out[58][0][0]=0.0; 358 out[58][1][0]=45.0*dxl4_x*l2_y-45.0*dxl4_x*l4_y; 359 out[59][0][0]=0.0; 360 out[59][1][0]=63.0*dxl4_x*l3_y-63.0*dxl4_x*l5_y; 361 362 // y-component 363 out[0][0][1]=0.0; 364 out[0][1][1]=0.0; 365 out[1][0][1]=-(1.5)*l4_x*dyl1_y+1.5*l5_x*dyl1_y; 366 out[1][1][1]=0.0; 367 out[2][0][1]=sign0*(-(2.5)*l4_x*dyl2_y+2.5*l5_x*dyl2_y); 368 out[2][1][1]=0.0; 369 out[3][0][1]=-(3.5)*l4_x*dyl3_y+3.5*l5_x*dyl3_y; 370 out[3][1][1]=0.0; 371 out[4][0][1]=sign0*(-(4.5)*l4_x*dyl4_y+4.5*l5_x*dyl4_y); 372 out[4][1][1]=0.0; 373 374 out[5][0][1]=0.0; 375 out[5][1][1]=0.0; 376 out[6][0][1]=-(1.5)*l4_x*dyl1_y-1.5*l5_x*dyl1_y; 377 out[6][1][1]=0.0; 378 out[7][0][1]=sign1*(2.5*l4_x*dyl2_y+2.5*l5_x*dyl2_y); 379 out[7][1][1]=0.0; 380 out[8][0][1]=-(3.5)*l4_x*dyl3_y-3.5*l5_x*dyl3_y; 381 out[8][1][1]=0.0; 382 out[9][0][1]=sign1*(4.5*l4_x*dyl4_y+4.5*l5_x*dyl4_y); 383 out[9][1][1]=0.0; 384 385 out[10][0][1]=0.0; 386 out[10][1][1]=sign2*(0.5*(-dyl4_y)+0.5*dyl5_y); 387 out[11][0][1]=0.0; 388 out[11][1][1]=1.5*l1_x*dyl4_y-1.5*l1_x*dyl5_y; 389 out[12][0][1]=0.0; 390 out[12][1][1]=sign2*(-(2.5)*l2_x*dyl4_y+2.5*l2_x*dyl5_y); 391 out[13][0][1]=0.0; 392 out[13][1][1]=3.5*l3_x*dyl4_y-3.5*l3_x*dyl5_y; 393 out[14][0][1]=0.0; 394 out[14][1][1]=sign2*(-(4.5)*l4_x*dyl4_y+4.5*l4_x*dyl5_y); 395 396 out[15][0][1]=0.0; 397 out[15][1][1]=sign3*(0.5*dyl4_y+0.5*dyl5_y); 398 out[16][0][1]=0.0; 399 out[16][1][1]=1.5*l1_x*dyl4_y+1.5*l1_x*dyl5_y; 400 out[17][0][1]=0.0; 401 out[17][1][1]=sign3*(2.5*l2_x*dyl4_y+2.5*l2_x*dyl5_y); 402 out[18][0][1]=0.0; 403 out[18][1][1]=3.5*l3_x*dyl4_y+3.5*l3_x*dyl5_y; 404 out[19][0][1]=0.0; 405 out[19][1][1]=sign3*(4.5*l4_x*dyl4_y+4.5*l4_x*dyl5_y); 406 407 out[20][0][1]=0.0; 408 out[20][1][1]=0.0; 409 out[21][0][1]=3.0*dyl1_y-3.0*l4_x*dyl1_y; 410 out[21][1][1]=0.0; 411 out[22][0][1]=5.0*dyl2_y-5.0*l4_x*dyl2_y; 412 out[22][1][1]=0.0; 413 out[23][0][1]=7.0*dyl3_y-7.0*l4_x*dyl3_y; 414 out[23][1][1]=0.0; 415 out[24][0][1]=9.0*dyl4_y-9.0*l4_x*dyl4_y; 416 out[24][1][1]=0.0; 417 out[25][0][1]=0.0; 418 out[25][1][1]=0.0; 419 out[26][0][1]=9.0*l1_x*dyl1_y-9.0*l5_x*dyl1_y; 420 out[26][1][1]=0.0; 421 out[27][0][1]=15.0*l1_x*dyl2_y-15.0*l5_x*dyl2_y; 422 out[27][1][1]=0.0; 423 out[28][0][1]=21.0*l1_x*dyl3_y-21.0*l5_x*dyl3_y; 424 out[28][1][1]=0.0; 425 out[29][0][1]=27.0*l1_x*dyl4_y-27.0*l5_x*dyl4_y; 426 out[29][1][1]=0.0; 427 out[30][0][1]=0.0; 428 out[30][1][1]=0.0; 429 out[31][0][1]=15.0*l2_x*dyl1_y-15.0*l4_x*dyl1_y; 430 out[31][1][1]=0.0; 431 out[32][0][1]=25.0*l2_x*dyl2_y-25.0*l4_x*dyl2_y; 432 out[32][1][1]=0.0; 433 out[33][0][1]=35.0*l2_x*dyl3_y-35.0*l4_x*dyl3_y; 434 out[33][1][1]=0.0; 435 out[34][0][1]=45.0*l2_x*dyl4_y-45.0*l4_x*dyl4_y; 436 out[34][1][1]=0.0; 437 out[35][0][1]=0.0; 438 out[35][1][1]=0.0; 439 out[36][0][1]=21.0*l3_x*dyl1_y-21.0*l5_x*dyl1_y; 440 out[36][1][1]=0.0; 441 out[37][0][1]=35.0*l3_x*dyl2_y-35.0*l5_x*dyl2_y; 442 out[37][1][1]=0.0; 443 out[38][0][1]=49.0*l3_x*dyl3_y-49.0*l5_x*dyl3_y; 444 out[38][1][1]=0.0; 445 out[39][0][1]=63.0*l3_x*dyl4_y-63.0*l5_x*dyl4_y; 446 out[39][1][1]=0.0; 447 out[40][0][1]=0.0; 448 out[40][1][1]=-dyl4_y; 449 out[41][0][1]=0.0; 450 out[41][1][1]=3.0*dyl1_y-3.0*dyl5_y; 451 out[42][0][1]=0.0; 452 out[42][1][1]=5.0*dyl2_y-5.0*dyl4_y; 453 out[43][0][1]=0.0; 454 out[43][1][1]=7.0*dyl3_y-7.0*dyl5_y; 455 out[44][0][1]=0.0; 456 out[44][1][1]=-3.0*l1_x*dyl4_y; 457 out[45][0][1]=0.0; 458 out[45][1][1]=9.0*l1_x*dyl1_y-9.0*l1_x*dyl5_y; 459 out[46][0][1]=0.0; 460 out[46][1][1]=15.0*l1_x*dyl2_y-15.0*l1_x*dyl4_y; 461 out[47][0][1]=0.0; 462 out[47][1][1]=21.0*l1_x*dyl3_y-21.0*l1_x*dyl5_y; 463 out[48][0][1]=0.0; 464 out[48][1][1]=-5.0*l2_x*dyl4_y; 465 out[49][0][1]=0.0; 466 out[49][1][1]=15.0*l2_x*dyl1_y-15.0*l2_x*dyl5_y; 467 out[50][0][1]=0.0; 468 out[50][1][1]=25.0*l2_x*dyl2_y-25.0*l2_x*dyl4_y; 469 out[51][0][1]=0.0; 470 out[51][1][1]=35.0*l2_x*dyl3_y-35.0*l2_x*dyl5_y; 471 out[52][0][1]=0.0; 472 out[52][1][1]=-7.0*l3_x*dyl4_y; 473 out[53][0][1]=0.0; 474 out[53][1][1]=21.0*l3_x*dyl1_y-21.0*l3_x*dyl5_y; 475 out[54][0][1]=0.0; 476 out[54][1][1]=35.0*l3_x*dyl2_y-35.0*l3_x*dyl4_y; 477 out[55][0][1]=0.0; 478 out[55][1][1]=49.0*l3_x*dyl3_y-49.0*l3_x*dyl5_y; 479 out[56][0][1]=0.0; 480 out[56][1][1]=-9.0*l4_x*dyl4_y; 481 out[57][0][1]=0.0; 482 out[57][1][1]=27.0*l4_x*dyl1_y-27.0*l4_x*dyl5_y; 483 out[58][0][1]=0.0; 484 out[58][1][1]=45.0*l4_x*dyl2_y-45.0*l4_x*dyl4_y; 485 out[59][0][1]=0.0; 486 out[59][1][1]=63.0*l4_x*dyl3_y-63.0*l4_x*dyl5_y; 487 } 488 489 //! \brief Evaluate partial derivatives of all shape functions partial(const std::array<unsigned int,2> & order,const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const490 void partial (const std::array<unsigned int, 2>& order, 491 const typename Traits::DomainType& in, // position 492 std::vector<typename Traits::RangeType>& out) const // return value 493 { 494 auto totalOrder = std::accumulate(order.begin(), order.end(), 0); 495 if (totalOrder == 0) { 496 evaluateFunction(in, out); 497 } else if (totalOrder == 1) { 498 out.resize(size()); 499 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1)); 500 auto const& x = in[0], y = in[1]; 501 502 auto l1_x = 2*x - 1; 503 auto l2_x = x*(6*x - 6) + 1; 504 auto l3_x = x*(x*(20*x - 30) + 12) - 1; 505 auto l4_x = x*(x*(x*(70*x - 140) + 90) - 20) + 1; 506 auto l5_x = x*(x*(x*(x*(252*x - 630) + 560) - 210) + 30) - 1; 507 auto l1_y = 2*y - 1; 508 auto l2_y = y*(6*y - 6) + 1; 509 auto l3_y = y*(y*(20*y - 30) + 12) - 1; 510 auto l4_y = y*(y*(y*(70*y - 140) + 90) - 20) + 1; 511 auto l5_y = y*(y*(y*(y*(252*y - 630) + 560) - 210) + 30) - 1; 512 513 if (direction == 0) { 514 auto dxl1_x = 2.0; 515 auto dxl2_x = 12*x - 6; 516 auto dxl3_x = x*(60*x - 60) + 12; 517 auto dxl4_x = x*(x*(280*x - 420) + 180) - 20; 518 auto dxl5_x = x*(x*(x*(1260*x - 2520) + 1680) - 420) + 30; 519 520 out[0][0]=sign0*(0.5*(-dxl4_x)+0.5*dxl5_x); 521 out[0][1]=0.0; 522 out[1][0]=-(1.5)*dxl4_x*l1_y+1.5*dxl5_x*l1_y; 523 out[1][1]=0.0; 524 out[2][0]=sign0*(-(2.5)*dxl4_x*l2_y+2.5*dxl5_x*l2_y); 525 out[2][1]=0.0; 526 out[3][0]=-(3.5)*dxl4_x*l3_y+3.5*dxl5_x*l3_y; 527 out[3][1]=0.0; 528 out[4][0]=sign0*(-(4.5)*dxl4_x*l4_y+4.5*dxl5_x*l4_y); 529 out[4][1]=0.0; 530 531 out[5][0]=sign1*(0.5*dxl4_x+0.5*dxl5_x); 532 out[5][1]=0.0; 533 out[6][0]=-(1.5)*dxl4_x*l1_y-1.5*dxl5_x*l1_y; 534 out[6][1]=0.0; 535 out[7][0]=sign1*(2.5*dxl4_x*l2_y+2.5*dxl5_x*l2_y); 536 out[7][1]=0.0; 537 out[8][0]=-(3.5)*dxl4_x*l3_y-3.5*dxl5_x*l3_y; 538 out[8][1]=0.0; 539 out[9][0]=sign1*(4.5*dxl4_x*l4_y+4.5*dxl5_x*l4_y); 540 out[9][1]=0.0; 541 542 out[10][0]=0.0; 543 out[10][1]=0.0; 544 out[11][0]=0.0; 545 out[11][1]=1.5*dxl1_x*l4_y-1.5*dxl1_x*l5_y; 546 out[12][0]=0.0; 547 out[12][1]=sign2*(-(2.5)*dxl2_x*l4_y+2.5*dxl2_x*l5_y); 548 out[13][0]=0.0; 549 out[13][1]=3.5*dxl3_x*l4_y-3.5*dxl3_x*l5_y; 550 out[14][0]=0.0; 551 out[14][1]=sign2*(-(4.5)*dxl4_x*l4_y+4.5*dxl4_x*l5_y); 552 553 out[15][0]=0.0; 554 out[15][1]=0.0; 555 out[16][0]=0.0; 556 out[16][1]=1.5*dxl1_x*l4_y+1.5*dxl1_x*l5_y; 557 out[17][0]=0.0; 558 out[17][1]=sign3*(2.5*dxl2_x*l4_y+2.5*dxl2_x*l5_y); 559 out[18][0]=0.0; 560 out[18][1]=3.5*dxl3_x*l4_y+3.5*dxl3_x*l5_y; 561 out[19][0]=0.0; 562 out[19][1]=sign3*(4.5*dxl4_x*l4_y+4.5*dxl4_x*l5_y); 563 564 out[20][0]=-dxl4_x; 565 out[20][1]=0.0; 566 out[21][0]=-3.0*dxl4_x*l1_y; 567 out[21][1]=0.0; 568 out[22][0]=-5.0*dxl4_x*l2_y; 569 out[22][1]=0.0; 570 out[23][0]=-7.0*dxl4_x*l3_y; 571 out[23][1]=0.0; 572 out[24][0]=-9.0*dxl4_x*l4_y; 573 out[24][1]=0.0; 574 out[25][0]=3.0*dxl1_x-3.0*dxl5_x; 575 out[25][1]=0.0; 576 out[26][0]=9.0*dxl1_x*l1_y-9.0*dxl5_x*l1_y; 577 out[26][1]=0.0; 578 out[27][0]=15.0*dxl1_x*l2_y-15.0*dxl5_x*l2_y; 579 out[27][1]=0.0; 580 out[28][0]=21.0*dxl1_x*l3_y-21.0*dxl5_x*l3_y; 581 out[28][1]=0.0; 582 out[29][0]=27.0*dxl1_x*l4_y-27.0*dxl5_x*l4_y; 583 out[29][1]=0.0; 584 out[30][0]=5.0*dxl2_x-5.0*dxl4_x; 585 out[30][1]=0.0; 586 out[31][0]=15.0*dxl2_x*l1_y-15.0*dxl4_x*l1_y; 587 out[31][1]=0.0; 588 out[32][0]=25.0*dxl2_x*l2_y-25.0*dxl4_x*l2_y; 589 out[32][1]=0.0; 590 out[33][0]=35.0*dxl2_x*l3_y-35.0*dxl4_x*l3_y; 591 out[33][1]=0.0; 592 out[34][0]=45.0*dxl2_x*l4_y-45.0*dxl4_x*l4_y; 593 out[34][1]=0.0; 594 out[35][0]=7.0*dxl3_x-7.0*dxl5_x; 595 out[35][1]=0.0; 596 out[36][0]=21.0*dxl3_x*l1_y-21.0*dxl5_x*l1_y; 597 out[36][1]=0.0; 598 out[37][0]=35.0*dxl3_x*l2_y-35.0*dxl5_x*l2_y; 599 out[37][1]=0.0; 600 out[38][0]=49.0*dxl3_x*l3_y-49.0*dxl5_x*l3_y; 601 out[38][1]=0.0; 602 out[39][0]=63.0*dxl3_x*l4_y-63.0*dxl5_x*l4_y; 603 out[39][1]=0.0; 604 out[40][0]=0.0; 605 out[40][1]=0.0; 606 out[41][0]=0.0; 607 out[41][1]=0.0; 608 out[42][0]=0.0; 609 out[42][1]=0.0; 610 out[43][0]=0.0; 611 out[43][1]=0.0; 612 out[44][0]=0.0; 613 out[44][1]=3.0*dxl1_x-3.0*dxl1_x*l4_y; 614 out[45][0]=0.0; 615 out[45][1]=9.0*dxl1_x*l1_y-9.0*dxl1_x*l5_y; 616 out[46][0]=0.0; 617 out[46][1]=15.0*dxl1_x*l2_y-15.0*dxl1_x*l4_y; 618 out[47][0]=0.0; 619 out[47][1]=21.0*dxl1_x*l3_y-21.0*dxl1_x*l5_y; 620 out[48][0]=0.0; 621 out[48][1]=5.0*dxl2_x-5.0*dxl2_x*l4_y; 622 out[49][0]=0.0; 623 out[49][1]=15.0*dxl2_x*l1_y-15.0*dxl2_x*l5_y; 624 out[50][0]=0.0; 625 out[50][1]=25.0*dxl2_x*l2_y-25.0*dxl2_x*l4_y; 626 out[51][0]=0.0; 627 out[51][1]=35.0*dxl2_x*l3_y-35.0*dxl2_x*l5_y; 628 out[52][0]=0.0; 629 out[52][1]=7.0*dxl3_x-7.0*dxl3_x*l4_y; 630 out[53][0]=0.0; 631 out[53][1]=21.0*dxl3_x*l1_y-21.0*dxl3_x*l5_y; 632 out[54][0]=0.0; 633 out[54][1]=35.0*dxl3_x*l2_y-35.0*dxl3_x*l4_y; 634 out[55][0]=0.0; 635 out[55][1]=49.0*dxl3_x*l3_y-49.0*dxl3_x*l5_y; 636 out[56][0]=0.0; 637 out[56][1]=9.0*dxl4_x-9.0*dxl4_x*l4_y; 638 out[57][0]=0.0; 639 out[57][1]=27.0*dxl4_x*l1_y-27.0*dxl4_x*l5_y; 640 out[58][0]=0.0; 641 out[58][1]=45.0*dxl4_x*l2_y-45.0*dxl4_x*l4_y; 642 out[59][0]=0.0; 643 out[59][1]=63.0*dxl4_x*l3_y-63.0*dxl4_x*l5_y; 644 645 } else if (direction == 1) { 646 auto dyl1_y = 2.0; 647 auto dyl2_y = 12*y - 6; 648 auto dyl3_y = y*(60*y - 60) + 12; 649 auto dyl4_y = y*(y*(280*y - 420) + 180) - 20; 650 auto dyl5_y = y*(y*(y*(1260*y - 2520) + 1680) - 420) + 30; 651 652 out[0][0]=0.0; 653 out[0][1]=0.0; 654 out[1][0]=-(1.5)*l4_x*dyl1_y+1.5*l5_x*dyl1_y; 655 out[1][1]=0.0; 656 out[2][0]=sign0*(-(2.5)*l4_x*dyl2_y+2.5*l5_x*dyl2_y); 657 out[2][1]=0.0; 658 out[3][0]=-(3.5)*l4_x*dyl3_y+3.5*l5_x*dyl3_y; 659 out[3][1]=0.0; 660 out[4][0]=sign0*(-(4.5)*l4_x*dyl4_y+4.5*l5_x*dyl4_y); 661 out[4][1]=0.0; 662 663 out[5][0]=0.0; 664 out[5][1]=0.0; 665 out[6][0]=-(1.5)*l4_x*dyl1_y-1.5*l5_x*dyl1_y; 666 out[6][1]=0.0; 667 out[7][0]=sign1*(2.5*l4_x*dyl2_y+2.5*l5_x*dyl2_y); 668 out[7][1]=0.0; 669 out[8][0]=-(3.5)*l4_x*dyl3_y-3.5*l5_x*dyl3_y; 670 out[8][1]=0.0; 671 out[9][0]=sign1*(4.5*l4_x*dyl4_y+4.5*l5_x*dyl4_y); 672 out[9][1]=0.0; 673 674 out[10][0]=0.0; 675 out[10][1]=sign2*(0.5*(-dyl4_y)+0.5*dyl5_y); 676 out[11][0]=0.0; 677 out[11][1]=1.5*l1_x*dyl4_y-1.5*l1_x*dyl5_y; 678 out[12][0]=0.0; 679 out[12][1]=sign2*(-(2.5)*l2_x*dyl4_y+2.5*l2_x*dyl5_y); 680 out[13][0]=0.0; 681 out[13][1]=3.5*l3_x*dyl4_y-3.5*l3_x*dyl5_y; 682 out[14][0]=0.0; 683 out[14][1]=sign2*(-(4.5)*l4_x*dyl4_y+4.5*l4_x*dyl5_y); 684 685 out[15][0]=0.0; 686 out[15][1]=sign3*(0.5*dyl4_y+0.5*dyl5_y); 687 out[16][0]=0.0; 688 out[16][1]=1.5*l1_x*dyl4_y+1.5*l1_x*dyl5_y; 689 out[17][0]=0.0; 690 out[17][1]=sign3*(2.5*l2_x*dyl4_y+2.5*l2_x*dyl5_y); 691 out[18][0]=0.0; 692 out[18][1]=3.5*l3_x*dyl4_y+3.5*l3_x*dyl5_y; 693 out[19][0]=0.0; 694 out[19][1]=sign3*(4.5*l4_x*dyl4_y+4.5*l4_x*dyl5_y); 695 696 out[20][0]=0.0; 697 out[20][1]=0.0; 698 out[21][0]=3.0*dyl1_y-3.0*l4_x*dyl1_y; 699 out[21][1]=0.0; 700 out[22][0]=5.0*dyl2_y-5.0*l4_x*dyl2_y; 701 out[22][1]=0.0; 702 out[23][0]=7.0*dyl3_y-7.0*l4_x*dyl3_y; 703 out[23][1]=0.0; 704 out[24][0]=9.0*dyl4_y-9.0*l4_x*dyl4_y; 705 out[24][1]=0.0; 706 out[25][0]=0.0; 707 out[25][1]=0.0; 708 out[26][0]=9.0*l1_x*dyl1_y-9.0*l5_x*dyl1_y; 709 out[26][1]=0.0; 710 out[27][0]=15.0*l1_x*dyl2_y-15.0*l5_x*dyl2_y; 711 out[27][1]=0.0; 712 out[28][0]=21.0*l1_x*dyl3_y-21.0*l5_x*dyl3_y; 713 out[28][1]=0.0; 714 out[29][0]=27.0*l1_x*dyl4_y-27.0*l5_x*dyl4_y; 715 out[29][1]=0.0; 716 out[30][0]=0.0; 717 out[30][1]=0.0; 718 out[31][0]=15.0*l2_x*dyl1_y-15.0*l4_x*dyl1_y; 719 out[31][1]=0.0; 720 out[32][0]=25.0*l2_x*dyl2_y-25.0*l4_x*dyl2_y; 721 out[32][1]=0.0; 722 out[33][0]=35.0*l2_x*dyl3_y-35.0*l4_x*dyl3_y; 723 out[33][1]=0.0; 724 out[34][0]=45.0*l2_x*dyl4_y-45.0*l4_x*dyl4_y; 725 out[34][1]=0.0; 726 out[35][0]=0.0; 727 out[35][1]=0.0; 728 out[36][0]=21.0*l3_x*dyl1_y-21.0*l5_x*dyl1_y; 729 out[36][1]=0.0; 730 out[37][0]=35.0*l3_x*dyl2_y-35.0*l5_x*dyl2_y; 731 out[37][1]=0.0; 732 out[38][0]=49.0*l3_x*dyl3_y-49.0*l5_x*dyl3_y; 733 out[38][1]=0.0; 734 out[39][0]=63.0*l3_x*dyl4_y-63.0*l5_x*dyl4_y; 735 out[39][1]=0.0; 736 out[40][0]=0.0; 737 out[40][1]=-dyl4_y; 738 out[41][0]=0.0; 739 out[41][1]=3.0*dyl1_y-3.0*dyl5_y; 740 out[42][0]=0.0; 741 out[42][1]=5.0*dyl2_y-5.0*dyl4_y; 742 out[43][0]=0.0; 743 out[43][1]=7.0*dyl3_y-7.0*dyl5_y; 744 out[44][0]=0.0; 745 out[44][1]=-3.0*l1_x*dyl4_y; 746 out[45][0]=0.0; 747 out[45][1]=9.0*l1_x*dyl1_y-9.0*l1_x*dyl5_y; 748 out[46][0]=0.0; 749 out[46][1]=15.0*l1_x*dyl2_y-15.0*l1_x*dyl4_y; 750 out[47][0]=0.0; 751 out[47][1]=21.0*l1_x*dyl3_y-21.0*l1_x*dyl5_y; 752 out[48][0]=0.0; 753 out[48][1]=-5.0*l2_x*dyl4_y; 754 out[49][0]=0.0; 755 out[49][1]=15.0*l2_x*dyl1_y-15.0*l2_x*dyl5_y; 756 out[50][0]=0.0; 757 out[50][1]=25.0*l2_x*dyl2_y-25.0*l2_x*dyl4_y; 758 out[51][0]=0.0; 759 out[51][1]=35.0*l2_x*dyl3_y-35.0*l2_x*dyl5_y; 760 out[52][0]=0.0; 761 out[52][1]=-7.0*l3_x*dyl4_y; 762 out[53][0]=0.0; 763 out[53][1]=21.0*l3_x*dyl1_y-21.0*l3_x*dyl5_y; 764 out[54][0]=0.0; 765 out[54][1]=35.0*l3_x*dyl2_y-35.0*l3_x*dyl4_y; 766 out[55][0]=0.0; 767 out[55][1]=49.0*l3_x*dyl3_y-49.0*l3_x*dyl5_y; 768 out[56][0]=0.0; 769 out[56][1]=-9.0*l4_x*dyl4_y; 770 out[57][0]=0.0; 771 out[57][1]=27.0*l4_x*dyl1_y-27.0*l4_x*dyl5_y; 772 out[58][0]=0.0; 773 out[58][1]=45.0*l4_x*dyl2_y-45.0*l4_x*dyl4_y; 774 out[59][0]=0.0; 775 out[59][1]=63.0*l4_x*dyl3_y-63.0*l4_x*dyl5_y; 776 } else { 777 DUNE_THROW(RangeError, "Component out of range."); 778 } 779 } else { 780 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented"); 781 } 782 } 783 784 //! \brief Polynomial order of the shape functions order() const785 unsigned int order () const 786 { 787 return 9; 788 } 789 790 private: 791 R sign0, sign1, sign2, sign3; 792 }; 793 } 794 795 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS3_CUBE2D_LOCALBASIS_HH 796