New problem class for nonlinear elasticity problems

dune/elasticity/common/CMakeLists.txt

 install(FILES
     elasticityhelpers.hh
+    nonlinearelasticityproblem.hh
+    nonlinearelasticityproblem.cc
     DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/elasticity/common)

dune/elasticity/common/Makefile.am

@@ -2,7 +2,7 @@ SUBDIRS =
 
 commondir = $(includedir)/dune/elasticity/common
 
-common_HEADERS =  elasticityhelpers.hh
+common_HEADERS =  elasticityhelpers.hh nonlinearelasticityproblem.hh nonlinearelasticityproblem.cc
 
 include $(top_srcdir)/am/global-rules
 

dune/elasticity/common/nonlinearelasticityproblem.cc

+#include <dune/fufem/assemblers/functionalassembler.hh>
+#include <dune/fufem/assemblers/operatorassembler.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+
+template <class VectorType, class MatrixType, class Basis>
+void NonlinearElasticityProblem<VectorType, MatrixType, Basis>::assembleQP(const VectorType& iterate)
+{
+
+    const Basis& basis = material_->basis();
+
+    GridFunctionPtr displace = std::make_shared<BasisGridFunction<Basis,VectorType> >(basis,iterate);A
+
+    // assemble quadratic term
+    OperatorAssembler<Basis,Basis> globalAssembler(basis,basis);
+    globalAssembler.assemble(material_->secondDerivative(displace), A);
+
+    // assemble linear term
+    FunctionalAssembler<Basis> funcAssembler(basis);
+    funcAssembler.assemble(material_->firstDerivative(displace), f);
+    f -= *extForces_;
+    f *= -1;
+}
+
+template <class VectorType, class MatrixType, class Basis>
+template <class DiagonalMatrixType>
+void NonlinearElasticityProblem<VectorType,MatrixType, Basis>::assembleDefectQP(
+                          const VectorType& iterate,
+                          VectorType& residual,
+                          DiagonalMatrixType& hessian)
+{
+
+    const Basis& basis = material_->basis();
+
+    GridFunctionPtr displace = std::make_shared<BasisGridFunction<Basis,VectorType> >(basis,iterate);
+
+    // assemble quadratic term
+    hessian.setSize(basis.size(),basis.size());
+    hessian = 0;
+
+    const auto& locHessianAssembler =  material_->secondDerivative(displace);
+
+    // ////////////////////////////////////////////////////////////////////////
+    //   Loop over all elements and assemble diagonal entries of the hessian
+    // ////////////////////////////////////////////////////////////////////////
+
+    const typename Basis::GridView& gridView = basis.getGridView();
+    typename Basis::GridView::template Codim<0>::Iterator eIt    = gridView.template begin<0>();
+    typename Basis::GridView::template Codim<0>::Iterator eEndIt = gridView.template end<0>();
+
+    for (; eIt!=eEndIt; ++eIt) {
+
+        // Get set of shape functions
+        const auto& lfe = basis.getLocalFiniteElement(*eIt);
+        const auto& localBasis = lfe.localBasis();
+        typename MaterialType::LocalHessian::LocalMatrix localMatrix(localBasis.size(),localBasis.size());
+        locHessianAssembler.assemble(*eIt,localMatrix, lfe, lfe);
+
+        // Add to residual vector
+        for (size_t j=0; j<localBasis.size(); j++) {
+
+            int globalRow = basis.index(*eIt, j);
+
+            // Assemble diagonal part of second order stiffness matrix
+            hessian[globalRow][globalRow] += localMatrix[j][j];
+
+        }
+    }
+
+    // assemble linear term
+    FunctionalAssembler<Basis> funcAssembler(basis);
+    funcAssembler.assemble(material_->firstDerivative(displace),residual);
+    residual -= *extForces_;
+    residual *= -1;
+}
+
+template <class VectorType, class MatrixType, class Basis>
+typename VectorType::field_type NonlinearElasticityProblem<VectorType,MatrixType, Basis>::energy(const VectorType& iterate) const
+{
+    // The energy is given by the material energy and the external forces
+    field_type energy(0);
+
+    // potential energy
+    GridFunctionPtr displace = std::make_shared<BasisGridFunction<Basis,VectorType> >(material_->basis(),iterate);
+    energy += material_->energy(displace);
+
+    // external terms
+    energy -= (*extForces_)*iterate;
+
+    return energy;
+}
+
+template <class VectorType, class MatrixType, class Basis>
+typename VectorType::field_type NonlinearElasticityProblem<VectorType,MatrixType, Basis>::modelDecrease(const VectorType& correction) const
+{
+    // the model decrease is simply <f,corr> - 0.5 <corr, H corr>
+    VectorType tmp(correction.size());
+    A.mv(correction, tmp);
+
+    return (f*correction) - 0.5 * (correction*tmp);
+}
+
+
+

dune/elasticity/common/nonlinearelasticityproblem.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_ELASTICITY_COMMON_NONLINEAR_ELASTICITY_CONTACT_PROBLEM_HH
+#define DUNE_ELASTICITY_COMMON_NONLINEAR_ELASTICITY_CONTACT_PROBLEM_HH
+
+#include <dune/elasticity/materials/material.hh>
+
+template <class VectorTypeTEMPLATE, class MatrixTypeTEMPLATE, class Basis>
+class NonlinearElasticityProblem
+{
+public:
+    typedef Material<Basis> MaterialType;
+    typedef typename MaterialType::GridType GridType;
+    typedef VectorTypeTEMPLATE VectorType;
+    typedef MatrixTypeTEMPLATE MatrixType;
+    typedef typename VectorType::field_type field_type;
+    typedef std::shared_ptr< <BasisGridFunction<Basis,VectorType> > GridFunctionPtr;
+
+    NonlinearElasticityProblem(const MaterialType& material,
+                         VectorType& extForces):
+        material_(Dune::stackobject_to_shared_ptr(material)),
+        extForces_(Dune::stackobject_to_shared_ptr(extForces))
+    {}
+
+     //! Set external forces
+    void setExternalForces(const VectorType& extForces)
+    {
+        extForces_ = Dune::stackobject_to_shared_ptr(extForces);
+    }
+
+    //! Assemble the quadratic problem
+    void assembleQP(const VectorType& iterate);
+
+    //! Assemble the preconditioned defect problem
+    template <class DiagonalMatrixType>
+    void assembleDefectQP(const VectorType& iterate,
+                          VectorType& residual,
+                          DiagonalMatrixType& hessian);
+
+    //! Compute the energy of the displacement minimization problem
+    field_type energy(const VectorType& iterate) const;
+
+    //! Compute the predicted decrease
+    field_type modelDecrease(const VectorType& correction) const;
+
+    //! The quadratic part
+    MatrixType A;
+    //! The linear functional
+    VectorType f;
+
+private:
+    //! The involved materials
+    std::shared_ptr<MaterialType> material_;
+    //! The external forces
+    std::shared_ptr<VectorType> extForces_;
+};
+
+#include "nonlinearelasticityproblem.cc"
+
+#endif