diff --git a/dune/tectonic/spatial-solving/fixedpointiterator.cc b/dune/tectonic/spatial-solving/fixedpointiterator.cc
index bef1fea2ba3666075dbe5e1d0fbe6323737977a6..5b5c4c294d5e33aea18b6367bdba82e4027cf230 100644
--- a/dune/tectonic/spatial-solving/fixedpointiterator.cc
+++ b/dune/tectonic/spatial-solving/fixedpointiterator.cc
@@ -5,6 +5,7 @@
#include <dune/common/exceptions.hh>
#include <dune/matrix-vector/axpy.hh>
+#include <dune/solvers/norms/sumnorm.hh>
#include <dune/tectonic/utils/reductionfactors.hh>
@@ -40,6 +41,49 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIte
solverParset_(parset.sub("v.solver")),
errorNorms_(errorNorms) {}
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+bool FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::displacementCriterion(const Updaters& updaters, std::vector<Vector>& last_u) const {
+ bool criterion = true;
+
+ std::vector<Vector> u;
+ updaters.rate_->extractDisplacement(u);
+
+ const auto& mats = updaters.rate_->getMatrices();
+
+ for (int i=0; i<u.size(); i++) {
+ const EnergyNorm<Matrix, Vector> ANorm(*mats.elasticity[i]), MNorm(*mats.mass[i]);
+ const SumNorm<Vector> AMNorm(1.0, ANorm, 1.0, MNorm);
+
+ if (u[i].size()==0 || last_u[i].size()==0)
+ continue;
+
+ if (AMNorm.diff(u[i], last_u[i]) >= fixedPointTolerance_) {
+ criterion = false;
+ break;
+ }
+ }
+
+ last_u = u;
+
+ return criterion;
+}
+
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+bool FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::stateCriterion(const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& newAlpha) const {
+ bool criterion = true;
+ for (int i=0; i<alpha.size(); i++) {
+ if (alpha[i].size()==0 || newAlpha[i].size()==0)
+ continue;
+
+ if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
+ criterion = false;
+ break;
+ }
+ }
+
+ return criterion;
+}
+
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
template <class LinearSolver>
FixedPointIterationCounter
@@ -48,16 +92,14 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
- //std::cout << "FixedPointIterator::run()" << std::endl;
+ // debugging
+ /*const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
+ for (size_t i=0; i<contactCouplings.size(); i++) {
+ print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
+ }*/
const auto nBodies = nBodyAssembler_.nGrids();
- // debugging
- /*const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
- for (size_t i=0; i<contactCouplings.size(); i++) {
- print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
- }*/
-
FunctionalFactory<std::decay_t<decltype(nBodyAssembler_)>, decltype(globalFriction_), Matrix, Vector> functionalFactory(nBodyAssembler_, globalFriction_, velocityMatrices, velocityRHSs);
functionalFactory.build();
auto functional = functionalFactory.functional();
@@ -66,91 +108,60 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
size_t fixedPointIteration;
size_t multigridIterations = 0;
+
+ std::vector<Vector> last_u;
+ updaters.rate_->extractOldDisplacement(last_u);
+
std::vector<ScalarVector> alpha(nBodies);
updaters.state_->extractAlpha(alpha);
- Vector totalVelocityIterate, old_v;
- nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
- old_v = totalVelocityIterate;
+ std::vector<Vector> v_o;
+ updaters.rate_->extractOldVelocity(v_o);
+
+ Vector total_v_o;
+ nBodyAssembler_.nodalToTransformed(v_o, total_v_o);
+
+ Vector total_v;
+ nBodyAssembler_.nodalToTransformed(velocityIterates, total_v);
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) {
- //print(alpha, "alpha:");
-
- // contribution from nonlinearity
+ // update friction
globalFriction_.updateAlpha(alpha);
- auto res = solver.solve(solverParset_, totalVelocityIterate);
-
- nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
- //nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
-
- //print(totalVelocityIterate, "totalVelocityIterate:");
- //print(velocityIterates, "velocityIterates:");
-
- //DUNE_THROW(Dune::Exception, "Just need to stop here!");
-
+ // solve rate problem
+ auto res = solver.solve(solverParset_, total_v);
multigridIterations += res.iterations;
- Vector v_m = old_v;
+ Vector v_m = total_v_o;
v_m *= 1.0 - lambda_;
- Dune::MatrixVector::addProduct(v_m, lambda_, totalVelocityIterate);
+ Dune::MatrixVector::addProduct(v_m, lambda_, total_v);
// extract relative velocities in mortar basis
std::vector<Vector> v_rel;
split(v_m, v_rel);
- //print(v_m, "v_m: ");
-
- //print(v_rel, "v_rel");
-
- //std::cout << "- State problem set" << std::endl;
-
- // solve a state problem
+ // solve state problem
updaters.state_->solve(v_rel);
-
- //std::cout << "-- Solved" << std::endl;
-
std::vector<ScalarVector> newAlpha(nBodies);
updaters.state_->extractAlpha(newAlpha);
- //print(newAlpha, "new alpha:");
-
- bool breakCriterion = true;
- for (int i=0; i<nBodies; i++) {
- if (alpha[i].size()==0 || newAlpha[i].size()==0)
- continue;
-
- //print(alpha[i], "alpha i:");
- //print(newAlpha[i], "new alpha i:");
- if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
- breakCriterion = false;
- //std::cout << "fixedPoint error: " << errorNorms_[i]->diff(alpha[i], newAlpha[i]) << std::endl;
- break;
- }
- }
-
- //std::cout << fixedPointIteration << std::endl;
+ nBodyAssembler_.postprocess(total_v, velocityIterates);
+ updaters.rate_->postProcess(velocityIterates);
+ bool breakCriterion = displacementCriterion(updaters, last_u); //stateCriterion(alpha, newAlpha);
if (lambda_ < 1e-12 or breakCriterion) {
- //std::cout << "-FixedPointIteration finished! " << (lambda_ < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
fixedPointIteration++;
break;
}
alpha = newAlpha;
}
- //TODO: recently added, might be wrong or superfluous
- globalFriction_.updateAlpha(alpha);
-
- //print(alpha, "alpha: ");
-
- //std::cout << "-FixedPointIteration finished with " << fixedPointIteration << " iterations, lambda " << lambda_ << "! " << std::endl;
-
if (fixedPointIteration == fixedPointMaxIterations_)
DUNE_THROW(Dune::Exception, "FPI failed to converge");
+ nBodyAssembler_.postprocess(total_v, velocityIterates);
updaters.rate_->postProcess(velocityIterates);
// Cannot use return { fixedPointIteration, multigridIterations };
diff --git a/dune/tectonic/spatial-solving/fixedpointiterator.hh b/dune/tectonic/spatial-solving/fixedpointiterator.hh
index 3742ed987aefad6bd569a10611e6bbf3e7892106..d9a844e7db52ad7326935c5b55357d8532304f05 100644
--- a/dune/tectonic/spatial-solving/fixedpointiterator.hh
+++ b/dune/tectonic/spatial-solving/fixedpointiterator.hh
@@ -45,6 +45,8 @@ class FixedPointIterator {
private:
void split(const Vector& v, std::vector<Vector>& v_rel) const;
+ bool displacementCriterion(const Updaters& updaters, std::vector<Vector>& last_u) const;
+ bool stateCriterion(const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& newAlpha) const;
public:
FixedPointIterator(const Dune::ParameterTree& parset,
diff --git a/dune/tectonic/time-stepping/rate/rateupdater.cc b/dune/tectonic/time-stepping/rate/rateupdater.cc
index 3b2d9272ec5a769d6ace38006a1fff547d0ba17c..d6cde5688423e7bf6e09f60cd96ba3ce62d80692 100644
--- a/dune/tectonic/time-stepping/rate/rateupdater.cc
+++ b/dune/tectonic/time-stepping/rate/rateupdater.cc
@@ -59,6 +59,12 @@ void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractAccel
acceleration = a;
}
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+const Matrices<Matrix,2>& RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::getMatrices() const {
+ return matrices;
+}
+
+
#include "backward_euler.cc"
#include "newmark.cc"
#include "rateupdater_tmpl.cc"
diff --git a/dune/tectonic/time-stepping/rate/rateupdater.hh b/dune/tectonic/time-stepping/rate/rateupdater.hh
index 8a332206f1de86e96b6b47b65fb82ae3ff46fa55..7020b546a3a538c7eab07160afebe5784c559544 100644
--- a/dune/tectonic/time-stepping/rate/rateupdater.hh
+++ b/dune/tectonic/time-stepping/rate/rateupdater.hh
@@ -36,6 +36,8 @@ class RateUpdater {
void extractOldVelocity(std::vector<Vector>& velocity) const;
void extractAcceleration(std::vector<Vector>& acceleration) const;
+ const Matrices<Matrix,2>& getMatrices() const;
+
std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> virtual clone() const = 0;
protected: