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_HIERARCHICAL_PRISM_P2_LOCALBASIS_HH 4 #define DUNE_HIERARCHICAL_PRISM_P2_LOCALBASIS_HH 5 6 /** \file 7 \brief Hierarchical prism p2 shape functions for the simplex 8 */ 9 10 #include <numeric> 11 12 #include <dune/common/fvector.hh> 13 #include <dune/common/fmatrix.hh> 14 15 #include <dune/localfunctions/common/localbasis.hh> 16 17 namespace Dune 18 { 19 template<class D, class R> 20 class HierarchicalPrismP2LocalBasis 21 { 22 public: 23 //! \brief export type traits for function signature 24 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>, Dune::FieldMatrix<R,1,3> > Traits; 25 26 //! \brief number of shape functions size() const27 unsigned int size () const 28 { 29 return 18; 30 } 31 32 //! \brief Evaluate all shape functions evaluateFunction(const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const33 void evaluateFunction (const typename Traits::DomainType& in, 34 std::vector<typename Traits::RangeType> & out) const 35 { 36 out.resize(18); 37 38 out[0]=(1.0-in[0]-in[1])*(1.0-in[2]); 39 out[1]= in[0]*(1-in[2]); 40 out[2]=in[1]*(1-in[2]); 41 out[3]=in[2]*(1.0-in[0]-in[1]); 42 out[4]=in[0]*in[2]; 43 out[5]=in[1]*in[2]; 44 45 //edges 46 out[6]=2*(1.0-in[0]-in[1])*(0.5-in[0]-in[1])*(4*in[2]-4*in[2]*in[2]); 47 out[7]=2*in[0]*(-0.5+in[0])*(4*in[2]-4*in[2]*in[2]); 48 out[8]=2*in[1]*(-0.5+in[1])*(4*in[2]-4*in[2]*in[2]); 49 out[9]=4*in[0]*(1-in[0]-in[1])*(1-3*in[2]+2*in[2]*in[2]); 50 out[10]=4*in[1]*(1-in[0]-in[1])*(1-3*in[2]+2*in[2]*in[2]); 51 out[11]=4*in[0]*in[1]*(1-3*in[2]+2*in[2]*in[2]); 52 out[12]=4*in[0]*(1-in[0]-in[1])*(-in[2]+2*in[2]*in[2]); 53 out[13]=4*in[1]*(1-in[0]-in[1])*(-in[2]+2*in[2]*in[2]); 54 out[14]=4*in[0]*in[1]*(-in[2]+2*in[2]*in[2]); 55 56 //faces 57 out[15]=4*in[0]*(1-in[0]-in[1])*(4*in[2]-4*in[2]*in[2]); 58 out[16]=4*in[1]*(1-in[0]-in[1])*(4*in[2]-4*in[2]*in[2]); 59 out[17]=4*in[0]*in[1]*(4*in[2]-4*in[2]*in[2]); 60 } 61 62 63 64 //! \brief Evaluate Jacobian of all shape functions evaluateJacobian(const typename Traits::DomainType & in,std::vector<typename Traits::JacobianType> & out) const65 void evaluateJacobian (const typename Traits::DomainType& in, //position 66 std::vector<typename Traits::JacobianType>& out) const //return value 67 { 68 out.resize(18); 69 70 //vertices 71 out[0][0][0] = in[2]-1; 72 out[0][0][1] = in[2]-1; 73 out[0][0][2] = in[0]+in[1]-1; 74 75 out[1][0][0] = 1-in[2]; 76 out[1][0][1] = 0; 77 out[1][0][2] =-in[0]; 78 79 out[2][0][0] = 0; 80 out[2][0][1] = 1-in[2]; 81 out[2][0][2] = -in[1]; 82 83 out[3][0][0] = -in[2]; 84 out[3][0][1] = -in[2]; 85 out[3][0][2] = 1-in[0]-in[1]; 86 87 out[4][0][0] = in[2]; 88 out[4][0][1] = 0; 89 out[4][0][2] = in[0]; 90 91 out[5][0][0] = 0; 92 out[5][0][1] = in[2]; 93 out[5][0][2] = in[1]; 94 95 //edges 96 out[6][0][0] = (-3+4*in[0]+4*in[1])*(4*in[2]-4*in[2]*in[2]); 97 out[6][0][1] = (-3+4*in[0]+4*in[1])*(4*in[2]-4*in[2]*in[2]); 98 out[6][0][2] = 2*(1-in[0]-in[1])*(0.5-in[0]-in[1])*(4-8*in[2]); 99 100 out[7][0][0] = (-1+4*in[0])*(4*in[2]-4*in[2]*in[2]); 101 out[7][0][1] = 0; 102 out[7][0][2] = 2*in[0]*(-0.5+in[0])*(4-8*in[2]); 103 104 out[8][0][0] = 0; 105 out[8][0][1] = (-1+4*in[1])*(4*in[2]-4*in[2]*in[2]); 106 out[8][0][2] = 2*in[1]*(-0.5+in[1])*(4-8*in[2]); 107 108 out[9][0][0] = (4-8*in[0]-4*in[1])*(1-3*in[2]+2*in[2]*in[2]); 109 out[9][0][1] = -4*in[0]*(1-3*in[2]+2*in[2]*in[2]); 110 out[9][0][2] = 4*in[0]*(1-in[0]-in[1])*(-3+4*in[2]); 111 112 out[10][0][0] = (-4*in[1])*(1-3*in[2]+2*in[2]*in[2]); 113 out[10][0][1] = (4-4*in[0]-8*in[1])*(1-3*in[2]+2*in[2]*in[2]); 114 out[10][0][2] = 4*in[1]*(1-in[0]-in[1])*(-3+4*in[2]); 115 116 out[11][0][0] = 4*in[1]*(1-3*in[2]+2*in[2]*in[2]); 117 out[11][0][1] = 4*in[0]*(1-3*in[2]+2*in[2]*in[2]); 118 out[11][0][2] = 4*in[0]*in[1]*(-3+4*in[2]); 119 120 out[12][0][0] = (4-8*in[0]-4*in[1])*(-in[2]+2*in[2]*in[2]); 121 out[12][0][1] = (-4*in[0])*(-in[2]+2*in[2]*in[2]); 122 out[12][0][2] = 4*in[0]*(1-in[0]-in[1])*(-1+4*in[2]); 123 124 out[13][0][0] = -4*in[1]*(-in[2]+2*in[2]*in[2]); 125 out[13][0][1] = (4-4*in[0]-8*in[1])*(-in[2]+2*in[2]*in[2]); 126 out[13][0][2] = 4*in[1]*(1-in[0]-in[1])*(-1+4*in[2]); 127 128 out[14][0][0] = 4*in[1]*(-in[2]+2*in[2]*in[2]); 129 out[14][0][1] = 4*in[0]*(-in[2]+2*in[2]*in[2]); 130 out[14][0][2] = 4*in[0]*in[1]*(-1+4*in[2]); 131 132 //faces 133 out[15][0][0] = (4-8*in[0]-4*in[1])*(4*in[2]-4*in[2]*in[2]); 134 out[15][0][1] = -4*in[0]*(4*in[2]-4*in[2]*in[2]); 135 out[15][0][2] = 4*in[0]*(1-in[0]-in[1])*(4-8*in[2]); 136 137 out[16][0][0] = -4*in[1]*(4*in[2]-4*in[2]*in[2]); 138 out[16][0][1] = (4-4*in[0]-8*in[1])*(4*in[2]-4*in[2]*in[2]); 139 out[16][0][2] = 4*in[1]*(1-in[0]-in[1])*(4-8*in[2]); 140 141 out[17][0][0] = 4*in[1]*(4*in[2]-4*in[2]*in[2]); 142 out[17][0][1] = 4*in[0]*(4*in[2]-4*in[2]*in[2]); 143 out[17][0][2] = 4*in[0]*in[1]*(4-8*in[2]); 144 } 145 146 //! \brief Evaluate partial derivatives of all shape functions partial(const std::array<unsigned int,3> & order,const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const147 void partial (const std::array<unsigned int, 3>& order, 148 const typename Traits::DomainType& in, // position 149 std::vector<typename Traits::RangeType>& out) const // return value 150 { 151 auto totalOrder = std::accumulate(order.begin(), order.end(), 0); 152 if (totalOrder == 0) { 153 evaluateFunction(in, out); 154 } else if (totalOrder == 1) { 155 out.resize(size()); 156 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1)); 157 158 switch (direction) { 159 case 0: 160 out[0] = in[2]-1; 161 out[1] = 1-in[2]; 162 out[2] = 0; 163 out[3] = -in[2]; 164 out[4] = in[2]; 165 out[5] = 0; 166 out[6] = (-3+4*in[0]+4*in[1])*(4*in[2]-4*in[2]*in[2]); 167 out[7] = (-1+4*in[0])*(4*in[2]-4*in[2]*in[2]); 168 out[8] = 0; 169 out[9] = (4-8*in[0]-4*in[1])*(1-3*in[2]+2*in[2]*in[2]); 170 out[10] = (-4*in[1])*(1-3*in[2]+2*in[2]*in[2]); 171 out[11] = 4*in[1]*(1-3*in[2]+2*in[2]*in[2]); 172 out[12] = (4-8*in[0]-4*in[1])*(-in[2]+2*in[2]*in[2]); 173 out[13] = -4*in[1]*(-in[2]+2*in[2]*in[2]); 174 out[14] = 4*in[1]*(-in[2]+2*in[2]*in[2]); 175 out[15] = (4-8*in[0]-4*in[1])*(4*in[2]-4*in[2]*in[2]); 176 out[16] = -4*in[1]*(4*in[2]-4*in[2]*in[2]); 177 out[17] = 4*in[1]*(4*in[2]-4*in[2]*in[2]); 178 break; 179 case 1: 180 out[0] = in[2]-1; 181 out[1] = 0; 182 out[2] = 1-in[2]; 183 out[3] = -in[2]; 184 out[4] = 0; 185 out[5] = in[2]; 186 out[6] = (-3+4*in[0]+4*in[1])*(4*in[2]-4*in[2]*in[2]); 187 out[7] = 0; 188 out[8] = (-1+4*in[1])*(4*in[2]-4*in[2]*in[2]); 189 out[9] = -4*in[0]*(1-3*in[2]+2*in[2]*in[2]); 190 out[10] = (4-4*in[0]-8*in[1])*(1-3*in[2]+2*in[2]*in[2]); 191 out[11] = 4*in[0]*(1-3*in[2]+2*in[2]*in[2]); 192 out[12] = (-4*in[0])*(-in[2]+2*in[2]*in[2]); 193 out[13] = (4-4*in[0]-8*in[1])*(-in[2]+2*in[2]*in[2]); 194 out[14] = 4*in[0]*(-in[2]+2*in[2]*in[2]); 195 out[15] = -4*in[0]*(4*in[2]-4*in[2]*in[2]); 196 out[16] = (4-4*in[0]-8*in[1])*(4*in[2]-4*in[2]*in[2]); 197 out[17] = 4*in[0]*(4*in[2]-4*in[2]*in[2]); 198 break; 199 case 2: 200 out[0] = in[0]+in[1]-1; 201 out[1] =-in[0]; 202 out[2] = -in[1]; 203 out[3] = 1-in[0]-in[1]; 204 out[4] = in[0]; 205 out[5] = in[1]; 206 out[6] = 2*(1-in[0]-in[1])*(0.5-in[0]-in[1])*(4-8*in[2]); 207 out[7] = 2*in[0]*(-0.5+in[0])*(4-8*in[2]); 208 out[8] = 2*in[1]*(-0.5+in[1])*(4-8*in[2]); 209 out[9] = 4*in[0]*(1-in[0]-in[1])*(-3+4*in[2]); 210 out[10] = 4*in[1]*(1-in[0]-in[1])*(-3+4*in[2]); 211 out[11] = 4*in[0]*in[1]*(-3+4*in[2]); 212 out[12] = 4*in[0]*(1-in[0]-in[1])*(-1+4*in[2]); 213 out[13] = 4*in[1]*(1-in[0]-in[1])*(-1+4*in[2]); 214 out[14] = 4*in[0]*in[1]*(-1+4*in[2]); 215 out[15] = 4*in[0]*(1-in[0]-in[1])*(4-8*in[2]); 216 out[16] = 4*in[1]*(1-in[0]-in[1])*(4-8*in[2]); 217 out[17] = 4*in[0]*in[1]*(4-8*in[2]); 218 break; 219 default: 220 DUNE_THROW(RangeError, "Component out of range."); 221 } 222 } else { 223 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented"); 224 } 225 } 226 227 /** \brief Polynomial order of the shape functions 228 */ order() const229 unsigned int order() const 230 { 231 return 2; 232 } 233 234 }; 235 } 236 #endif 237