From 282c8f13f8fb4cb670602469300e251c5d895b7b Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Tue, 20 Mar 2018 17:17:50 +0100
Subject: [PATCH] .
---
src/spatial-solving/fixedpointiterator.cc | 88 ++++++++++++++++-------
src/spatial-solving/fixedpointiterator.hh | 4 +-
src/time-stepping/adaptivetimestepper.cc | 2 +-
src/time-stepping/adaptivetimestepper.hh | 4 +-
src/time-stepping/coupledtimestepper.cc | 2 +-
src/time-stepping/coupledtimestepper.hh | 4 +-
6 files changed, 70 insertions(+), 34 deletions(-)
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 81be4536..521c513b 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -52,28 +52,22 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterationCounter
FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
- Updaters updaters, const std::vector<const Matrix*>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+ Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
- EnergyNorm<Matrix, Vector> energyNorm(velocityMatrices[0]);
- LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
- velocityTolerance_, &energyNorm,
- verbosity_, false); // absolute error
-
- // assemble full global contact problem
auto contactAssembler = factory_.getNBodyAssembler();
- Matrix bilinearForm;
- contactAssembler.assembleJacobian(velocityMatrices, bilinearForm);
-
- Vector totalRhs;
- contactAssembler.assembleRightHandSide(velocityRHSs, totalRhs);
-
- Vector totalVelocityIterate;
- contactAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+ Matrix mat;
+ std::vector<const Matrix*> matrices_ptr(velocityMatrices.size());
+ for (size_t i=0; i<matrices_ptr.size(); i++) {
+ matrices_ptr[i] = &velocityMatrices[i];
+ }
+ contactAssembler.assembleJacobian(matrices_ptr, mat);
- // contribution from nonlinearity
- // add gradient to rhs, hessian to matrix!?
+ EnergyNorm<Matrix, Vector> energyNorm(mat);
+ LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
+ velocityTolerance_, &energyNorm,
+ verbosity_, false); // absolute error
size_t fixedPointIteration;
size_t multigridIterations = 0;
@@ -82,15 +76,43 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) {
- // solve a velocity problem
+ // compute relative velocities
+ std::vector<Vector> v_rel;
+ relativeVelocities(velocityIterates, v_rel);
+
+ // contribution from nonlinearity
for (size_t i=0; i<alpha.size(); i++) {
globalFriction_[i]->updateAlpha(alpha[i]);
}
+ std::vector<Matrix> matrices(velocityMatrices.size());
+ std::vector<Vector> rhs(velocityRHSs.size());
+ for (size_t i=0; i<globalFriction_.size(); i++) {
+ matrices[i] = velocityMatrices[i];
+ rhs[i] = velocityRHSs[i];
+
+ globalFriction_[i]->addHessian(v_rel[i], matrices[i]);
+ globalFriction_[i]->addGradient(v_rel[i], rhs[i]);
+
+ matrices_ptr[i] = &matrices[i];
+ }
+
+ // assemble full global contact problem
+ Matrix bilinearForm;
+ contactAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+ Vector totalRhs;
+ contactAssembler.assembleRightHandSide(rhs, totalRhs);
+
+ Vector totalVelocityIterate;
+ contactAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+
+ // solve a velocity problem
step_->setProblem(bilinearForm, totalVelocityIterate, totalRhs);
velocityProblemSolver.preprocess();
velocityProblemSolver.solve();
+ contactAssembler.postprocess(totalVelocityIterate, velocityIterates);
multigridIterations += velocityProblemSolver.getResult().iterations;
@@ -103,19 +125,26 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
}
// compute relative velocities on contact boundaries
- relativeVelocities(v_m);
+ relativeVelocities(v_m, v_rel);
// solve a state problem
- updaters.state_->solve(v_m);
- ScalarVector newAlpha;
- /* updaters.state_->extractAlpha(newAlpha);
+ updaters.state_->solve(v_rel);
+ std::vector<ScalarVector> newAlpha;
+ updaters.state_->extractAlpha(newAlpha);
- if (lambda_ < 1e-12 or
- errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
+ bool breakCriterion = true;
+ for (size_t i=0; i<alpha.size(); i++) {
+ if (errorNorm_.diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
+ breakCriterion = false;
+ break;
+ }
+ }
+
+ if (lambda_ < 1e-12 or breakCriterion) {
fixedPointIteration++;
break;
}
- alpha = newAlpha;*/
+ alpha = newAlpha;
}
if (fixedPointIteration == fixedPointMaxIterations_)
DUNE_THROW(Dune::Exception, "FPI failed to converge");
@@ -137,7 +166,14 @@ std::ostream &operator<<(std::ostream &stream,
}
template <class Factory, class Updaters, class ErrorNorm>
-void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::vector<Vector>& v_m) const {
+void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const {
+ // init result
+ v_rel.resize(v.size());
+ for (size_t i=0; i<v_rel.size(); i++) {
+ v_rel[i].resize(v[i].size());
+ v_rel[i] = 0;
+ }
+
// needs assemblers to obtain basis
/*
std::vector<std::shared_ptr<MyAssembler>> assemblers(bodyCount);
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index ed1fe9b4..e579645b 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -32,7 +32,7 @@ class FixedPointIterator {
using DeformedGrid = typename Factory::DeformedGrid;
private:
- void relativeVelocities(std::vector<Vector>& v_m) const;
+ void relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const;
public:
FixedPointIterator(Factory& factory, const Dune::ParameterTree& parset,
@@ -40,7 +40,7 @@ class FixedPointIterator {
const ErrorNorm& errorNorm);
FixedPointIterationCounter run(Updaters updaters,
- const std::vector<const Matrix*>& velocityMatrices,
+ const std::vector<Matrix>& velocityMatrices,
const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates);
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index 2598ecaf..0808750f 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -22,7 +22,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
Factory &factory, Dune::ParameterTree const &parset,
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters ¤t,
double relativeTime, double relativeTau,
- std::function<void(double, Vector &)> externalForces,
+ std::vector<std::function<void(double, Vector &)>>& externalForces,
ErrorNorm const &errorNorm,
std::function<bool(Updaters &, Updaters &)> mustRefine)
: relativeTime_(relativeTime),
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index 5d559a8b..a0bf6113 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -33,7 +33,7 @@ class AdaptiveTimeStepper {
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
Updaters ¤t, double relativeTime,
double relativeTau,
- std::function<void(double, Vector &)> externalForces,
+ std::vector<std::function<void(double, Vector &)>>& externalForces,
ErrorNorm const &errorNorm,
std::function<bool(Updaters &, Updaters &)> mustRefine);
@@ -53,7 +53,7 @@ class AdaptiveTimeStepper {
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
Updaters ¤t_;
UpdatersWithCount R1_;
- std::function<void(double, Vector &)> externalForces_;
+ std::vector<std::function<void(double, Vector &)>>& externalForces_;
std::function<bool(Updaters &, Updaters &)> mustRefine_;
ErrorNorm const &errorNorm_;
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index 410a48f8..7501f68a 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -9,7 +9,7 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
double finalTime, Factory &factory, Dune::ParameterTree const &parset,
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters updaters,
ErrorNorm const &errorNorm,
- std::function<void(double, Vector &)> externalForces)
+ std::vector<std::function<void(double, Vector &)>>& externalForces)
: finalTime_(finalTime),
factory_(factory),
parset_(parset),
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 52543d1e..92e3cd4a 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -19,7 +19,7 @@ class CoupledTimeStepper {
Dune::ParameterTree const &parset,
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
Updaters updaters, ErrorNorm const &errorNorm,
- std::function<void(double, Vector &)> externalForces);
+ std::vector<std::function<void(double, Vector &)>>& externalForces);
FixedPointIterationCounter step(double relativeTime, double relativeTau);
@@ -29,7 +29,7 @@ class CoupledTimeStepper {
Dune::ParameterTree const &parset_;
std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
Updaters updaters_;
- std::function<void(double, Vector &)> externalForces_;
+ std::vector<std::function<void(double, Vector &)>>& externalForces_;
ErrorNorm const &errorNorm_;
};
#endif
--
GitLab