From 396b7ddbcb3a85ea24cd81b2e49abce83b836f5a Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 23 Mar 2015 08:30:41 +0000
Subject: [PATCH] Get rid of the mechanism to set point loads
[[Imported from SVN: r10104]]
---
dune/gfe/feassembler.hh | 22 +++++++---------------
dune/gfe/localadolcstiffness.hh | 12 ++++--------
dune/gfe/localfestiffness.hh | 14 +-------------
dune/gfe/stvenantkirchhoffenergy.hh | 10 ++--------
dune/gfe/trustregionsolver.cc | 9 +++------
dune/gfe/trustregionsolver.hh | 5 +----
src/finite-strain-elasticity.cc | 8 +-------
7 files changed, 19 insertions(+), 61 deletions(-)
diff --git a/dune/gfe/feassembler.hh b/dune/gfe/feassembler.hh
index 700d93b7..540ba64e 100644
--- a/dune/gfe/feassembler.hh
+++ b/dune/gfe/feassembler.hh
@@ -50,14 +50,12 @@ public:
* anyway to compute the Riemannian Hessian.
*/
virtual void assembleGradientAndHessian(const VectorType& sol,
- const VectorType& pointLoads,
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 VectorType& pointLoads) const;
+ virtual double computeEnergy(const VectorType& sol) const;
//protected:
void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
@@ -101,7 +99,6 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
template <class Basis, class VectorType>
void FEAssembler<Basis,VectorType>::
assembleGradientAndHessian(const VectorType& sol,
- const VectorType& pointLoads,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern) const
@@ -130,15 +127,13 @@ assembleGradientAndHessian(const VectorType& sol,
VectorType localSolution(numOfBaseFct);
VectorType localPointLoads(numOfBaseFct);
- for (int i=0; i<numOfBaseFct; i++) {
+ for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[basis_.index(*it,i)];
- localPointLoads[i] = pointLoads[basis_.index(*it,i)];
- }
std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
// setup local matrix and gradient
- localStiffness_->assembleGradientAndHessian(*it, basis_.getLocalFiniteElement(*it), localSolution, localPointLoads, localGradient);
+ localStiffness_->assembleGradientAndHessian(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
@@ -164,12 +159,12 @@ assembleGradientAndHessian(const VectorType& sol,
template <class Basis, class VectorType>
double FEAssembler<Basis, VectorType>::
-computeEnergy(const VectorType& sol, const VectorType& pointLoads) const
+computeEnergy(const VectorType& sol) const
{
double energy = 0;
if (sol.size()!=basis_.size())
- DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
+ 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>();
@@ -181,14 +176,11 @@ computeEnergy(const VectorType& sol, const VectorType& pointLoads) const
size_t nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
VectorType localSolution(nDofs);
- VectorType localPointLoads(nDofs);
- for (size_t i=0; i<nDofs; i++) {
+ for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[basis_.index(*it,i)];
- localPointLoads[i] = pointLoads[basis_.index(*it,i)];
- }
- energy += localStiffness_->energy(*it, basis_.getLocalFiniteElement(*it), localSolution, localPointLoads);
+ energy += localStiffness_->energy(*it, basis_.getLocalFiniteElement(*it), localSolution);
}
diff --git a/dune/gfe/localadolcstiffness.hh b/dune/gfe/localadolcstiffness.hh
index 97b9acbe..ef25d653 100644
--- a/dune/gfe/localadolcstiffness.hh
+++ b/dune/gfe/localadolcstiffness.hh
@@ -45,8 +45,7 @@ public:
/** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& e,
const LocalFiniteElement& localFiniteElement,
- const VectorType& localConfiguration,
- const std::vector<Dune::FieldVector<double,3> >& localPointLoads) const;
+ const VectorType& localConfiguration) const;
/** \brief Assemble the local stiffness matrix at the current position
@@ -55,7 +54,6 @@ public:
virtual void assembleGradientAndHessian(const Entity& e,
const LocalFiniteElement& localFiniteElement,
const VectorType& localConfiguration,
- const std::vector<Dune::FieldVector<double,3> >& localPointLoads,
VectorType& localGradient);
const LocalFEStiffness<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
@@ -68,8 +66,7 @@ typename LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::RT
LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
- const VectorType& localSolution,
- const std::vector<Dune::FieldVector<double,3> >& localPointLoads) const
+ const VectorType& localSolution) const
{
double pureEnergy;
@@ -83,7 +80,7 @@ energy(const Entity& element,
for (size_t j=0; j<localSolution[i].size(); j++)
localASolution[i][j] <<= localSolution[i][j];
- energy = localEnergy_->energy(element,localFiniteElement,localASolution,localPointLoads);
+ energy = localEnergy_->energy(element,localFiniteElement,localASolution);
energy >>= pureEnergy;
@@ -104,11 +101,10 @@ void LocalADOLCStiffness<GridType, LocalFiniteElement, VectorType>::
assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const VectorType& localSolution,
- const std::vector<Dune::FieldVector<double,3> >& localPointLoads,
VectorType& localGradient)
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
- energy(element, localFiniteElement, localSolution, localPointLoads);
+ energy(element, localFiniteElement, localSolution);
/////////////////////////////////////////////////////////////////
// Compute the energy gradient
diff --git a/dune/gfe/localfestiffness.hh b/dune/gfe/localfestiffness.hh
index c6ea0110..413d0eb2 100644
--- a/dune/gfe/localfestiffness.hh
+++ b/dune/gfe/localfestiffness.hh
@@ -24,27 +24,16 @@ public:
/** \brief Assemble the local stiffness matrix at the current position
This default implementation used finite-difference approximations to compute the second derivatives
-
- The formula for the Riemannian Hessian has been taken from Absil, Mahony, Sepulchre:
- 'Optimization algorithms on matrix manifolds', page 107. There it says that
- \f[
- \langle Hess f(x)[\xi], \eta \rangle
- = \frac 12 \frac{d^2}{dt^2} \Big(f(\exp_x(t(\xi + \eta))) - f(\exp_x(t\xi)) - f(\exp_x(t\eta))\Big)\Big|_{t=0}.
- \f]
- We compute that using a finite difference approximation.
-
*/
virtual void assembleGradientAndHessian(const Entity& e,
const LocalFiniteElement& localFiniteElement,
const VectorType& localConfiguration,
- const std::vector<Dune::FieldVector<double,3> >& localPointLoads,
VectorType& localGradient);
/** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& e,
const LocalFiniteElement& localFiniteElement,
- const VectorType& localConfiguration,
- const std::vector<Dune::FieldVector<double,3> >& localPointLoads) const = 0;
+ const VectorType& localConfiguration) const = 0;
// assembled data
Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
@@ -60,7 +49,6 @@ void LocalFEStiffness<GridType, LocalFiniteElement, VectorType>::
assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const VectorType& localConfiguration,
- const std::vector<Dune::FieldVector<double,3> >& localPointLoads,
VectorType& localGradient)
{
DUNE_THROW(Dune::NotImplemented, "!");
diff --git a/dune/gfe/stvenantkirchhoffenergy.hh b/dune/gfe/stvenantkirchhoffenergy.hh
index 223e5152..e54b53ae 100644
--- a/dune/gfe/stvenantkirchhoffenergy.hh
+++ b/dune/gfe/stvenantkirchhoffenergy.hh
@@ -45,8 +45,7 @@ public:
/** \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 std::vector<Dune::FieldVector<double,gridDim> >& localPointLoads) const;
+ const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
/** \brief Lame constants */
double mu_, lambda_;
@@ -63,8 +62,7 @@ field_type
StVenantKirchhoffEnergy<GridView,LocalFiniteElement,field_type>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
- const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration,
- const std::vector<Dune::FieldVector<double,gridDim> >& localPointLoads) const
+ const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
{
assert(element.type() == localFiniteElement.type());
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
@@ -146,10 +144,6 @@ energy(const Entity& element,
// Assemble boundary contributions
//////////////////////////////////////////////////////////////////////////////
- for (size_t i=0; i<localPointLoads.size(); i++)
- for (size_t j=0; j<dim; j++)
- energy -= localConfiguration[i][j] * localPointLoads[i][j];
-
if (not neumannFunction_)
return energy;
diff --git a/dune/gfe/trustregionsolver.cc b/dune/gfe/trustregionsolver.cc
index eb2bb37d..50adc5c5 100644
--- a/dune/gfe/trustregionsolver.cc
+++ b/dune/gfe/trustregionsolver.cc
@@ -37,8 +37,7 @@ setup(const typename BasisType::GridView::Grid& grid,
int nu1,
int nu2,
int baseIterations,
- double baseTolerance,
- const SolutionType& pointLoads)
+ double baseTolerance)
{
grid_ = &grid;
assembler_ = assembler;
@@ -49,7 +48,6 @@ setup(const typename BasisType::GridView::Grid& grid,
innerIterations_ = multigridIterations;
innerTolerance_ = mgTolerance;
ignoreNodes_ = &dirichletNodes;
- pointLoads_ = pointLoads;
int numLevels = grid_->maxLevel()+1;
@@ -206,7 +204,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
// Trust-Region Solver
// /////////////////////////////////////////////////////
- double oldEnergy = assembler_->computeEnergy(x_, pointLoads_);
+ double oldEnergy = assembler_->computeEnergy(x_);
bool recomputeGradientHessian = true;
CorrectionType rhs;
@@ -228,7 +226,6 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (recomputeGradientHessian) {
assembler_->assembleGradientAndHessian(x_,
- pointLoads_,
rhs,
*hessianMatrix_,
i==0 // assemble occupation pattern only for the first call
@@ -326,7 +323,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
for (size_t j=0; j<newIterate.size(); j++)
newIterate[j] += corr[j];
- double energy = assembler_->computeEnergy(newIterate, pointLoads_);
+ double energy = assembler_->computeEnergy(newIterate);
// compute the model decrease
// It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
diff --git a/dune/gfe/trustregionsolver.hh b/dune/gfe/trustregionsolver.hh
index 3fc5c871..8ce1df7c 100644
--- a/dune/gfe/trustregionsolver.hh
+++ b/dune/gfe/trustregionsolver.hh
@@ -60,8 +60,7 @@ public:
int nu1,
int nu2,
int baseIterations,
- double baseTolerance,
- const SolutionType& pointLoads);
+ double baseTolerance);
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
@@ -124,8 +123,6 @@ protected:
/** \brief An L2-norm, really. The H1SemiNorm class is badly named */
std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
- SolutionType pointLoads_;
-
public:
VectorType identity_;
};
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index c45b2a58..67d15999 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -273,11 +273,6 @@ int main (int argc, char *argv[]) try
FEAssembler<FEBasis,SolutionType> assembler(gridView, &localADOLCStiffness);
- std::vector<FieldVector<double,3> > pointLoads(x.size());
- std::fill(pointLoads.begin(), pointLoads.end(), 0);
-// pointLoads[1372] = parameterSet.get<FieldVector<double,3> >("neumannValues");
-// pointLoads[1372] *= 0.5;
-
// /////////////////////////////////////////////////
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
@@ -294,8 +289,7 @@ int main (int argc, char *argv[]) try
mgTolerance,
mu, nu1, nu2,
baseIterations,
- baseTolerance,
- pointLoads
+ baseTolerance
);
solver.identity_ = identity;
--
GitLab