Skip to content
Snippets Groups Projects
Commit 396b7ddb authored by Oliver Sander's avatar Oliver Sander Committed by sander
Browse files

Get rid of the mechanism to set point loads 

[[Imported from SVN: r10104]]
parent ec5ea5e6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/feassembler.hh
+ 7
− 15
View file @ 396b7ddb
...@@ -50,14 +50,12 @@ public: ...@@ -50,14 +50,12 @@ public:
* anyway to compute the Riemannian Hessian. * anyway to compute the Riemannian Hessian.
*/ */
virtual void assembleGradientAndHessian(const VectorType& sol, virtual void assembleGradientAndHessian(const VectorType& sol,
const VectorType& pointLoads,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient, Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian, Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern=true) const; bool computeOccupationPattern=true) const;
/** \brief Compute the energy of a deformation state */ /** \brief Compute the energy of a deformation state */
virtual double computeEnergy(const VectorType& sol, virtual double computeEnergy(const VectorType& sol) const;
const VectorType& pointLoads) const;
//protected: //protected:
void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const; void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
...@@ -101,7 +99,6 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const ...@@ -101,7 +99,6 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
template <class Basis, class VectorType> template <class Basis, class VectorType>
void FEAssembler<Basis,VectorType>:: void FEAssembler<Basis,VectorType>::
assembleGradientAndHessian(const VectorType& sol, assembleGradientAndHessian(const VectorType& sol,
const VectorType& pointLoads,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient, Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian, Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern) const bool computeOccupationPattern) const
...@@ -130,15 +127,13 @@ assembleGradientAndHessian(const VectorType& sol, ...@@ -130,15 +127,13 @@ assembleGradientAndHessian(const VectorType& sol,
VectorType localSolution(numOfBaseFct); VectorType localSolution(numOfBaseFct);
VectorType localPointLoads(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)]; localSolution[i] = sol[basis_.index(*it,i)];
localPointLoads[i] = pointLoads[basis_.index(*it,i)];
}
std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct); std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
// setup local matrix and gradient // 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 // Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) { for(int i=0; i<numOfBaseFct; i++) {
...@@ -164,12 +159,12 @@ assembleGradientAndHessian(const VectorType& sol, ...@@ -164,12 +159,12 @@ assembleGradientAndHessian(const VectorType& sol,
template <class Basis, class VectorType> template <class Basis, class VectorType>
double FEAssembler<Basis, VectorType>:: double FEAssembler<Basis, VectorType>::
computeEnergy(const VectorType& sol, const VectorType& pointLoads) const computeEnergy(const VectorType& sol) const
{ {
double energy = 0; double energy = 0;
if (sol.size()!=basis_.size()) 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 it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.getGridView().template end<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 ...@@ -181,14 +176,11 @@ computeEnergy(const VectorType& sol, const VectorType& pointLoads) const
size_t nDofs = basis_.getLocalFiniteElement(*it).localBasis().size(); size_t nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
VectorType localSolution(nDofs); 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)]; 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);
} }
......
dune/gfe/localadolcstiffness.hh
+ 4
− 8
View file @ 396b7ddb
...@@ -45,8 +45,7 @@ public: ...@@ -45,8 +45,7 @@ public:
/** \brief Compute the energy at the current configuration */ /** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& e, virtual RT energy (const Entity& e,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const VectorType& localConfiguration, const VectorType& localConfiguration) const;
const std::vector<Dune::FieldVector<double,3> >& localPointLoads) const;
/** \brief Assemble the local stiffness matrix at the current position /** \brief Assemble the local stiffness matrix at the current position
...@@ -55,7 +54,6 @@ public: ...@@ -55,7 +54,6 @@ public:
virtual void assembleGradientAndHessian(const Entity& e, virtual void assembleGradientAndHessian(const Entity& e,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const VectorType& localConfiguration, const VectorType& localConfiguration,
const std::vector<Dune::FieldVector<double,3> >& localPointLoads,
VectorType& localGradient); VectorType& localGradient);
const LocalFEStiffness<GridView, LocalFiniteElement, AVectorType>* localEnergy_; const LocalFEStiffness<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
...@@ -68,8 +66,7 @@ typename LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::RT ...@@ -68,8 +66,7 @@ typename LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::RT
LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>:: LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::
energy(const Entity& element, energy(const Entity& element,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const VectorType& localSolution, const VectorType& localSolution) const
const std::vector<Dune::FieldVector<double,3> >& localPointLoads) const
{ {
double pureEnergy; double pureEnergy;
...@@ -83,7 +80,7 @@ energy(const Entity& element, ...@@ -83,7 +80,7 @@ energy(const Entity& element,
for (size_t j=0; j<localSolution[i].size(); j++) for (size_t j=0; j<localSolution[i].size(); j++)
localASolution[i][j] <<= localSolution[i][j]; localASolution[i][j] <<= localSolution[i][j];
energy = localEnergy_->energy(element,localFiniteElement,localASolution,localPointLoads); energy = localEnergy_->energy(element,localFiniteElement,localASolution);
energy >>= pureEnergy; energy >>= pureEnergy;
...@@ -104,11 +101,10 @@ void LocalADOLCStiffness<GridType, LocalFiniteElement, VectorType>:: ...@@ -104,11 +101,10 @@ void LocalADOLCStiffness<GridType, LocalFiniteElement, VectorType>::
assembleGradientAndHessian(const Entity& element, assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const VectorType& localSolution, const VectorType& localSolution,
const std::vector<Dune::FieldVector<double,3> >& localPointLoads,
VectorType& localGradient) VectorType& localGradient)
{ {
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap. // 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 // Compute the energy gradient
......
dune/gfe/localfestiffness.hh
+ 1
− 13
View file @ 396b7ddb
...@@ -24,27 +24,16 @@ public: ...@@ -24,27 +24,16 @@ public:
/** \brief Assemble the local stiffness matrix at the current position /** \brief Assemble the local stiffness matrix at the current position
This default implementation used finite-difference approximations to compute the second derivatives 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, virtual void assembleGradientAndHessian(const Entity& e,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const VectorType& localConfiguration, const VectorType& localConfiguration,
const std::vector<Dune::FieldVector<double,3> >& localPointLoads,
VectorType& localGradient); VectorType& localGradient);
/** \brief Compute the energy at the current configuration */ /** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& e, virtual RT energy (const Entity& e,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const VectorType& localConfiguration, const VectorType& localConfiguration) const = 0;
const std::vector<Dune::FieldVector<double,3> >& localPointLoads) const = 0;
// assembled data // assembled data
Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_; Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
...@@ -60,7 +49,6 @@ void LocalFEStiffness<GridType, LocalFiniteElement, VectorType>:: ...@@ -60,7 +49,6 @@ void LocalFEStiffness<GridType, LocalFiniteElement, VectorType>::
assembleGradientAndHessian(const Entity& element, assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const VectorType& localConfiguration, const VectorType& localConfiguration,
const std::vector<Dune::FieldVector<double,3> >& localPointLoads,
VectorType& localGradient) VectorType& localGradient)
{ {
DUNE_THROW(Dune::NotImplemented, "!"); DUNE_THROW(Dune::NotImplemented, "!");
......
dune/gfe/stvenantkirchhoffenergy.hh
+ 2
− 8
View file @ 396b7ddb
...@@ -45,8 +45,7 @@ public: ...@@ -45,8 +45,7 @@ public:
/** \brief Assemble the energy for a single element */ /** \brief Assemble the energy for a single element */
field_type energy (const Entity& e, field_type energy (const Entity& e,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration, const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
const std::vector<Dune::FieldVector<double,gridDim> >& localPointLoads) const;
/** \brief Lame constants */ /** \brief Lame constants */
double mu_, lambda_; double mu_, lambda_;
...@@ -63,8 +62,7 @@ field_type ...@@ -63,8 +62,7 @@ field_type
StVenantKirchhoffEnergy<GridView,LocalFiniteElement,field_type>:: StVenantKirchhoffEnergy<GridView,LocalFiniteElement,field_type>::
energy(const Entity& element, energy(const Entity& element,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration, const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
const std::vector<Dune::FieldVector<double,gridDim> >& localPointLoads) const
{ {
assert(element.type() == localFiniteElement.type()); assert(element.type() == localFiniteElement.type());
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry; typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
...@@ -146,10 +144,6 @@ energy(const Entity& element, ...@@ -146,10 +144,6 @@ energy(const Entity& element,
// Assemble boundary contributions // 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_) if (not neumannFunction_)
return energy; return energy;
......
dune/gfe/trustregionsolver.cc
+ 3
− 6
View file @ 396b7ddb
...@@ -37,8 +37,7 @@ setup(const typename BasisType::GridView::Grid& grid, ...@@ -37,8 +37,7 @@ setup(const typename BasisType::GridView::Grid& grid,
int nu1, int nu1,
int nu2, int nu2,
int baseIterations, int baseIterations,
double baseTolerance, double baseTolerance)
const SolutionType& pointLoads)
{ {
grid_ = &grid; grid_ = &grid;
assembler_ = assembler; assembler_ = assembler;
...@@ -49,7 +48,6 @@ setup(const typename BasisType::GridView::Grid& grid, ...@@ -49,7 +48,6 @@ setup(const typename BasisType::GridView::Grid& grid,
innerIterations_ = multigridIterations; innerIterations_ = multigridIterations;
innerTolerance_ = mgTolerance; innerTolerance_ = mgTolerance;
ignoreNodes_ = &dirichletNodes; ignoreNodes_ = &dirichletNodes;
pointLoads_ = pointLoads;
int numLevels = grid_->maxLevel()+1; int numLevels = grid_->maxLevel()+1;
...@@ -206,7 +204,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve() ...@@ -206,7 +204,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
// Trust-Region Solver // Trust-Region Solver
// ///////////////////////////////////////////////////// // /////////////////////////////////////////////////////
double oldEnergy = assembler_->computeEnergy(x_, pointLoads_); double oldEnergy = assembler_->computeEnergy(x_);
bool recomputeGradientHessian = true; bool recomputeGradientHessian = true;
CorrectionType rhs; CorrectionType rhs;
...@@ -228,7 +226,6 @@ void TrustRegionSolver<BasisType,VectorType>::solve() ...@@ -228,7 +226,6 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (recomputeGradientHessian) { if (recomputeGradientHessian) {
assembler_->assembleGradientAndHessian(x_, assembler_->assembleGradientAndHessian(x_,
pointLoads_,
rhs, rhs,
*hessianMatrix_, *hessianMatrix_,
i==0 // assemble occupation pattern only for the first call i==0 // assemble occupation pattern only for the first call
...@@ -326,7 +323,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve() ...@@ -326,7 +323,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
for (size_t j=0; j<newIterate.size(); j++) for (size_t j=0; j<newIterate.size(); j++)
newIterate[j] += corr[j]; newIterate[j] += corr[j];
double energy = assembler_->computeEnergy(newIterate, pointLoads_); double energy = assembler_->computeEnergy(newIterate);
// compute the model decrease // compute the model decrease
// It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs> // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
......
dune/gfe/trustregionsolver.hh
+ 1
− 4
View file @ 396b7ddb
...@@ -60,8 +60,7 @@ public: ...@@ -60,8 +60,7 @@ public:
int nu1, int nu1,
int nu2, int nu2,
int baseIterations, int baseIterations,
double baseTolerance, double baseTolerance);
const SolutionType& pointLoads);
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes) void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{ {
...@@ -124,8 +123,6 @@ protected: ...@@ -124,8 +123,6 @@ protected:
/** \brief An L2-norm, really. The H1SemiNorm class is badly named */ /** \brief An L2-norm, really. The H1SemiNorm class is badly named */
std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_; std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
SolutionType pointLoads_;
public: public:
VectorType identity_; VectorType identity_;
}; };
......
src/finite-strain-elasticity.cc
+ 1
− 7
View file @ 396b7ddb
...@@ -273,11 +273,6 @@ int main (int argc, char *argv[]) try ...@@ -273,11 +273,6 @@ int main (int argc, char *argv[]) try
FEAssembler<FEBasis,SolutionType> assembler(gridView, &localADOLCStiffness); 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 // Create a Riemannian trust-region solver
// ///////////////////////////////////////////////// // /////////////////////////////////////////////////
...@@ -294,8 +289,7 @@ int main (int argc, char *argv[]) try ...@@ -294,8 +289,7 @@ int main (int argc, char *argv[]) try
mgTolerance, mgTolerance,
mu, nu1, nu2, mu, nu1, nu2,
baseIterations, baseIterations,
baseTolerance, baseTolerance
pointLoads
); );
solver.identity_ = identity; solver.identity_ = identity;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment