diff --git a/dune/elasticity/common/CMakeLists.txt b/dune/elasticity/common/CMakeLists.txt index e090d8c6feb90303482645033675c02f79d6f56f..09932dd0daf3bf0d89c4d0ace0b736e2c7f1ffd6 100644 --- a/dune/elasticity/common/CMakeLists.txt +++ b/dune/elasticity/common/CMakeLists.txt @@ -1,3 +1,5 @@ install(FILES elasticityhelpers.hh + nonlinearelasticityproblem.hh + nonlinearelasticityproblem.cc DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/elasticity/common) diff --git a/dune/elasticity/common/Makefile.am b/dune/elasticity/common/Makefile.am index f8fce82d01436da91d062190654b90acd343a7a4..327f3192a9d3956bdaee557778120be99f3036c6 100644 --- a/dune/elasticity/common/Makefile.am +++ b/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 diff --git a/dune/elasticity/common/nonlinearelasticityproblem.cc b/dune/elasticity/common/nonlinearelasticityproblem.cc new file mode 100644 index 0000000000000000000000000000000000000000..39fc4b5f740ffe3345656340906892ad6f218c4a --- /dev/null +++ b/dune/elasticity/common/nonlinearelasticityproblem.cc @@ -0,0 +1,105 @@ +#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); +} + + + diff --git a/dune/elasticity/common/nonlinearelasticityproblem.hh b/dune/elasticity/common/nonlinearelasticityproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..2a80781011b6a40f9577be989e40d6a8224f3581 --- /dev/null +++ b/dune/elasticity/common/nonlinearelasticityproblem.hh @@ -0,0 +1,60 @@ +// -*- 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