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_BREZZIDOUGLASMARINI2_SIMPLEX2D_LOCALBASIS_HH 4 #define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI2_SIMPLEX2D_LOCALBASIS_HH 5 6 #include <array> 7 #include <bitset> 8 #include <numeric> 9 #include <vector> 10 11 #include <dune/common/fmatrix.hh> 12 13 #include "../../common/localbasis.hh" 14 15 namespace Dune 16 { 17 /** 18 * \ingroup LocalBasisImplementation 19 * \brief First order Brezzi-Douglas-Marini shape functions on quadrilaterals. 20 * 21 * \tparam D Type to represent the field in the domain. 22 * \tparam R Type to represent the field in the range. 23 * 24 * \nosubgrouping 25 */ 26 template<class D, class R> 27 class BDM2Simplex2DLocalBasis 28 { 29 30 public: 31 typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>, 32 R,2,Dune::FieldVector<R,2>, 33 Dune::FieldMatrix<R,2,2> > Traits; 34 35 //! \brief Standard constructor BDM2Simplex2DLocalBasis()36 BDM2Simplex2DLocalBasis() 37 { 38 for (size_t i=0; i<3; i++) 39 sign_[i] = 1.0; 40 } 41 42 /** 43 * \brief Make set number s, where 0 <= s < 8 44 * 45 * \param s Edge orientation indicator 46 */ BDM2Simplex2DLocalBasis(std::bitset<3> s)47 BDM2Simplex2DLocalBasis(std::bitset<3> s) 48 { 49 for (size_t i=0; i<3; i++) 50 sign_[i] = s[i] ? -1.0 : 1.0; 51 } 52 53 //! \brief number of shape functions size() const54 unsigned int size() const 55 { 56 return 12; 57 } 58 59 /** 60 * \brief Evaluate all shape functions 61 * 62 * \param in Position 63 * \param out return value 64 */ evaluateFunction(const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const65 inline void evaluateFunction(const typename Traits::DomainType& in, 66 std::vector<typename Traits::RangeType>& out) const 67 { 68 out.resize(size()); 69 70 out[0][0] = sign_[0]*(-2*in[0]*in[1] + in[0]*in[0]); 71 out[0][1] = sign_[0]*(-1 + 6*in[1] -2*in[0]*in[1] - 5*in[1]*in[1]); 72 73 out[1][0] = 1.5*in[0] + 3*in[0]*in[1] - 4.5*in[0]*in[0]; 74 out[1][1] = -3 + 6*in[0] + 10.5*in[1] - 15*in[0]*in[1] - 7.5*in[1]*in[1]; 75 76 out[2][0] = sign_[0]*(-7.5*in[0] + 5*in[0]*in[1] + 12.5*in[0]*in[0]); 77 out[2][1] = sign_[0]*(-5 + 30*in[0] + 7.5*in[1] - 25*in[0]*in[1] - 30*in[0]*in[0] - 2.5*in[1]*in[1]); 78 79 80 81 out[3][0] = sign_[1]*(-1 + 6*in[0] - 2*in[0]*in[1] - 5*in[0]*in[0]); 82 out[3][1] = sign_[1]*(-2*in[0]*in[1] + in[1]*in[1]); 83 84 out[4][0] = 3 - 10.5*in[0] - 6*in[1] + 15*in[0]*in[1] + 7.5*in[0]*in[0]; 85 out[4][1] = -1.5*in[1] - 3*in[0]*in[1] + 4.5*in[1]*in[1]; 86 87 out[5][0] = sign_[1]*(-5 + 7.5*in[0] + 30*in[1] - 25*in[0]*in[1] - 2.5*in[0]*in[0] - 30*in[1]*in[1]); 88 out[5][1] = sign_[1]*(-7.5*in[1] + 5*in[0]*in[1] + 12.5*in[1]*in[1]); 89 90 91 92 out[6][0] = sign_[2]*(-3*in[0] + 4*in[0]*in[1] + 4*in[0]*in[0]); 93 out[6][1] = sign_[2]*(-3*in[1] + 4*in[0]*in[1] + 4*in[1]*in[1]); 94 95 out[7][0] = -3*in[0] + 6*in[0]*in[0]; 96 out[7][1] = 3*in[1] - 6*in[1]*in[1]; 97 98 out[8][0] = sign_[2]*(-10*in[0]*in[1] + 5*in[0]*in[0]); 99 out[8][1] = sign_[2]*(-10*in[0]*in[1] + 5*in[1]*in[1]); 100 101 102 103 out[9][0] = 18*in[0] - 12*in[0]*in[1] - 18*in[0]*in[0]; 104 out[9][1] = 6*in[1] - 12*in[0]*in[1] - 6*in[1]*in[1]; 105 106 out[10][0] = 6*in[0] - 12*in[0]*in[1] - 6*in[0]*in[0]; 107 out[10][1] = 18*in[1] - 12*in[0]*in[1] - 18*in[1]*in[1]; 108 109 out[11][0] = 90*in[0] - 180*in[0]*in[1] - 90*in[0]*in[0]; 110 out[11][1] = -90*in[1] + 180*in[0]*in[1] + 90*in[1]*in[1]; 111 } 112 113 /** 114 * \brief Evaluate Jacobian of all shape functions 115 * 116 * \param in Position 117 * \param out return value 118 */ evaluateJacobian(const typename Traits::DomainType & in,std::vector<typename Traits::JacobianType> & out) const119 inline void evaluateJacobian(const typename Traits::DomainType& in, 120 std::vector<typename Traits::JacobianType>& out) const 121 { 122 out.resize(size()); 123 124 out[0][0][0] = sign_[0]*(-2*in[1] + 2*in[0]); 125 out[0][0][1] = sign_[0]*(-2*in[0]); 126 127 out[0][1][0] = sign_[0]*(-2*in[1]); 128 out[0][1][1] = sign_[0]*(6 -2*in[0] - 10*in[1]); 129 130 131 out[1][0][0] = 1.5 + 3*in[1] - 9*in[0]; 132 out[1][0][1] = 3*in[0]; 133 134 out[1][1][0] = 6 - 15*in[1]; 135 out[1][1][1] = 10.5 - 15*in[0] - 15*in[1]; 136 137 138 out[2][0][0] = sign_[0]*(-7.5 + 5*in[1] + 25*in[0]); 139 out[2][0][1] = sign_[0]*(5*in[0]); 140 141 out[2][1][0] = sign_[0]*(30 - 25*in[1] - 60*in[0]); 142 out[2][1][1] = sign_[0]*(7.5 - 25*in[0] - 5*in[1]); 143 144 145 146 out[3][0][0] = sign_[1]*(6 - 2*in[1] - 10*in[0]); 147 out[3][0][1] = sign_[1]*(-2*in[0]); 148 149 out[3][1][0] = sign_[1]*(-2*in[1]); 150 out[3][1][1] = sign_[1]*(-2*in[0] + 2*in[1]); 151 152 153 out[4][0][0] = -10.5 + 15*in[1] + 15*in[0]; 154 out[4][0][1] = -6 + 15*in[0]; 155 156 out[4][1][0] = -3*in[1]; 157 out[4][1][1] = -1.5 - 3*in[0] + 9*in[1]; 158 159 160 out[5][0][0] = sign_[1]*(7.5 - 25*in[1] - 5*in[0]); 161 out[5][0][1] = sign_[1]*(30 - 25*in[0] - 60*in[1]); 162 163 out[5][1][0] = sign_[1]*(5*in[1]); 164 out[5][1][1] = sign_[1]*(-7.5 + 5*in[0] + 25*in[1]); 165 166 167 168 out[6][0][0] = sign_[2]*(-3 + 4*in[1] + 8*in[0]); 169 out[6][0][1] = sign_[2]*(4*in[0]); 170 171 out[6][1][0] = sign_[2]*(4*in[1]); 172 out[6][1][1] = sign_[2]*(-3 + 4*in[0] + 8*in[1]); 173 174 175 out[7][0][0] = -3 + 12*in[0]; 176 out[7][0][1] = 0; 177 178 out[7][1][0] = 0; 179 out[7][1][1] = 3 - 12*in[1]; 180 181 182 out[8][0][0] = sign_[2]*(-10*in[1] + 10*in[0]); 183 out[8][0][1] = sign_[2]*(-10*in[0]); 184 185 out[8][1][0] = sign_[2]*(-10*in[1]); 186 out[8][1][1] = sign_[2]*(-10*in[0] + 10*in[1]); 187 188 189 out[9][0][0] = 18 - 12*in[1] - 36*in[0]; 190 out[9][0][1] = -12*in[0]; 191 192 out[9][1][0] = -12*in[1]; 193 out[9][1][1] = 6 - 12*in[0] - 12*in[1]; 194 195 out[10][0][0] = 6 - 12*in[1] - 12*in[0]; 196 out[10][0][1] = -12*in[0]; 197 198 out[10][1][0] = -12*in[1]; 199 out[10][1][1] = 18 - 12*in[0] - 36*in[1]; 200 201 out[11][0][0] = 90 - 180*in[1] - 180*in[0]; 202 out[11][0][1] = -180*in[0]; 203 204 out[11][1][0] = 180*in[1]; 205 out[11][1][1] = -90 + 180*in[0] + 180*in[1]; 206 } 207 208 //! \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) const209 void partial (const std::array<unsigned int, 2>& order, 210 const typename Traits::DomainType& in, // position 211 std::vector<typename Traits::RangeType>& out) const // return value 212 { 213 auto totalOrder = std::accumulate(order.begin(), order.end(), 0); 214 if (totalOrder == 0) { 215 evaluateFunction(in, out); 216 } else if (totalOrder == 1) { 217 out.resize(size()); 218 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1)); 219 220 switch (direction) { 221 case 0: 222 out[0][0] = sign_[0]*(-2*in[1] + 2*in[0]); 223 out[0][1] = sign_[0]*(-2*in[1]); 224 225 out[1][0] = 1.5 + 3*in[1] - 9*in[0]; 226 out[1][1] = 6 - 15*in[1]; 227 228 out[2][0] = sign_[0]*(-7.5 + 5*in[1] + 25*in[0]); 229 out[2][1] = sign_[0]*(30 - 25*in[1] - 60*in[0]); 230 231 out[3][0] = sign_[1]*(6 - 2*in[1] - 10*in[0]); 232 out[3][1] = sign_[1]*(-2*in[1]); 233 234 out[4][0] = -10.5 + 15*in[1] + 15*in[0]; 235 out[4][1] = -3*in[1]; 236 237 out[5][0] = sign_[1]*(7.5 - 25*in[1] - 5*in[0]); 238 out[5][1] = sign_[1]*(5*in[1]); 239 240 out[6][0] = sign_[2]*(-3 + 4*in[1] + 8*in[0]); 241 out[6][1] = sign_[2]*(4*in[1]); 242 243 out[7][0] = -3 + 12*in[0]; 244 out[7][1] = 0; 245 246 out[8][0] = sign_[2]*(-10*in[1] + 10*in[0]); 247 out[8][1] = sign_[2]*(-10*in[1]); 248 249 out[9][0] = 18 - 12*in[1] - 36*in[0]; 250 out[9][1] = -12*in[1]; 251 252 out[10][0] = 6 - 12*in[1] - 12*in[0]; 253 out[10][1] = -12*in[1]; 254 255 out[11][0] = 90 - 180*in[1] - 180*in[0]; 256 out[11][1] = 180*in[1]; 257 break; 258 case 1: 259 out[0][0] = sign_[0]*(-2*in[0]); 260 out[0][1] = sign_[0]*(6 -2*in[0] - 10*in[1]); 261 262 out[1][0] = 3*in[0]; 263 out[1][1] = 10.5 - 15*in[0] - 15*in[1]; 264 265 out[2][0] = sign_[0]*(5*in[0]); 266 out[2][1] = sign_[0]*(7.5 - 25*in[0] - 5*in[1]); 267 268 out[3][0] = sign_[1]*(-2*in[0]); 269 out[3][1] = sign_[1]*(-2*in[0] + 2*in[1]); 270 271 out[4][0] = -6 + 15*in[0]; 272 out[4][1] = -1.5 - 3*in[0] + 9*in[1]; 273 274 out[5][0] = sign_[1]*(30 - 25*in[0] - 60*in[1]); 275 out[5][1] = sign_[1]*(-7.5 + 5*in[0] + 25*in[1]); 276 277 out[6][0] = sign_[2]*(4*in[0]); 278 out[6][1] = sign_[2]*(-3 + 4*in[0] + 8*in[1]); 279 280 out[7][0] = 0; 281 out[7][1] = 3 - 12*in[1]; 282 283 out[8][0] = sign_[2]*(-10*in[0]); 284 out[8][1] = sign_[2]*(-10*in[0] + 10*in[1]); 285 286 out[9][0] = -12*in[0]; 287 out[9][1] = 6 - 12*in[0] - 12*in[1]; 288 289 out[10][0] = -12*in[0]; 290 out[10][1] = 18 - 12*in[0] - 36*in[1]; 291 292 out[11][0] = -180*in[0]; 293 out[11][1] = -90 + 180*in[0] + 180*in[1]; 294 break; 295 default: 296 DUNE_THROW(RangeError, "Component out of range."); 297 } 298 } else { 299 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented"); 300 } 301 } 302 303 //! \brief Polynomial order of the shape functions order() const304 unsigned int order() const 305 { 306 return 2; // TODO: check whether this is not order 3 307 } 308 309 private: 310 std::array<R,3> sign_; 311 }; 312 } // end namespace Dune 313 #endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI2_SIMPLEX2D_LOCALBASIS_HH 314