From b40abd4330487b226490c61afbbcd0ac85289734 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Thu, 4 Jul 2019 18:53:42 +0200
Subject: [PATCH] .
---
src/CMakeLists.txt | 2 +-
src/data-structures/contactnetwork.cc | 2 +-
src/data-structures/program_state.hh | 10 +-
src/factories/stackedblocksfactory.cc | 208 +++++++-----
src/factories/stackedblocksfactory.hh | 49 +--
src/factories/stackedblocksfactory_tmpl.cc | 3 +-
src/multi-body-problem-2D.cfg | 3 +
src/multi-body-problem.cc | 51 ++-
src/multi-body-problem.cfg | 10 +-
src/spatial-solving/fixedpointiterator.cc | 6 +-
src/spatial-solving/fixedpointiterator.hh | 2 +
.../preconditioners/CMakeLists.txt | 1 +
.../levelglobalpreconditioner.hh | 132 ++------
.../levelpatchpreconditioner.hh | 78 ++---
.../multilevelpatchpreconditioner.hh | 209 ++++++------
.../preconditioners/nbodycontacttransfer.cc | 297 ++++++++++++++++++
.../preconditioners/nbodycontacttransfer.hh | 88 ++++++
.../preconditioners/supportpatchfactory.hh | 10 +-
src/spatial-solving/solverfactory.cc | 16 +-
src/spatial-solving/solverfactory.hh | 22 +-
src/time-stepping/adaptivetimestepper.cc | 18 +-
src/time-stepping/adaptivetimestepper.hh | 7 +-
src/time-stepping/coupledtimestepper.cc | 7 +-
src/time-stepping/coupledtimestepper.hh | 3 +-
24 files changed, 788 insertions(+), 446 deletions(-)
create mode 100644 src/spatial-solving/preconditioners/nbodycontacttransfer.cc
create mode 100644 src/spatial-solving/preconditioners/nbodycontacttransfer.hh
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 6bf43fe7..26276d41 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -33,7 +33,7 @@ set(MSW_SOURCE_FILES
data-structures/enumparser.cc
#factories/cantorfactory.cc
factories/threeblocksfactory.cc
- #factories/stackedblocksfactory.cc
+ factories/stackedblocksfactory.cc
io/vtk.cc
io/hdf5/frictionalboundary-writer.cc
io/hdf5/iteration-writer.cc
diff --git a/src/data-structures/contactnetwork.cc b/src/data-structures/contactnetwork.cc
index 81d9f8a9..4e66cf76 100644
--- a/src/data-structures/contactnetwork.cc
+++ b/src/data-structures/contactnetwork.cc
@@ -24,7 +24,7 @@ ContactNetwork<HostGridType, VectorType>::ContactNetwork(
levelContactNetworks_(nLevels),
bodies_(nBodies),
couplings_(nCouplings),
- frictionBoundaries_(nCouplings),
+ frictionBoundaries_(nBodies),
stateEnergyNorms_(nBodies),
nBodyAssembler_(nBodies, nCouplings)
{}
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 9c5dd4af..1e8f679d 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -151,10 +151,10 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// set up functional and TNMMG solver
using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
- Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper);
+ Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper); //TODO
- using Factory = SolverFactory<Functional, Dune::BitSetVector<dims>>;
- Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
+ // using Factory = SolverFactory<Functional, Dune::BitSetVector<dims>>;
+ // Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
/* std::vector<Dune::BitSetVector<dims>> bodyDirichletNodes;
nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
@@ -166,7 +166,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 +184,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;
diff --git a/src/factories/stackedblocksfactory.cc b/src/factories/stackedblocksfactory.cc
index 5a2aaaff..9f7fc32a 100644
--- a/src/factories/stackedblocksfactory.cc
+++ b/src/factories/stackedblocksfactory.cc
@@ -14,75 +14,97 @@
#include "stackedblocksfactory.hh"
-template <class GridType, int dims>
-void StackedBlocksFactory<GridType, dims>::setBodies() {
- // set up cuboid geometries
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
+
+ // set up cuboid geometries
double const length = 1.00;
- double const width = 0.27;
- double const weakLength = 0.20;
+ double const height = 0.27;
+ double const weakLength = 0.60;
- std::vector<LocalVector> origins(this->bodyCount_);
- std::vector<LocalVector> lowerWeakOrigins(this->bodyCount_);
- std::vector<LocalVector> upperWeakOrigins(this->bodyCount_);
+ std::vector<GlobalCoords> origins(this->bodyCount_);
+
+ const auto& subParset = this->parset_.sub("boundary.friction.weakening");
#if MY_DIM == 3
double const depth = 0.60;
+ auto weakBound = [&] (const GlobalCoords& origin) {
+ GlobalCoords h = {origin[0] + weakLength, origin[1], origin[2]};
+ return h;
+ };
+
origins[0] = {0,0,0};
- lowerWeakOrigins[0] = {0.2,0,0};
- upperWeakOrigins[0] = {0.2,width,0};
+ GlobalCoords lowerWeakOrigin = {0.2, 0, 0};
+ GlobalCoords upperWeakOrigin = {0.2, height, 0};
+
+ cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height, depth);
+ cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
- cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width, depth);
for (size_t i=1; i<this->bodyCount_-1; i++) {
- origins[i] = cuboidGeometries_[i-1]->D;
- lowerWeakOrigins[i] = {0.6,i*width,0};
- upperWeakOrigins[i] = {0.6,(i+1)*width,0};
+ origins[i] = cuboidGeometries_[i-1]->upperLeft();
+ GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
+ GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];
- cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width, depth);
+ cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height, depth);
+ cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
+ cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
}
+
const size_t idx = this->bodyCount_-1;
- origins[idx] = cuboidGeometries_[idx-1]->D;
- lowerWeakOrigins[idx] = {0.6,idx*width,0};
- upperWeakOrigins[idx] = {0.6,(idx+1)*width,0};
+ origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
+ lowerWeakOrigin += origins[idx];
- cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, length, width, depth);
+ cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height, depth);
+ cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
#elif MY_DIM == 2
+ auto weakBound = [&] (const GlobalCoords& origin) {
+ GlobalCoords h = {origin[0] + weakLength, origin[1]};
+ return h;
+ };
+
origins[0] = {0,0};
- lowerWeakOrigins[0] = {0.2,0};
- upperWeakOrigins[0] = {0.6,width};
+ GlobalCoords lowerWeakOrigin = {0.2, 0};
+ GlobalCoords upperWeakOrigin = {0.2, height};
+
+ cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height);
+ cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
- cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width);
for (size_t i=1; i<this->bodyCount_-1; i++) {
- origins[i] = cuboidGeometries_[i-1]->D;
- lowerWeakOrigins[i] = {0.6,i*width};
- upperWeakOrigins[i] = {0.6,(i+1)*width};
+ origins[i] = cuboidGeometries_[i-1]->upperLeft();
+ GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
+ GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];
- cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width);
+ cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height);
+
+
+ cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
+ cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
}
+
const size_t idx = this->bodyCount_-1;
- origins[idx] = cuboidGeometries_[idx-1]->D + LocalVector({0.25,0});
- lowerWeakOrigins[idx] = {0.6,idx*width};
- upperWeakOrigins[idx] = {0.6,(idx+1)*width};
+ origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
+ lowerWeakOrigin += origins[idx];
- cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, 5.0/8*length, width);
+ cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height);
+ cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif
// set up reference grids
- gridConstructor_ = new GridsConstructor<GridType>(cuboidGeometries_);
+ gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(cuboidGeometries_);
auto& grids = gridConstructor_->getGrids();
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
// define weak patch and refine grid
- auto lowerWeakPatch = lowerWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
- auto upperWeakPatch = upperWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
- getWeakPatch<LocalVector>(this->parset_.sub("boundary.friction.weakening"), cuboidGeometry, *lowerWeakPatch, *upperWeakPatch);
- refine(*grids[i], *lowerWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
- refine(*grids[i], *upperWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
+ const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
+ for (size_t j=0; j<weakeningRegions.size(); j++) {
+ refine(*grids[i], weakeningRegions[j], this->parset_.template get<double>("boundary.friction.smallestDiameter"), CuboidGeometry::lengthScale);
+ }
// determine minDiameter and maxDiameter
double minDiameter = std::numeric_limits<double>::infinity();
@@ -98,96 +120,116 @@ void StackedBlocksFactory<GridType, dims>::setBodies() {
}
for (size_t i=0; i<this->bodyCount_; i++) {
- this->bodies_[i] = std::make_shared<typename Base::Body>(bodyData_, grids[i]);
+ this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_, grids[i]);
}
}
-template <class GridType, int dims>
-void StackedBlocksFactory<GridType, dims>::setCouplings() {
- const auto deformedGrids = this->levelContactNetwork_.deformedGrids();
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setCouplings() {
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
- leafFaces_[i] = std::make_shared<LeafFaces>(this->levelContactNetwork_.leafView(i), cuboidGeometry);
- levelFaces_[i] = std::make_shared<LevelFaces>(this->levelContactNetwork_.levelView(i), cuboidGeometry);
+
+ leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
+ levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
}
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
+
+
for (size_t i=0; i<this->couplings_.size(); i++) {
auto& coupling = this->couplings_[i];
coupling = std::make_shared<typename Base::FrictionCouplingPair>();
nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
+ weakPatches_[i] = std::make_shared<WeakeningRegion>(cuboidGeometries_[i]->weakeningRegions()[0]);
coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
- coupling->setWeakeningPatch(upperWeakPatches_[i]);
- coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *upperWeakPatches_[i], CuboidGeometry::lengthScale));
+ coupling->setWeakeningPatch(weakPatches_[i]);
+ coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale));
}
}
-template <class GridType, int dims>
-void StackedBlocksFactory<GridType, dims>::setBoundaryConditions() {
- using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setLevelBodies() {
+ const size_t nBodies = this->bodyCount_;
- using Function = Dune::VirtualFunction<double, double>;
- std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
- std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();
-
- const double lengthScale = CuboidGeometry::lengthScale;
+ size_t maxLevel = 0;
+ std::vector<size_t> maxLevels(nBodies);
+ for (size_t i=0; i<nBodies; i++) {
+ const size_t bodyMaxLevel = this->bodies_[i]->grid()->maxLevel();
+ maxLevels[i] = bodyMaxLevel;
+ maxLevel = std::max(maxLevel, bodyMaxLevel);
+ }
- for (size_t i=0; i<this->bodyCount_; i++) {
- const auto& body = this->levelContactNetwork_.body(i);
- const auto& leafVertexCount = body->leafVertexCount();
+ for (int l=maxLevel; l>=0; l--) {
+ this->contactNetwork_.addLevel(maxLevels, l);
- std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;
+ for (size_t i=0; i<nBodies; i++) {
+ maxLevels[i]--;
+ }
+ }
+}
- // Neumann boundary
- std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->leafView()), neumannFunction, "neumann");
- body->addBoundaryCondition(neumannBoundary);
- // right Dirichlet Boundary
- // identify Dirichlet nodes on leaf level
- std::shared_ptr<Dune::BitSetVector<dims>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(leafVertexCount);
- for (int j=0; j<leafVertexCount; j++) {
- if (leafFaces_[i]->upper.containsVertex(j))
- (*velocityDirichletNodes)[j][0] = true;
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setBoundaryConditions() {
+ using LeafBoundaryCondition = BoundaryCondition<LeafGridView, dim>;
- #if MY_DIM == 3 //TODO: wrong, needs revision
- if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
- zeroDirichletNodes->at(j)[2] = true;
- #endif
- }
+ using Function = Dune::VirtualFunction<double, double>;
+ std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+ std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();
- if (i>0) {
- std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ const double lengthScale = CuboidGeometry::lengthScale;
- velocityDirichletBoundary->setBoundaryPatch(body->leafView(), velocityDirichletNodes);
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ const auto& body = this->contactNetwork_.body(i);
+ const auto& leafVertexCount = body->nVertices();
+
+ // Neumann boundary
+ std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
+ body->addBoundaryCondition(neumannBoundary);
+
+ // upper Dirichlet Boundary
+ // identify Dirichlet nodes on leaf level
+ if (i==this->bodyCount_-1) {
+ std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
+ for (int j=0; j<leafVertexCount; j++) {
+ if (leafFaces_[i]->upper.containsVertex(j))
+ (*velocityDirichletNodes)[j][0] = true;
+
+ #if MY_DIM == 3 //TODO: wrong, needs revision
+ if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
+ zeroDirichletNodes->at(j)[2] = true;
+ #endif
+ }
- velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
- body->addBoundaryCondition(velocityDirichletBoundary);
- } else {
- //velocityDirichletBoundary->setBoundaryFunction(neumannFunction);
- }
+ std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ velocityDirichletBoundary->setBoundaryPatch(body->gridView(), velocityDirichletNodes);
+ velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+ body->addBoundaryCondition(velocityDirichletBoundary);
+ } else {
+ }
}
// lower Dirichlet Boundary
- const auto& firstBody = this->levelContactNetwork_.body(0);
- const auto& firstLeafVertexCount = firstBody->leafVertexCount();
- std::shared_ptr<Dune::BitSetVector<dims>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(firstLeafVertexCount);
+ const auto& firstBody = this->contactNetwork_.body(0);
+ const auto& firstLeafVertexCount = firstBody->nVertices();
+ std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(firstLeafVertexCount);
for (int j=0; j<firstLeafVertexCount; j++) {
if (leafFaces_[0]->lower.containsVertex(j)) {
- for (size_t d=0; d<dims; d++) {
+ for (size_t d=0; d<dim; d++) {
(*zeroDirichletNodes)[j][d] = true;
}
}
}
std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
- zeroDirichletBoundary->setBoundaryPatch(firstBody->leafView(), zeroDirichletNodes);
+ zeroDirichletBoundary->setBoundaryPatch(firstBody->gridView(), zeroDirichletNodes);
std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
diff --git a/src/factories/stackedblocksfactory.hh b/src/factories/stackedblocksfactory.hh
index ed38c0ff..b9471853 100644
--- a/src/factories/stackedblocksfactory.hh
+++ b/src/factories/stackedblocksfactory.hh
@@ -7,38 +7,46 @@
#include <dune/fufem/boundarypatch.hh>
-#include "levelcontactnetworkfactory.hh"
+#include "contactnetworkfactory.hh"
#include "../multi-body-problem-data/mybody.hh"
#include "../multi-body-problem-data/grid/mygrids.hh"
+#include "../multi-body-problem-data/grid/cuboidgeometry.hh"
-template <class GridType, int dims> class StackedBlocksFactory : public LevelContactNetworkFactory<GridType, dims>{
+template <class HostGridType, class VectorType> class StackedBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
private:
- using Base = LevelContactNetworkFactory<GridType, dims>;
+ using Base = ContactNetworkFactory<HostGridType, VectorType>;
- using LocalVector = typename Base::LevelContactNetwork::LocalVector;
+public:
+ using ContactNetwork = typename Base::ContactNetwork;
+
+private:
+ using GlobalCoords = typename ContactNetwork::LocalVector;
- using DeformedLeafGridView = typename Base::LevelContactNetwork::DeformedLeafGridView;
- using DeformedLevelGridView = typename Base::LevelContactNetwork::DeformedLevelGridView;
+ using LeafGridView = typename ContactNetwork::GridView;
+ using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
- using LevelBoundaryPatch = BoundaryPatch<DeformedLevelGridView>;
- using LeafBoundaryPatch = BoundaryPatch<DeformedLeafGridView>;
+ using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+ using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
- using LeafFaces = MyFaces<DeformedLeafGridView>;
- using LevelFaces = MyFaces<DeformedLevelGridView>;
+ using LeafFaces = MyFaces<LeafGridView>;
+ using LevelFaces = MyFaces<LevelGridView>;
-private:
- const std::shared_ptr<MyBodyData<dims>> bodyData_; // material properties of bodies
+ using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
+ using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
- GridsConstructor<GridType>* gridConstructor_;
+ static const int dim = ContactNetwork::dim;
+
+ const std::shared_ptr<MyBodyData<dim>> bodyData_; // material properties of bodies
+
+ std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
- std::vector<std::shared_ptr<ConvexPolyhedron<LocalVector>>> lowerWeakPatches_;
- std::vector<std::shared_ptr<ConvexPolyhedron<LocalVector>>> upperWeakPatches_;
+ std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
@@ -46,22 +54,19 @@ template <class GridType, int dims> class StackedBlocksFactory : public LevelCon
public:
StackedBlocksFactory(const Dune::ParameterTree& parset) :
Base(parset, parset.get<size_t>("problem.bodyCount"), parset.get<size_t>("problem.bodyCount")-1),
- bodyData_(std::make_shared<MyBodyData<dims>>(this->parset_, zenith_())),
+ bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_, zenith_())),
cuboidGeometries_(this->bodyCount_),
leafFaces_(this->bodyCount_),
levelFaces_(this->bodyCount_),
- lowerWeakPatches_(this->bodyCount_),
- upperWeakPatches_(this->bodyCount_),
+ weakPatches_(this->bodyCount_),
nonmortarPatch_(this->couplingCount_),
mortarPatch_(this->couplingCount_)
{}
- ~StackedBlocksFactory() {
- delete gridConstructor_;
- }
-
void setBodies();
+ void setLevelBodies();
void setCouplings();
+ void setLevelCouplings() {}
void setBoundaryConditions();
private:
diff --git a/src/factories/stackedblocksfactory_tmpl.cc b/src/factories/stackedblocksfactory_tmpl.cc
index c72f5199..ea0c3039 100644
--- a/src/factories/stackedblocksfactory_tmpl.cc
+++ b/src/factories/stackedblocksfactory_tmpl.cc
@@ -3,5 +3,6 @@
#endif
#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
-template class StackedBlocksFactory<Grid, MY_DIM>;
+template class StackedBlocksFactory<Grid, Vector>;
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index a308c80b..c0b758a8 100644
--- a/src/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem-2D.cfg
@@ -19,3 +19,6 @@ tolerance = 1e-5
[solver.tnnmg.linear]
tolerance = 1e-10
+
+[solver.tnnmg.linear.preconditioner]
+tolerance = 1e-10
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 51f691ad..6682374d 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -33,12 +33,8 @@
#include <dune/fufem/hdf5/file.hh>
#include <dune/solvers/norms/energynorm.hh>
-
-
-/*
#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/solvers/solver.hh>*/
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include <dune/contact/common/deformedcontinuacomplex.hh>
#include <dune/contact/common/couplingpair.hh>
@@ -47,7 +43,6 @@
#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/globalfriction.hh>
-
#include "assemblers.hh"
#include "gridselector.hh"
#include "explicitgrid.hh"
@@ -59,6 +54,7 @@
#include "data-structures/matrices.hh"
#include "data-structures/program_state.hh"
+#include "factories/stackedblocksfactory.hh"
#include "factories/threeblocksfactory.hh"
//#include "io/hdf5-levelwriter.hh"
@@ -70,7 +66,7 @@
#include "multi-body-problem-data/grid/mygrids.hh"
#include "spatial-solving/tnnmg/functional.hh"
-#include "spatial-solving/preconditioners/supportpatchfactory.hh"
+//#include "spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
#include "spatial-solving/tnnmg/localbisectionsolver.hh"
#include "spatial-solving/contact/nbodyassembler.hh"
#include "spatial-solving/solverfactory.hh"
@@ -151,15 +147,20 @@ int main(int argc, char *argv[]) {
using Assembler = MyAssembler<DefLeafGridView, dims>;
using field_type = Matrix::field_type;
-
// ----------------------
// set up contact network
// ----------------------
- ThreeBlocksFactory<Grid, Vector> threeBlocksFactory(parset);
+ StackedBlocksFactory<Grid, Vector> stackedBlocksFactory(parset);
+ using ContactNetwork = typename StackedBlocksFactory<Grid, Vector>::ContactNetwork;
+ stackedBlocksFactory.build();
+
+ ContactNetwork& contactNetwork = stackedBlocksFactory.contactNetwork();
+
+ /*ThreeBlocksFactory<Grid, Vector> threeBlocksFactory(parset);
using ContactNetwork = typename ThreeBlocksFactory<Grid, Vector>::ContactNetwork;
threeBlocksFactory.build();
- ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork();
+ ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork(); */
const size_t bodyCount = contactNetwork.nBodies();
for (size_t i=0; i<bodyCount; i++) {
@@ -191,7 +192,7 @@ int main(int argc, char *argv[]) {
const auto & coarseContactNetwork = *contactNetwork.level(3);
const auto & fineContactNetwork = *contactNetwork.level(4);
- SupportPatchFactory<typename ContactNetwork::LevelContactNetwork> supportPatchFactory(coarseContactNetwork, fineContactNetwork);
+ // SupportPatchFactory<typename ContactNetwork::LevelContactNetwork> supportPatchFactory(coarseContactNetwork, fineContactNetwork);
size_t bodyID = 2;
size_t patchDepth = 3;
@@ -223,6 +224,7 @@ int main(int argc, char *argv[]) {
writeToVTK(gv, dofs, "", "body_" + std::to_string(i) + "_fine");
}
+ /*
using Patch = typename SupportPatchFactory<typename ContactNetwork::LevelContactNetwork>::Patch;
Patch patch;
@@ -263,7 +265,7 @@ int main(int argc, char *argv[]) {
}
}
}
- }
+ } */
std::cout << "Done! :)" << std::endl;
return 1;
@@ -506,6 +508,27 @@ int main(int argc, char *argv[]) {
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
+ /*
+ // 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<Preconditioner::Mode>("mode"));
+ preconditioner.setPatchSolver(patchSolver);
+ preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+
+ // set adaptive time stepper
typename ContactNetwork::ExternalForces externalForces;
contactNetwork.externalForces(externalForces);
@@ -517,7 +540,8 @@ int main(int argc, char *argv[]) {
while (!adaptiveTimeStepper.reachedEnd()) {
programState.timeStep++;
- iterationCount = adaptiveTimeStepper.advance();
+ preconditioner.build();
+ iterationCount = adaptiveTimeStepper.advance(preconditioner);
programState.relativeTime = adaptiveTimeStepper.relativeTime_;
programState.relativeTau = adaptiveTimeStepper.relativeTau_;
@@ -545,6 +569,7 @@ int main(int argc, char *argv[]) {
break;
}
}
+ */
std::cout.rdbuf(coutbuf); //reset to standard output again
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index 5fa9847b..8678a5dd 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -11,7 +11,7 @@ vtk.write = true
[problem]
finalTime = 100 # [s] #1000
-bodyCount = 2
+bodyCount = 4
[body]
bulkModulus = 0.5e5 # [Pa]
@@ -69,8 +69,16 @@ maximumIterations = 100
pre = 3
cycle = 1 # 1 = V, 2 = W, etc.
post = 3
+[solver.tnnmg.linear.preconditioner]
+mode = additive
+patchDepth = 0
+maximumIterations = 10000
+verbosity = quiet
[solver.tnnmg.main]
pre = 1
multi = 5 # number of multigrid steps
post = 0
+
+
+
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 2c741b77..99456a3b 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -59,9 +59,11 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIte
errorNorms_(errorNorms) {}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Preconditioner>
FixedPointIterationCounter
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
- Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+ Updaters updaters, Preconditioner&& preconditioner,
+ const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
std::cout << "FixedPointIterator::run()" << std::endl;
@@ -124,7 +126,7 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
// set up functional and TNMMG solver
Functional J(bilinearForm, totalRhs, globalFriction_, vLower, vUpper);
- Factory solverFactory(parset_.sub("solver.tnnmg"), J, ignoreNodes_);
+ Factory solverFactory(parset_.sub("solver.tnnmg"), J, std::forward<Preconditioner>(preconditioner), ignoreNodes_);
auto step = solverFactory.step();
std::cout << "- Functional and TNNMG step setup: success" << std::endl;
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index 867e8058..c7c74189 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -52,7 +52,9 @@ class FixedPointIterator {
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
const ErrorNorms& errorNorms);
+ template <class Preconditioner>
FixedPointIterationCounter run(Updaters updaters,
+ Preconditioner&& preconditioner,
const std::vector<Matrix>& velocityMatrices,
const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates);
diff --git a/src/spatial-solving/preconditioners/CMakeLists.txt b/src/spatial-solving/preconditioners/CMakeLists.txt
index abb4ff76..c5f8a1b5 100644
--- a/src/spatial-solving/preconditioners/CMakeLists.txt
+++ b/src/spatial-solving/preconditioners/CMakeLists.txt
@@ -3,5 +3,6 @@ add_custom_target(dune-tectonic_spatial-solving_preconditioners_src SOURCES
levelglobalpreconditioner.hh
levelpatchpreconditioner.hh
multilevelpatchpreconditioner.hh
+ nbodycontacttransfer.hh
supportpatchfactory.hh
)
diff --git a/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh b/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh
index ac9fd3fc..e593c60d 100644
--- a/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh
@@ -3,6 +3,7 @@
#include <string>
+#include <dune/common/timer.hh>
#include <dune/common/fvector.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
@@ -10,105 +11,61 @@
//#include <dune/istl/superlu.hh>
#include <dune/istl/umfpack.hh>
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
+#include "../../data-structures/levelcontactnetwork.hh"
-#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+template <class LevelContactNetwork, class Solver, class MatrixType, class VectorType>
class LevelGlobalPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-public:
- enum BoundaryMode {homogeneous, fromIterate};
-
private:
- typedef typename BasisType::GridView GridView;
- typedef typename GridView::Grid GridType;
+ using Base = LinearIterationStep<MatrixType, VectorType>;
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
- const LocalAssembler& localAssembler_;
- const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
+ using ctype = typename LevelContactNetwork::ctype;
- const GridType& grid_;
+ const LevelContactNetwork& levelContactNetwork_;
const int level_;
- const BasisType basis_;
-
- Dune::BitSetVector<1> boundaryDofs_;
- MatrixType matrix_;
- VectorType rhs_;
- bool requireBuild_;
+ std::shared_ptr<Solver> solver_;
public:
- // for each active fault in levelInterfaceNetwork: set up local problem, compute local correction
- LevelGlobalPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
- const LocalAssembler& localAssembler,
- const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers) :
- levelInterfaceNetwork_(levelInterfaceNetwork),
- localAssembler_(localAssembler),
- localInterfaceAssemblers_(localInterfaceAssemblers),
- grid_(levelInterfaceNetwork_.grid()),
- level_(levelInterfaceNetwork_.level()),
- basis_(levelInterfaceNetwork_)
- {
- assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
- }
-
- void build() {
- //printBasisDofLocation(basis_);
+ LevelGlobalPreconditioner(const LevelContactNetwork& levelContactNetwork, const int level = 0) :
+ levelContactNetwork_(levelContactNetwork),
+ level_(level) {
- GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
- globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
+ this->verbosity_ = QUIET;
+ }
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
+ void setSolver(Solver&& patchSolver) {
+ solver_ = Dune::Solvers::wrap_own_share<PatchSolver>(std::forward<PatchSolver>(patchSolver));
+ }
- const GridView& gridView = levelInterfaceNetwork_.levelGridView();
+ void iterate() override {
+ Dune::Timer timer;
- // set boundary conditions
- BoundaryPatch<GridView> boundaryPatch(gridView, true);
- constructBoundaryDofs(boundaryPatch, basis_, boundaryDofs_);
+ if (this->verbosity_ == FULL) {
+ timer.reset();
+ timer.start();
+ }
- for(size_t i=0; i<boundaryDofs_.size(); i++) {
- if(!boundaryDofs_[i][0])
- continue;
+ *(this->x_) = 0;
- RowType& row = matrix_[i];
+ auto& step = solver_->getIterationStep();
+ step.setProblem(this->mat_, this->x_, this->rhs_);
+ step.setIgnore(this->ignoreNodes_);
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
+ solver_->check();
+ solver_->preprocess();
+ solver_->solve();
- for(; cIt!=cEndIt; ++cIt) {
- row[cIt.index()] = 0;
- }
+ *(this->x_) = solver_->getSol();
- row[i] = 1;
+ if (this->verbosity_ == FULL) {
+ std::cout << "LevelGlobalPreconditioner::iterate() Solving global problem took: " << timer.stop() << " seconds" << std::endl;
}
-
- requireBuild_ = false;
- }
-
-
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- updateRhs(rhs);
- this->mat_ = Dune::stackobject_to_shared_ptr(mat); // mat will be set but not used
}
- virtual void iterate() {
- //print(matrix_, "coarse matrix: ");
- //print(rhs_, "coarse rhs: ");
-
- // compute solution directly
-
- if (requireBuild_)
- DUNE_THROW(Dune::Exception, "LevelGlobalPreconditioner::iterate() Call build() before solving the global problem!");
-
+ /*void iterate() override {
Dune::Timer timer;
timer.reset();
timer.start();
@@ -129,31 +86,8 @@ class LevelGlobalPreconditioner : public LinearIterationStep<MatrixType, VectorT
//std::cout << "LevelGlobalPreconditioner::iterate() Solving global problem took: " << timer.elapsed() << " seconds" << std::endl;
- }
-
- const BasisType& basis() const {
- return basis_;
- }
-
- const GridView& gridView() const {
- return basis_.getGridView();
- }
+ }*/
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
- return levelInterfaceNetwork_;
- }
-
-private:
- void updateRhs(const VectorType& rhs){
- rhs_.resize(matrix_.M());
- rhs_ = 0;
-
- for (size_t i=0; i<rhs.size(); i++) {
- if (!boundaryDofs_[i][0])
- rhs_[i] = rhs[i];
- }
- }
};
-
#endif
diff --git a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
index 1c8257c1..52d5f1d1 100644
--- a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
@@ -43,15 +43,10 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
PatchFactory patchFactory_;
std::vector<Patch> patches_;
- PatchSolver patchSolver_;
-
-
+ std::shared_ptr<PatchSolver> patchSolver_;
size_t patchDepth_;
- MatrixType matrix_;
- std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
-
public:
// for each coarse patch given by levelContactNetwork: set up local problem, compute local correction
@@ -78,12 +73,13 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
patches_.resize(totalNVertices);
// init local fine level corrections
+ Dune::Timer timer;
if (this->verbosity_ == FULL) {
std::cout << std::endl;
std::cout << "---------------------------------------------" << std::endl;
std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
- Dune::Timer timer;
+
timer.reset();
timer.start();
}
@@ -91,6 +87,8 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
Dune::BitSetVector<1> vertexVisited(totalNVertices);
vertexVisited.unsetAll();
+ const auto& levelIndices = patchFactory_.coarseIndices();
+
for (size_t bodyIdx=0; bodyIdx<levelContactNetwork_.nBodies(); bodyIdx++) {
const auto& gridView = levelContactNetwork_.body(bodyIdx)->gridView();
@@ -98,16 +96,10 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
const auto& refElement = Dune::ReferenceElements<double, dim>::general(e.type());
for (size_t i=0; i<refElement.size(dim); i++) {
- auto globalIdx = patchFactory_.coarseIndices().vertexIndex(bodyIdx, e, i);
+ auto globalIdx = levelIndices.vertexIndex(bodyIdx, e, i);
if (!vertexVisited[globalIdx][0]) {
vertexVisited[globalIdx][0] = true;
-
- /* if (this->verbosity_ != QUIET) {
- std::cout << i << std::endl;
- std::cout << "---------------" << std::endl;
- }*/
-
patchFactory_.build(bodyIdx, e, i, patches_[globalIdx], patchDepth_);
}
}
@@ -121,11 +113,15 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
timer.stop();
}
- if (this->verbosity_ != REDUCED) {
+ if (this->verbosity_ != QUIET) {
std::cout << "Level: " << level_ << " LevelPatchPreconditioner::build() - SUCCESS" << std::endl;
}
}
+ void setPatchSolver(PatchSolver&& patchSolver) {
+ patchSolver_ = Dune::Solvers::wrap_own_share<PatchSolver>(std::forward<PatchSolver>(patchSolver));
+ }
+
void setPatchDepth(const size_t patchDepth = 0) {
patchDepth_ = patchDepth;
}
@@ -134,24 +130,6 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
this->verbosity_ = verbosity;
}
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
-
- /* VectorType rhs;
- rhs.resize(matrix_.N());
- rhs = 0;
-
- //print(localToGlobal, "localToGlobal: ");
- //print(boundaryDofs, "boundaryDofs: ");
-
- localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs);*/
-
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- patchSolver_.setProblem(mat, x, rhs);
- }
-
virtual void iterate() {
if (mode_ == additive)
iterateAdd();
@@ -161,43 +139,37 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
void iterateAdd() {
*(this->x_) = 0;
-
VectorType x = *(this->x_);
- for (size_t i=0; i<patches_.size(); i++) {
- patchSolver_.setIgnore(patches_[i]);
- patchSolver_.solve();
- *(this->x_) += x;
+ for (const auto& p : patches_) {
+ x = 0;
+
+ auto& step = patchSolver_->getIterationStep();
+ step.setProblem(this->mat_, x, this->rhs_);
+ step.setIgnore(p);
+
+ patchSolver_->check();
+ patchSolver_->preprocess();
+ patchSolver_->solve();
+
+ *(this->x_) += patchSolver_->getSol();
}
}
void iterateMult() {
*(this->x_) = 0;
-
+/*
VectorType x = *(this->x_);
for (size_t i=0; i<patches_.size(); i++) {
VectorType updatedRhs(*(this->rhs_));
this->mat_->mmv(*(this->x_), updatedRhs);
- //step(i);
- //print(updatedRhs, "localRhs: ");
- //writeToVTK(basis_, updatedRhs, "/storage/mi/podlesny/data/faultnetworks/rhs/local", "exactvertexdata_step_"+std::to_string(i));
-
patchSolver_.setProblem(*this->mat_, *(this->x_), updatedRhs);
patchSolver_.setIgnore(patches_[i]);
patchSolver_.solve();
- //print(it, "it: ");
- //print(x, "x: ");
-
- //writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/sol/local", "exactvertexdata_step_"+std::to_string(i));
-
- /*if (i==5) {
- writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
- }*/
-
*(this->x_) += x;
- }
+ } */
}
size_t size() const {
diff --git a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index 93529434..3d1f9afc 100644
--- a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -9,60 +9,50 @@
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-#include <dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh>
-#include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/interfacenetwork.hh>
+#include "nbodytransfer.hh"
+#include "levelpatchpreconditioner.hh"
#include <dune/faultnetworks/dgmgtransfer.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-template <class BasisType, class LocalOperatorAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+template <class ContactNetwork, class PatchSolver, class MatrixType, class VectorType>
class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
public:
enum Mode {additive, multiplicative};
- typedef LevelPatchPreconditioner<BasisType, LocalOperatorAssembler, LocalInterfaceAssembler, MatrixType, VectorType> LevelPatchPreconditionerType;
- typedef LevelGlobalPreconditioner<BasisType, LocalOperatorAssembler, LocalInterfaceAssembler, MatrixType, VectorType> LevelGlobalPreconditionerType;
private:
- typedef typename BasisType::GridView GridView;
- typedef typename GridView::Grid GridType;
-
+ 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;
- typedef typename Dune::BitSetVector<1> BitVector;
+ using MGTransfer = ;
- const ContactNetwork<GridType>& interfaceNetwork_;
- const BitVector& activeLevels_;
+ const ContactNetwork& contactNetwork_;
+ const Dune::BitSetVector<1>& activeLevels_;
const Mode mode_;
- int itCounter_;
+ std::vector<std::shared_ptr<LevelPatchPreconditioner>> levelPatchPreconditioners_;
- std::shared_ptr<LevelGlobalPreconditionerType> coarseGlobalPreconditioner_;
- std::vector<std::shared_ptr<LevelPatchPreconditionerType>> levelPatchPreconditioners_;
+ std::vector<MatrixType> levelMat_;
std::vector<VectorType> levelX_;
std::vector<VectorType> levelRhs_;
+ std::vector<BitVector> levelIgnore_;
- std::vector<std::shared_ptr<DGMGTransfer<BasisType>>> mgTransfer_;
+ PatchSolver coarseSolver_;
+
+ std::vector<std::shared_ptr<MGTransfer>> mgTransfer_; //TODO
public:
- MultilevelPatchPreconditioner(const ContactNetwork<GridType>& interfaceNetwork,
- const BitVector& activeLevels,
+ MultilevelPatchPreconditioner(const ContactNetwork& contactNetwork,
+ const Dune::BitSetVector<1>& activeLevels,
const Mode mode = MultilevelPatchPreconditioner::Mode::additive) :
- interfaceNetwork_(interfaceNetwork),
+ contactNetwork_(contactNetwork),
activeLevels_(activeLevels),
- localAssemblers_(localAssemblers),
- localInterfaceAssemblers_(localInterfaceAssemblers),
- mode_(mode)
- {
-
- if (activeLevels_.size() > (size_t) interfaceNetwork.maxLevel() +1)
- DUNE_THROW(Dune::Exception, "MultilevelFaultPreconditioner: too many active levels; preconditioner supports at most (grid.maxLevel + 1) levels!");
+ mode_(mode) {
- assert(activeLevels.size() == localAssemblers_.size() && activeLevels.size() == localInterfaceAssemblers_.size());
+ assert(activeLevels.size() == contactNetwork_.nLevels());
- // init level fault preconditioners and multigrid transfer
+ // init level patch preconditioners and multigrid transfer
levelPatchPreconditioners_.resize(0);
mgTransfer_.resize(0);
@@ -70,74 +60,73 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
for (size_t i=0; i<activeLevels_.size(); i++) {
if (activeLevels_[i][0]) {
// init global problem on coarsest level
- const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork.levelInterfaceNetwork(i);
- coarseGlobalPreconditioner_ = std::make_shared<LevelGlobalPreconditionerType>(levelNetwork, *localAssemblers_[i], localInterfaceAssemblers_[i]);
+ const auto& levelNetwork = *contactNetwork_.level(i);
+ coarseGlobalPreconditioner_ = std::make_shared<LevelGlobalPreconditioner>(levelNetwork);
maxLevel = i;
break;
}
}
- typename LevelPatchPreconditionerType::Mode levelMode = LevelPatchPreconditionerType::Mode::additive;
+ const int coarsestLevel = maxLevel;
+
+ typename LevelPatchPreconditioner::Mode levelMode = LevelPatchPreconditioner::Mode::additive;
if (mode_ != additive){
- levelMode = LevelPatchPreconditionerType::Mode::multiplicative;
+ levelMode = LevelPatchPreconditioner::Mode::multiplicative;
}
for (size_t i=maxLevel+1; i<activeLevels_.size(); i++) {
if (activeLevels_[i][0]) {
// init local patch preconditioners on each level
- const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork.levelInterfaceNetwork(i);
-
- //levelPatchPreconditioners_.push_back(std::make_shared<LevelGlobalPreconditionerType>(levelNetwork, *localAssemblers_[i], localInterfaceAssemblers_[i]));
-
- if (levelPatchPreconditioners_.size() == 0) {
- levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditionerType>(levelNetwork, coarseGlobalPreconditioner_->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], levelMode));
- } else {
- //levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditionerType>(levelNetwork, coarseGlobalPreconditioner_->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], levelMode));
-
- levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditionerType>(levelNetwork, levelPatchPreconditioners_.back()->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], levelMode));
- }
+ const auto& fineNetwork = contactNetwork_.level(i);
+ levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditioner(*contactNetwork_.level(maxLevel), fineNetwork, levelMode));
maxLevel = i;
}
}
- levelX_.resize(levelPatchPreconditioners_.size()+1);
- levelRhs_.resize(levelPatchPreconditioners_.size()+1);
-
- allFaultBasis_ = std::make_shared<BasisType>(interfaceNetwork_.levelInterfaceNetwork(maxLevel));
+ levelX_.resize(size());
+ levelRhs_.resize(size());
// init multigrid transfer
for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- mgTransfer_.push_back(std::make_shared<DGMGTransfer<BasisType>>(levelPatchPreconditioners_[i]->basis(), *allFaultBasis_));
+ mgTransfer_.push_back(std::make_shared<MGTransfer>(...));
}
- mgTransfer_.push_back(std::make_shared<DGMGTransfer<BasisType>>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
-
- itCounter_ = 0;
+ mgTransfer_.push_back(std::make_shared<MGTransfer>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
}
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
+ 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++) {
- mgTransfer_[i]->restrict(x, levelX_[i]);
- mgTransfer_[i]->restrict(rhs, levelRhs_[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);
- levelPatchPreconditioners_[i]->setProblem(mat, levelX_[i], levelRhs_[i]);
+ levelPatchPreconditioners_[i]->setProblem(_mat, _x, _rhs);
}
- size_t j = levelPatchPreconditioners_.size();
- mgTransfer_[j]->restrict(x, levelX_[j]);
- mgTransfer_[j]->restrict(rhs, levelRhs_[j]);
+ // set coarse global problem
+ const auto& coarseTransfer = *mgTransfer_.back();
+ auto& coarseX = levelX_.back();
+ auto& coarseRhs = levelRhs_.back();
+ auto& coarseMat = levelMat_.back();
+
+ coarseTransfer.restrict(x, coarseX);
+ coarseTransfer.restrict(rhs, coarseRhs);
+ coarseTransfer.galerkinRestrict(mat, coarseMat);
- coarseGlobalPreconditioner_->setProblem(mat, levelX_[j], levelRhs_[j]);
+ coarseSolver_.getIterationStep().setProblem(coarseMat, coarseX, coarseRhs);
}
- virtual void iterate() {
+ void iterate() {
if (mode_ == additive)
iterateAdd();
else
@@ -149,67 +138,52 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
VectorType x;
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- levelPatchPreconditioners_[i]->iterate();
- const VectorType& it = levelPatchPreconditioners_[i]->getSol();
+ // solve coarse global problem
+ auto& coarseIgnore = levelIgnore_.back();
+ mgTransfer_.back()->restrict(this->ignore(), coarseIgnore);
+ coarseSolver_.getIterationStep().setIgnore(coarseIgnore);
- mgTransfer_[i]->prolong(it, x);
+ coarseSolver_.check();
+ coarseSolver_.preprocess();
+ coarseSolver_.solve();
- //writeToVTK(*allFaultBasis_, x, "/home/mi/podlesny/data/faultnetworks/sol/", "preconditioner_step_"+std::to_string(i));
+ mgTransfer_.back()->prolong(coarseSolver_.getSol(), x);
- *(this->x_) += x;
- }
+ for (size_t i=levelPatchPreconditioners_.size()-2; i>=0; i--) {
+ const auto& transfer = *mgTransfer_[i];
+ auto& preconditioner = *levelPatchPreconditioners_[i];
- coarseGlobalPreconditioner_->iterate();
- const VectorType& it = coarseGlobalPreconditioner_->getSol();
+ auto& _ignore = levelIgnore_[i];
+ transfer.restrict(this->ignore(), _ignore);
+ preconditioner.setIgnore(_ignore);
- mgTransfer_[levelPatchPreconditioners_.size()]->prolong(it, x);
- *(this->x_) += x;
+ preconditioner.iterate();
+ x += preconditioner.getSol();
- /*
- levelPatchPreconditioners_[itCounter_]->iterate();
- const VectorType& it = levelPatchPreconditioners_[itCounter_]->getSol();
+ VectorType newX;
+ transfer.prolong(x, newX);
- mgTransfer_[itCounter_]->prolong(it, x);
-
- //writeToVTK(*allFaultBasis_, x, "/home/mi/podlesny/data/faultnetworks/sol/", "preconditioner_step_"+std::to_string(i));
-
- *(this->x_) += x;
-
- coarseGlobalPreconditioner_->iterate();
- const VectorType& it2 = coarseGlobalPreconditioner_->getSol();
-
- mgTransfer_[levelPatchPreconditioners_.size()]->prolong(it2, x);
- *(this->x_) += x;
-
- itCounter_++;
- itCounter_ = itCounter_ % (levelPatchPreconditioners_.size());
-
- if (itCounter_==0)
- std::cout << "Multilevel pass complete!" << std::endl;
+ x = newX;
+ }
- */
+ *(this->x_) = x;
}
void iterateMult() {
*(this->x_) = 0;
+ /*
VectorType x;
for (int i=(levelPatchPreconditioners_.size()-1); i>=0; i--) {
VectorType updatedRhs(*(this->rhs_));
this->mat_->mmv(*(this->x_), updatedRhs);
- //print(updatedRhs, "globalRhs: ");
-
mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
mgTransfer_[i]->restrict(updatedRhs, levelRhs_[i]);
- //print(levelRhs_[i], "levelLocalCoarseRhs: ");
- //writeToVTK(levelPatchPreconditioners_[i]->basis(), levelRhs_[i], "/storage/mi/podlesny/data/faultnetworks/rhs/fine", "exactvertexdata_step_"+std::to_string(itCounter_));
-
levelPatchPreconditioners_[i]->setProblem(*(this->mat_), levelX_[i], levelRhs_[i]);
levelPatchPreconditioners_[i]->iterate();
@@ -217,9 +191,6 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
mgTransfer_[i]->prolong(it, x);
- //writeToVTK(*allFaultBasis_, x, "/storage/mi/podlesny/data/faultnetworks/fineIterates/multilevel", "exactvertexdata_step_"+std::to_string(itCounter_));
-
-
*(this->x_) += x;
}
@@ -230,20 +201,13 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
mgTransfer_[j]->restrict(*(this->x_), levelX_[j]);
mgTransfer_[j]->restrict(updatedCoarseRhs, levelRhs_[j]);
- //print(levelRhs_[j], "localCoarseRhs: ");
- //writeToVTK(coarseGlobalPreconditioner_->basis(), levelRhs_[j], "/storage/mi/podlesny/data/faultnetworks/rhs/coarse", "exactvertexdata_step_"+std::to_string(itCounter_));
-
coarseGlobalPreconditioner_->setProblem(*(this->mat_), levelX_[j], levelRhs_[j]);
coarseGlobalPreconditioner_->iterate();
const VectorType& it = coarseGlobalPreconditioner_->getSol();
mgTransfer_[j]->prolong(it, x);
- //writeToVTK(*allFaultBasis_, x, "/storage/mi/podlesny/data/faultnetworks/coarseIterates/multilevel", "exactvertexdata_step_"+std::to_string(itCounter_));
-
- *(this->x_) += x;
-
- itCounter_++;
+ *(this->x_) += x; */
}
void build() {
@@ -254,12 +218,17 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
}
}
- std::shared_ptr<LevelPatchPreconditionerType> levelPatchPreconditioner(const int level) {
- return levelPatchPreconditioners_[level];
- }
+ void setPatchSolver(PatchSolver&& patchSolver) {
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ levelPatchPreconditioners_[i]->setPatchSolver(std::forward<PatchSolver>(patchSolver));
+ }
+ coarseGlobalPreconditioner_->setPatchSolver(std::forward<PatchSolver>(patchSolver));
+ };
- const BasisType& basis() const {
- return *allFaultBasis_;
+ void setPatchDepth(const size_t patchDepth = 0) {
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ levelPatchPreconditioners_[i]->setPatchDepth(patchDepth);
+ }
}
size_t size() const {
diff --git a/src/spatial-solving/preconditioners/nbodycontacttransfer.cc b/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
new file mode 100644
index 00000000..961ee3f5
--- /dev/null
+++ b/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
@@ -0,0 +1,297 @@
+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/istl/bdmatrix.hh>
+
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+#include <dune/matrix-vector/blockmatrixview.hh>
+
+#include <dune/matrix-vector/addtodiagonal.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/solvers/common/numproc.hh>
+#include <dune/solvers/transferoperators/densemultigridtransfer.hh>
+
+template<class ContactNetwork, class VectorType>
+void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwork& contactNetwork, const int coarseLevel, const int fineLevel)
+
+ (const std::vector<const GridType*>& grids, int colevel,
+ const std::vector<const MatrixType*>& mortarTransferOperator,
+ const CoordSystemVector& fineLocalCoordSystems,
+ const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+ const std::vector<std::array<int,2> >& gridIdx)
+
+transfer->setup(*grids[0], *grids[1], i,
+ contactAssembler.contactCoupling_[0]->mortarLagrangeMatrix(),
+ contactAssembler.localCoordSystems_,
+ *contactAssembler.contactCoupling_[0]->nonmortarBoundary().getVertices());
+
+{
+ const size_t nBodies = contactNetwork.nBodies();
+ const size_t nCouplings = contactNetwork.nCouplings();
+
+ const auto& coarseContactNetwork = contactNetwork.level(coarseLevel);
+ const auto& fineContactNetwork = contactNetwork.level(fineLevel);
+
+ std::vector<size_t> maxL(nBodies), cL(nBodies), fL(nBodies);
+
+ for (size_t i=0; i<nBodies; i++) {
+ maxL[i] = contactNetwork.body(i)->grid()->maxLevel();
+ cL[i] = coarseContactNetwork.body(i)->level();
+ fL[i] = fineContactNetwork.body(i)->level();
+ }
+
+ // ////////////////////////////////////////////////////////////
+ // Create the standard prolongation matrices for each grid
+ // ////////////////////////////////////////////////////////////
+
+ std::vector<TruncatedDenseMGTransfer<VectorType>* > gridTransfer(nBodies);
+ std::vector<const MatrixType*> submat(nBodies);
+
+ for (size_t i=0; i<nBodies; i++) {
+ if (fL[i] > cL[i]) {
+ gridTransfer[i] = new TruncatedDenseMGTransfer<VectorType>;
+ gridTransfer[i]->setup(*contactNetwork.body(i)->grid(), cL[i], fL[i]);
+ submat[i] = &gridTransfer[i]->getMatrix();
+
+ } else {
+ // gridTransfer is identity if coarse and fine level coincide for body
+ BDMatrix<MatrixBlock>* newMatrix = new BDMatrix<MatrixBlock>(coarseContactNetwork.body(i)->nVertices());
+
+ for (size_t j=0; j<newMatrix->N(); j++)
+ for (int k=0; k<blocksize; k++)
+ for (int l=0; l<blocksize; l++)
+ (*newMatrix)[j][j][k][l] = (k==l);
+
+ submat[i] = newMatrix;
+ }
+ }
+
+ // ///////////////////////////////////
+ // Combine the submatrices
+ // ///////////////////////////////////
+
+ if (fineLevel == contactNetwork.maxLevel())
+ combineSubMatrices(submat, mortarTransferOperator, fineLocalCoordSystems, fineHasObstacle, gridIdx);
+ else
+ Dune::MatrixVector::BlockMatrixView<MatrixType>::setupBlockMatrix(submat, this->matrix_);
+
+ for (size_t i=0; i<nBodies; i++) {
+ if (fL[i]==0) {
+ delete(submat[i]);
+ }
+ }
+
+ // Delete separate transfer objects
+ for (size_t i=0; i<nBodies; i++) {
+ delete(gridTransfer[i]);
+ }
+}
+
+template<class VectorType>
+template<class GridType>
+void NBodyContactTransfer<VectorType>::setupHierarchy(std::vector<std::shared_ptr<NonSmoothNewtonContactTransfer<VectorType> > >& mgTransfers,
+ const std::vector<const GridType*> grids,
+ const std::vector<const MatrixType*> mortarTransferOperator,
+ const CoordSystemVector& fineLocalCoordSystems,
+ const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+ const std::vector<std::array<int,2> >& gridIdx)
+{
+ const size_t nGrids = grids.size();
+
+ std::vector<std::vector<TruncatedDenseMGTransfer<VectorType>* > > gridTransfer(nGrids);
+ std::vector<const MatrixType*> submat(nGrids);
+
+ // ////////////////////////////////////////////////////////////
+ // Create the standard prolongation matrices for each grid
+ // ////////////////////////////////////////////////////////////
+
+ for (size_t i=0; i<nGrids; i++) {
+
+ gridTransfer[i].resize(grids[i]->maxLevel());
+ for (size_t j=0; j<gridTransfer[i].size(); j++)
+ gridTransfer[i][j] = new TruncatedDenseMGTransfer<VectorType>;
+
+ // Assemble standard transfer operator
+ TransferOperatorAssembler<GridType> transferOperatorAssembler(*grids[i]);
+ transferOperatorAssembler.assembleOperatorPointerHierarchy(gridTransfer[i]);
+ }
+
+ // ////////////////////////////////////////////////////////////////////////////
+ // Combine matrices in one matrix and add mortar entries on the fine level
+ // ///////////////////////////////////////////////////////////////////////////
+ int toplevel = mgTransfers.size();
+
+ for (size_t colevel=0; colevel<mgTransfers.size(); colevel++) {
+
+ std::vector<int> fL(nGrids);
+
+ for (size_t i=0; i<nGrids; i++) {
+ fL[i] = std::max(size_t(0), grids[i]->maxLevel() - colevel);
+
+ // If the prolongation matrix exists, take it
+ if (fL[i] > 0) {
+
+ submat[i] =&gridTransfer[i][fL[i]-1]->getMatrix();
+ } else {
+ // when the maxLevels of the grids differ then we add copys of the coarsest level to the "smaller" grid
+ BDMatrix<MatrixBlock>* newMatrix = new BDMatrix<MatrixBlock>(grids[i]->size(0,GridType::dimension));
+ *newMatrix = 0;
+
+ for (size_t j=0; j<newMatrix->N(); j++)
+ for (int k=0; k<blocksize; k++)
+ (*newMatrix)[j][j][k][k] = 1.0;
+
+ submat[i] = newMatrix;
+ }
+ }
+
+ if (colevel == 0)
+ mgTransfers[toplevel-colevel-1]->combineSubMatrices(submat, mortarTransferOperator, fineLocalCoordSystems, fineHasObstacle, gridIdx);
+ else
+ Dune::MatrixVector::BlockMatrixView<MatrixType>::setupBlockMatrix(submat, mgTransfers[toplevel-colevel-1]->matrix_);
+
+ for (size_t i=0; i<nGrids; i++)
+ if (fL[i]==0)
+ delete(submat[i]);
+ }
+
+ // Delete separate transfer objects
+ for (size_t i=0; i<nGrids; i++)
+ for (size_t j=0; j<gridTransfer[i].size(); j++)
+ delete(gridTransfer[i][j]);
+}
+
+
+template<class VectorType>
+void NBodyContactTransfer<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<std::array<int,2> >& gridIdx)
+{
+ // ///////////////////////////////////
+ // Combine the submatrices
+ // ///////////////////////////////////
+
+ const size_t nGrids = submat.size();
+ const size_t nCouplings = mortarTransferOperator.size();
+
+ Dune::MatrixVector::BlockMatrixView<MatrixType> view(submat);
+
+ MatrixIndexSet totalIndexSet(view.nRows(), view.nCols());
+
+ // import indices of canonical transfer operator
+ for (size_t i=0; i<nGrids; i++)
+ totalIndexSet.import(*submat[i], view.row(i, 0), view.col(i, 0));
+
+ // ///////////////////////////////////////////////////////////////
+ // Add additional matrix entries $ -D^{-1}M I_{mm} $
+ // ///////////////////////////////////////////////////////////////
+
+ typedef typename OperatorType::row_type RowType;
+ typedef typename RowType::ConstIterator ConstIterator;
+
+ std::vector<std::vector<int> > localToGlobal(nCouplings);
+
+ for (size_t i=0; i<nCouplings; i++) {
+
+ if (fineHasObstacle[i]->size() != submat[gridIdx[i][0]]->N())
+ DUNE_THROW(Dune::Exception,
+ "fineHasObstacle[" << i << "] doesn't have the proper length!");
+
+ localToGlobal[i].resize(mortarTransferOperator[i]->N());
+ size_t idx = 0;
+ for (size_t j=0; j<fineHasObstacle[i]->size(); j++)
+ if ((*fineHasObstacle[i])[j][0])
+ localToGlobal[i][idx++] = j;
+
+ assert(idx==localToGlobal[i].size());
+
+ for (size_t j=0; j<mortarTransferOperator[i]->N(); j++) {
+
+ ConstIterator cTIt = (*mortarTransferOperator[i])[j].begin();
+ ConstIterator cEndIt = (*mortarTransferOperator[i])[j].end();
+ for (; cTIt != cEndIt; ++cTIt) {
+
+ int k = cTIt.index();
+
+ ConstIterator cMIt = (*submat[gridIdx[i][1]])[k].begin();
+ ConstIterator cMEndIt = (*submat[gridIdx[i][1]])[k].end();
+ for (; cMIt != cMEndIt; ++cMIt)
+ totalIndexSet.add(view.row(gridIdx[i][0], localToGlobal[i][j]),
+ view.col(gridIdx[i][1], cMIt.index()));
+
+ }
+
+ }
+
+ }
+
+ totalIndexSet.exportIdx(this->matrix_);
+ this->matrix_ = 0;
+
+ // Copy matrices
+ for (size_t i=0; i<nGrids; i++) {
+
+ for(size_t rowIdx=0; rowIdx<submat[i]->N(); rowIdx++) {
+
+ const RowType& row = (*submat[i])[rowIdx];
+
+ ConstIterator cIt = row.begin();
+ ConstIterator cEndIt = row.end();
+
+ for(; cIt!=cEndIt; ++cIt)
+ this->matrix_[view.row(i, rowIdx)][view.col(i, cIt.index())] = *cIt;
+
+ }
+
+ }
+
+ // ///////////////////////////////////////////////////////////////
+ // Add additional matrix entries $ -D^{-1}M I_{mm} $
+ // ///////////////////////////////////////////////////////////////
+
+ for (size_t i=0; i<nCouplings; i++) {
+
+ for (size_t j=0; j<mortarTransferOperator[i]->N(); j++) {
+
+ ConstIterator cTIt = (*mortarTransferOperator[i])[j].begin();
+ ConstIterator cTEndIt = (*mortarTransferOperator[i])[j].end();
+
+ for (; cTIt != cTEndIt; ++cTIt) {
+
+ int k = cTIt.index();
+
+ ConstIterator cMIt = (*submat[gridIdx[i][1]])[k].begin();
+ ConstIterator cMEndIt = (*submat[gridIdx[i][1]])[k].end();
+
+ for (; cMIt != cMEndIt; ++cMIt) {
+
+ auto& currentMatrixBlock = this->matrix_[view.row(gridIdx[i][0], localToGlobal[i][j])][view.col(gridIdx[i][1], cMIt.index())];
+
+ // the entry in the prolongation matrix of the mortar grid
+ const auto& subMatBlock = this->matrix_[view.row(gridIdx[i][1], k)][view.col(gridIdx[i][1], cMIt.index())];
+
+ // the - is due to the negation formula for BT
+ Dune::MatrixVector::subtractProduct(currentMatrixBlock,(*cTIt), subMatBlock);
+ }
+ }
+ }
+ }
+
+ // ///////////////////////////////////////////////////////////////
+ // Transform matrix to account for the local coordinate systems
+ // ///////////////////////////////////////////////////////////////
+
+ /** \todo Hack. Those should really be equal */
+ assert(fineLocalCoordSystems.size() <= this->matrix_.N());
+
+ for(size_t rowIdx=0; rowIdx<fineLocalCoordSystems.size(); rowIdx++) {
+
+ auto& row = this->matrix_[rowIdx];
+
+ for(auto& col : row)
+ col.leftmultiply(fineLocalCoordSystems[rowIdx]);
+
+ }
+
+}
diff --git a/src/spatial-solving/preconditioners/nbodycontacttransfer.hh b/src/spatial-solving/preconditioners/nbodycontacttransfer.hh
new file mode 100644
index 00000000..7cc8ae01
--- /dev/null
+++ b/src/spatial-solving/preconditioners/nbodycontacttransfer.hh
@@ -0,0 +1,88 @@
+#ifndef N_BODY_CONTACT_TRANSFER_HH
+#define N_BODY_CONTACT_TRANSFER_HH
+
+#include <vector>
+#include <dune/common/bitsetvector.hh>
+#include <dune/solvers/transferoperators/truncateddensemgtransfer.hh>
+#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
+
+
+/** \brief Multigrid restriction and prolongation operator for contact problems using the TNNMG method.
+ *
+ * Currently only works for first-order Lagrangian elements!
+ *
+ * This class sets up the multigrid transfer operator for n-body contact
+ * problems. It first constructs the standard prolongation matrices for the
+ * grids and then combines them into a single matrix. Then, this matrix
+ * is modified using the mortar transformation basis of Wohlmuth, Krause: "Monotone
+ * Methods on Nonmatching Grids for Nonlinear Contact Problems". In the TNNMG method
+ * the nonlinear obstacle problem is only solved on the finest level using a smoother.
+ * As coarse grid correction a linear multigrid step is applied which is later projected
+ * back onto the admissible set. Consequently the mortar basis is only needed on the finest
+ * level. This transfer operator thus additionally to the restriction and truncation,
+ * transformes the mortar basis back to the nodal basis on the finest level and otherwise
+ * just does standard restriction/prolongation. On the finest level the prolongation matrix
+ * is of the form
+
+ * \f[ \begin{pmatrix}
+ * (I_l^{l+1})_{ii} & (I_l^{l+1})_{im} & (I_l^{l+1})_{is} \\
+ * 0 & (I_l^{l+1})_{mm} & 0 \\
+ * 0 & -M_{l+1}(I_l^{l+1})_{mm} & (I_l^{l+1})_{ss}
+ * \end{pmatrix} \f].
+ *
+ * where \f[ I_l^{l+1} \f] denote the standard transfer operators
+ * and \f[ M_{l+1} \f] denotes the weighted mortar matrix
+ */
+template<class ContactNetwork, class VectorType>
+class NBodyContactTransfer : public TruncatedDenseMGTransfer<VectorType> {
+
+ enum {blocksize = VectorType::block_type::dimension};
+
+ using field_type = typename VectorType::field_type;
+
+ using MatrixBlock = Dune::FieldMatrix<field_type, blocksize, blocksize>;
+ using MatrixType = Dune::BCRSMatrix<MatrixBlock>;
+
+ using CoordSystemVector = std::vector<MatrixBlock>;
+
+public:
+
+ typedef typename MultigridTransfer<VectorType>::OperatorType OperatorType;
+
+ /** \brief Setup transfer operator for a N-body contact problem.
+ *
+ * \param grids The involved grids
+ * \param colevel The co-level of the fine grid level to which the transfer operator belongs
+ * \param mortarTransferOperator The weighted mortar matrices of the couplings
+ * \param fineLocalCoordSystems A vector containing the local coordinate systems (householder transformations) of all grids
+ * \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);
+
+protected:
+ /** \brief Combine all the standard prolongation and mortar operators to one big matrix.
+ *
+ * \param submat The standard transfer operators
+ * \param mortarTransferOperator The weighted mortar matrices of the couplings
+ * \param fineLocalCoordSystems A vector containing the local coordinate systems (householder transformations) of all grids
+ * \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.
+
+ */
+ 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<std::array<int,2> >& gridIdx);
+
+};
+
+#include "nbodycontacttransfer.cc"
+
+#endif
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh b/src/spatial-solving/preconditioners/supportpatchfactory.hh
index 07091fee..15348097 100644
--- a/src/spatial-solving/preconditioners/supportpatchfactory.hh
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh
@@ -102,12 +102,12 @@ template <class LevelContactNetwork>
class SupportPatchFactory
{
public:
- using Patch = Dune::BitSetVector<1>;
-
-private:
static const int dim = LevelContactNetwork::dim; //TODO
using ctype = typename LevelContactNetwork::ctype;
+ using Patch = Dune::BitSetVector<dim>;
+
+private:
using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;
struct RemoteIntersection {
@@ -555,14 +555,14 @@ class SupportPatchFactory
void addToPatchBoundary(const std::set<size_t>& dofs, std::set<size_t>& patchBoundary, Patch& patchDofs) {
for (auto& dof : dofs) {
patchBoundary.insert(dof);
- patchDofs[dof][0] = true;
+ patchDofs[dof] = true;
}
}
void addToPatch(const std::set<size_t>& dofs, const std::set<size_t>& patchBoundary, Patch& patchDofs) {
for (auto& dof : dofs) {
if (!patchBoundary.count(dof)) {
- patchDofs[dof][0] = false;
+ patchDofs[dof] = false;
}
}
}
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index b40d02d1..6196d3ea 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -11,22 +11,20 @@
#include "../utils/debugutils.hh"
template <class Functional, class BitVector>
+template <class Preconditioner>
SolverFactory<Functional, BitVector>::SolverFactory(
const Dune::ParameterTree& parset,
Functional& J,
+ Preconditioner&& preconditioner,
const BitVector& ignoreNodes) :
- J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J)))
- //linear solver
- //preconditioner_(),
- //linearSolverStep_(),
- {
+ J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))) {
+
nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
- // linearSolverStep_.setPreconditioner(preconditioner_);
- linearSolverStep_ = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::ldlt());
+ 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);
- linearSolver_ = std::make_shared<Dune::Solvers::UMFPackSolver<Matrix, Vector>>();
+ 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());
tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 607ac365..458b54d7 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -8,9 +8,7 @@
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-//#include <dune/solvers/iterationsteps/cgstep.hh>
-//#include <dune/solvers/iterationsteps/amgstep.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
@@ -22,8 +20,6 @@
#include "tnnmg/linesearchsolver.hh"
#include "tnnmg/localbisectionsolver.hh"
-//#include "preconditioners/supportpatchfactory.hh"
-
template <class FunctionalTEMPLATE, class BitVectorType>
class SolverFactory {
public:
@@ -34,25 +30,18 @@ class SolverFactory {
using LocalSolver = LocalBisectionSolver;
using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalSolver>;
-
using Linearization = Linearization<Functional, BitVector>;
- //using Preconditioner = Dune::Solvers::Preconditioner; //TODO Platzhalter für PatchPreconditioner
- //using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
- //using LinearSolverStep = typename Dune::Solvers::AMGStep<Matrix, Vector>;
- using LinearSolverStep = LinearIterationStep<Matrix, Vector>;
- //using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
- using LinearSolver = typename Dune::Solvers::LinearSolver<Matrix, Vector>;
-
- //using LinearSolver = typename Dune::Solvers::UMFPackSolver<Matrix, Vector>;
- //using LinearSolver = typename Dune::TNNMG::FastAMGSolver<Matrix, Vector>;
-
using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
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>
SolverFactory(Dune::ParameterTree const &,
Functional&,
+ Preconditioner&& preconditioner,
const BitVector&);
void setProblem(Vector& x);
@@ -67,7 +56,6 @@ class SolverFactory {
std::shared_ptr<NonlinearSmoother> nonlinearSmoother_;
// linear solver
- //Preconditioner preconditioner_;
std::shared_ptr<LinearSolverStep> linearSolverStep_;
std::shared_ptr<EnergyNorm<Matrix, Vector>> energyNorm_;
std::shared_ptr<LinearSolver> linearSolver_;
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index e2608f9b..fae06ed6 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -50,7 +50,8 @@ bool AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::reached
}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::advance() {
+template <class Preconditioner>
+IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::advance(Preconditioner&& preconditioner) {
/*
| C | We check here if making the step R1 of size tau is a
| R1 | R2 | good idea. To check if we can coarsen, we compare
@@ -63,7 +64,7 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
std::cout << "AdaptiveTimeStepper::advance()" << std::endl;
if (R1_.updaters == Updaters())
- R1_ = step(current_, relativeTime_, relativeTau_);
+ R1_ = step(current_, relativeTime_, relativeTau_, std::forward<Preconditioner>(preconditioner));
std::cout << "AdaptiveTimeStepper Step 1" << std::endl;
@@ -73,9 +74,9 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
UpdatersWithCount R2;
UpdatersWithCount C;
while (relativeTime_ + relativeTau_ <= 1.0) {
- R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
+ R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_, std::forward<Preconditioner>(preconditioner));
std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
- C = step(current_, relativeTime_, 2 * relativeTau_);
+ C = step(current_, relativeTime_, 2 * relativeTau_, std::forward<Preconditioner>(preconditioner));
std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
if (mustRefine_(R2.updaters, C.updaters))
break;
@@ -90,10 +91,10 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
UpdatersWithCount F2;
if (!didCoarsen) {
while (true) {
- F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
+ F1 = step(current_, relativeTime_, relativeTau_ / 2.0, std::forward<Preconditioner>(preconditioner));
std::cout << "AdaptiveTimeStepper F1 computed!" << std::endl << std::endl;
F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
- relativeTau_ / 2.0);
+ relativeTau_ / 2.0, std::forward<Preconditioner>(preconditioner));
std::cout << "AdaptiveTimeStepper F2 computed!" << std::endl << std::endl;
if (!mustRefine_(F2.updaters, R1_.updaters)) {
std::cout << "Sufficiently refined!" << std::endl;
@@ -122,15 +123,16 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Preconditioner>
typename AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::UpdatersWithCount
AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(
- Updaters const &oldUpdaters, double rTime, double rTau) {
+ Updaters const &oldUpdaters, double rTime, double rTau, Preconditioner&& preconditioner) {
UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
newUpdatersAndCount.count =
MyCoupledTimeStepper(finalTime_, parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_,
newUpdatersAndCount.updaters, errorNorms_,
externalForces_)
- .step(rTime, rTau);
+ .step(rTime, rTau, std::forward<Preconditioner>(preconditioner));
iterationRegister_.registerCount(newUpdatersAndCount.count);
return newUpdatersAndCount;
}
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index af9fb0b1..3e48b221 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -48,14 +48,17 @@ class AdaptiveTimeStepper {
std::function<bool(Updaters &, Updaters &)> mustRefine);
bool reachedEnd();
- IterationRegister advance();
+
+ template <class Preconditioner>
+ IterationRegister advance(Preconditioner&& preconditioner);
double relativeTime_;
double relativeTau_;
private:
+ template <class Preconditioner>
UpdatersWithCount step(Updaters const &oldUpdaters, double rTime,
- double rTau);
+ double rTau, Preconditioner&& preconditioner);
double finalTime_;
Dune::ParameterTree const &parset_;
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index d3b9b1e2..6b9daac7 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -25,9 +25,10 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::CoupledTimeSt
errorNorms_(errorNorms) {}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Preconditioner>
FixedPointIterationCounter
CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(double relativeTime,
- double relativeTau) {
+ double relativeTau, Preconditioner&& preconditioner) {
std::cout << "CoupledTimeStepper::step()" << std::endl;
@@ -50,8 +51,8 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(double r
FixedPointIterator fixedPointIterator(
parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
- auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
- velocityRHS, velocityIterate);
+ auto const iterations = fixedPointIterator.run(updaters_, std::forward<Preconditioner>(preconditioner),
+ velocityMatrix, velocityRHS, velocityIterate);
return iterations;
}
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 8b017775..0b2ac468 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -31,7 +31,8 @@ class CoupledTimeStepper {
const ErrorNorms& errorNorms,
ExternalForces& externalForces);
- FixedPointIterationCounter step(double relativeTime, double relativeTau);
+ template <class Preconditioner>
+ FixedPointIterationCounter step(double relativeTime, double relativeTau, Preconditioner&& preconditioner);
private:
double finalTime_;
--
GitLab