diff --git a/dune/elasticity/materials/mooneyrivlinenergy.hh b/dune/elasticity/materials/mooneyrivlinenergy.hh new file mode 100644 index 0000000000000000000000000000000000000000..1acc797fde3257a9ab8329853a9f91cd0298c202 --- /dev/null +++ b/dune/elasticity/materials/mooneyrivlinenergy.hh @@ -0,0 +1,192 @@ +#ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH +#define DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH + +#include <dune/common/fmatrix.hh> +#include <dune/common/fmatrixev.hh> + +#include <dune/geometry/quadraturerules.hh> + +#include <dune/elasticity/assemblers/localfestiffness.hh> + +namespace Dune { + +template<class GridView, class LocalFiniteElement, class field_type=double> +class MooneyRivlinEnergy + : public LocalFEStiffness<GridView, + LocalFiniteElement, + std::vector<Dune::FieldVector<field_type,GridView::dimension> > > +{ + // grid types + typedef typename GridView::Grid::ctype DT; + typedef typename GridView::template Codim<0>::Entity Entity; + + // some other sizes + enum {gridDim=GridView::dimension}; + enum {dim=GridView::dimension}; + +public: + + /** \brief Constructor with a set of material parameters + * \param parameters The material parameters + */ + MooneyRivlinEnergy(const Dune::ParameterTree& parameters) + { + // Mooneyrivlin constants + mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a"); + mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b"); + mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c"); + + mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10"); + mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01"); + mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20"); + mooneyrivlin_02 = parameters.get<double>("mooneyrivlin_02"); + mooneyrivlin_11 = parameters.get<double>("mooneyrivlin_11"); + mooneyrivlin_30 = parameters.get<double>("mooneyrivlin_30"); + mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03"); + mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21"); + mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12"); + mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k"); + + mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy"); + + } + + /** \brief Assemble the energy for a single element */ + field_type energy (const Entity& e, + const LocalFiniteElement& localFiniteElement, + const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const; + + /** \brief Lame constants */ + field_type mooneyrivlin_a, + mooneyrivlin_b, + mooneyrivlin_c, + mooneyrivlin_10, + mooneyrivlin_01, + mooneyrivlin_20, + mooneyrivlin_02, + mooneyrivlin_11, + mooneyrivlin_30, + mooneyrivlin_21, + mooneyrivlin_12, + mooneyrivlin_03, + mooneyrivlin_k; + std::string mooneyrivlin_energy; +}; + +template <class GridView, class LocalFiniteElement, class field_type> +field_type +MooneyRivlinEnergy<GridView, LocalFiniteElement, field_type>:: +energy(const Entity& element, + const LocalFiniteElement& localFiniteElement, + const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const +{ + + assert(element.type() == localFiniteElement.type()); + typedef typename GridView::template Codim<0>::Entity::Geometry Geometry; + + field_type strainEnergyCiarlet = 0; + field_type strainEnergy = 0; + field_type strainEnergyWithLog = 0; + field_type strainEnergyWithSquare = 0; + + // store gradients of shape functions and base functions + std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size()); + std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size()); + + int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order() + : localFiniteElement.localBasis().order() * gridDim; + + const Dune::QuadratureRule<DT, gridDim>& quad + = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder); + + for (size_t pt=0; pt<quad.size(); pt++) + { + // Local position of the quadrature point + const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position(); + + const DT integrationElement = element.geometry().integrationElement(quadPos); + + const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos); + + DT weight = quad[pt].weight() * integrationElement; + + // get gradients of shape functions + localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients); + + // compute gradients of base functions + for (size_t i=0; i<gradients.size(); ++i) + jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); + + Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0); + for (size_t i=0; i<gradients.size(); i++) + for (int j=0; j<gridDim; j++) + derivative[j].axpy(localConfiguration[i][j], gradients[i]); + + + Dune::FieldMatrix<field_type,gridDim,gridDim> C(0); + + for (int i=0; i<gridDim; i++) { + for (int j=0; j<gridDim; j++) { + for (int k=0; k<gridDim; k++) + C[i][j] += derivative[k][i] * derivative[k][j]; + } + } + + ////////////////////////////////////////////////////////// + // Eigenvalues of FTF + ////////////////////////////////////////////////////////// + + // eigenvalues of C, i.e., squared singular values \sigma_i of F + Dune::FieldVector<field_type, dim> sigmaSquared; + FMatrixHelp::eigenValues(C, sigmaSquared); + + field_type normFSquared = derivative.frobenius_norm2(); + field_type detF = derivative.determinant(); + + field_type normFinvSquared = 0; + + field_type c2Tilde = 0; + for (int i = 0; i < dim; i++) { + normFinvSquared += 1/sigmaSquared[i]; + // compute D, which is the sum of the squared eigenvalues + for (int j = i+1; j < dim; j++) + c2Tilde += sigmaSquared[j]*sigmaSquared[i]; + } + + field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3; + // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3) + c2Tilde /= pow(detF, 4.0/dim); + field_type c2TildeMinus3 = c2Tilde - 3; + field_type detFMinus1 = detF - 1; + + // Add the local energy density + strainEnergyCiarlet += weight * (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF)); + strainEnergy = 0; + strainEnergy = mooneyrivlin_10 * trCTildeMinus3 + + mooneyrivlin_01 * c2TildeMinus3 + + mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 + + mooneyrivlin_02 * c2TildeMinus3 * c2TildeMinus3 + + mooneyrivlin_11 * trCTildeMinus3 * c2TildeMinus3 + + mooneyrivlin_30 * trCTildeMinus3 * trCTildeMinus3 * trCTildeMinus3 + + mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3 + + mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3 * c2TildeMinus3 + + mooneyrivlin_03 * c2TildeMinus3 * c2TildeMinus3 * c2TildeMinus3; + field_type logDetF = std::log(detF); + strainEnergyWithLog += weight * ( strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF ); + strainEnergyWithSquare += weight * ( strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1); + } + + std::cout << std::scientific; + if (mooneyrivlin_energy == "log") { + return strainEnergyWithLog; + } else if (mooneyrivlin_energy == "square") { + return strainEnergyWithSquare; + } else { + return strainEnergyCiarlet; + } + std::cout << std::fixed; +} + +} // namespace Dune + +#endif //#ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH diff --git a/src/finite-strain-elasticity-bending.parset b/src/finite-strain-elasticity-bending.parset index c05f7828ba35399306471e2b6a575116f0dcb599..a8e3d15661a36284c2cce19a71157d6241fe41e6 100644 --- a/src/finite-strain-elasticity-bending.parset +++ b/src/finite-strain-elasticity-bending.parset @@ -62,6 +62,22 @@ energy = stvenantkirchhoff mu = 2.7191e+6 lambda = 4.4364e+6 +mooneyrivlin_a = 2.7191e+6 +mooneyrivlin_b = 2.7191e+6 +mooneyrivlin_c = 2.7191e+6 +mooneyrivlin_10 = -7.28e+5 +mooneyrivlin_01 = 9.17e+5 +mooneyrivlin_20 = 1.23e+5 +mooneyrivlin_02 = 8.15e+5 +mooneyrivlin_11 = -5.14e+5 +mooneyrivlin_30 = 0 +mooneyrivlin_21 = 0 +mooneyrivlin_12 = 0 +mooneyrivlin_03 = 0 + +mooneyrivlin_k = 1e+6 # volume-preserving parameter +mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy + [] ############################################# diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc index d40c095f48ce2bb6a1e011d3dce028f482076638..62586d2a65ec317886d3cef12c798def076a3b8a 100644 --- a/src/finite-strain-elasticity.cc +++ b/src/finite-strain-elasticity.cc @@ -38,6 +38,7 @@ #include <dune/elasticity/materials/stvenantkirchhoffenergy.hh> #include <dune/elasticity/materials/henckyenergy.hh> #include <dune/elasticity/materials/exphenckyenergy.hh> +#include <dune/elasticity/materials/mooneyrivlinenergy.hh> #include <dune/elasticity/materials/neohookeenergy.hh> #include <dune/elasticity/materials/neumannenergy.hh> #include <dune/elasticity/materials/sumenergy.hh> @@ -285,11 +286,14 @@ int main (int argc, char *argv[]) try elasticEnergy = std::make_shared<HenckyEnergy<GridView, FEBasis::LocalView::Tree::FiniteElement, adouble> >(materialParameters); - if (parameterSet.get<std::string>("energy") == "exphencky") elasticEnergy = std::make_shared<ExpHenckyEnergy<GridView, FEBasis::LocalView::Tree::FiniteElement, adouble> >(materialParameters); + if (parameterSet.get<std::string>("energy") == "mooneyrivlin") + elasticEnergy = std::make_shared<MooneyRivlinEnergy<GridView, + FEBasis::LocalView::Tree::FiniteElement, + adouble> >(materialParameters); if(!elasticEnergy) DUNE_THROW(Exception, "Error: Selected energy not available!");