diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh index 1ffe631f27a56470fc1bb55917bd1ae25b8515bf..eff72c112f659fb8b3a4a548803932dcf25bc65e 100644 --- a/dune/gfe/localgeodesicfestiffness.hh +++ b/dune/gfe/localgeodesicfestiffness.hh @@ -10,18 +10,18 @@ template<class GridView, class LocalFiniteElement, class TargetSpace> -class LocalGeodesicFEStiffness +class LocalGeodesicFEStiffness { // grid types typedef typename GridView::Grid::ctype DT; typedef typename TargetSpace::ctype RT; typedef typename GridView::template Codim<0>::Entity Entity; - + // some other sizes enum {gridDim=GridView::dimension}; public: - + //! Dimension of a tangent space enum { blocksize = TargetSpace::TangentVector::dimension }; @@ -50,17 +50,17 @@ public: const LocalFiniteElement& localFiniteElement, const std::vector<TargetSpace>& localSolution) const = 0; - /** \brief Assemble the element gradient of the energy functional + /** \brief Assemble the element gradient of the energy functional The default implementation in this class uses a finite difference approximation */ virtual void assembleGradient(const Entity& element, const LocalFiniteElement& localFiniteElement, const std::vector<TargetSpace>& solution, std::vector<typename TargetSpace::TangentVector>& gradient) const; - + // assembled data Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > A_; - + }; @@ -84,18 +84,18 @@ assembleGradient(const Entity& element, std::vector<TargetSpace> backwardSolution = localSolution; for (size_t i=0; i<localSolution.size(); i++) { - + // basis vectors of the tangent space of the i-th entry of localSolution const Dune::FieldMatrix<double,blocksize,embeddedBlocksize> B = localSolution[i].orthonormalFrame(); for (int j=0; j<blocksize; j++) { - + typename TargetSpace::EmbeddedTangentVector forwardCorrection = B[j]; forwardCorrection *= eps; - + typename TargetSpace::EmbeddedTangentVector backwardCorrection = B[j]; backwardCorrection *= -eps; - + forwardSolution[i] = TargetSpace::exp(localSolution[i], forwardCorrection); backwardSolution[i] = TargetSpace::exp(localSolution[i], backwardCorrection); @@ -105,7 +105,7 @@ assembleGradient(const Entity& element, forwardSolution[i] = localSolution[i]; backwardSolution[i] = localSolution[i]; - + } } @@ -129,7 +129,7 @@ assembleHessian(const Entity& element, A_ = 0; const double eps = 1e-4; - + std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size()); for (size_t i=0; i<B.size(); i++) B[i] = localSolution[i].orthonormalFrame(); @@ -137,7 +137,7 @@ assembleHessian(const Entity& element, // Precompute negative energy at the current configuration // (negative because that is how we need it as part of the 2nd-order fd formula) double centerValue = -energy(element, localFiniteElement, localSolution); - + // Precompute energy infinitesimal corrections in the directions of the local basis vectors std::vector<Dune::array<double,blocksize> > forwardEnergy(nDofs); std::vector<Dune::array<double,blocksize> > backwardEnergy(nDofs); @@ -149,18 +149,18 @@ assembleHessian(const Entity& element, epsXi *= eps; Dune::FieldVector<double,embeddedBlocksize> minusEpsXi = epsXi; minusEpsXi *= -1; - + std::vector<TargetSpace> forwardSolution = localSolution; std::vector<TargetSpace> backwardSolution = localSolution; - + forwardSolution[i] = TargetSpace::exp(localSolution[i],epsXi); backwardSolution[i] = TargetSpace::exp(localSolution[i],minusEpsXi); - + forwardEnergy[i][i2] = energy(element, localFiniteElement, forwardSolution); backwardEnergy[i][i2] = energy(element, localFiniteElement, backwardSolution); - + } - + } // finite-difference approximation @@ -174,20 +174,20 @@ assembleHessian(const Entity& element, std::vector<TargetSpace> forwardSolutionXiEta = localSolution; std::vector<TargetSpace> backwardSolutionXiEta = localSolution; - + Dune::FieldVector<double,embeddedBlocksize> epsXi = B[i][i2]; epsXi *= eps; Dune::FieldVector<double,embeddedBlocksize> epsEta = B[j][j2]; epsEta *= eps; - + Dune::FieldVector<double,embeddedBlocksize> minusEpsXi = epsXi; minusEpsXi *= -1; Dune::FieldVector<double,embeddedBlocksize> minusEpsEta = epsEta; minusEpsEta *= -1; - + if (i==j) forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi+epsEta); else { forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi); forwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],epsEta); } - + if (i==j) backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi+minusEpsEta); else { @@ -197,9 +197,9 @@ assembleHessian(const Entity& element, double forwardValue = energy(element, localFiniteElement, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2]; double backwardValue = energy(element, localFiniteElement, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2]; - + A_[i][j][i2][j2] = A_[j][i][j2][i2] = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps); - + } } }