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_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALBASIS_HH 4 #define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALBASIS_HH 5 6 #include <algorithm> 7 #include <array> 8 #include <bitset> 9 #include <numeric> 10 #include <vector> 11 12 #include <dune/common/fmatrix.hh> 13 #include <dune/common/fvector.hh> 14 #include <dune/common/math.hh> 15 #include <dune/common/rangeutilities.hh> 16 #include <dune/common/typetraits.hh> 17 18 #include <dune/localfunctions/common/localbasis.hh> 19 20 namespace Dune 21 { 22 /** 23 * \ingroup LocalBasisImplementation 24 * \brief Brezzi-Douglas-Fortin-Marini shape functions on a reference cube. 25 * 26 * \tparam D Type to represent the field in the domain. 27 * \tparam R Type to represent the field in the range. 28 * \tparam dim dimension of the reference element, must be >= 2. 29 * \tparam order order of the element, must be >= 1. 30 * 31 * \nosubgroup 32 */ 33 template<class D, class R, unsigned int dim, unsigned int order> 34 class BDFMCubeLocalBasis 35 { 36 static_assert( AlwaysFalse<D>::value, 37 "`BDFMCubeLocalBasis` not implemented for chosen `dim` and `order`." ); 38 }; 39 40 41 #ifndef DOXYGEN 42 template<class D, class R, unsigned int dim> 43 class BDFMCubeLocalBasis<D, R, dim, 0> 44 { 45 static_assert( AlwaysFalse<D>::value, 46 "`BDFMCubeLocalBasis` not defined for order 0." ); 47 }; 48 #endif // #ifndef DOXYGEN 49 50 51 /** 52 * \brief First order Brezzi-Douglas-Fortin-Marini shape functions on the refrence quadrialteral. 53 * 54 * \nosubgrouping 55 */ 56 template<class D, class R > 57 class BDFMCubeLocalBasis<D, R, 2, 1> 58 { 59 using DomainType = FieldVector<D, 2>; 60 using RangeType = FieldVector<R, 2>; 61 using JacobianType = FieldMatrix<R, 2, 2>; 62 63 public: 64 using Traits = LocalBasisTraits<D, 2, DomainType, R, 2, RangeType, JacobianType>; 65 66 //! \brief Standard constructor BDFMCubeLocalBasis()67 BDFMCubeLocalBasis () 68 { 69 std::fill(s_.begin(), s_.end(), 1); 70 } 71 72 /** 73 * \brief Make set number s, where 0<= s < 16 74 * 75 * \param s Edge orientation indicator 76 */ BDFMCubeLocalBasis(std::bitset<4> s)77 BDFMCubeLocalBasis (std::bitset<4> s) 78 { 79 for (auto i : range(4)) 80 s_[i] = s[i] ? -1 : 1; 81 } 82 83 //! \brief number of shape functions size() const84 unsigned int size () const { return 4; } 85 86 /** 87 * \brief Evaluate all shape functions 88 * 89 * \param in Position 90 * \param out return value 91 */ evaluateFunction(const DomainType & in,std::vector<RangeType> & out) const92 inline void evaluateFunction (const DomainType& in, std::vector<RangeType>& out) const 93 { 94 out.resize(4); 95 96 out[0] = {s_[0]*(in[0]-1), 0 }; 97 out[1] = {s_[1]*(in[0]) , 0 }; 98 out[2] = {0, s_[2]*(in[1]-1)}; 99 out[3] = {0, s_[3]*(in[1]) }; 100 } 101 102 /** 103 * \brief Evaluate Jacobian of all shape functions 104 * 105 * \param in Position 106 * \param out return value 107 */ evaluateJacobian(const DomainType & in,std::vector<JacobianType> & out) const108 inline void evaluateJacobian (const DomainType& in, std::vector<JacobianType>& out) const 109 { 110 out.resize(4); 111 112 out[0] = {{s_[0], 0}, {0, 0}}; 113 out[1] = {{s_[1], 0}, {0, 0}}; 114 out[2] = {{0, 0}, {0, s_[2]}}; 115 out[3] = {{0, 0}, {0, s_[3]}}; 116 } 117 118 /** 119 * \brief Evaluate all partial derivatives of all shape functions 120 * 121 * \param order order the partial derivative 122 * \param in Position 123 * \param out return value 124 */ partial(const std::array<unsigned int,2> & order,const DomainType & in,std::vector<RangeType> & out) const125 void partial (const std::array<unsigned int, 2>& order, 126 const DomainType& in, 127 std::vector<RangeType>& out) const 128 { 129 if (std::accumulate(order.begin(), order.end(), 0) == 0) 130 evaluateFunction(in, out); 131 else 132 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented"); 133 } 134 135 //! \brief Polynomial order of the shape functions order() const136 unsigned int order () const { return 1; } 137 138 private: 139 std::array<R, 4> s_; 140 }; 141 142 143 /** 144 * \brief Second order Brezzi-Douglas-Fortin-Marini shape functions on the refrence quadrialteral. 145 * 146 * \nosubgrouping 147 */ 148 template<class D, class R > 149 class BDFMCubeLocalBasis<D, R, 2, 2> 150 { 151 using DomainType = FieldVector<D, 2>; 152 using RangeType = FieldVector<R, 2>; 153 using JacobianType = FieldMatrix<R, 2, 2>; 154 155 public: 156 using Traits = LocalBasisTraits<D, 2, DomainType, R, 2, RangeType, JacobianType>; 157 158 //! \brief Standard constructor BDFMCubeLocalBasis()159 BDFMCubeLocalBasis () 160 { 161 std::fill(s_.begin(), s_.end(), 1); 162 } 163 164 /** 165 * \brief Make set number s, where 0<= s < 16 166 * 167 * \param s Edge orientation indicator 168 */ BDFMCubeLocalBasis(std::bitset<4> s)169 BDFMCubeLocalBasis (std::bitset<4> s) 170 { 171 for (auto i : range(4)) 172 s_[i] = s[i] ? -1 : 1; 173 } 174 175 //! \brief number of shape functions size() const176 unsigned int size () const { return 10; } 177 178 /** 179 * \brief Evaluate all shape functions 180 * 181 * \param in Position 182 * \param out return value 183 */ evaluateFunction(const DomainType & in,std::vector<RangeType> & out) const184 inline void evaluateFunction (const DomainType& in, std::vector<RangeType>& out) const 185 { 186 out.resize(10); 187 188 const auto& x = in[0]; 189 const auto& y = in[1]; 190 191 auto pre = 1-x; 192 out[0] = {pre*s_[0]*(3*x-1), 0}; 193 out[1] = { pre*3*(2*y-1), 0}; 194 195 pre = x; 196 out[2] = {pre*s_[1]*(3*x-2), 0}; 197 out[3] = {pre* 3*(2*y-1), 0}; 198 199 pre = 1-y; 200 out[4] = {0, pre*s_[2]*(3*y-1)}; 201 out[5] = {0, pre*3*(2*x-1)}; 202 203 pre = y; 204 out[6] = {0, pre*s_[3]*(3*y-2)}; 205 out[7] = {0, pre*3*(2*x-1)}; 206 207 out[8] = {6*x*(1-x), 0}; 208 209 out[9] = {0, 6*y*(1-y)}; 210 } 211 212 /** 213 * \brief Evaluate Jacobian of all shape functions 214 * 215 * \param in Position 216 * \param out return value 217 */ evaluateJacobian(const DomainType & in,std::vector<JacobianType> & out) const218 inline void evaluateJacobian (const DomainType& in, std::vector<JacobianType>& out) const 219 { 220 out.resize(10); 221 222 const auto& x = in[0]; 223 const auto& y = in[1]; 224 225 out[0] = {{s_[0]*(4-6*x), 0}, {0, 0}}; 226 out[1] = {{3-6*y, 6-6*x}, {0, 0}}; 227 228 out[2] = {{s_[1]*(6*x-2), 0}, {0, 0}}; 229 out[3] = {{6*y-3, 6*x}, {0, 0}}; 230 231 out[4] = {{0, 0}, {0, s_[2]*(4-6*y)}}; 232 out[5] = {{0, 0}, {6-6*y, 3-6*x}}; 233 234 out[6] = {{0, 0}, {0, s_[3]*(6*y-2)}}; 235 out[7] = {{0, 0}, {6*y, 6*x-3}}; 236 237 out[8] = {{6-12*x, 0}, {0, 0}}; 238 239 out[9] = {{0, 0}, {0, 6-12*y}}; 240 } 241 242 /** 243 * \brief Evaluate all partial derivatives of all shape functions 244 * 245 * \param order order the partial derivative 246 * \param in Position 247 * \param out return value 248 */ partial(const std::array<unsigned int,2> & order,const DomainType & in,std::vector<RangeType> & out) const249 void partial (const std::array<unsigned int, 2>& order, 250 const DomainType& in, 251 std::vector<RangeType>& out) const 252 { 253 if (std::accumulate(order.begin(), order.end(), 0) == 0) 254 evaluateFunction(in, out); 255 else 256 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented"); 257 } 258 259 //! \brief Polynomial order of the shape functions order() const260 unsigned int order () const { return 2; } 261 262 private: 263 std::array<R, 4> s_; 264 }; 265 266 267 /** 268 * \brief Third order Brezzi-Douglas-Fortin-Marini shape functions on the refrence quadrialteral. 269 * 270 * \nosubgrouping 271 */ 272 template<class D, class R > 273 class BDFMCubeLocalBasis<D, R, 2, 3> 274 { 275 using DomainType = FieldVector<D, 2>; 276 using RangeType = FieldVector<R, 2>; 277 using JacobianType = FieldMatrix<R, 2, 2>; 278 279 public: 280 using Traits = LocalBasisTraits<D, 2, DomainType, R, 2, RangeType, JacobianType>; 281 282 //! \brief Standard constructor BDFMCubeLocalBasis()283 BDFMCubeLocalBasis () 284 { 285 std::fill(s_.begin(), s_.end(), 1); 286 } 287 288 /** 289 * \brief Make set number s, where 0<= s < 16 290 * 291 * \param s Edge orientation indicator 292 */ BDFMCubeLocalBasis(std::bitset<4> s)293 BDFMCubeLocalBasis (std::bitset<4> s) 294 { 295 for (auto i : range(4)) 296 s_[i] = s[i] ? -1 : 1; 297 } 298 299 //! \brief number of shape functions size() const300 unsigned int size () const { return 18; } 301 302 /** 303 * \brief Evaluate all shape functions 304 * 305 * \param in Position 306 * \param out return value 307 */ evaluateFunction(const DomainType & in,std::vector<RangeType> & out) const308 inline void evaluateFunction (const DomainType& in, std::vector<RangeType>& out) const 309 { 310 out.resize(18); 311 312 const auto& x = in[0]; 313 const auto& y = in[1]; 314 315 auto pre = 1-x; 316 out[0] = {pre*s_[0]*-1*(10*x*x-8*x+1), 0}; 317 out[1] = {pre* 3*(1-3*x)*(2*y-1), 0}; 318 out[2] = {pre* s_[0]*-5*(6*y*y-6*y+1), 0}; 319 320 pre = x; 321 out[3] = {pre*s_[1]*(10*x*x-12*x+3), 0}; 322 out[4] = {pre* 3*(3*x-2)*(2*y-1), 0}; 323 out[5] = {pre*s_[1]*5*(6*y*y-6*y+1), 0}; 324 325 pre = 1-y; 326 out[6] = {0, pre*s_[2]*-1*(10*y*y-8*y+1)}; 327 out[7] = {0, pre* 3*(1-3*y)*(2*x-1)}; 328 out[8] = {0, pre*s_[2]*-5*( 6*x*x-6*x+1)}; 329 330 pre = y; 331 out[9] = {0, pre*s_[3]*(10*y*y-12*y+3)}; 332 out[10] = {0, pre* 3*(3*y-2)*(2*x-1)}; 333 out[11] = {0, pre*s_[3]*5*(6*x*x-6*x+1)}; 334 335 pre = x*(1-x); 336 out[12] = {pre* 6 , 0}; 337 out[14] = {pre*30*(2*x-1), 0}; 338 out[16] = {pre*18*(2*y-1), 0}; 339 340 pre = y*(1-y); 341 out[13] = {0, pre* 6 }; 342 out[15] = {0, pre*18*(2*x-1)}; 343 out[17] = {0, pre*30*(2*y-1)}; 344 } 345 346 /** 347 * \brief Evaluate Jacobian of all shape functions 348 * 349 * \param in Position 350 * \param out return value 351 */ evaluateJacobian(const DomainType & in,std::vector<JacobianType> & out) const352 inline void evaluateJacobian (const DomainType& in, std::vector<JacobianType>& out) const 353 { 354 out.resize(18); 355 356 const auto& x = in[0]; 357 const auto& y = in[1]; 358 359 out[0] = {{s_[0]*(30*x*x-36*x+9), 0}, {0, 0}}; 360 out[1] = {{ 6*(6*x*y-3*x-4*y+2), 6*(3*x*x-4*x+1)}, {0, 0}}; 361 out[2] = {{s_[0]*5*(6*y*y-6*y+1), s_[0]*30*(x-1)*(2*y-1)}, {0, 0}}; 362 363 out[3] = {{ s_[1]*30*x*x-24*x+3, 0}, {0, 0}}; 364 out[4] = {{ 6*(3*x-1)*(2*y-1), 6*x*(3*x-2)}, {0, 0}}; 365 out[5] = {{s_[1]*5*(6*y*y-6*y+1), s_[1]*30*x*(2*y-1)}, {0, 0}}; 366 367 out[6] = {{0, 0}, { 0, s_[2]*(30*y*y-36*y+9)}}; 368 out[7] = {{0, 0}, { 6*(3*y*y-4*y+1), 6*(6*y*x-3*y-4*x+2)}}; 369 out[8] = {{0, 0}, {s_[2]*30*(y-1)*(2*x-1), s_[2]*5*(6*x*x-6*x+1)}}; 370 371 out[9] = {{0, 0}, { 0, s_[3]*30*y*y-24*y+3}}; 372 out[10] = {{0, 0}, { 6*y*(3*y-2), 6*(3*y-1)*(2*x-1)}}; 373 out[11] = {{0, 0}, {s_[3]*30*y*(2*x-1), s_[3]*5*(6*x*x-6*x+1)}}; 374 375 out[12] = {{ -6*(2*x-1), 0}, {0, 0}}; 376 out[14] = {{ -30*(6*x*x-6*x+1), 0}, {0, 0}}; 377 out[16] = {{-18*(2*x-1)*(2*y-1), 36*x*(1-x)}, {0, 0}}; 378 379 out[13] = {{0, 0}, { 0, -6*(2*y-1)}}; 380 out[15] = {{0, 0}, {36*y*(1-y), -18*(2*y-1)*(2*x-1)}}; 381 out[17] = {{0, 0}, { 0, -30*(6*y*y-6*y+1)}}; 382 } 383 384 /** 385 * \brief Evaluate all partial derivatives of all shape functions 386 * 387 * \param order order the partial derivative 388 * \param in Position 389 * \param out return value 390 */ partial(const std::array<unsigned int,2> & order,const DomainType & in,std::vector<RangeType> & out) const391 void partial (const std::array<unsigned int, 2>& order, 392 const DomainType& in, 393 std::vector<RangeType>& out) const 394 { 395 if (std::accumulate(order.begin(), order.end(), 0) == 0) 396 evaluateFunction(in, out); 397 else 398 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented"); 399 } 400 401 //! \brief Polynomial order of the shape functions order() const402 unsigned int order () const { return 3; } 403 404 private: 405 std::array<R, 4> s_; 406 }; 407 408 } // namespace Dune 409 410 #endif // #ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALBASIS_HH 411