diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 1e8f679d15ac84017312e38b4f7a4df33d800834..22a2ddfe452e08650f608f86f744e04d51ab3ed6 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -10,11 +10,17 @@
#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+
#include <dune/tectonic/bodydata.hh>
#include "../assemblers.hh"
#include "contactnetwork.hh"
#include "matrices.hh"
+#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
#include "../spatial-solving/solverfactory.hh"
#include "../spatial-solving/tnnmg/functional.hh"
#include "../spatial-solving/tnnmg/zerononlinearity.hh"
@@ -51,6 +57,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
using LocalVector = typename Vector::block_type;
//using LocalMatrix = typename Matrix::block_type;
const static int dims = LocalVector::dimension;
+ using BitVector = Dune::BitSetVector<dims>;
public:
ProgramState(const std::vector<size_t>& leafVertexCounts)
@@ -85,14 +92,43 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
}
// Set up initial conditions
- template <class GridType>
- void setupInitialConditions(const Dune::ParameterTree& parset, const ContactNetwork<GridType, Vector>& contactNetwork) {
- using Matrix = typename ContactNetwork<GridType, Vector>::Matrix;
+ template <class ContactNetwork>
+ void setupInitialConditions(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
+ using Matrix = typename ContactNetwork::Matrix;
const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+ // set patch preconditioner for linear correction in TNNMG method
+ using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;
+
+ const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
+
+ auto gsStep = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+ PatchSolver patchSolver(gsStep,
+ preconditionerParset.get<size_t>("maximumIterations"),
+ preconditionerParset.get<double>("tolerance"),
+ nullptr,
+ preconditionerParset.get<Solver::VerbosityMode>("verbosity"),
+ false); // absolute error
+
+ Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+ Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<MPPMode>("mode"));
+ preconditioner.setPatchSolver(std::forward<PatchSolver>(patchSolver));
+ preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+ preconditioner.build();
+
+ using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
+
+ LinearSolverStep linearSolverStep;
+ linearSolverStep.setPreconditioner(preconditioner);
+
+ EnergyNorm<Matrix, Vector> energyNorm(linearSolverStep);
+ LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), energyNorm, Solver::QUIET);
+
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
- const Dune::BitSetVector<dims>& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+ const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
Dune::ParameterTree const &_localParset) {
@@ -149,14 +185,15 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
print(lower, "lower");
print(upper, "upper");*/
- // set up functional and TNMMG solver
+ // set up functional
using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper); //TODO
- // using Factory = SolverFactory<Functional, Dune::BitSetVector<dims>>;
- // Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
+ // set up TNMMG solver
+ using Factory = SolverFactory<Functional, BitVector>;
+ Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
- /* std::vector<Dune::BitSetVector<dims>> bodyDirichletNodes;
+ /* std::vector<BitVector> bodyDirichletNodes;
nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
@@ -166,7 +203,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
print(totalX, "totalX: ");
print(totalRhs, "totalRhs: ");*/
- /* auto tnnmgStep = factory.step();
+ auto tnnmgStep = factory.step();
factory.setProblem(totalX);
const EnergyNorm<Matrix, Vector> norm(bilinearForm);
@@ -184,7 +221,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
Vector res = totalRhs;
bilinearForm.mmv(tnnmgStep->getSol(), res);
- std::cout << "TNNMG Res - energy norm: " << norm.operator ()(res) << std::endl; */
+ std::cout << "TNNMG Res - energy norm: " << norm.operator ()(res) << std::endl;
// TODO: remove after debugging
/* using DeformedGridType = typename LevelContactNetwork<GridType, dims>::DeformedGridType;
@@ -233,7 +270,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Initial displacement: Start from a situation of minimal stress,
// which is automatically attained in the case [v = 0 = a].
// Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
- Dune::BitSetVector<dims> dirichletNodes;
+ BitVector dirichletNodes;
contactNetwork.totalNodes("dirichlet", dirichletNodes);
/*for (size_t i=0; i<dirichletNodes.size(); i++) {
bool val = false;
@@ -262,7 +299,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Initial normal stress
const auto& body = contactNetwork.body(i);
- std::vector<std::shared_ptr<typename ContactNetwork<GridType, Vector>::LeafBody::BoundaryCondition>> frictionBoundaryConditions;
+ std::vector<std::shared_ptr<typename ContactNetwork::LeafBody::BoundaryCondition>> frictionBoundaryConditions;
body->boundaryConditions("friction", frictionBoundaryConditions);
for (size_t j=0; j<frictionBoundaryConditions.size(); j++) {
ScalarVector frictionBoundaryStress(weightedNormalStress[i].size());
@@ -277,7 +314,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
}
- Dune::BitSetVector<dims> noNodes(dirichletNodes.size(), false);
+ BitVector noNodes(dirichletNodes.size(), false);
solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a,
parset.sub("a0.solver"));
diff --git a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
index a8216abd0169da350d749f166e868b4c709f7fb6..0b0f6b59cb7de6490be85202f89dc28dc7052335 100644
--- a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
@@ -8,6 +8,7 @@
#include <dune/common/bitsetvector.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+#include <dune/solvers/common/numproc.hh>
#include "../../data-structures/levelcontactnetwork.hh"
@@ -60,7 +61,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
patchFactory_(levelContactNetwork_, fineContactNetwork_) {
setPatchDepth();
- this->verbosity_ = QUIET;
+ this->verbosity_ = NumProc::QUIET;
}
// build support patches
@@ -74,7 +75,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
// init local fine level corrections
Dune::Timer timer;
- if (this->verbosity_ == FULL) {
+ if (this->verbosity_ == NumProc::FULL) {
std::cout << std::endl;
std::cout << "---------------------------------------------" << std::endl;
std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
@@ -106,14 +107,14 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
}
}
- if (this->verbosity_ == FULL) {
+ if (this->verbosity_ == NumProc::FULL) {
std::cout << std::endl;
std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
std::cout << "---------------------------------------------" << std::endl;
timer.stop();
}
- if (this->verbosity_ != QUIET) {
+ if (this->verbosity_ != NumProc::QUIET) {
std::cout << "Level: " << level_ << " LevelPatchPreconditioner::build() - SUCCESS" << std::endl;
}
}
@@ -126,7 +127,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
patchDepth_ = patchDepth;
}
- void setVerbosity(VerbosityMode verbosity) {
+ void setVerbosity(NumProc::VerbosityMode verbosity) {
this->verbosity_ = verbosity;
}
diff --git a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index 9bab0af77022e1af2c08cc6cac40d027ca00cf74..6e4ef3bf8ea2c3b3c50e67e37580e50b16849c3b 100644
--- a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -12,15 +12,16 @@
#include "nbodycontacttransfer.hh"
#include "levelpatchpreconditioner.hh"
+enum MPPMode {additive, multiplicative};
+
template <class ContactNetwork, class PatchSolver, class MatrixType, class VectorType>
class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
public:
- enum Mode {additive, multiplicative};
+ //struct Mode { enum ModeType {additive, multiplicative}; };
private:
using LevelContactNetwork = typename ContactNetwork::LevelContactNetwork;
using LevelPatchPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
- using LevelGlobalPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
using BitVector = typename LinearIterationStep<MatrixType, VectorType>::BitVector;
using MGTransfer = NBodyContactTransfer<ContactNetwork, VectorType>;
@@ -30,7 +31,7 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
const ContactNetwork& contactNetwork_;
const Dune::BitSetVector<1>& activeLevels_;
- const Mode mode_;
+ const MPPMode mode_;
std::vector<std::shared_ptr<LevelPatchPreconditioner>> levelPatchPreconditioners_;
@@ -39,15 +40,16 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
std::vector<VectorType> levelRhs_;
std::vector<BitVector> levelIgnore_;
- PatchSolver coarseSolver_;
+ std::unique_ptr<PatchSolver> coarseSolver_;
std::vector<BitVector> recompute_;
+ std::vector<BitVector> truncation_;
std::vector<std::shared_ptr<MGTransfer>> mgTransfer_; //TODO
public:
MultilevelPatchPreconditioner(const ContactNetwork& contactNetwork,
const Dune::BitSetVector<1>& activeLevels,
- const Mode mode = MultilevelPatchPreconditioner::Mode::additive) :
+ const MPPMode mode = additive) :
contactNetwork_(contactNetwork),
activeLevels_(activeLevels),
mode_(mode) {
@@ -55,15 +57,12 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
assert(activeLevels.size() == contactNetwork_.nLevels());
// init level patch preconditioners and multigrid transfer
- levelPatchPreconditioners_.resize(0);
+ levelPatchPreconditioners_.resize(1);
+ levelPatchPreconditioners_[0] = nullptr; // blank level representing global coarse problem to make further indexing easier
int maxLevel = -1;
for (size_t i=0; i<activeLevels_.size(); i++) {
if (activeLevels_[i][0]) {
- // init global problem on coarsest level
- const auto& levelNetwork = *contactNetwork_.level(i);
- coarseGlobalPreconditioner_ = std::make_shared<LevelGlobalPreconditioner>(levelNetwork);
-
maxLevel = i;
break;
}
@@ -71,16 +70,16 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
const int coarsestLevel = maxLevel;
- typename LevelPatchPreconditioner::Mode levelMode = LevelPatchPreconditioner::Mode::additive;
+ MPPMode levelMode = additive;
if (mode_ != additive){
- levelMode = LevelPatchPreconditioner::Mode::multiplicative;
+ levelMode = multiplicative;
}
for (size_t i=maxLevel+1; i<activeLevels_.size(); i++) {
if (activeLevels_[i][0]) {
// init local patch preconditioners on each level
const auto& fineNetwork = contactNetwork_.level(i);
- levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditioner(*contactNetwork_.level(maxLevel), fineNetwork, levelMode));
+ levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditioner>(*contactNetwork_.level(maxLevel), fineNetwork, levelMode));
maxLevel = i;
}
@@ -90,14 +89,14 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
levelRhs_.resize(size());
// init multigrid transfer
- mgTransfer_.resize(levelPatchPreconditioners_.size());
+ mgTransfer_.resize(levelPatchPreconditioners_.size()-1);
for (size_t i=0; i<mgTransfer_.size(); i++) {
mgTransfer_[i] = std::make_shared<MGTransfer>();
- mgTransfer_.setup(contactNetwork_, levelPatchPreconditioners_[i]->level(), levelPatchPreconditioners_[i]->fineLevel());
+ mgTransfer_[i]->setup(contactNetwork_, levelPatchPreconditioners_[i+1]->level(), levelPatchPreconditioners_[i+1]->fineLevel());
}
//TODO
- mgTransfer_.push_back(std::make_shared<MGTransfer>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
+ //mgTransfer_.push_back(std::make_shared<MGTransfer>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
// Remove any recompute filed so that initially the full transferoperator is assembled
for (size_t i=0; i<mgTransfer_.size(); i++)
@@ -106,49 +105,51 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
// specify the subset of entries to be reassembled at each iteration
recompute_.resize(size());
recompute_[recompute_.size()-1] = contactNetwork_.nBodyAssembler().totalHasObstacle_;
- for (size_t i=mgTransfer_.size(); i>=1; i--)
- mgTransfer_[i-1]->restrict(recompute_[i], recompute_[i-1]);
-
- for (size_t i=0; i<this->mgTransfer_.size(); i++)
+ for (int i=mgTransfer_.size()-1; i>=0; i--) {
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
+ // std::dynamic_pointer_cast<DenseMultigridTransfer<VectorType, BitVector, MatrixType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(&recompute_[i]);
-
+ }
}
- void setIgnore(const BitVector& i) override {
- this->setIgnore(i);
+ void setIgnore(const BitVector& ignore) {
+ this->setIgnore(ignore);
- auto mgTransfer = std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1]);
- mgTransfer->setCriticalBitField(&critical);
+ truncation_.resize(size());
+ truncation_[size()-1] = ignore;
+ for (int i=mgTransfer_.size()-1; i>=0; i--) {
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(truncation_[i+1], truncation_[i]);
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setCriticalBitField(&truncation_[i]);
+ }
}
void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override {
this->setProblem(mat, x, rhs);
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ size_t maxLevel = levelPatchPreconditioners_.size()-1;
+ levelX_[maxLevel] = x;
+ levelRhs_[maxLevel] = rhs;
+ levelPatchPreconditioners_[maxLevel]->setProblem(mat, levelX_[maxLevel], levelRhs_[maxLevel]);
+
+ for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
const auto& transfer = *mgTransfer_[i];
- auto& _x = levelX_[i];
- auto& _rhs = levelRhs_[i];
- auto& _mat = levelMat_[i];
- transfer.restrict(x, _x);
- transfer.restrict(rhs, _rhs);
- transfer.galerkinRestrict(mat, _mat);
+ transfer.restrict(levelX_[i+1], levelX_[i]);
+ transfer.restrict(levelRhs_[i+1], levelRhs_[i]);
+ transfer.galerkinRestrict(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
- levelPatchPreconditioners_[i]->setProblem(_mat, _x, _rhs);
+ levelPatchPreconditioners_[i]->setProblem(levelMat_[i], levelX_[i], levelRhs_[i]);
}
// set coarse global problem
- const auto& coarseTransfer = *mgTransfer_.back();
- auto& coarseX = levelX_.back();
- auto& coarseRhs = levelRhs_.back();
- auto& coarseMat = levelMat_.back();
+ const auto& coarseTransfer = *mgTransfer_[0];
- coarseTransfer.restrict(x, coarseX);
- coarseTransfer.restrict(rhs, coarseRhs);
- coarseTransfer.galerkinRestrict(mat, coarseMat);
+ coarseTransfer.restrict(levelX_[1], levelX_[0]);
+ coarseTransfer.restrict(levelRhs_[1], levelRhs_[0]);
+ coarseTransfer.galerkinRestrict(levelMat_[1], levelMat_[0]);
- coarseSolver_.getIterationStep().setProblem(coarseMat, coarseX, coarseRhs);
+ coarseSolver_.getIterationStep().setProblem(levelMat_[0], levelX_[0], levelRhs_[0]);
}
@@ -237,19 +238,18 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
}
void build() {
- coarseGlobalPreconditioner_->build();
-
for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
levelPatchPreconditioners_[i]->build();
}
}
- void setPatchSolver(PatchSolver&& patchSolver) {
+ void setPatchSolver(PatchSolver&& patchSolver) { //TODO create copies of patch solver for each level
for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
levelPatchPreconditioners_[i]->setPatchSolver(std::forward<PatchSolver>(patchSolver));
}
- coarseGlobalPreconditioner_->setPatchSolver(std::forward<PatchSolver>(patchSolver));
- };
+
+ coarseSolver_ = Dune::Solvers::wrap_own_share<PatchSolver>(std::forward<PatchSolver>(patchSolver));
+ }
void setPatchDepth(const size_t patchDepth = 0) {
for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
@@ -258,7 +258,7 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
}
size_t size() const {
- return levelPatchPreconditioners_.size()+1;
+ return levelPatchPreconditioners_.size();
}
};
diff --git a/src/spatial-solving/preconditioners/nbodycontacttransfer.cc b/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
index 0f0a58c5236c55a7b1b64d53fa37d9e0645f376b..8fba1420da21e8fabda0c1459b29a0a4ece5121e 100644
--- a/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
+++ b/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
@@ -1,6 +1,7 @@
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/istl/bdmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
#include <dune/fufem/assemblers/transferoperatorassembler.hh>
#include <dune/matrix-vector/blockmatrixview.hh>
@@ -41,7 +42,7 @@ void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwor
} else {
// gridTransfer is identity if coarse and fine level coincide for body
- BDMatrix<MatrixBlock>* newMatrix = new BDMatrix<MatrixBlock>(coarseContactNetwork.body(i)->nVertices());
+ Dune::BDMatrix<MatrixBlock>* newMatrix = new Dune::BDMatrix<MatrixBlock>(coarseContactNetwork.body(i)->nVertices());
for (size_t j=0; j<newMatrix->N(); j++)
for (int k=0; k<blocksize; k++)
@@ -75,7 +76,7 @@ void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwor
const auto& contactCouplings = nBodyAssembler.getContactCouplings();
const std::vector<const MatrixType*> mortarTransferOperators(nBodyAssembler.nCouplings());
- const std::vector<const BitSetVector<1>*> fineHasObstacle(nBodyAssembler.nCouplings());
+ const std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nBodyAssembler.nCouplings());
const std::vector<std::array<int,2> > gridIdx(nBodyAssembler.nCouplings());
for (size_t i=0; i<nBodyAssembler.nCouplings(); i++) {
@@ -173,11 +174,11 @@ void NBodyContactTransfer<VectorType>::setupHierarchy(std::vector<std::shared_pt
}
*/
-template<class VectorType>
-void NBodyContactTransfer<VectorType>::combineSubMatrices(const std::vector<const MatrixType*>& submat,
+template<class ContactNetwork, class VectorType>
+void NBodyContactTransfer<ContactNetwork, VectorType>::combineSubMatrices(const std::vector<const MatrixType*>& submat,
const std::vector<const MatrixType*>& mortarTransferOperator,
const CoordSystemVector& fineLocalCoordSystems,
- const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+ const std::vector<const Dune::BitSetVector<1>*>& fineHasObstacle,
const std::vector<std::array<int,2> >& gridIdx)
{
// ///////////////////////////////////
@@ -189,7 +190,7 @@ void NBodyContactTransfer<VectorType>::combineSubMatrices(const std::vector<cons
Dune::MatrixVector::BlockMatrixView<MatrixType> view(submat);
- MatrixIndexSet totalIndexSet(view.nRows(), view.nCols());
+ Dune::MatrixIndexSet totalIndexSet(view.nRows(), view.nCols());
// import indices of canonical transfer operator
for (size_t i=0; i<nGrids; i++)
diff --git a/src/spatial-solving/preconditioners/nbodycontacttransfer.hh b/src/spatial-solving/preconditioners/nbodycontacttransfer.hh
index 7cc8ae010efccb6b58026845c1e30126f0133414..61d4e7c06d5c157786fccf01964ad95018284910 100644
--- a/src/spatial-solving/preconditioners/nbodycontacttransfer.hh
+++ b/src/spatial-solving/preconditioners/nbodycontacttransfer.hh
@@ -58,12 +58,8 @@ class NBodyContactTransfer : public TruncatedDenseMGTransfer<VectorType> {
* \param fineHasObstacle Bitfields determining for each coupling which fine grid nodes belong to the nonmortar boundary.
* \param gridIdx For each coupling store the indices of the nonmortar and mortar grid.
*/
- template <class GridType>
- void setup(const std::vector<const GridType*>& grids, int colevel,
- const std::vector<const MatrixType*>& mortarTransferOper,
- const CoordSystemVector& fineLocalCoordSystems,
- const std::vector<const BitSetVector<1>*>& fineHasObstacle,
- const std::vector<std::array<int,2> >& gridIdx);
+ void setup(const ContactNetwork& contactNetwork, const int coarseLevel, const int fineLevel);
+
protected:
/** \brief Combine all the standard prolongation and mortar operators to one big matrix.
@@ -78,7 +74,7 @@ class NBodyContactTransfer : public TruncatedDenseMGTransfer<VectorType> {
void combineSubMatrices(const std::vector<const MatrixType*>& submat,
const std::vector<const MatrixType*>& mortarTransferOperator,
const CoordSystemVector& fineLocalCoordSystems,
- const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+ const std::vector<const Dune::BitSetVector<1>*>& fineHasObstacle,
const std::vector<std::array<int,2> >& gridIdx);
};
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 6196d3ea69f3207149285cbfc357a7add0793995..7b6875ce48eb1f4cfec7303629cf0e4d236261fd 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -11,22 +11,18 @@
#include "../utils/debugutils.hh"
template <class Functional, class BitVector>
-template <class Preconditioner>
+template <class LinearSolver>
SolverFactory<Functional, BitVector>::SolverFactory(
const Dune::ParameterTree& parset,
Functional& J,
- Preconditioner&& preconditioner,
+ LinearSolver&& linearSolver,
const BitVector& ignoreNodes) :
J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))) {
nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
- linearSolverStep_ = std::make_shared<LinearSolverStep>();
- linearSolverStep_.setPreconditioner(std::forward<Preconditioner>(preconditioner));
- energyNorm_ = std::make_shared<EnergyNorm<Matrix, Vector>>(*linearSolverStep_);
- linearSolver_ = std::make_shared<LinearSolver>(*linearSolverStep_, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), *energyNorm_, Solver::QUIET);
-
- tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_, DefectProjection(), LineSearchSolver());
+ auto linearSolver = Dune::Solvers::wrap_own_share<LinearSolver>(std::forward<LinearSolver>(linearSolver));
+ tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver, DefectProjection(), LineSearchSolver());
tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
tnnmgStep_->setIgnore(ignoreNodes);
}
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 458b54d7b98331bf1a01fad804e96c23020fbab8..05235798fd608925bded68104f37e11f66c0fe19 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -13,7 +13,6 @@
//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
#include <dune/tnnmg/projections/obstacledefectprojection.hh>
-#include <dune/tnnmg/linearsolvers/fastamgsolver.hh>
#include "tnnmg/tnnmgstep.hh"
#include "tnnmg/linearization.hh"
@@ -35,13 +34,11 @@ class SolverFactory {
using Step = typename Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
- using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
- using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
- template <class Preconditioner>
+ template <class LinearSolver>
SolverFactory(Dune::ParameterTree const &,
Functional&,
- Preconditioner&& preconditioner,
+ LinearSolver&&,
const BitVector&);
void setProblem(Vector& x);
@@ -55,11 +52,6 @@ class SolverFactory {
// nonlinear smoother
std::shared_ptr<NonlinearSmoother> nonlinearSmoother_;
- // linear solver
- std::shared_ptr<LinearSolverStep> linearSolverStep_;
- std::shared_ptr<EnergyNorm<Matrix, Vector>> energyNorm_;
- std::shared_ptr<LinearSolver> linearSolver_;
-
// TNNMGStep
std::shared_ptr<Step> tnnmgStep_;
};
diff --git a/src/spatial-solving/solverfactory_tmpl.cc b/src/spatial-solving/solverfactory_tmpl.cc
index 1355a1fe359e43beebf60264207a0f9e10977222..8e3b2a153e3ec91877c0c6d2b8f78f4e19d601cc 100644
--- a/src/spatial-solving/solverfactory_tmpl.cc
+++ b/src/spatial-solving/solverfactory_tmpl.cc
@@ -5,6 +5,8 @@
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
+#include <dune/solvers/solvers/loopsolver.hh>
+
#include "../../dune/tectonic/globalfriction.hh"
#include "tnnmg/functional.hh"
#include "tnnmg/zerononlinearity.hh"
@@ -15,5 +17,16 @@ using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Ve
using MyZeroFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, double>;
using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
+using MyLinearSolver = Dune::Solvers::LoopSolver<Vector, BitVector>;
+
template class SolverFactory<MyFunctional, BitVector>;
+template SolverFactory<MyFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
+ MyFunctional&,
+ MyLinearSolver&&,
+ const BitVector&);
+
template class SolverFactory<MyZeroFunctional, BitVector>;
+template SolverFactory<MyZeroFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
+ MyZeroFunctional&,
+ MyLinearSolver&&,
+ const BitVector&);