elasticityhelpers.hh 9.71 KB
 Oliver Sander committed Jan 12, 2006 1 2 3 #ifndef DUNE_ELASTICITY_HELPERS_HH #define DUNE_ELASTICITY_HELPERS_HH  Oliver Sander committed Apr 16, 2009 4 5 6 7 #include #include  akbib committed Apr 17, 2013 8 9 #include  Oliver Sander committed Jan 12, 2006 10 11 12 13 namespace Dune { namespace Elasticity {  akbib committed Mar 29, 2013 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68  /** \brief Compute nonlinear strain tensor from a displacement gradient. * * \param grad The gradient of the direction in which the linearisation is computed. * \param strain The tensor to store the strain in. */ template void computeNonlinearStrain(const Dune::FieldMatrix& grad, SymmetricTensor& strain) { strain = 0; for (int i=0; i void computeLinearisedStrain(const Dune::FieldMatrix& grad, SymmetricTensor& strain) { for (int i=0; i void computeLinearisedStrain(const Dune::FieldMatrix& grad, const Dune::FieldMatrix& conf, SymmetricTensor& strain) { strain = 0; for (int i=0;i double E_du(const Dune::FieldMatrix& u, int i, int j, int p, int q) { double result = delta(i,p)*delta(j,q) + delta(j,p)*delta(i,q); result += ( delta(i,q)*u[p][j] + u[p][i]*delta(j,q) ); return 0.5*result; } // template // static inline double e_dudu() // { // return 0.5 * Delta::delta // * ( Delta::delta*Delta::delta + Delta::delta*Delta::delta ); // } static inline double e_dudu(int i, int j, int p, int q, int r, int s) { return 0.5 * delta(p,r)* ( delta(i,q)*delta(j,s) + delta(i,s)*delta(j,q) ); } // The Kronecker delta function template struct Delta { static int const delta = 0; }; template struct Delta { static int const delta = 1; }; // \parder {det I + u}{ u_pq } template double det_du_tmpl(const Dune::FieldMatrix& u) { int const q = qq+3; return (Delta<(p+1)%3,(q+1)%3>::delta + u[(p+1)%3][(q+1)%3]) * (Delta<(p+2)%3,(q+2)%3>::delta+u[(p+2)%3][(q+2)%3]) - (Delta<(p+1)%3,(q-1)%3>::delta+u[(p+1)%3][(q-1)%3]) * (Delta<(p+2)%3,(q-2)%3>::delta+u[(p+2)%3][(q-2)%3]); }  akbib committed Apr 16, 2013 128   akbib committed Apr 16, 2013 129  /** \brief Compute linearization of the determinant of a deformation gradient. Multiplication with a test function(scalar) gradient, gives the linearization in the direction of the test function.  akbib committed Apr 16, 2013 130 131 132  * * \param u The displacement gradient(!) at which the determinant is evaluated */  akbib committed Apr 16, 2013 133  void linearisedDefDet(const Dune::FieldMatrix& u, Dune::FieldMatrix& linDet) {  akbib committed Apr 16, 2013 134 135 136 137  linDet = 0; for (int i=0; i<2; i++) for (int j=(i+1)%3; j<3; j++) {  akbib committed Apr 17, 2013 138  int k=(-(i+j))%3;  akbib committed Apr 16, 2013 139 140 141 142 143 144 145 146  linDet[i][j] = u[j][k]*u[k][i] - u[j][i]*(1+u[k][k]); linDet[j][i] = u[k][j]*u[i][k] - u[i][j]*(1+u[k][k]); } // the diagonal parts for (int i=0; i<3; i++) linDet[i][i] = (u[(i+1)%3][(i+1)%3]+1)*(u[(i+2)%3][(i+2)%3]+1)-u[(i+1)%3][(i+2)%3]*u[(i+2)%3][(i+1)%3] }  Oliver Sander committed Jan 12, 2006 147   akbib committed Apr 16, 2013 148  /** \brief Compute linearization of the determinant of a deformation gradient. Multiplication with a test function(scalar) gradient, gives the linearization in the direction of the test function.  akbib committed Apr 16, 2013 149 150 151  * * \param u The displacement gradient(!) at which the determinant is evaluated */  akbib committed Apr 16, 2013 152  void linearisedDefDet(const Dune::FieldMatrix& u, Dune::FieldMatrix& linDet) {  akbib committed Apr 16, 2013 153 154 155 156 157 158 159 160  linDet = 0; for (int i=0; i<2; i++) { linDet[i][i] = 1+u[(i+1)%2][(i+1)%2]; linDet[i][(i+1)%2] = -u[(i+1)%2][i]; } }  Oliver Sander committed Jan 12, 2006 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217  // \parder {det I + u}{ u_pq } template double det_du_tmpl(const Dune::FieldMatrix& u) { return Delta::delta * (1+u[(p+1)%2][(q+1)%2]) - Delta::delta * u[(p+1)%2][(q+1)%2]; // p = 0, q = 0; (1+u[1][1]) // p = 0, q = 1; -u[1][0] // p = 1, q = 0; -u[0][1] // p = 1, q = 1; (1+u[0][0]) // } // template // double det_dudu_tmpl(const Dune::FieldMatrix& u) { // int const q = qq+3; // return // Delta<(p+1)%3,r>::delta*Delta<(q+1)%3,s>::delta * (Delta<(p+2)%3,(q+2)%3>::delta+u[(p+2)%3][(q+2)%3]) // + Delta<(p+2)%3,r>::delta*Delta<(q+2)%3,s>::delta * (Delta<(p+1)%3,(q+1)%3>::delta+u[(p+1)%3][(q+1)%3]) // - Delta<(p+1)%3,r>::delta*Delta<(q-1)%3,s>::delta * (Delta<(p+2)%3,(q-2)%3>::delta+u[(p+2)%3][(q-2)%3]) // - Delta<(p+2)%3,r>::delta*Delta<(q-2)%3,s>::delta * (Delta<(p+1)%3,(q-1)%3>::delta+u[(p+1)%3][(q-1)%3]); // } double det_dudu_tmpl(const Dune::FieldMatrix& u, int p, int qq, int r, int s) { int const q = qq+3; return delta((p+1)%3,r) * delta((q+1)%3,s) * (delta((p+2)%3,(q+2)%3)+u[(p+2)%3][(q+2)%3]) + delta((p+2)%3,r)*delta((q+2)%3,s) * (delta((p+1)%3,(q+1)%3)+u[(p+1)%3][(q+1)%3]) - delta((p+1)%3,r)*delta((q-1)%3,s) * (delta((p+2)%3,(q-2)%3)+u[(p+2)%3][(q-2)%3]) - delta((p+2)%3,r)*delta((q-2)%3,s) * (delta((p+1)%3,(q-1)%3)+u[(p+1)%3][(q-1)%3]); } // parder det (I + u) partial u_pq \partial u_rs //template double det_dudu_tmpl(const Dune::FieldMatrix& u, int p, int q, int r, int s) { return delta(p,q)*delta(r,(p+1)%2)*delta(s,(q+1)%2) - delta(p,(q+1)%2)*delta(r,(p+1)%2)*delta(s,(q+1)%2); // 0 0 0 0 --> 0 - 0 // 0 0 0 1 --> 0 - 0 // 0 0 1 0 --> 0 - 0 // 0 0 1 1 --> 1 - 0 = 1 // 0 1 0 0 --> 0 - 0 // 0 1 0 1 --> 0 - 0 // 0 1 1 0 --> 0 - 1 = -1 // 0 1 1 1 --> 0 - 0 // 1 0 0 0 --> 0 - 0 // 1 0 0 1 --> 0 - 1 = -1 // 1 0 1 0 --> 0 - 0 // 1 0 1 1 --> 0 - 0 // 1 1 0 0 --> 1 - 0 = 1 // 1 1 0 1 --> 0 - 0 // 1 1 1 0 --> 0 - 0 // 1 1 1 1 --> 0 - 0 }  akbib committed Apr 17, 2013 218   Oliver Sander committed Jan 12, 2006 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248  double det_val(const Dune::FieldMatrix& u) { return (1+u[0][0])*(1+u[1][1])*(1+u[2][2]) + u[0][1]*u[1][2]*u[2][0] + u[0][2]*u[1][0]*u[2][1] - u[2][0]*(1+u[1][1])*u[0][2] - u[2][1]*u[1][2]*(1+u[0][0]) - (1+u[2][2])*u[1][0]*u[0][1]; } double det_val(const Dune::FieldMatrix& u) { return (1+u[0][0])*(1+u[1][1]) - u[1][0]*u[0][1]; } template double E_val(const Dune::FieldMatrix& u, int i, int j) { double result = u[i][j] + u[j][i]; for (int k=0; k