Add membrane-energy and parset files for Patrizios experiment

@osander: can you have a look?

dune/gfe/cosseratmembranestiffness.hh

+#ifndef COSSERAT_MEMBRANE_STIFFNESS_HH
+#define COSSERAT_MEMBRANE_STIFFNESS_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/gfe/localenergy.hh>
+#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
+#ifdef PROJECTED_INTERPOLATION
+#include <dune/gfe/localprojectedfefunction.hh>
+#else
+#include "localgeodesicfefunction.hh"
+#endif
+#include <dune/gfe/rigidbodymotion.hh>
+#include <dune/gfe/tensor3.hh>
+#include <dune/gfe/orthogonalmatrix.hh>
+#include <dune/gfe/cosseratstrain.hh>
+#include <dune/gfe/cosseratenergystiffness.hh>
+
+template<class Basis, int dim, class field_type=double>
+class CosseratMembraneStiffness
+    : public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >,
+      public MixedLocalGeodesicFEStiffness<Basis,
+                                           RealTuple<field_type,dim>,
+                                           Rotation<field_type,dim> >
+{
+    // grid types
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::ctype DT;
+    typedef RigidBodyMotion<field_type,dim> TargetSpace;
+    typedef typename TargetSpace::ctype RT;
+    typedef typename GridView::template Codim<0>::Entity Entity;
+
+    // some other sizes
+    enum {gridDim=GridView::dimension};
+    enum {dimworld=GridView::dimensionworld};
+    // use this with: gridDim = 2, dimworld = 2, dim = 3 (the shell moves out of the 2D-setting and bends/wrinkles/moves into 3D)
+
+public:
+
+    /** \brief Constructor with a set of material parameters
+     * \param parameters The material parameters
+     */
+    CosseratMembraneStiffness(const Dune::ParameterTree& parameters,
+                                 const BoundaryPatch<GridView>* neumannBoundary,
+                                 const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
+                                 const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
+    : neumannBoundary_(neumannBoundary),
+      neumannFunction_(neumannFunction)
+    {
+        // The shell thickness
+        thickness_ = parameters.template get<double>("thickness");
+
+        // Lame constants
+        mu_ = parameters.template get<double>("mu");
+        lambda_ = parameters.template get<double>("lambda");
+
+        // Cosserat couple modulus
+        mu_c_ = parameters.template get<double>("mu_c");
+
+        // Length scale parameter
+        L_c_ = parameters.template get<double>("L_c");
+
+        // Curvature exponent
+        q_ = parameters.template get<double>("q");
+
+        // Flag for using geometric mean or not
+        useGeometricMean_ = parameters.template get<bool>("useGeometricMean", false);
+
+    }
+
+    /** \brief Assemble the energy for a single element */
+    RT energy (const typename Basis::LocalView& localView,
+               const std::vector<TargetSpace>& localSolution) const override;
+
+    /** \brief Assemble the energy for a single element */
+    RT energy (const typename Basis::LocalView& localView,
+               const std::vector<RealTuple<field_type,dim> >& localDisplacementConfiguration,
+               const std::vector<Rotation<field_type,dim> >& localOrientationConfiguration) const override;
+
+    // !!! dimworld = 2, dim = 3 !!!!
+    /** \brief Part 1 of the Membrane energy, everything but the term with the harmonic mean of mu and mu_c
+     */
+    RT membraneEnergyPart1(const Dune::FieldMatrix<field_type,dim,dimworld>& deformationGradient,
+                           const Dune::FieldMatrix<field_type,dim,dim>& R) const // R - Rotation matrix, in the format (R1|R2|R3)
+    {
+        Dune::FieldMatrix<field_type,dim,dimworld> R1R2Transposed(0);
+        for (int i=0; i<dimworld; i++)
+            for (int j=0; j<dim; j++)
+                R1R2Transposed[i][j] = R[j][i]; //check all actions on R!!! does that not need to be the other way round?
+
+        // Multiply R1R2Transposed with the deformationGradient, this returns a 2x2 matrix
+        // the "elastic non-symmetric stretch tensor in 2D" - U2D
+        // And directly subtract 1 on the diagonal ;)
+        Dune::FieldMatrix<field_type,dimworld,dimworld> U2DMinus1(0);
+        for (int i=0; i<dimworld; i++) {
+            for (int j=0; j<dimworld; j++) {
+                for (int k=0; k<dim; k++) {
+                    U2DMinus1[i][j] += R1R2Transposed[i][k]*deformationGradient[k][j];
+                }
+            }
+            U2DMinus1[i][i] -= 1;
+        }
+
+        return mu_ * Dune::GFE::sym(U2DMinus1).frobenius_norm2()
+                + mu_c_ * Dune::GFE::skew(U2DMinus1).frobenius_norm2()
+                + (mu_*lambda_)/(2*mu_ + lambda_) * Dune::GFE::traceSquared(U2DMinus1);
+    }
+
+    /** \brief Part 2 of the Membrane energy, only the term with the harmonic mean of mu and mu_c
+     */
+    RT membraneEnergyPart2(const Dune::FieldMatrix<field_type,dim,dimworld>& deformationGradient,
+                           const Dune::FieldMatrix<field_type,dim,dim>& R) const
+                           // R - Rotation matrix, in the format (R1|R2|R3), so the column R3 is (R[0][2],R[1][2],R[2][2])^T
+    {
+        //(R_3,m_x)^2 + (R_3,m_y)
+        field_type R3mx = 0;
+        field_type R3my = 0;
+        for (int i=0; i<dim; i++) {
+            R3mx += R[i][2]*deformationGradient[i][0];
+            R3my += R[i][2]*deformationGradient[i][1];
+        }
+        if (useGeometricMean_)
+            return (R3mx*R3mx + R3my*R3my)*(mu_ + mu_c_)/2;
+        else
+            return (R3mx*R3mx + R3my*R3my)*(2* mu_*mu_c_)/(mu_ + mu_c_);
+    }
+
+    RT curvatureEnergy(const Tensor3<field_type,3,3,gridDim>& DR) const
+    {
+        using std::pow;
+        return mu_ * pow(L_c_ * L_c_ * DR.frobenius_norm2(),q_/2.0)/2.0;
+    }
+
+    /** \brief The shell thickness */
+    double thickness_;
+
+    /** \brief Lame constants */
+    double mu_, lambda_;
+
+    /** \brief Cosserat couple modulus, preferably 0 */
+    double mu_c_;
+
+    /** \brief Length scale parameter */
+    double L_c_;
+
+    /** \brief Curvature exponent */
+    double q_;
+
+    /** \brief Flag for using geometric mean or not */
+    bool useGeometricMean_;
+
+    /** \brief The Neumann boundary */
+    const BoundaryPatch<GridView>* neumannBoundary_;
+
+    /** \brief The function implementing the Neumann data */
+    const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction_;
+
+    /** \brief The function implementing a volume load */
+    const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
+};
+
+template <class Basis, int dim, class field_type>
+typename CosseratMembraneStiffness<Basis,dim,field_type>::RT
+CosseratMembraneStiffness<Basis,dim,field_type>::
+energy(const typename Basis::LocalView& localView,
+       const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
+{
+    RT energy = 0;
+    DUNE_THROW(Dune::NotImplemented, "CosseratMembraneStiffness for usage with RigidBodyMotion");
+    return energy;
+}
+
+template <class Basis, int dim, class field_type>
+typename CosseratMembraneStiffness<Basis,dim,field_type>::RT
+CosseratMembraneStiffness<Basis,dim,field_type>::
+energy(const typename Basis::LocalView& localView,
+       const std::vector<RealTuple<field_type,dim> >& localDeformationConfiguration,
+       const std::vector<Rotation<field_type,dim> >& localOrientationConfiguration) const
+{
+    auto element = localView.element();
+
+    RT energy = 0;
+
+    using namespace Dune::TypeTree::Indices;
+    const auto& deformationLocalFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
+    const auto& orientationLocalFiniteElement = LocalFiniteElementFactory<Basis,1>::get(localView,_1);
+
+#ifdef PROJECTED_INTERPOLATION
+    typedef Dune::GFE::LocalProjectedFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> >
+    LocalDeformationGFEFunctionType;
+#else
+    typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> >
+    LocalDeformationGFEFunctionType;
+#endif
+    LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localDeformationConfiguration);
+
+#ifdef PROJECTED_INTERPOLATION
+    typedef Dune::GFE::LocalProjectedFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
+#else
+    typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
+#endif
+    LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localOrientationConfiguration);
+
+    // \todo Implement smarter quadrature rule selection for more efficiency, i.e., less evaluations of the Rotation GFE function
+    int quadOrder = deformationLocalFiniteElement.localBasis().order() * ((element.type().isSimplex()) ? 1 : gridDim);
+
+    const auto& 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 auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+
+        DT weight = quad[pt].weight() * integrationElement;
+
+        // The value of the local deformation
+        RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
+        Rotation<field_type,dim>  orientationValue = localOrientationGFEFunction.evaluate(quadPos);
+
+        // The derivative of the local function defined on the reference element
+        typename LocalDeformationGFEFunctionType::DerivativeType deformationReferenceDerivative = localDeformationGFEFunction.evaluateDerivative(quadPos,deformationValue);
+        typename LocalOrientationGFEFunctionType::DerivativeType orientationReferenceDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
+
+        // The derivative of the function defined on the actual element
+        typename LocalDeformationGFEFunctionType::DerivativeType deformationDerivative;
+        typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
+
+        for (size_t comp=0; comp<deformationReferenceDerivative.N(); comp++)
+            jacobianInverseTransposed.mv(deformationReferenceDerivative[comp], deformationDerivative[comp]);
+
+        for (size_t comp=0; comp<orientationReferenceDerivative.N(); comp++)
+            jacobianInverseTransposed.mv(orientationReferenceDerivative[comp], orientationDerivative[comp]);
+
+        /////////////////////////////////////////////////////////
+        // compute U, the Cosserat strain
+        /////////////////////////////////////////////////////////
+        static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
+
+        //
+        Dune::FieldMatrix<field_type,dim,dim> R;
+        orientationValue.matrix(R);
+
+        //////////////////////////////////////////////////////////
+        //  Compute the derivative of the rotation
+        //  Note: we need it in matrix coordinates
+        //////////////////////////////////////////////////////////
+
+        Tensor3<field_type,3,3,gridDim> DR = orientationValue.quaternionTangentToMatrixTangent(orientationDerivative);
+
+        // Add the local energy density
+        if (gridDim==2) {
+            energy += weight * thickness_ * membraneEnergyPart1(deformationDerivative,R);
+            energy += weight * thickness_ * membraneEnergyPart2(deformationDerivative,R);
+            energy += weight * thickness_ * curvatureEnergy(DR);
+        } else if (gridDim==3) {
+            DUNE_THROW(Dune::NotImplemented, "CosseratMembraneStiffness for 3d grids");
+            //energy += weight * quadraticMembraneEnergy(U);
+            //energy += weight * curvatureEnergy(DR);
+        } else
+            DUNE_THROW(Dune::NotImplemented, "CosseratMembraneStiffness for 1d grids");
+
+    }
+
+    //////////////////////////////////////////////////////////////////////////////
+    //   Assemble boundary contributions
+    //////////////////////////////////////////////////////////////////////////////
+
+    if (not neumannFunction_)
+        return energy;
+
+    for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
+    {
+        if (not neumannBoundary_ or not neumannBoundary_->contains(it))
+            continue;
+
+        const auto& quad = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
+
+        for (size_t pt=0; pt<quad.size(); pt++) {
+
+            // Local position of the quadrature point
+            const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
+
+            const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
+
+            // The value of the local function
+            RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
+
+            // Value of the Neumann data at the current position
+            auto neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
+
+            // Only translational dofs are affected by the Neumann force
+            for (size_t i=0; i<neumannValue.size(); i++)
+                energy += thickness_ * (neumannValue[i] * deformationValue.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
+
+        }
+
+    }
+
+    return energy;
+}
+
+#endif   //#ifndef COSSERAT_MEMBRANE_STIFFNESS_HH
+