From 90d7a0c77ea171a227b1de3724f8ecba0f342e7b Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Mon, 19 Mar 2018 17:16:11 +0100
Subject: [PATCH] .
---
src/multi-body-problem.cc | 2 +-
src/spatial-solving/fixedpointiterator.cc | 37 ++++++++++++++-----
src/spatial-solving/fixedpointiterator.hh | 11 +++---
.../fixedpointiterator_tmpl.cc | 2 +-
src/spatial-solving/solverfactory.cc | 32 +++++++++-------
src/spatial-solving/solverfactory.hh | 7 +++-
src/spatial-solving/solverfactory_tmpl.cc | 10 ++---
src/time-stepping/adaptivetimestepper.cc | 2 +-
src/time-stepping/adaptivetimestepper.hh | 8 ++--
src/time-stepping/adaptivetimestepper_tmpl.cc | 2 +-
src/time-stepping/coupledtimestepper.cc | 2 +-
src/time-stepping/coupledtimestepper.hh | 7 ++--
src/time-stepping/coupledtimestepper_tmpl.cc | 6 +--
src/time-stepping/rate/rateupdater.hh | 2 +-
14 files changed, 77 insertions(+), 53 deletions(-)
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index fc4b2f35..d66d8cc1 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -480,7 +480,7 @@ int main(int argc, char *argv[]) {
// Set up TNNMG solver
- using NonlinearFactory = SolverFactory<DeformedGrid, Matrix, Vector>;
+ using NonlinearFactory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, dirichletNodes);
using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index c479c32d..81be4536 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -36,8 +36,9 @@ void FixedPointIterationCounter::operator+=(
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
Factory &factory, Dune::ParameterTree const &parset,
- std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
- : step_(factory.getStep()),
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, const ErrorNorm& errorNorm)
+ : factory_(factory),
+ step_(factory_.getStep()),
parset_(parset),
globalFriction_(globalFriction),
fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
@@ -51,27 +52,43 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterationCounter
FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
- Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+ Updaters updaters, const std::vector<const Matrix*>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
- //EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
+ 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);
+
+ // contribution from nonlinearity
+ // add gradient to rhs, hessian to matrix!?
+
size_t fixedPointIteration;
size_t multigridIterations = 0;
std::vector<ScalarVector> alpha;
updaters.state_->extractAlpha(alpha);
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) {
+
// solve a velocity problem
- // globalFriction_->updateAlpha(alpha);
- //ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
- // velocityRHSs, velocityIterate);
- //BlockProblem velocityProblem(parset_, convexProblem);
- //step_->setProblem(velocityIterate, velocityProblem);
- //step_->setProblem(velocityIterate);
+ for (size_t i=0; i<alpha.size(); i++) {
+ globalFriction_[i]->updateAlpha(alpha[i]);
+ }
+
+ step_->setProblem(bilinearForm, totalVelocityIterate, totalRhs);
+
velocityProblemSolver.preprocess();
velocityProblemSolver.solve();
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index 2e004e1a..ed1fe9b4 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -25,6 +25,7 @@ class FixedPointIterator {
using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
+ using Nonlinearity = typename Factory::Nonlinearity;
// using Nonlinearity = typename ConvexProblem::NonlinearityType;
@@ -34,20 +35,20 @@ class FixedPointIterator {
void relativeVelocities(std::vector<Vector>& v_m) const;
public:
- FixedPointIterator(Factory &factory, const Dune::ParameterTree& parset,
- std::shared_ptr<Nonlinearity> globalFriction,
+ FixedPointIterator(Factory& factory, const Dune::ParameterTree& parset,
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
const ErrorNorm& errorNorm);
FixedPointIterationCounter run(Updaters updaters,
- const std::vector<Matrix>& velocityMatrices,
+ const std::vector<const Matrix*>& velocityMatrices,
const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates);
private:
-
+ Factory& factory_;
std::shared_ptr<typename Factory::Step> step_;
Dune::ParameterTree const &parset_;
- std::shared_ptr<Nonlinearity> globalFriction_;
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
size_t fixedPointMaxIterations_;
double fixedPointTolerance_;
diff --git a/src/spatial-solving/fixedpointiterator_tmpl.cc b/src/spatial-solving/fixedpointiterator_tmpl.cc
index ca39e910..97dfe834 100644
--- a/src/spatial-solving/fixedpointiterator_tmpl.cc
+++ b/src/spatial-solving/fixedpointiterator_tmpl.cc
@@ -21,7 +21,7 @@
#include "solverfactory.hh"
using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<DeformedGrid, Matrix, Vector>;
+using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index c68ba188..faf14f1c 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -13,8 +13,8 @@
#include "solverfactory.hh"
-template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
-SolverFactory<DeformedGrid, Matrix, Vector>::SolverFactory(
+template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
+SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
const std::vector<Dune::BitSetVector<DeformedGrid::dimension>>& ignoreNodes)
: nBodyAssembler_(nBodyAssembler),
@@ -27,11 +27,11 @@ SolverFactory<DeformedGrid, Matrix, Vector>::SolverFactory(
// tnnmg iteration step
multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
- multigridStep_->setIgnore(ignoreNodes);
+ //multigridStep_->setIgnore(ignoreNodes);
multigridStep_->setBaseSolver(baseSolver_);
multigridStep_->setSmoother(&presmoother_, &postsmoother_);
- multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
- multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
+ // multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
+ // multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
// create the transfer operators
const int nCouplings = nBodyAssembler_.nCouplings();
@@ -43,11 +43,12 @@ SolverFactory<DeformedGrid, Matrix, Vector>::SolverFactory(
std::vector<const Matrix*> mortarTransfer(nCouplings);
std::vector<std::array<int,2> > gridIdx(nCouplings);
+ /*
for (int j=0; j<nCouplings; j++) {
- coarseHasObstacle[j] = nBodyAssembler_.nonmortarBoundary_[j][i].getVertices();
+ coarseHasObstacle[j] = nBodyAssembler_.contactCoupling_[j]->nonmortarBoundary_[i].getVertices();
fineHasObstacle[j] = nBodyAssembler_.nonmortarBoundary_[j][i+1].getVertices();
- mortarTransfer[j] = &nBodyAssembler_.contactCoupling_[j].mortarLagrangeMatrix(i);
+ mortarTransfer[j] = &nBodyAssembler_.contactCoupling_[j]->mortarLagrangeMatrix(i);
gridIdx[j] = nBodyAssembler_.coupling_[j].gridIdx_;
}
@@ -57,20 +58,25 @@ SolverFactory<DeformedGrid, Matrix, Vector>::SolverFactory(
nBodyAssembler_.localCoordSystems_[i+1],
coarseHasObstacle, fineHasObstacle,
mortarTransfer,
- gridIdx);
+ gridIdx); */
}
-
+/*
+ grids[0], *grids[1], i,
+ contactAssembler.contactCoupling_[0]->mortarLagrangeMatrix(),
+ contactAssembler.localCoordSystems_,
+ *contactAssembler.contactCoupling_[0]->nonmortarBoundary().getVertices());
+*/
multigridStep_->setTransferOperators(transferOperators_);
}
-template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
-SolverFactory<DeformedGrid, Matrix, Vector>::~SolverFactory() {
+template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
+SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::~SolverFactory() {
for (auto &&x : transferOperators_)
delete x;
}
-template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
-auto SolverFactory<DeformedGrid, Matrix, Vector>::getStep()
+template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
+auto SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::getStep()
-> std::shared_ptr<Step> {
return multigridStep_;
}
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 2c4bc74e..932ad8db 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -23,7 +23,7 @@
#define USE_OLD_TNNMG //needed for old bisection.hh in tnnmg
-template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
+template <class DeformedGridTEMPLATE, class NonlinearityTEMPLATE, class MatrixType, class VectorType>
class SolverFactory {
//protected:
// const size_t dim = DeformedGrid::dimension;
@@ -33,6 +33,7 @@ class SolverFactory {
using Matrix = MatrixType;
using DeformedGrid = DeformedGridTEMPLATE;
+ using Nonlinearity = NonlinearityTEMPLATE;
public:
using Step = Dune::Contact::NonSmoothNewtonMGStep<Matrix, Vector>;
@@ -44,6 +45,10 @@ class SolverFactory {
std::shared_ptr<Step> getStep();
+ const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& getNBodyAssembler() const {
+ return nBodyAssembler_;
+ }
+
private:
const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler_;
diff --git a/src/spatial-solving/solverfactory_tmpl.cc b/src/spatial-solving/solverfactory_tmpl.cc
index 506d9042..7e06390b 100644
--- a/src/spatial-solving/solverfactory_tmpl.cc
+++ b/src/spatial-solving/solverfactory_tmpl.cc
@@ -12,11 +12,9 @@
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/myblockproblem.hh>
-template class SolverFactory<
- MY_DIM,
- MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
- Grid>;
-template class SolverFactory<
+template class SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+
+/*template class SolverFactory<
MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
- Grid>;
+ Grid>;*/
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index d6bce3d5..2598ecaf 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -20,7 +20,7 @@ void IterationRegister::reset() {
template <class Factory, class Updaters, class ErrorNorm>
AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
Factory &factory, Dune::ParameterTree const &parset,
- std::shared_ptr<Nonlinearity> globalFriction, Updaters ¤t,
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters ¤t,
double relativeTime, double relativeTau,
std::function<void(double, Vector &)> externalForces,
ErrorNorm const &errorNorm,
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index 73f75ba6..5d559a8b 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -23,14 +23,14 @@ class AdaptiveTimeStepper {
};
using Vector = typename Factory::Vector;
- using ConvexProblem = typename Factory::ConvexProblem;
- using Nonlinearity = typename ConvexProblem::NonlinearityType;
+ //using ConvexProblem = typename Factory::ConvexProblem;
+ using Nonlinearity = typename Factory::Nonlinearity;
using MyCoupledTimeStepper = CoupledTimeStepper<Factory, Updaters, ErrorNorm>;
public:
AdaptiveTimeStepper(Factory &factory, Dune::ParameterTree const &parset,
- std::shared_ptr<Nonlinearity> globalFriction,
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
Updaters ¤t, double relativeTime,
double relativeTau,
std::function<void(double, Vector &)> externalForces,
@@ -50,7 +50,7 @@ class AdaptiveTimeStepper {
double finalTime_;
Factory &factory_;
Dune::ParameterTree const &parset_;
- std::shared_ptr<Nonlinearity> globalFriction_;
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
Updaters ¤t_;
UpdatersWithCount R1_;
std::function<void(double, Vector &)> externalForces_;
diff --git a/src/time-stepping/adaptivetimestepper_tmpl.cc b/src/time-stepping/adaptivetimestepper_tmpl.cc
index e369a8f5..97c383ed 100644
--- a/src/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/src/time-stepping/adaptivetimestepper_tmpl.cc
@@ -19,7 +19,7 @@
#include "updaters.hh"
using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<DeformedGrid, Matrix, Vector>;
+using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index b8a4129a..410a48f8 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -7,7 +7,7 @@
template <class Factory, class Updaters, class ErrorNorm>
CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
double finalTime, Factory &factory, Dune::ParameterTree const &parset,
- std::shared_ptr<Nonlinearity> globalFriction, Updaters updaters,
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters updaters,
ErrorNorm const &errorNorm,
std::function<void(double, Vector &)> externalForces)
: finalTime_(finalTime),
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index d99f191c..52543d1e 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -12,13 +12,12 @@ template <class Factory, class Updaters, class ErrorNorm>
class CoupledTimeStepper {
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
- using ConvexProblem = typename Factory::ConvexProblem;
- using Nonlinearity = typename ConvexProblem::NonlinearityType;
+ using Nonlinearity = typename Factory::Nonlinearity;
public:
CoupledTimeStepper(double finalTime, Factory &factory,
Dune::ParameterTree const &parset,
- std::shared_ptr<Nonlinearity> globalFriction,
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
Updaters updaters, ErrorNorm const &errorNorm,
std::function<void(double, Vector &)> externalForces);
@@ -28,7 +27,7 @@ class CoupledTimeStepper {
double finalTime_;
Factory &factory_;
Dune::ParameterTree const &parset_;
- std::shared_ptr<Nonlinearity> globalFriction_;
+ std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
Updaters updaters_;
std::function<void(double, Vector &)> externalForces_;
ErrorNorm const &errorNorm_;
diff --git a/src/time-stepping/coupledtimestepper_tmpl.cc b/src/time-stepping/coupledtimestepper_tmpl.cc
index 946c282d..a2bd2dfa 100644
--- a/src/time-stepping/coupledtimestepper_tmpl.cc
+++ b/src/time-stepping/coupledtimestepper_tmpl.cc
@@ -19,10 +19,8 @@
#include "updaters.hh"
using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
- MY_DIM,
- MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
- Grid>;
+using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
diff --git a/src/time-stepping/rate/rateupdater.hh b/src/time-stepping/rate/rateupdater.hh
index 6e2a9863..5266259d 100644
--- a/src/time-stepping/rate/rateupdater.hh
+++ b/src/time-stepping/rate/rateupdater.hh
@@ -20,7 +20,7 @@ class RateUpdater {
void virtual setup(Vector const &ell, double _tau, double relativeTime,
Vector &rhs, Vector &iterate, Matrix &AB) = 0;
- void virtual postProcess(Vector const &iterate) = 0;
+ void virtual postProcess(const std::vector<Vector>& iterate) = 0;
void extractDisplacement(std::vector<Vector>& displacements) const;
void extractVelocity(std::vector<Vector>& velocity) const;
void extractOldVelocity(std::vector<Vector>& velocity) const;
--
GitLab