#ifndef DUNE_ELASTICITY_HELPERS_HH #define DUNE_ELASTICITY_HELPERS_HH #include #include namespace Dune { namespace Elasticity { /** \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]); } /** \brief Compute linearization of the determinant of a deformation gradient * * \param u The displacement gradient(!) at which the determinant is evaluated */ void linearizedDefDet(const Dune::FieldMatrix& u, Dune::FieldMatrix& linDet) { linDet = 0; for (int i=0; i<2; i++) for (int j=(i+1)%3; j<3; j++) { int k=(j+1)%3 + (j+1)/3; 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] } /** \brief Compute linearization of the determinant of a deformation gradient * * \param u The displacement gradient(!) at which the determinant is evaluated */ void linearizedDefDet(const Dune::FieldMatrix& u, Dune::FieldMatrix& linDet) { 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]; } } // \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 } 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