From 9dcba1903abc563124215c87d7cf1c1728b3038d Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 15 Apr 2015 15:54:25 +0200
Subject: [PATCH] Move finite-strain-elasticity into the dune-elasticity module
It was only in the dune-gfe module for historic reasons.
---
dune/gfe/exphenckyenergy.hh | 150 -----------
dune/gfe/feassembler.hh | 194 --------------
dune/gfe/henckyenergy.hh | 139 ----------
dune/gfe/localadolcstiffness.hh | 202 --------------
dune/gfe/localfestiffness.hh | 58 ----
dune/gfe/neohookeenergy.hh | 140 ----------
dune/gfe/neumannenergy.hh | 106 --------
dune/gfe/stvenantkirchhoffenergy.hh | 134 ----------
dune/gfe/sumenergy.hh | 60 -----
dune/gfe/trustregionsolver.cc | 396 ----------------------------
dune/gfe/trustregionsolver.hh | 130 ---------
src/CMakeLists.txt | 1 -
src/finite-strain-elasticity.cc | 362 -------------------------
13 files changed, 2072 deletions(-)
delete mode 100644 dune/gfe/exphenckyenergy.hh
delete mode 100644 dune/gfe/feassembler.hh
delete mode 100644 dune/gfe/henckyenergy.hh
delete mode 100644 dune/gfe/localadolcstiffness.hh
delete mode 100644 dune/gfe/localfestiffness.hh
delete mode 100644 dune/gfe/neohookeenergy.hh
delete mode 100644 dune/gfe/neumannenergy.hh
delete mode 100644 dune/gfe/stvenantkirchhoffenergy.hh
delete mode 100644 dune/gfe/sumenergy.hh
delete mode 100644 dune/gfe/trustregionsolver.cc
delete mode 100644 dune/gfe/trustregionsolver.hh
delete mode 100644 src/finite-strain-elasticity.cc
diff --git a/dune/gfe/exphenckyenergy.hh b/dune/gfe/exphenckyenergy.hh
deleted file mode 100644
index 64eecfac..00000000
--- a/dune/gfe/exphenckyenergy.hh
+++ /dev/null
@@ -1,150 +0,0 @@
-#ifndef DUNE_GFE_EXPHENCKYENERGY_HH
-#define DUNE_GFE_EXPHENCKYENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/gfe/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class ExpHenckyEnergy
- : public LocalFEStiffness<GridView,
- LocalFiniteElement,
- std::vector<Dune::FieldVector<field_type,3> > >
-{
- // 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
- */
- ExpHenckyEnergy(const Dune::ParameterTree& parameters)
- {
- // Lame constants
- mu_ = parameters.get<double>("mu");
- lambda_ = parameters.get<double>("lambda");
- kappa_ = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
-
- // Exponential Hencky constants with naming convention from Neff, Ghiba and Lankeit,
- // "The exponentiated Hencky-logarithmic strain energy. Part I: Constitutive issues
- // and rank one convexity"
- k_ = parameters.get<double>("k");
- khat_ = parameters.get<double>("khat");
-
- // weights for distortion and dilatation parts of the energy
- alpha_ = mu_ / k_;
- beta_ = kappa_ / (2.0 * khat_);
- }
-
- /** \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 mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-ExpHenckyEnergy<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;
-
- // 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);
-
- field_type distortionEnergy = 0, dilatationEnergy = 0;
- 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]);
-
- /////////////////////////////////////////////////////////
- // compute C = F^T F
- /////////////////////////////////////////////////////////
-
- 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
- //////////////////////////////////////////////////////////
-
- // compute eigenvalues of C, i.e., squared singular values \sigma_i of F
- Dune::FieldVector<field_type, dim> sigmaSquared;
- FMatrixHelp::eigenValues(C, sigmaSquared);
-
- // logarithms of the singular values of F, i.e., eigenvalues of U
- std::array<field_type, dim> lnSigma;
- for (int i = 0; i < dim; i++)
- lnSigma[i] = 0.5 * log(sigmaSquared[i]);
-
- field_type trLogU = 0;
- for (int i = 0; i < dim; i++)
- trLogU += lnSigma[i];
-
- field_type normDevLogUSquared = 0;
- for (int i = 0; i < dim; i++)
- {
- // the deviator shifts the
- field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
- normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
- }
-
- // Add the local energy density
- distortionEnergy += weight * alpha_ * exp(khat_ * normDevLogUSquared);
- dilatationEnergy += weight * beta_ * exp(k_ * (trLogU * trLogU));
- }
-
- return distortionEnergy + dilatationEnergy;
-}
-
-} // namespace Dune
-
-#endif //#ifndef DUNE_GFE_EXPHENCKYENERGY_HH
-
-
diff --git a/dune/gfe/feassembler.hh b/dune/gfe/feassembler.hh
deleted file mode 100644
index 43644144..00000000
--- a/dune/gfe/feassembler.hh
+++ /dev/null
@@ -1,194 +0,0 @@
-#ifndef DUNE_GFE_FE_ASSEMBLER_HH
-#define DUNE_GFE_FE_ASSEMBLER_HH
-
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrixindexset.hh>
-#include <dune/istl/matrix.hh>
-
-//#include "localgeodesicfestiffness.hh"
-
-
-/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
- */
-template <class Basis, class VectorType>
-class FEAssembler {
-
- typedef typename Basis::GridView GridView;
- typedef typename GridView::template Codim<0>::template Partition<Dune::Interior_Partition>::Iterator ElementIterator;
-
- //! Dimension of the grid.
- enum { gridDim = GridView::dimension };
-
- //! Dimension of a tangent space
- enum { blocksize = VectorType::value_type::dimension };
-
- //!
- typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-
-public:
- const Basis basis_;
-
-protected:
-
- LocalFEStiffness<GridView,
- typename Basis::LocalFiniteElement,
- VectorType>* localStiffness_;
-
-public:
-
- /** \brief Constructor for a given grid */
- FEAssembler(const Basis& basis,
- LocalFEStiffness<GridView,typename Basis::LocalFiniteElement, VectorType>* localStiffness)
- : basis_(basis),
- localStiffness_(localStiffness)
- {}
-
- /** \brief Assemble the tangent stiffness matrix and the functional gradient together
- *
- * This is more efficient than computing them separately, because you need the gradient
- * anyway to compute the Riemannian Hessian.
- */
- virtual void assembleGradientAndHessian(const VectorType& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
- Dune::BCRSMatrix<MatrixBlock>& hessian,
- bool computeOccupationPattern=true) const;
-
- /** \brief Compute the energy of a deformation state */
- virtual double computeEnergy(const VectorType& sol) const;
-
- //protected:
- void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
-
-}; // end class
-
-
-
-template <class Basis, class TargetSpace>
-void FEAssembler<Basis,TargetSpace>::
-getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
-{
- int n = basis_.size();
-
- nb.resize(n, n);
-
- ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
- ElementIterator endit = basis_.getGridView().template end<0,Dune::Interior_Partition> ();
-
- for (; it!=endit; ++it) {
-
- const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(*it);
-
- for (size_t i=0; i<lfe.localBasis().size(); i++) {
-
- for (size_t j=0; j<lfe.localBasis().size(); j++) {
-
- int iIdx = basis_.index(*it,i);
- int jIdx = basis_.index(*it,j);
-
- nb.add(iIdx, jIdx);
-
- }
-
- }
-
- }
-
-}
-
-template <class Basis, class VectorType>
-void FEAssembler<Basis,VectorType>::
-assembleGradientAndHessian(const VectorType& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
- Dune::BCRSMatrix<MatrixBlock>& hessian,
- bool computeOccupationPattern) const
-{
- if (computeOccupationPattern) {
-
- Dune::MatrixIndexSet neighborsPerVertex;
- getNeighborsPerVertex(neighborsPerVertex);
- neighborsPerVertex.exportIdx(hessian);
-
- }
-
- hessian = 0;
-
- gradient.resize(sol.size());
- gradient = 0;
-
- ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
- ElementIterator endit = basis_.getGridView().template end<0,Dune::Interior_Partition> ();
-
- for( ; it != endit; ++it ) {
-
- const int numOfBaseFct = basis_.getLocalFiniteElement(*it).localBasis().size();
-
- // Extract local solution
- VectorType localSolution(numOfBaseFct);
- VectorType localPointLoads(numOfBaseFct);
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[basis_.index(*it,i)];
-
- VectorType localGradient(numOfBaseFct);
-
- // setup local matrix and gradient
- localStiffness_->assembleGradientAndHessian(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
-
- // Add element matrix to global stiffness matrix
- for(int i=0; i<numOfBaseFct; i++) {
-
- int row = basis_.index(*it,i);
-
- for (int j=0; j<numOfBaseFct; j++ ) {
-
- int col = basis_.index(*it,j);
- hessian[row][col] += localStiffness_->A_[i][j];
-
- }
- }
-
- // Add local gradient to global gradient
- for (int i=0; i<numOfBaseFct; i++)
- gradient[basis_.index(*it,i)] += localGradient[i];
-
- }
-
-}
-
-
-template <class Basis, class VectorType>
-double FEAssembler<Basis, VectorType>::
-computeEnergy(const VectorType& sol) const
-{
- double energy = 0;
-
- if (sol.size()!=basis_.size())
- DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
-
- ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
- ElementIterator endIt = basis_.getGridView().template end<0,Dune::Interior_Partition>();
-
- // Loop over all elements
- for (; it!=endIt; ++it) {
-
- // Number of degrees of freedom on this element
- size_t nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
-
- VectorType localSolution(nDofs);
-
- for (size_t i=0; i<nDofs; i++)
- localSolution[i] = sol[basis_.index(*it,i)];
-
- energy += localStiffness_->energy(*it, basis_.getLocalFiniteElement(*it), localSolution);
-
- }
-
- return energy;
-
-}
-
-
-
-#endif
-
diff --git a/dune/gfe/henckyenergy.hh b/dune/gfe/henckyenergy.hh
deleted file mode 100644
index 51a6ca49..00000000
--- a/dune/gfe/henckyenergy.hh
+++ /dev/null
@@ -1,139 +0,0 @@
-#ifndef DUNE_GFE_HENCKYENERGY_HH
-#define DUNE_GFE_HENCKYENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/gfe/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class HenckyEnergy
- : public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > >
-{
- // 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: // for testing
-
- /** \brief Constructor with a set of material parameters
- * \param parameters The material parameters
- */
- HenckyEnergy(const Dune::ParameterTree& parameters)
- {
- // Lame constants
- mu_ = parameters.template get<double>("mu");
- lambda_ = parameters.template get<double>("lambda");
- kappa_ = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
- }
-
- /** \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 */
- double mu_, lambda_, kappa_;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-HenckyEnergy<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 energy = 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);
-
- field_type distortionEnergy = 0, dilatationEnergy = 0;
- 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]);
-
- /////////////////////////////////////////////////////////
- // compute C = F^TF
- /////////////////////////////////////////////////////////
-
- 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 C = F^TF
- //////////////////////////////////////////////////////////
-
- Dune::FieldVector<field_type,dim> sigmaSquared;
- FMatrixHelp::eigenValues(C, sigmaSquared);
-
- // logarithms of the singular values of F, i.e., eigenvalues of U
- std::array<field_type, dim> lnSigma;
- for (int i = 0; i < dim; i++)
- lnSigma[i] = 0.5 * log(sigmaSquared[i]);
-
- field_type trLogU = 0;
- for (int i = 0; i < dim; i++)
- trLogU += lnSigma[i];
-
- field_type normDevLogUSquared = 0;
- for (int i = 0; i < dim; i++)
- {
- // the deviator shifts the spectrum by the trace
- field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
- normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
- }
-
- // Add the local energy density
- distortionEnergy += weight * mu_ * normDevLogUSquared;
- dilatationEnergy += weight * 0.5 * kappa_ * (trLogU * trLogU);
- }
-
- return distortionEnergy + dilatationEnergy;
-}
-
-} // namespace Dune
-
-#endif //#ifndef DUNE_GFE_HENCKYENERGY_HH
-
-
diff --git a/dune/gfe/localadolcstiffness.hh b/dune/gfe/localadolcstiffness.hh
deleted file mode 100644
index ef25d653..00000000
--- a/dune/gfe/localadolcstiffness.hh
+++ /dev/null
@@ -1,202 +0,0 @@
-#ifndef DUNE_GFE_LOCAL_ADOL_C_STIFFNESS_HH
-#define DUNE_GFE_LOCAL_ADOL_C_STIFFNESS_HH
-
-#include <adolc/adouble.h> // use of active doubles
-#include <adolc/drivers/drivers.h> // use of "Easy to Use" drivers
-// gradient(.) and hessian(.)
-#include <adolc/interfaces.h> // use of "Easy to Use" drivers
-#include <adolc/taping.h> // use of taping
-
-#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
-
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrix.hh>
-
-#include <dune/gfe/localfestiffness.hh>
-
-//#define ADOLC_VECTOR_MODE
-
-/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
- */
-template<class GridView, class LocalFiniteElement, class VectorType>
-class LocalADOLCStiffness
- : public LocalFEStiffness<GridView,LocalFiniteElement,VectorType>
-{
- // grid types
- typedef typename GridView::Grid::ctype DT;
- typedef typename VectorType::value_type::field_type RT;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- // some other sizes
- enum {gridDim=GridView::dimension};
-
- // Hack
- typedef std::vector<Dune::FieldVector<adouble,gridDim> > AVectorType;
-
-public:
-
- //! Dimension of a tangent space
- enum { blocksize = VectorType::value_type::dimension };
-
- LocalADOLCStiffness(const LocalFEStiffness<GridView, LocalFiniteElement, AVectorType>* energy)
- : localEnergy_(energy)
- {}
-
- /** \brief Compute the energy at the current configuration */
- virtual RT energy (const Entity& e,
- const LocalFiniteElement& localFiniteElement,
- const VectorType& localConfiguration) const;
-
- /** \brief Assemble the local stiffness matrix at the current position
-
- This uses the automatic differentiation toolbox ADOL_C.
- */
- virtual void assembleGradientAndHessian(const Entity& e,
- const LocalFiniteElement& localFiniteElement,
- const VectorType& localConfiguration,
- VectorType& localGradient);
-
- const LocalFEStiffness<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
-
-};
-
-
-template <class GridView, class LocalFiniteElement, class VectorType>
-typename LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::RT
-LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::
-energy(const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const VectorType& localSolution) const
-{
- double pureEnergy;
-
- std::vector<Dune::FieldVector<adouble,blocksize> > localASolution(localSolution.size());
-
- trace_on(1);
-
- adouble energy = 0;
-
- for (size_t i=0; i<localSolution.size(); i++)
- for (size_t j=0; j<localSolution[i].size(); j++)
- localASolution[i][j] <<= localSolution[i][j];
-
- energy = localEnergy_->energy(element,localFiniteElement,localASolution);
-
- energy >>= pureEnergy;
-
- trace_off();
-
- return pureEnergy;
-}
-
-
-
-// ///////////////////////////////////////////////////////////
-// Compute gradient and Hessian together
-// To compute the Hessian we need to compute the gradient anyway, so we may
-// as well return it. This saves assembly time.
-// ///////////////////////////////////////////////////////////
-template <class GridType, class LocalFiniteElement, class VectorType>
-void LocalADOLCStiffness<GridType, LocalFiniteElement, VectorType>::
-assembleGradientAndHessian(const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const VectorType& localSolution,
- VectorType& localGradient)
-{
- // Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
- energy(element, localFiniteElement, localSolution);
-
- /////////////////////////////////////////////////////////////////
- // Compute the energy gradient
- /////////////////////////////////////////////////////////////////
-
- // Compute the actual gradient
- size_t nDofs = localSolution.size();
- size_t nDoubles = nDofs*blocksize;
- std::vector<double> xp(nDoubles);
- int idx=0;
- for (size_t i=0; i<nDofs; i++)
- for (size_t j=0; j<blocksize; j++)
- xp[idx++] = localSolution[i][j];
-
- // Compute gradient
- std::vector<double> g(nDoubles);
- gradient(1,nDoubles,xp.data(),g.data()); // gradient evaluation
-
- // Copy into Dune type
- localGradient.resize(localSolution.size());
-
- idx=0;
- for (size_t i=0; i<nDofs; i++)
- for (size_t j=0; j<blocksize; j++)
- localGradient[i][j] = g[idx++];
-
- /////////////////////////////////////////////////////////////////
- // Compute Hessian
- /////////////////////////////////////////////////////////////////
-
- this->A_.setSize(nDofs,nDofs);
-
-#ifndef ADOLC_VECTOR_MODE
- std::vector<double> v(nDoubles);
- std::vector<double> w(nDoubles);
-
- std::fill(v.begin(), v.end(), 0.0);
-
- for (size_t i=0; i<nDofs; i++)
- for (size_t ii=0; ii<blocksize; ii++)
- {
- // Evaluate Hessian in the direction of each vector of the orthonormal frame
- for (size_t k=0; k<blocksize; k++)
- v[i*blocksize + k] = (k==ii);
-
- int rc= 3;
- MINDEC(rc, hess_vec(1, nDoubles, xp.data(), v.data(), w.data()));
- if (rc < 0)
- DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
-
- for (size_t j=0; j<nDoubles; j++)
- this->A_[i][j/blocksize][ii][j%blocksize] = w[j];
-
- // Make v the null vector again
- std::fill(&v[i*blocksize], &v[(i+1)*blocksize], 0.0);
- }
-#else
- int n = nDoubles;
- int nDirections = nDofs * blocksize;
- double* tangent[nDoubles];
- for(size_t i=0; i<nDoubles; i++)
- tangent[i] = (double*)malloc(nDirections*sizeof(double));
-
- double* rawHessian[nDoubles];
- for(size_t i=0; i<nDoubles; i++)
- rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
-
- for (int j=0; j<nDirections; j++)
- {
- for (int i=0; i<n; i++)
- tangent[i][j] = 0.0;
-
- for (int i=0; i<embeddedBlocksize; i++)
- tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
- }
-
- hess_mat(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
-
- // Copy Hessian into Dune data type
- for(size_t i=0; i<nDoubles; i++)
- for (size_t j=0; j<nDirections; j++)
- this->A_[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
-
- for(size_t i=0; i<nDoubles; i++) {
- free(rawHessian[i]);
- free(tangent[i]);
- }
-#endif
-
-// std::cout << "ADOL-C stiffness:\n";
-// printmatrix(std::cout, this->A_, "foo", "--");
-}
-
-#endif
-
diff --git a/dune/gfe/localfestiffness.hh b/dune/gfe/localfestiffness.hh
deleted file mode 100644
index 413d0eb2..00000000
--- a/dune/gfe/localfestiffness.hh
+++ /dev/null
@@ -1,58 +0,0 @@
-#ifndef LOCAL_FE_STIFFNESS_HH
-#define LOCAL_FE_STIFFNESS_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrix.hh>
-
-
-template<class GridView, class LocalFiniteElement, class VectorType>
-class LocalFEStiffness
-{
- // grid types
- typedef typename GridView::Grid::ctype DT;
- typedef typename VectorType::value_type::field_type RT;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- // some other sizes
- enum {gridDim=GridView::dimension};
-
-public:
-
- //! Dimension of a tangent space
- enum { blocksize = VectorType::value_type::dimension };
-
- /** \brief Assemble the local stiffness matrix at the current position
-
- This default implementation used finite-difference approximations to compute the second derivatives
- */
- virtual void assembleGradientAndHessian(const Entity& e,
- const LocalFiniteElement& localFiniteElement,
- const VectorType& localConfiguration,
- VectorType& localGradient);
-
- /** \brief Compute the energy at the current configuration */
- virtual RT energy (const Entity& e,
- const LocalFiniteElement& localFiniteElement,
- const VectorType& localConfiguration) const = 0;
-
- // assembled data
- Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
-
-};
-
-
-// ///////////////////////////////////////////////////////////
-// Compute gradient by finite-difference approximation
-// ///////////////////////////////////////////////////////////
-template <class GridType, class LocalFiniteElement, class VectorType>
-void LocalFEStiffness<GridType, LocalFiniteElement, VectorType>::
-assembleGradientAndHessian(const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const VectorType& localConfiguration,
- VectorType& localGradient)
-{
- DUNE_THROW(Dune::NotImplemented, "!");
-}
-
-#endif
-
diff --git a/dune/gfe/neohookeenergy.hh b/dune/gfe/neohookeenergy.hh
deleted file mode 100644
index 975935da..00000000
--- a/dune/gfe/neohookeenergy.hh
+++ /dev/null
@@ -1,140 +0,0 @@
-#ifndef DUNE_GFE_NEOHOOKEENERGY_HH
-#define DUNE_GFE_NEOHOOKEENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/gfe/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeoHookeEnergy
- : public LocalFEStiffness<GridView,
- LocalFiniteElement,
- std::vector<Dune::FieldVector<field_type,3> > >
-{
- // 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
- */
- NeoHookeEnergy(const Dune::ParameterTree& parameters)
- {
- // Lame constants
- mu_ = parameters.get<double>("mu");
- lambda_ = parameters.get<double>("lambda");
- kappa_ = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
- }
-
- /** \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 mu_, lambda_, kappa_;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeoHookeEnergy<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 distortionEnergy = 0, dilatationEnergy = 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]);
-
- /////////////////////////////////////////////////////////
- // compute C = F^T F
- /////////////////////////////////////////////////////////
-
- 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);
-
- // singular values of F, i.e., eigenvalues of U
- std::array<field_type, dim> sigma;
- for (int i = 0; i < dim; i++)
- sigma[i] = std::sqrt(sigmaSquared[i]);
-
- field_type detC = 1.0;
- for (int i = 0; i < dim; i++)
- detC *= sigmaSquared[i];
- field_type detF = std::sqrt(detC);
-
- // \tilde{C} = \tilde{F}^T\tilde{F} = \frac{1}{\det{F}^{2/3}}C
- field_type trCTilde = 0;
- for (int i = 0; i < dim; i++)
- trCTilde += sigmaSquared[i];
- trCTilde /= pow(detC, 1.0/3.0);
-
- // Add the local energy density
- distortionEnergy += weight * (0.5 * mu_) * (trCTilde - 3);
- dilatationEnergy += weight * (0.5 * kappa_) * ((detF - 1) * (detF - 1));
- }
-
- return distortionEnergy + dilatationEnergy;
-}
-
-} // namespace Dune
-
-#endif //#ifndef DUNE_GFE_NEOHOOKEENERGY_HH
-
-
diff --git a/dune/gfe/neumannenergy.hh b/dune/gfe/neumannenergy.hh
deleted file mode 100644
index 72b095e0..00000000
--- a/dune/gfe/neumannenergy.hh
+++ /dev/null
@@ -1,106 +0,0 @@
-#ifndef DUNE_GFE_NEUMANNENERGY_HH
-#define DUNE_GFE_NEUMANNENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/fufem/functions/virtualgridfunction.hh>
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/gfe/localgeodesicfestiffness.hh>
-#include <dune/gfe/localfestiffness.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/realtuple.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeumannEnergy
-: public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > >
-{
- // grid types
- typedef typename GridView::ctype ctype;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- enum {dim=GridView::dimension};
-
-public:
-
- /** \brief Constructor with a set of material parameters
- * \param parameters The material parameters
- */
- NeumannEnergy(const BoundaryPatch<GridView>* neumannBoundary,
- const Dune::VirtualFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,3> >* neumannFunction)
- : neumannBoundary_(neumannBoundary),
- neumannFunction_(neumannFunction)
- {}
-
- /** \brief Assemble the energy for a single element */
- field_type energy (const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const std::vector<Dune::FieldVector<field_type,dim> >& localConfiguration) const
- {
- assert(element.type() == localFiniteElement.type());
-
- field_type energy = 0;
-
- for (auto&& it : intersections(neumannBoundary_->gridView(), element)) {
-
- if (not neumannBoundary_ or not neumannBoundary_->contains(it))
- continue;
-
- int quadOrder = localFiniteElement.localBasis().order();
-
- const Dune::QuadratureRule<ctype, dim-1>& quad
- = Dune::QuadratureRules<ctype, dim-1>::rule(it.type(), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<ctype,dim>& quadPos = it.geometryInInside().global(quad[pt].position());
-
- const auto integrationElement = it.geometry().integrationElement(quad[pt].position());
-
- // The value of the local function
- std::vector<Dune::FieldVector<ctype,1> > shapeFunctionValues;
- localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
-
- Dune::FieldVector<field_type,dim> value(0);
- for (size_t i=0; i<localFiniteElement.size(); i++)
- for (int j=0; j<dim; j++)
- value[j] += shapeFunctionValues[i] * localConfiguration[i][j];
-
- // Value of the Neumann data at the current position
- Dune::FieldVector<double,3> neumannValue;
-
- if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_))
- dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
- else
- neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
-
- // Only translational dofs are affected by the Neumann force
- for (size_t i=0; i<neumannValue.size(); i++)
- energy += (neumannValue[i] * value[i]) * quad[pt].weight() * integrationElement;
-
- }
-
- }
-
- return energy;
- }
-
-private:
- /** \brief The Neumann boundary */
- const BoundaryPatch<GridView>* neumannBoundary_;
-
- /** \brief The function implementing the Neumann data */
- const Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,3> >* neumannFunction_;
-};
-
-}
-
-#endif //#ifndef DUNE_GFE_NEUMANNENERGY_HH
-
-
diff --git a/dune/gfe/stvenantkirchhoffenergy.hh b/dune/gfe/stvenantkirchhoffenergy.hh
deleted file mode 100644
index 96b95c20..00000000
--- a/dune/gfe/stvenantkirchhoffenergy.hh
+++ /dev/null
@@ -1,134 +0,0 @@
-#ifndef DUNE_GFE_STVENANTKIRCHHOFFENERGY_HH
-#define DUNE_GFE_STVENANTKIRCHHOFFENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/gfe/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class StVenantKirchhoffEnergy
- : 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
- */
- StVenantKirchhoffEnergy(const Dune::ParameterTree& parameters)
- {
- // Lame constants
- mu_ = parameters.template get<double>("mu");
- lambda_ = parameters.template get<double>("lambda");
- }
-
- /** \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 */
- double mu_, lambda_;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-StVenantKirchhoffEnergy<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 energy = 0;
-
- // store gradients of shape functions and base functions
- std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.localBasis().size());
- std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.localBasis().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) {
-
- // transform gradients
- 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]);
-
- /////////////////////////////////////////////////////////
- // compute strain E = 1/2 *( F^T F - I)
- /////////////////////////////////////////////////////////
-
- Dune::FieldMatrix<field_type,gridDim,gridDim> FTF(0);
- for (int i=0; i<gridDim; i++)
- for (int j=0; j<gridDim; j++)
- for (int k=0; k<gridDim; k++)
- FTF[i][j] += derivative[k][i] * derivative[k][j];
-
- Dune::FieldMatrix<field_type,dim,dim> E = FTF;
- for (int i=0; i<dim; i++)
- E[i][i] -= 1.0;
- E *= 0.5;
-
- /////////////////////////////////////////////////////////
- // Compute energy
- /////////////////////////////////////////////////////////
-
- field_type trE = field_type(0);
- for (int i=0; i<dim; i++)
- trE += E[i][i];
-
- // TODO Wasteful, we only need the trace, not the full product
- Dune::FieldMatrix<field_type,dim,dim> ESquared = E*E;
-
- field_type trESquared = field_type(0);
- for (int i=0; i<dim; i++)
- trESquared += ESquared[i][i];
-
- energy += weight * mu_ * trESquared + weight * 0.5 * lambda_ * trE * trE;
-
- }
-
- return energy;
-}
-
-} // namespace Dune
-
-#endif //#ifndef DUNE_GFE_STVENANTKIRCHHOFFENERGY_HH
-
-
diff --git a/dune/gfe/sumenergy.hh b/dune/gfe/sumenergy.hh
deleted file mode 100644
index ed7e6516..00000000
--- a/dune/gfe/sumenergy.hh
+++ /dev/null
@@ -1,60 +0,0 @@
-#ifndef DUNE_GFE_SUMENERGY_HH
-#define DUNE_GFE_SUMENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/fufem/functions/virtualgridfunction.hh>
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/gfe/localgeodesicfestiffness.hh>
-#include <dune/gfe/localfestiffness.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/realtuple.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class SumEnergy
-: public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > >
-{
- // grid types
- typedef typename GridView::ctype ctype;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- enum {dim=GridView::dimension};
-
-public:
-
- /** \brief Constructor with a set of material parameters
- * \param parameters The material parameters
- */
- SumEnergy(std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > > a,
- std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > > b)
- : a_(a),
- b_(b)
- {}
-
- /** \brief Assemble the energy for a single element */
- field_type energy (const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const std::vector<Dune::FieldVector<field_type,dim> >& localConfiguration) const
- {
- return a_->energy(element, localFiniteElement, localConfiguration)
- + b_->energy(element, localFiniteElement, localConfiguration);
- }
-
-private:
-
- std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > > a_;
-
- std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > > b_;
-};
-
-}
-
-#endif //#ifndef DUNE_GFE_SUMENERGY_HH
-
-
diff --git a/dune/gfe/trustregionsolver.cc b/dune/gfe/trustregionsolver.cc
deleted file mode 100644
index 16326368..00000000
--- a/dune/gfe/trustregionsolver.cc
+++ /dev/null
@@ -1,396 +0,0 @@
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/timer.hh>
-
-#include <dune/istl/io.hh>
-
-#include <dune/functions/functionspacebases/pq1nodalbasis.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
-
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/assemblers/operatorassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
-#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
-
-// Using a monotone multigrid as the inner solver
-#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
-#include <dune/solvers/iterationsteps/mmgstep.hh>
-#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
-#if defined THIRD_ORDER || defined SECOND_ORDER
-#include <dune/gfe/pktop1mgtransfer.hh>
-#endif
-#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include "maxnormtrustregion.hh"
-
-#include <dune/solvers/norms/twonorm.hh>
-#include <dune/solvers/norms/h1seminorm.hh>
-
-template <class BasisType, class VectorType>
-void TrustRegionSolver<BasisType,VectorType>::
-setup(const typename BasisType::GridView::Grid& grid,
- const FEAssembler<BasisType, VectorType>* assembler,
- const SolutionType& x,
- const Dune::BitSetVector<blocksize>& dirichletNodes,
- double tolerance,
- int maxTrustRegionSteps,
- double initialTrustRegionRadius,
- int multigridIterations,
- double mgTolerance,
- int mu,
- int nu1,
- int nu2,
- int baseIterations,
- double baseTolerance)
-{
- grid_ = &grid;
- assembler_ = assembler;
- x_ = x;
- this->tolerance_ = tolerance;
- maxTrustRegionSteps_ = maxTrustRegionSteps;
- initialTrustRegionRadius_ = initialTrustRegionRadius;
- innerIterations_ = multigridIterations;
- innerTolerance_ = mgTolerance;
- ignoreNodes_ = &dirichletNodes;
-
- int numLevels = grid_->maxLevel()+1;
-
- // ////////////////////////////////
- // Create a multigrid solver
- // ////////////////////////////////
-
-#ifdef HAVE_IPOPT
- // First create an IPOpt base solver
- QuadraticIPOptSolver<MatrixType, CorrectionType>* baseSolver = new QuadraticIPOptSolver<MatrixType,CorrectionType>;
- baseSolver->verbosity_ = NumProc::QUIET;
- baseSolver->tolerance_ = baseTolerance;
- baseSolver->linearSolverType_ = "mumps";
-#else
- // First create a Gauss-seidel base solver
- TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
-
- // Hack: the two-norm may not scale all that well, but it is fast!
- TwoNorm<CorrectionType>* baseNorm = new TwoNorm<CorrectionType>;
-
- ::LoopSolver<CorrectionType>* baseSolver = new ::LoopSolver<CorrectionType>(baseSolverStep,
- baseIterations,
- baseTolerance,
- baseNorm,
- Solver::QUIET);
-#endif
-
- // Make pre and postsmoothers
- TrustRegionGSStep<MatrixType, CorrectionType>* presmoother = new TrustRegionGSStep<MatrixType, CorrectionType>;
- TrustRegionGSStep<MatrixType, CorrectionType>* postsmoother = new TrustRegionGSStep<MatrixType, CorrectionType>;
-
- MonotoneMGStep<MatrixType, CorrectionType>* mmgStep = new MonotoneMGStep<MatrixType, CorrectionType>;
-
- mmgStep->setMGType(mu, nu1, nu2);
- mmgStep->ignoreNodes_ = &dirichletNodes;
- mmgStep->basesolver_ = baseSolver;
- mmgStep->setSmoother(presmoother, postsmoother);
- mmgStep->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType>();
- mmgStep->verbosity_ = Solver::QUIET;
-
- // //////////////////////////////////////////////////////////////////////////////////////
- // Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
- // //////////////////////////////////////////////////////////////////////////////////////
-
- BasisType basis(grid.leafGridView());
- OperatorAssembler<BasisType,BasisType> operatorAssembler(basis, basis);
-
- LaplaceAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> laplaceStiffness;
- typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
- ScalarMatrixType localA;
-
- operatorAssembler.assemble(laplaceStiffness, localA);
-
- if (h1SemiNorm_)
- delete h1SemiNorm_;
-
- ScalarMatrixType* A = new ScalarMatrixType(localA);
-
- h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
-
- innerSolver_ = std::shared_ptr<LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(mmgStep,
- innerIterations_,
- innerTolerance_,
- h1SemiNorm_,
- Solver::REDUCED));
-
- // //////////////////////////////////////////////////////////////////////////////////////
- // Assemble a mass matrix to create a norm that's equivalent to the L2-norm
- // This will be used to monitor the gradient
- // //////////////////////////////////////////////////////////////////////////////////////
-
- MassAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> massStiffness;
- ScalarMatrixType localMassMatrix;
-
- operatorAssembler.assemble(massStiffness, localMassMatrix);
-
- ScalarMatrixType* massMatrix = new ScalarMatrixType(localMassMatrix);
- l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
-
- // ////////////////////////////////////////////////////////////
- // Create Hessian matrix and its occupation structure
- // ////////////////////////////////////////////////////////////
-
- hessianMatrix_ = std::auto_ptr<MatrixType>(new MatrixType);
- Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
- assembler_->getNeighborsPerVertex(indices);
- indices.exportIdx(*hessianMatrix_);
-
- // ////////////////////////////////////
- // Create the transfer operators
- // ////////////////////////////////////
-
- for (size_t k=0; k<mmgStep->mgTransfer_.size(); k++)
- delete(mmgStep->mgTransfer_[k]);
-
- ////////////////////////////////////////////////////////////////////////
- // The P1 space (actually P1/Q1, depending on the grid) is special:
- // If we work in such a space, then the multigrid hierarchy of spaces
- // is constructed in the usual way. For all other space, there is
- // an additional restriction operator on the top of the hierarchy, which
- // restricts the FE space to the P1/Q1 space on the same grid.
- // On the lower grid levels a hierarchy of P1/Q1 spaces is used again.
- ////////////////////////////////////////////////////////////////////////
-
- typedef BasisType Basis;
- bool isP1Basis = std::is_same<Basis,DuneFunctionsBasis<Dune::Functions::PQ1NodalBasis<typename Basis::GridView> > >::value
- || std::is_same<Basis,DuneFunctionsBasis<Dune::Functions::PQKNodalBasis<typename Basis::GridView, 1> > >::value;
-
- if (isP1Basis)
- mmgStep->mgTransfer_.resize(numLevels-1);
- else
- mmgStep->mgTransfer_.resize(numLevels);
-
- // Here we set up the restriction of the leaf grid space into the leaf grid P1/Q1 space
- if (not isP1Basis)
- {
- typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
- P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
-
- TransferOperatorType pkToP1TransferMatrix;
- assembleBasisInterpolationMatrix<TransferOperatorType,
- P1NodalBasis<typename GridType::LeafGridView,double>,
- BasisType>(pkToP1TransferMatrix,p1Basis,basis);
- mmgStep->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType>;
- Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
- dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType>*>(mmgStep->mgTransfer_.back())->setMatrix(topTransferOperator);
- }
-
- // Now the P1/Q1 restriction operators for the remaining levels
- for (int i=0; i<numLevels-1; i++) {
-
- // Construct the local multigrid transfer matrix
- TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
- newTransferOp->setup(*grid_,i,i+1);
-
- typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
- mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>;
- std::shared_ptr<TransferOperatorType> transferOperatorMatrix = Dune::make_shared<TransferOperatorType>(newTransferOp->getMatrix());
- dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType>*>(mmgStep->mgTransfer_[i])->setMatrix(transferOperatorMatrix);
- }
-
- // //////////////////////////////////////////////////////////
- // Create obstacles
- // //////////////////////////////////////////////////////////
-
- hasObstacle_.resize(basis.size(), true);
- mmgStep->hasObstacle_ = &hasObstacle_;
-
-}
-
-
-template <class BasisType, class VectorType>
-void TrustRegionSolver<BasisType,VectorType>::solve()
-{
- MonotoneMGStep<MatrixType,CorrectionType>* mgStep = NULL;
-
- // if the inner solver is a monotone multigrid set up a max-norm trust-region
- if (dynamic_cast<LoopSolver<CorrectionType>*>(innerSolver_.get())) {
- mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(dynamic_cast<LoopSolver<CorrectionType>*>(innerSolver_.get())->iterationStep_);
-
- }
-
- BasisType basis(grid_->leafGridView());
- MaxNormTrustRegion<blocksize> trustRegion(basis.size(), initialTrustRegionRadius_);
-
- std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles;
-
- // /////////////////////////////////////////////////////
- // Trust-Region Solver
- // /////////////////////////////////////////////////////
-
- double oldEnergy = assembler_->computeEnergy(x_);
-
- bool recomputeGradientHessian = true;
- CorrectionType rhs;
- MatrixType stiffnessMatrix;
-
- for (int i=0; i<maxTrustRegionSteps_; i++) {
-
- Dune::Timer totalTimer;
- if (this->verbosity_ == Solver::FULL) {
- std::cout << "----------------------------------------------------" << std::endl;
- std::cout << " Trust-Region Step Number: " << i
- << ", radius: " << trustRegion.radius()
- << ", energy: " << oldEnergy << std::endl;
- std::cout << "----------------------------------------------------" << std::endl;
- }
-
- Dune::Timer gradientTimer;
-
- if (recomputeGradientHessian) {
-
- assembler_->assembleGradientAndHessian(x_,
- rhs,
- *hessianMatrix_,
- i==0 // assemble occupation pattern only for the first call
- );
-
- rhs *= -1; // The right hand side is the _negative_ gradient
-
- // Compute gradient norm to monitor convergence
- CorrectionType gradient = rhs;
- for (size_t j=0; j<gradient.size(); j++)
- for (size_t k=0; k<gradient[j].size(); k++)
- if ((*ignoreNodes_)[j][k])
- gradient[j][k] = 0;
-
- if (this->verbosity_ == Solver::FULL)
- std::cout << "Gradient norm: " << l2Norm_->operator()(gradient) << std::endl;
-
- if (this->verbosity_ == Solver::FULL)
- std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
-
- // Transfer matrix data
- stiffnessMatrix = *hessianMatrix_;
-
- recomputeGradientHessian = false;
-
- }
-
- CorrectionType corr(rhs.size());
- corr = 0;
-
- mgStep->setProblem(stiffnessMatrix, corr, rhs);
-
- trustRegionObstacles = trustRegion.obstacles();
- mgStep->obstacles_ = &trustRegionObstacles;
-
- innerSolver_->preprocess();
-
- ///////////////////////////////
- // Solve !
- ///////////////////////////////
-
- std::cout << "Solve quadratic problem..." << std::endl;
-
- Dune::Timer solutionTimer;
- innerSolver_->solve();
- std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
-
- if (mgStep)
- corr = mgStep->getSol();
-
- //std::cout << "Correction: " << std::endl << corr << std::endl;
-
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
-
- if (corr.infinity_norm() < this->tolerance_) {
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
-
- if (this->verbosity_ != NumProc::QUIET)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
- break;
- }
-
- // ////////////////////////////////////////////////////
- // Check whether trust-region step can be accepted
- // ////////////////////////////////////////////////////
-
- SolutionType newIterate = x_;
- for (size_t j=0; j<newIterate.size(); j++)
- newIterate[j] += corr[j];
-
- double energy = assembler_->computeEnergy(newIterate);
-
- // compute the model decrease
- // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
- // Note that rhs = -g
- CorrectionType tmp(corr.size());
- tmp = 0;
- hessianMatrix_->umv(corr, tmp);
- double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
-
- double relativeModelDecrease = modelDecrease / std::fabs(energy);
-
- if (this->verbosity_ == NumProc::FULL) {
- std::cout << "Absolute model decrease: " << modelDecrease
- << ", functional decrease: " << oldEnergy - energy << std::endl;
- std::cout << "Relative model decrease: " << relativeModelDecrease
- << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
- }
-
- assert(modelDecrease >= 0);
-
- if (energy >= oldEnergy) {
- if (this->verbosity_ == NumProc::FULL)
- printf("Richtung ist keine Abstiegsrichtung!\n");
- }
-
- if (energy >= oldEnergy &&
- (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Suspecting rounding problems" << std::endl;
-
- if (this->verbosity_ != NumProc::QUIET)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
-
- x_ = newIterate;
- break;
- }
-
- // //////////////////////////////////////////////
- // Check for acceptance of the step
- // //////////////////////////////////////////////
- if ( (oldEnergy-energy) / modelDecrease > 0.9) {
- // very successful iteration
-
- x_ = newIterate;
- trustRegion.scale(2);
-
- // current energy becomes 'oldEnergy' for the next iteration
- oldEnergy = energy;
-
- recomputeGradientHessian = true;
-
- } else if ( (oldEnergy-energy) / modelDecrease > 0.01
- || std::abs(oldEnergy-energy) < 1e-12) {
- // successful iteration
- x_ = newIterate;
-
- // current energy becomes 'oldEnergy' for the next iteration
- oldEnergy = energy;
-
- recomputeGradientHessian = true;
-
- } else {
-
- // unsuccessful iteration
-
- // Decrease the trust-region radius
- trustRegion.scale(0.5);
-
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Unsuccessful iteration!" << std::endl;
- }
-
- std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
- }
-
-}
diff --git a/dune/gfe/trustregionsolver.hh b/dune/gfe/trustregionsolver.hh
deleted file mode 100644
index e371d7e2..00000000
--- a/dune/gfe/trustregionsolver.hh
+++ /dev/null
@@ -1,130 +0,0 @@
-#ifndef DUNE_GFE_TRUST_REGION_SOLVER_HH
-#define DUNE_GFE_TRUST_REGION_SOLVER_HH
-
-#include <vector>
-
-#include <dune/common/bitsetvector.hh>
-
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/solvers/common/boxconstraint.hh>
-#include <dune/solvers/norms/h1seminorm.hh>
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
-#include <dune/fufem/functionspacebases/p3nodalbasis.hh>
-
-#include <dune/gfe/feassembler.hh>
-
-/** \brief Trust-region solver */
-template <class BasisType, class VectorType>
-class TrustRegionSolver
- : public IterativeSolver<VectorType,
- Dune::BitSetVector<VectorType::value_type::dimension> >
-{
- typedef typename BasisType::GridView::Grid GridType;
-
- const static int blocksize = VectorType::value_type::dimension;
-
- const static int gridDim = GridType::dimension;
-
- // Centralize the field type here
- typedef double field_type;
-
- // Some types that I need
- typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
- typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
- typedef VectorType SolutionType;
-
-public:
-
- TrustRegionSolver()
- : IterativeSolver<VectorType, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
- hessianMatrix_(std::auto_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL)
- {}
-
- /** \brief Set up the solver using a monotone multigrid method as the inner solver */
- void setup(const GridType& grid,
- const FEAssembler<BasisType,VectorType>* assembler,
- const SolutionType& x,
- const Dune::BitSetVector<blocksize>& dirichletNodes,
- double tolerance,
- int maxTrustRegionSteps,
- double initialTrustRegionRadius,
- int multigridIterations,
- double mgTolerance,
- int mu,
- int nu1,
- int nu2,
- int baseIterations,
- double baseTolerance);
-
- void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
- {
- ignoreNodes_ = &ignoreNodes;
- Dune::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
- assert(loopSolver);
- loopSolver->iterationStep_->ignoreNodes_ = ignoreNodes_;
- }
-
- void solve();
-
- void setInitialSolution(const SolutionType& x) DUNE_DEPRECATED {
- x_ = x;
- }
-
- void setInitialIterate(const SolutionType& x) {
- x_ = x;
- }
-
- SolutionType getSol() const {return x_;}
-
-protected:
-
- /** \brief The grid */
- const GridType* grid_;
-
- /** \brief The solution vector */
- SolutionType x_;
-
- /** \brief The initial trust-region radius in the maximum-norm */
- double initialTrustRegionRadius_;
-
- /** \brief Maximum number of trust-region steps */
- int maxTrustRegionSteps_;
-
- /** \brief Maximum number of multigrid iterations */
- int innerIterations_;
-
- /** \brief Error tolerance of the multigrid QP solver */
- double innerTolerance_;
-
- /** \brief Hessian matrix */
- std::auto_ptr<MatrixType> hessianMatrix_;
-
- /** \brief The assembler for the material law */
- const FEAssembler<BasisType, VectorType>* assembler_;
-
- /** \brief The solver for the quadratic inner problems */
- std::shared_ptr<Solver> innerSolver_;
-
- /** \brief Contains 'true' everywhere -- the trust-region is bounded */
- Dune::BitSetVector<blocksize> hasObstacle_;
-
- /** \brief The Dirichlet nodes */
- const Dune::BitSetVector<blocksize>* ignoreNodes_;
-
- /** \brief The norm used to measure multigrid convergence */
- H1SemiNorm<CorrectionType>* h1SemiNorm_;
-
- /** \brief An L2-norm, really. The H1SemiNorm class is badly named */
- std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
-
-};
-
-#include "trustregionsolver.cc"
-
-#endif
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index c029adc8..d74091c0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,5 +1,4 @@
set(programs cosserat-continuum
- finite-strain-elasticity
harmonicmaps
mixed-cosserat-continuum
rod-eoc
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
deleted file mode 100644
index 1bd5fd14..00000000
--- a/src/finite-strain-elasticity.cc
+++ /dev/null
@@ -1,362 +0,0 @@
-#include <config.h>
-
-// Includes for the ADOL-C automatic differentiation library
-// Need to come before (almost) all others.
-#include <adolc/adouble.h>
-#include <adolc/drivers/drivers.h> // use of "Easy to Use" drivers
-#include <adolc/taping.h>
-
-#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/parametertreeparser.hh>
-
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-
-#include <dune/grid/io/file/gmshreader.hh>
-#include <dune/grid/io/file/vtk.hh>
-
-#include <dune/functions/functionspacebases/pq2nodalbasis.hh>
-#include <dune/functions/functionspacebases/interpolate.hh>
-#include <dune/functions/gridfunctions/discretescalarglobalbasisfunction.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
-#include <dune/fufem/dunepython.hh>
-
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include <dune/solvers/norms/energynorm.hh>
-
-#include <dune/gfe/localadolcstiffness.hh>
-#include <dune/gfe/stvenantkirchhoffenergy.hh>
-#include <dune/gfe/henckyenergy.hh>
-#include <dune/gfe/exphenckyenergy.hh>
-#include <dune/gfe/neohookeenergy.hh>
-#include <dune/gfe/neumannenergy.hh>
-#include <dune/gfe/sumenergy.hh>
-#include <dune/gfe/feassembler.hh>
-#include <dune/gfe/trustregionsolver.hh>
-
-// grid dimension
-const int dim = 3;
-
-using namespace Dune;
-
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
- : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
-{
- NeumannFunction(const FieldVector<double,3> values,
- double homotopyParameter)
- : values_(values),
- homotopyParameter_(homotopyParameter)
- {}
-
- void evaluate(const FieldVector<double, dim>& x, FieldVector<double,3>& out) const
- {
- out = 0;
- out.axpy(-homotopyParameter_, values_);
- }
-
- FieldVector<double,3> values_;
- double homotopyParameter_;
-};
-
-
-int main (int argc, char *argv[]) try
-{
- // initialize MPI, finalize is done automatically on exit
- Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
-
- // Start Python interpreter
- Python::start();
- Python::Reference main = Python::import("__main__");
- Python::run("import math");
-
- //feenableexcept(FE_INVALID);
- Python::runStream()
- << std::endl << "import sys"
- << std::endl << "sys.path.append('/home/sander/dune/dune-gfe/')"
- << std::endl;
-
- typedef BlockVector<FieldVector<double,dim> > SolutionType;
-
- // parse data file
- ParameterTree parameterSet;
-
- ParameterTreeParser::readINITree(argv[1], parameterSet);
-
- ParameterTreeParser::readOptions(argc, argv, parameterSet);
-
- // read solver settings
- const int numLevels = parameterSet.get<int>("numLevels");
- int numHomotopySteps = parameterSet.get<int>("numHomotopySteps");
- const double tolerance = parameterSet.get<double>("tolerance");
- const int maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
- const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
- const int multigridIterations = parameterSet.get<int>("numIt");
- const int nu1 = parameterSet.get<int>("nu1");
- const int nu2 = parameterSet.get<int>("nu2");
- const int mu = parameterSet.get<int>("mu");
- const int baseIterations = parameterSet.get<int>("baseIt");
- const double mgTolerance = parameterSet.get<double>("mgTolerance");
- const double baseTolerance = parameterSet.get<double>("baseTolerance");
- std::string resultPath = parameterSet.get("resultPath", "");
-
- // ///////////////////////////////////////
- // Create the grid
- // ///////////////////////////////////////
- typedef UGGrid<dim> GridType;
-
- shared_ptr<GridType> grid;
-
- FieldVector<double,dim> lower(0), upper(1);
-
- if (parameterSet.get<bool>("structuredGrid")) {
-
- lower = parameterSet.get<FieldVector<double,dim> >("lower");
- upper = parameterSet.get<FieldVector<double,dim> >("upper");
-
- array<unsigned int,dim> elements = parameterSet.get<array<unsigned int,dim> >("elements");
- grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
-
- } else {
- std::string path = parameterSet.get<std::string>("path");
- std::string gridFile = parameterSet.get<std::string>("gridFile");
- grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
- }
-
- grid->globalRefine(numLevels-1);
-
- grid->loadBalance();
-
- if (mpiHelper.rank()==0)
- std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
-
- typedef GridType::LeafGridView GridView;
- GridView gridView = grid->leafGridView();
-
- // FE basis spanning the FE space that we are working in
- typedef Dune::Functions::PQ2NodalBasis<GridView> FEBasis;
- FEBasis feBasis(gridView);
-
- // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
- typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
- FufemFEBasis fufemFEBasis(feBasis);
-
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
-
- BitSetVector<1> dirichletVertices(gridView.size(dim), false);
- BitSetVector<1> neumannVertices(gridView.size(dim), false);
-
- const GridView::IndexSet& indexSet = gridView.indexSet();
-
- // Make Python function that computes which vertices are on the Dirichlet boundary,
- // based on the vertex positions.
- std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
-
- // Same for the Neumann boundary
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
-
- for (auto&& v : vertices(gridView))
- {
- bool isDirichlet;
- pythonDirichletVertices.evaluate(v.geometry().corner(0), isDirichlet);
- dirichletVertices[indexSet.index(v)] = isDirichlet;
-
- bool isNeumann;
- pythonNeumannVertices.evaluate(v.geometry().corner(0), isNeumann);
- neumannVertices[indexSet.index(v)] = isNeumann;
- }
-
- BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
-
- if (mpiHelper.rank()==0)
- std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
-
-
- BitSetVector<1> dirichletNodes(feBasis.indexSet().size(), false);
- constructBoundaryDofs(dirichletBoundary,fufemFEBasis,dirichletNodes);
-
- BitSetVector<1> neumannNodes(feBasis.indexSet().size(), false);
- constructBoundaryDofs(neumannBoundary,fufemFEBasis,neumannNodes);
-
- BitSetVector<dim> dirichletDofs(feBasis.indexSet().size(), false);
- for (size_t i=0; i<feBasis.indexSet().size(); i++)
- if (dirichletNodes[i][0])
- for (int j=0; j<dim; j++)
- dirichletDofs[i][j] = true;
-
- // //////////////////////////
- // Initial iterate
- // //////////////////////////
-
- SolutionType x(feBasis.indexSet().size());
-
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
-
- ::Functions::interpolate(fufemFEBasis, x, pythonInitialDeformation);
-
- ////////////////////////////////////////////////////////
- // Main homotopy loop
- ////////////////////////////////////////////////////////
-
- // Output initial iterate (of homotopy loop)
- SolutionType identity;
- Dune::Functions::interpolate(feBasis, identity, [](FieldVector<double,dim> x){ return x; });
-
- SolutionType displacement = x;
- displacement -= identity;
-
- Dune::Functions::DiscreteScalarGlobalBasisFunction<FEBasis,SolutionType> displacementFunction(feBasis,displacement);
- auto localDisplacementFunction = localFunction(displacementFunction);
-
- // We need to subsample, because VTK cannot natively display real second-order functions
- SubsamplingVTKWriter<GridView> vtkWriter(gridView,2);
- vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
- vtkWriter.write(resultPath + "finite-strain_homotopy_0");
-
- for (int i=0; i<numHomotopySteps; i++)
- {
- double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
- if (mpiHelper.rank()==0)
- std::cout << "Homotopy step: " << i << ", parameter: " << homotopyParameter << std::endl;
-
-
- // ////////////////////////////////////////////////////////////
- // Create an assembler for the energy functional
- // ////////////////////////////////////////////////////////////
-
- const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
- shared_ptr<NeumannFunction> neumannFunction;
- if (parameterSet.hasKey("neumannValues"))
- neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
- homotopyParameter);
-
- std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,3> >("neumannValues") << std::endl;
-
- if (mpiHelper.rank() == 0)
- {
- std::cout << "Material parameters:" << std::endl;
- materialParameters.report();
- }
-
- // Assembler using ADOL-C
- std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
- std::shared_ptr<LocalFEStiffness<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- std::vector<Dune::FieldVector<adouble, 3> > > > elasticEnergy;
-
- if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
- elasticEnergy = std::make_shared<StVenantKirchhoffEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "neohooke")
- elasticEnergy = std::make_shared<NeoHookeEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "hencky")
- 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(!elasticEnergy)
- DUNE_THROW(Exception, "Error: Selected energy not available!");
-
- auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> >(&neumannBoundary,neumannFunction.get());
-
- SumEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> totalEnergy(elasticEnergy, neumannEnergy);
-
- LocalADOLCStiffness<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- SolutionType> localADOLCStiffness(&totalEnergy);
-
- FEAssembler<FufemFEBasis,SolutionType> assembler(fufemFEBasis, &localADOLCStiffness);
-
- // /////////////////////////////////////////////////
- // Create a Riemannian trust-region solver
- // /////////////////////////////////////////////////
-
- TrustRegionSolver<FufemFEBasis,SolutionType> solver;
- solver.setup(*grid,
- &assembler,
- x,
- dirichletDofs,
- tolerance,
- maxTrustRegionSteps,
- initialTrustRegionRadius,
- multigridIterations,
- mgTolerance,
- mu, nu1, nu2,
- baseIterations,
- baseTolerance
- );
-
- ////////////////////////////////////////////////////////
- // Set Dirichlet values
- ////////////////////////////////////////////////////////
-
- // The python class that implements the Dirichlet values
- Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
-
- // The member method 'DirichletValues' of that class
- Python::Callable C = dirichletValuesClass.get("DirichletValues");
-
- // Call a constructor.
- Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
-
- // Extract object member functions as Dune functions
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > dirichletValues(dirichletValuesPythonObject.get("dirichletValues"));
-
- ::Functions::interpolate(fufemFEBasis, x, dirichletValues, dirichletDofs);
-
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
-
- solver.setInitialIterate(x);
- solver.solve();
-
- x = solver.getSol();
-
- /////////////////////////////////
- // Output result
- /////////////////////////////////
-
- // Compute the displacement
- auto displacement = x;
- displacement -= identity;
-
- Dune::Functions::DiscreteScalarGlobalBasisFunction<FEBasis,SolutionType> displacementFunction(feBasis,displacement);
- auto localDisplacementFunction = localFunction(displacementFunction);
-
- // We need to subsample, because VTK cannot natively display real second-order functions
- SubsamplingVTKWriter<GridView> vtkWriter(gridView,2);
- vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
- vtkWriter.write(resultPath + "finite-strain_homotopy_" + std::to_string(i+1));
- }
-
-} catch (Exception e) {
- std::cout << e << std::endl;
-}
--
GitLab