neohookeoperatorassembler.hh 8.36 KB
 akbib committed Apr 17, 2013 1 2 3 4 5 6 7 8 9 10 11 12 13 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 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 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 #ifndef NEO_HOOKE_OPERATOR_ASSEMBLER_HH #define NEO_HOOKE_OPERATOR_ASSEMBLER_HH #include #include #include #include #include #include #include #include #include /** \brief Assemble the second derivative of a neo-hookean material. * * Laursen: * \f[ * W(u)= \frac{\lambda}4(J^2-1)-(\frac{\lambda}2+\mu)log(J)+\mu tr(E(u)) * \f] * * Fischer/Wriggers: * \f[ * W(u)= \frac{\lambda}2(J-1)^2-2\mu)log(J)+\mu tr(E(u)) * \f] * * * - \f$D\varphi\f$: the deformation gradient * - \f$E\f$: the nonlinear strain tensor * - \f$E'_u(v)\f$: the linearized strain tensor at the point u in direction v * - \f$\sigma(u) = C:E(u)\f$: the second Piola-Kirchhoff stress tensor * - \f$C\f$: the Hooke tensor * - \f$Dv\f$: the gradient of the test function * */ template class NeoHookeOperatorAssembler : public LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, Dune::FieldMatrix > { static const int dim = GridType::dimension; typedef typename GridType::ctype ctype; typedef VirtualGridFunction > GridFunction; public: typedef typename Dune::FieldMatrix T; typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, T >::Element Element; typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE,T >::BoolMatrix BoolMatrix; typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE,T >::LocalMatrix LocalMatrix; NeoHookeOperatorAssembler(ctype lambda, ctype mu, const Dune::shared_ptr displacement): lambda_(lambda), mu_(mu), displacement_(displacement) {} NeoHookeOperatorAssembler(ctype lambda, ctype mu): lambda_(lambda), mu_(mu) {} void indices(const Element& element, BoolMatrix& isNonZero, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const { isNonZero = true; } void assemble(const Element& element, LocalMatrix& localMatrix, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const { typedef Dune::FieldVector FVdim; typedef Dune::FieldMatrix FMdimdim; typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType; int rows = tFE.localBasis().size(); int cols = aFE.localBasis().size(); localMatrix.setSize(rows,cols); localMatrix = 0.0; // get quadrature rule const int order = (element.type().isSimplex()) ? 4*(tFE.localBasis().order()) : 4*(tFE.localBasis().order()*dim); const Dune::template QuadratureRule& quad = QuadratureRuleCache::rule(element.type(), order, IsRefinedLocalFiniteElement::value(tFE) ); // store gradients of shape functions and base functions std::vector referenceGradients(tFE.localBasis().size()); std::vector gradients(tFE.localBasis().size()); // store entity geometry mapping const typename Element::Geometry geometry = element.geometry(); // loop over quadrature points for (size_t pt=0; pt < quad.size(); ++pt) { // get quadrature point const FVdim& quadPos = quad[pt].position(); // get transposed inverse of Jacobian of transformation const FMdimdim& invJacobian = geometry.jacobianInverseTransposed(quadPos); // get integration factor const ctype integrationElement = geometry.integrationElement(quadPos); // get gradients of shape functions tFE.localBasis().evaluateJacobian(quadPos, referenceGradients); // compute gradients of base functions for (size_t i=0; iisDefinedOn(element)) displacement_->evaluateDerivativeLocal(element, quadPos, localConfGrad); else displacement_->evaluateDerivative(geometry.global(quadPos),localConfGrad); // compute linearization of the determinante of the deformation gradient FMdimdim linDefDet; Dune::Elasticity::linearisedDefDet(localConfGrad,linDefDet); // compute deformation gradient of the configuration for (int i=0;i newConfig) { displacement_ = newConfig; } private: /** \brief Lame constant */ ctype lambda_; /** \brief Lame constant */ ctype mu_; /** \brief The configuration at which the linearized operator is evaluated.*/ Dune::shared_ptr displacement_; /** \brief Compute the matrix entry that arises from the second derivative of the determinant. * Note that it is assumed the defGrad is the deformation gradient and not the displacement gradient */ Dune::FieldMatrix hessDefDet(const Dune::FieldMatrix& defGrad,  akbib committed Apr 17, 2013 195  const Dune::FieldVector& testGrad, const Dune::FieldVector& ansatzGrad) const  akbib committed Apr 17, 2013 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212  { Dune::FieldMatrix res(0); if (dim==2) { res[0][1] = testGrad[0]*ansatzGrad[1]-ansatzGrad[0]*testGrad[1]; res[1][0] = -res[0][1]; return res; } //compute the cross product Dune::FieldVector cross(0); for (int k=0; k for (int i=0; i