diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc index fc4b2f35adfe4a91147683a1791f7cd95af0d55e..d66d8cc19d984c492aa035f452c613fef274a1e2 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 c479c32dc6bc830eded15e6a59e7d80bc799a860..81be4536caf79567e61f0eeab7df9994518f183e 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 2e004e1a7da9c89e8d8149cb4446c96dbc4fac34..ed1fe9b4957fddfb26291f99042dedc708fac353 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 ca39e910c456518c7ffc3a110c4becf891d89194..97dfe8341850194d93d6cba4dca335331c3342d6 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 c68ba188e06f05eecdae9b1d639dc13527c38914..faf14f1c281a5509fee6c4d824c9e2ca7183fc70 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 2c4bc74ed0c5ec9c1e63dcf50f3a0befe1a7377c..932ad8db5256c66ecf2738df778f101e8a1df183 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 506d9042a6814dab7d4a4ad43647f81262b7091a..7e06390bb45510eb2a9b4b209a23db01d29a4b1c 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 d6bce3d59ca2d87354f247b8d6d608a100be02f5..2598ecaf14e377f1bcf6b8dab695e83a1403c894 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 73f75ba6b55bfcb2f5cee02111531c5e3a5206bf..5d559a8b097758f4f0ea1633eef1cf8b2a11545f 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 e369a8f5e2a77dfd50ff5ef85ca000fd84637832..97c383edff112abcd9c782cc6d50bc2860cf3e33 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 b8a4129a891806ee3130814ceb3670998a7c5fa1..410a48f8d65812124e3edbeceb0b5700a015ddae 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 d99f191c02a5adf02d17c3fb8654cc084692e9e1..52543d1e8dd2b9204dd6fd978ba26aa2f833be0b 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 946c282d006a71096347c627031699851a3d8629..a2bd2dfa43a144219f95c7793de0520f85fea1c5 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 6e2a986306b407d6bf1ef31c86b65baacf4d459b..5266259d256d027d3bd3ebe9ba6ad32b4ab0702d 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;