Add membrane-energy and parset files for Patrizios experiment

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
+

problems/circle-09-dirichlet-values.py

+import math
+
+class DirichletValues:
+    def __init__(self, homotopyParameter):
+        self.homotopyParameter = homotopyParameter
+
+    def deformation(self, x):
+        norm = math.sqrt(x[0]*x[0] + x[1]*x[1])
+        return [0.9*x[0]/norm, 0.9*x[1]/norm, 0.0]
+
+    def orientation(self, x):
+        rotation = [[1,0,0], [0,1,0], [0,0,1]]
+        return rotation

problems/circle2ndorder.msh

+$MeshFormat
+2 0 8
+$EndMeshFormat
+$Nodes
+20
+1 -1 0 0
+2 1 0 0
+3 0 0 0
+4 -0.5000000000004422 -0.8660254037841834 0
+5 0.5000000000016158 -0.8660254037835057 0
+6 -0.8660254037845663 -0.4999999999997789 0
+7 6.776189039135044e-13 -1 0
+8 0.8660254037849052 -0.4999999999991919 0
+9 0.5000000000004422 0.8660254037841834 0
+10 -0.5000000000016158 0.8660254037835057 0
+11 0.8660254037845663 0.4999999999997789 0
+12 -6.776189039135044e-13 1 0
+13 -0.8660254037849052 0.4999999999991919 0
+14 -0.02222222222229404 0.03849001794593358 0
+15 -0.261111111111955 0.4522577108647196 0
+16 0.2388888888890741 0.4522577108650584 0
+17 0.488888888888853 0.01924500897296679 0
+18 -0.2611111111113681 -0.4137676929191249 0
+19 -0.511111111111147 0.01924500897296679 0
+20 0.2388888888896609 -0.4137676929187861 0
+$EndNodes
+$Elements
+15
+1 15 2 0 1 1
+2 15 2 0 2 2
+3 15 2 0 3 3
+4 8 3 0 1 0 1 4 6
+5 8 3 0 1 0 4 5 7
+6 8 3 0 1 0 5 2 8
+7 8 3 0 2 0 2 9 11
+8 8 3 0 2 0 9 10 12
+9 8 3 0 2 0 10 1 13
+10 9 3 0 200 0 9 10 14 12 15 16
+11 9 3 0 200 0 14 2 9 17 11 16
+12 9 3 0 200 0 4 14 1 18 19 6
+13 9 3 0 200 0 1 14 10 19 15 13
+14 9 3 0 200 0 2 14 5 17 20 8
+15 9 3 0 200 0 5 14 4 20 18 7
+$EndElements

problems/cosserat-continuum-membrane.parset

+#############################################
+#  Grid parameters
+#############################################
+
+path = /home/lnebel/git_repos/dune/dune-gfe/problems
+gridFile = circle2ndorder.msh
+
+#structuredGrid = cube
+
+# bounding box
+#lower = 0 0
+#upper = 1 1
+
+#elements = 1 1
+
+# Number of grid levels
+numLevels = 5
+
+#############################################
+#  Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-3
+
+# Max number of steps of the trust region solver
+maxSolverSteps = 200
+
+trustRegionScaling = 1 1 1 0.01 0.01 0.01
+
+# Initial trust-region radius
+initialTrustRegionRadius = 3.125
+
+# Number of multigrid iterations per trust-region step
+numIt = 400
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 1
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-5
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+# Measure convergence
+instrumented = 0
+
+############################
+#   Material parameters
+############################
+
+
+## For the Wriggers L-shape example
+[materialParameters]
+
+# shell thickness
+thickness = 1
+
+# Lame parameters
+# corresponds to E = 71240 N/mm^2, nu=0.31
+# However, we use units N/m^2
+mu = 2.7191e+4
+lambda = 4.4364e+4
+
+# Cosserat couple modulus
+mu_c = 2.7191e+4
+
+# Length scale parameter
+L_c = 0.001
+
+# Curvature exponent
+q = 2
+
+# Shear correction factor
+kappa = 1
+
+[]
+
+#############################################
+#  Boundary values
+#############################################
+
+problem = circle-09
+#problem = square-09
+
+###  Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0]*x[0] + x[1]*x[1] >= 0.99"
+#dirichletVerticesPredicate = "x[0] >= 0.99 or x[1] >= 0.99 or x[0] < 0.01 or x[1] < 0.01"
+
+###  Python predicate specifying all Neumann grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "0"
+
+###  Neumann values
+neumannValues = 0 0 0
+
+# Initial deformation
+initialDeformation = "[x[0], x[1], 0]"
+
+#startFromFile = yes
+#initialIterateFilename = initial_iterate.vtu

src/cosserat-continuum.cc

@@ -45,6 +45,7 @@
 #include <dune/gfe/localgeodesicfeadolcstiffness.hh>
 #include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
 #include <dune/gfe/cosseratenergystiffness.hh>
+#include <dune/gfe/cosseratmembranestiffness.hh>
 #include <dune/gfe/nonplanarcosseratshellenergy.hh>
 #include <dune/gfe/cosseratvtkwriter.hh>
 #include <dune/gfe/cosseratvtkreader.hh>
@@ -447,16 +448,32 @@ int main (int argc, char *argv[]) try
             x[_1][i].set(dOV[i]);
 
         if (dim==dimworld) {
-            CosseratEnergyLocalStiffness<CompositeBasis, 3,adouble> localCosseratEnergy(materialParameters,
+            std::shared_ptr<MixedLocalGFEADOLCStiffness<CompositeBasis,
+                                    RealTuple<double,3>,
+                                    Rotation<double,3> >> mixedLocalGFEADOLCStiffness;
+
+            if (materialParameters.template get<bool>("useMembraneEnergy", false) && dim == 3) {
+                DUNE_THROW(Exception, "The Cosserat-Membrane-Energy is only implemented for dim = dimworld = 2.");
+            } else if (materialParameters.template get<bool>("useMembraneEnergy", false) && dim == 2) {
+                CosseratMembraneStiffness<CompositeBasis, 3,adouble> localCosseratMembraneEnergy(materialParameters,
                                                                                             &neumannBoundary,
                                                                                             neumannFunction,
                                                                                             volumeLoad);
-            MixedLocalGFEADOLCStiffness<CompositeBasis,
-                                RealTuple<double,3>,
-                                Rotation<double,3> > localGFEADOLCStiffness(&localCosseratEnergy);
+                mixedLocalGFEADOLCStiffness = std::make_shared<MixedLocalGFEADOLCStiffness<CompositeBasis,
+                                    RealTuple<double,3>,
+                                    Rotation<double,3> >>(&localCosseratMembraneEnergy);
+            } else {
+                CosseratEnergyLocalStiffness<CompositeBasis, 3,adouble> localCosseratEnergy(materialParameters,
+                                                                                            &neumannBoundary,
+                                                                                            neumannFunction,
+                                                                                            volumeLoad);
+                mixedLocalGFEADOLCStiffness = std::make_shared<MixedLocalGFEADOLCStiffness<CompositeBasis,
+                                    RealTuple<double,3>,
+                                    Rotation<double,3> >>(&localCosseratEnergy);
+            }
             MixedGFEAssembler<CompositeBasis,
                       RealTuple<double,3>,
-                      Rotation<double,3> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
+                      Rotation<double,3> > mixedAssembler(compositeBasis, mixedLocalGFEADOLCStiffness.get());
 #if MIXED_SPACE
             MixedRiemannianTrustRegionSolver<GridType,
                                          CompositeBasis,