diff --git a/dune/elasticity/assemblers/localhyperdualstiffness.hh b/dune/elasticity/assemblers/localhyperdualstiffness.hh index a91b68ef0d31a6c0bd7da8507a20c301d9d528e0..058935904e1618b3a544e6e959077790bf26b474 100644 --- a/dune/elasticity/assemblers/localhyperdualstiffness.hh +++ b/dune/elasticity/assemblers/localhyperdualstiffness.hh @@ -8,7 +8,7 @@ namespace Dune::Elasticity { -/** \brief Assembles energy gradient and Hessian using Hyperdual Numbers +/** \brief Assembles energy gradient and Hessian using Hyperdual Numbers */ template<class LocalView> class LocalHyperDualStiffness @@ -36,7 +36,7 @@ public: std::vector<double>& localGradient, Dune::Matrix<double>& localHessian); - std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> localEnergy_; + std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> localEnergy_; }; @@ -50,10 +50,10 @@ energy(const LocalView& localView, std::vector<hyperdual> localHyperDualConfiguration(localConfiguration.size()); hyperdual energy; - + for (size_t i=0; i<localConfiguration.size(); i++) - localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]); - + localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]); + energy = localEnergy_->energy(localView,localHyperDualConfiguration); pureEnergy = energy.real(); @@ -70,7 +70,7 @@ assembleGradientAndHessian(const LocalView& localView, Dune::Matrix<double>& localHessian ) { - size_t nDoubles = localConfiguration.size(); + size_t nDoubles = localConfiguration.size(); localGradient.resize(nDoubles); @@ -79,28 +79,28 @@ assembleGradientAndHessian(const LocalView& localView, std::vector<hyperdual> localHyperDualConfiguration(nDoubles); for (size_t i=0; i<nDoubles; i++) localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]); - + // Compute gradient and hessian (symmetry is used) using hyperdual function evaluation // localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives - for (size_t i=0; i<nDoubles; i++) + for (size_t i=0; i<nDoubles; i++) { localHyperDualConfiguration[i].seteps1(1.0); localHyperDualConfiguration[i].seteps2(1.0); - + hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration); - localGradient[i] = temp.eps1(); + localGradient[i] = temp.eps1(); localHessian[i][i] = temp.eps1eps2(); - + localHyperDualConfiguration[i].seteps2(0.0); - for (size_t j=i+1; j<nDoubles; j++) - { + for (size_t j=i+1; j<nDoubles; j++) + { localHyperDualConfiguration[j].seteps2(1.0); localHessian[i][j] = localEnergy_->energy(localView, localHyperDualConfiguration).eps1eps2(); - localHessian[j][i] = localHessian[i][j]; // Use symmetry of hessian + localHessian[j][i] = localHessian[i][j]; // Use symmetry of hessian localHyperDualConfiguration[j].seteps2(0.0); } - + localHyperDualConfiguration[i].seteps1(0.0); }