diff --git a/dune/tectonic/globalfriction.hh b/dune/tectonic/globalfriction.hh
index 41ddd18f892ea5b1130a63d9620ab5502d4f7421..e87fc139e3354111c8a6bfbb8e56284ff5da99f0 100644
--- a/dune/tectonic/globalfriction.hh
+++ b/dune/tectonic/globalfriction.hh
@@ -9,9 +9,12 @@
#include <dune/solvers/common/interval.hh>
+#include <dune/tnnmg/functionals/nonsmoothconvexfunctional.hh>
+
#include <dune/tectonic/localfriction.hh>
-template <class Matrix, class Vector> class GlobalFriction {
+template <class Matrix, class Vector>
+class GlobalFriction { //: public Dune::TNNMG::NonsmoothConvexFunctional<> {
protected:
using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
@@ -82,5 +85,6 @@ template <class Matrix, class Vector> class GlobalFriction {
}
void virtual updateAlpha(ScalarVector const &alpha) = 0;
+
};
#endif
diff --git a/dune/tectonic/mydirectionalconvexfunction.hh b/dune/tectonic/mydirectionalconvexfunction.hh
deleted file mode 100644
index 11f20c52df933d36842b7d3bf6f2a6db4534db29..0000000000000000000000000000000000000000
--- a/dune/tectonic/mydirectionalconvexfunction.hh
+++ /dev/null
@@ -1,68 +0,0 @@
-#ifndef DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
-#define DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
-
-// Copied from dune/tnnmg/problem-classes/directionalconvexfunction.hh
-// Allows phi to be const
-
-#include <dune/fufem/arithmetic.hh>
-#include <dune/solvers/common/interval.hh>
-
-template <class Nonlinearity> class MyDirectionalConvexFunction {
-public:
- using Vector = typename Nonlinearity::VectorType;
- using Matrix = typename Nonlinearity::MatrixType;
-
- MyDirectionalConvexFunction(double A, double b, Nonlinearity const &phi,
- Vector const &u, Vector const &v)
- : A(A), b(b), phi(phi), u(u), v(v) {
- phi.directionalDomain(u, v, dom);
- }
-
- double quadraticPart() const { return A; }
-
- double linearPart() const { return b; }
-
- void subDiff(double x, Dune::Solvers::Interval<double> &D) const {
- Vector uxv = u;
- Arithmetic::addProduct(uxv, x, v);
- phi.directionalSubDiff(uxv, v, D);
- auto const Axmb = A * x - b;
- D[0] += Axmb;
- D[1] += Axmb;
- }
-
- void domain(Dune::Solvers::Interval<double> &domain) const {
- domain[0] = this->dom[0];
- domain[1] = this->dom[1];
- }
-
- double A;
- double b;
-
-private:
- Nonlinearity const φ
- Vector const &u;
- Vector const &v;
-
- Dune::Solvers::Interval<double> dom;
-};
-
-/*
- 1/2 <A(u + hv),u + hv> - <b, u + hv>
- = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
-
- localA = <Av,v>
- localb = <b - Au, v>
-*/
-
-template <class Matrix, class Vector, class Nonlinearity>
-MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
- Vector const &b,
- Vector const &u,
- Vector const &v,
- Nonlinearity const &phi) {
- return MyDirectionalConvexFunction<Nonlinearity>(
- Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
-}
-
-#endif
diff --git a/dune/tectonic/tectonic.hh b/dune/tectonic/tectonic.hh
deleted file mode 100644
index bd921f8f7bdcf36d5cdc187e95f5660edeb4dcd0..0000000000000000000000000000000000000000
--- a/dune/tectonic/tectonic.hh
+++ /dev/null
@@ -1,5 +0,0 @@
-#ifndef DUNE_TECTONIC_TECTONIC_HH
-#define DUNE_TECTONIC_TECTONIC_HH
-
-// add your classes here
-#endif
diff --git a/program_structure.txt b/program_structure.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a2db29cdfb2d6791876342fd694022dc111f93fb
--- /dev/null
+++ b/program_structure.txt
@@ -0,0 +1,29 @@
+--------------------------
+-- Program structure: --
+--------------------------
+
+
+1. build n-body system
+ - contains grids, couplings, gridViews, assembler
+
+ Data structure: LevelContactNetwork
+ Factories: StackedBlocksFactory
+
+2. initialize/set up program state
+ - holds bodyStates, u, v, a, alpha for each body
+ - defines time, timeStep
+ - computes initial conditions
+
+ Data structure: ProgramState
+
+-- tested until here
+
+3. assemble RSD friction
+
+4. set up TNNMG solver
+ - rate updater
+ - state updater
+
+5. adaptive time stepper
+
+
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d86dde2a5b076bb28e0578c31c69df606c69ca8f..8a2187759e4f0a749c1225193916e479fbf7ff29 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -11,9 +11,9 @@ set(SW_SOURCE_FILES
io/hdf5/restart-io.cc
io/hdf5/surface-writer.cc
io/hdf5/time-writer.cc
- one-body-problem-data/mygeometry.cc
- one-body-problem-data/mygrid.cc
- one-body-problem.cc
+# one-body-problem-data/mygeometry.cc
+# one-body-problem-data/mygrid.cc
+# one-body-problem.cc
spatial-solving/fixedpointiterator.cc
spatial-solving/solverfactory.cc
time-stepping/adaptivetimestepper.cc
@@ -46,6 +46,7 @@ set(MSW_SOURCE_FILES
multi-body-problem.cc
spatial-solving/fixedpointiterator.cc
spatial-solving/solverfactory.cc
+ spatial-solving/solverfactory_old.cc
time-stepping/adaptivetimestepper.cc
time-stepping/coupledtimestepper.cc
time-stepping/rate.cc
diff --git a/src/data-structures/.cc b/src/data-structures/.cc
new file mode 100644
index 0000000000000000000000000000000000000000..32c9a280f0a4babfa053e3bf37aebf154527b705
--- /dev/null
+++ b/src/data-structures/.cc
@@ -0,0 +1,7 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+
+template class ContactNetwork<Grid, MY_DIM>;
diff --git a/src/data-structures/levelcontactnetwork.hh b/src/data-structures/levelcontactnetwork.hh
index f6c1bf194f3889ff2effa0a446f317df4e09621f..e7672cd53828ed7c2045a6ee0d1cc1c176a10981 100644
--- a/src/data-structures/levelcontactnetwork.hh
+++ b/src/data-structures/levelcontactnetwork.hh
@@ -56,6 +56,8 @@ template <class GridType, int dimension> class LevelContactNetwork {
using StateEnergyNorms = std::vector<const typename Body::StateEnergyNorm *>;
using ExternalForces = std::vector<const typename Body::ExternalForce *>;
+ using NBodyAssembler = Dune::Contact::NBodyAssembler<DeformedGridType, Vector>;
+
public:
LevelContactNetwork(int nBodies, int nCouplings, int level = 0);
@@ -133,7 +135,7 @@ template <class GridType, int dimension> class LevelContactNetwork {
std::vector<std::vector<std::unique_ptr<const DeformedLevelGridView>>> levelViews_;
std::vector<size_t> leafVertexCounts_;
- Dune::Contact::NBodyAssembler<DeformedGridType, Vector> nBodyAssembler_;
+ NBodyAssembler nBodyAssembler_;
Matrices<Matrix,2> matrices_;
};
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 275697da77effe110850f77c239a4d1981b9a41c..181e5c86c568f83f8d81239fd809fabea2250a2c 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -6,7 +6,6 @@
#include <dune/matrix-vector/axpy.hh>
#include <dune/fufem/boundarypatch.hh>
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
#include <dune/contact/assemblers/nbodyassembler.hh>
@@ -17,7 +16,9 @@
#include "levelcontactnetwork.hh"
#include "matrices.hh"
#include "../spatial-solving/solverfactory.hh"
-
+#include "../spatial-solving/solverfactory_old.hh"
+#include "../spatial-solving/tnnmg/functional.hh"
+#include "../spatial-solving/tnnmg/zerononlinearity.hh"
#include "../utils/debugutils.hh"
template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class BodyState {
@@ -118,15 +119,43 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
- Vector totalRhs;
+ Vector totalRhs, oldTotalRhs;
nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
+ oldTotalRhs = totalRhs;
- Vector totalX;
+ Vector totalX, oldTotalX;
nBodyAssembler.nodalToTransformed(_x, totalX);
+ oldTotalX = totalX;
+
+ // get lower and upper obstacles
+ const auto& totalObstacles = nBodyAssembler.totalObstacles_;
+ Vector lower(totalObstacles.size());
+ Vector upper(totalObstacles.size());
+
+ for (size_t j=0; j<totalObstacles.size(); ++j) {
+ const auto& totalObstaclesj = totalObstacles[j];
+ auto& lowerj = lower[j];
+ auto& upperj = upper[j];
+ for (size_t d=0; d<dims; ++d) {
+ lowerj[d] = totalObstaclesj[d][0];
+ upperj[d] = totalObstaclesj[d][1];
+ }
+ }
+
+ // print problem
+ print(bilinearForm, "bilinearForm");
+ print(totalRhs, "totalRhs");
+ print(_dirichletNodes, "ignore");
+ print(totalObstacles, "totalObstacles");
+ print(lower, "lower");
+ print(upper, "upper");
- using DeformedGridType = typename LevelContactNetwork<GridType, dims>::DeformedGridType;
- using LinearFactory = SolverFactory<DeformedGridType, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
- LinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, _dirichletNodes);
+ // 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);
+
+ using Factory = SolverFactory<Functional, Dune::BitSetVector<dims>>;
+ Factory factory(parset.sub("solver.tnnmg"), J);
/* std::vector<Dune::BitSetVector<dims>> bodyDirichletNodes;
nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
@@ -138,13 +167,14 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
print(totalX, "totalX: ");
print(totalRhs, "totalRhs: ");*/
- auto multigridStep = factory.getStep();
- multigridStep->setProblem(bilinearForm, totalX, totalRhs);
+ auto tnnmgStep = factory.step();
+ tnnmgStep->setIgnore(_dirichletNodes);
+ factory.setProblem(totalX);
const EnergyNorm<Matrix, Vector> norm(bilinearForm);
LoopSolver<Vector> solver(
- multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
+ tnnmgStep.get(), _localParset.get<size_t>("maximumIterations"),
_localParset.get<double>("tolerance"), &norm,
_localParset.get<Solver::VerbosityMode>("verbosity"),
false); // absolute error
@@ -152,7 +182,38 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
solver.preprocess();
solver.solve();
- nBodyAssembler.postprocess(multigridStep->getSol(), _x);
+ nBodyAssembler.postprocess(tnnmgStep->getSol(), _x);
+
+ Vector res = totalRhs;
+ bilinearForm.mmv(tnnmgStep->getSol(), res);
+ std::cout << "TNNMG Res - energy norm: " << norm.operator ()(res) << std::endl;
+
+ // TODO: remove after debugging
+ /* using DeformedGridType = typename LevelContactNetwork<GridType, dims>::DeformedGridType;
+ using OldLinearFactory = SolverFactoryOld<DeformedGridType, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+ OldLinearFactory oldFactory(parset.sub("solver.tnnmg"), nBodyAssembler, _dirichletNodes);
+
+ auto oldStep = oldFactory.getStep();
+ oldStep->setProblem(bilinearForm, oldTotalX, oldTotalRhs);
+
+ LoopSolver<Vector> oldSolver(
+ oldStep.get(), _localParset.get<size_t>("maximumIterations"),
+ _localParset.get<double>("tolerance"), &norm,
+ _localParset.get<Solver::VerbosityMode>("verbosity"),
+ false); // absolute error
+
+ oldSolver.preprocess();
+ oldSolver.solve();
+
+ Vector oldRes = totalRhs;
+ bilinearForm.mmv(oldStep->getSol(), oldRes);
+ std::cout << "Old Res - energy norm: " << norm.operator ()(oldRes) << std::endl;*/
+
+ print(tnnmgStep->getSol(), "TNNMG Solution: ");
+ /* print(oldStep->getSol(), "Old Solution: ");
+ auto diff = tnnmgStep->getSol();
+ diff -= oldStep->getSol();
+ std::cout << "Energy norm: " << norm.operator ()(diff) << std::endl;*/
};
timeStep = 0;
@@ -162,6 +223,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
std::vector<Vector> ell0(bodyCount_);
for (size_t i=0; i<bodyCount_; i++) {
// Initial velocity
+ u[i] = 0.0;
v[i] = 0.0;
ell0[i].resize(u[i].size());
@@ -175,7 +237,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
Dune::BitSetVector<dims> dirichletNodes;
levelContactNetwork.totalNodes("dirichlet", dirichletNodes);
- for (size_t i=0; i<dirichletNodes.size(); i++) {
+ /*for (size_t i=0; i<dirichletNodes.size(); i++) {
bool val = false;
for (size_t d=0; d<dims; d++) {
val = val || dirichletNodes[i][d];
@@ -185,7 +247,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
for (size_t d=0; d<dims; d++) {
dirichletNodes[i][d] = val;
}
- }
+ }*/
solveLinearProblem(dirichletNodes, levelContactNetwork.matrices().elasticity, ell0, u,
parset.sub("u0.solver"));
diff --git a/src/explicitgrid.hh b/src/explicitgrid.hh
index beefb64cc800f06604d3ac12484eb0be305e7902..fb010a25306c621807104a97869c02312c610b34 100644
--- a/src/explicitgrid.hh
+++ b/src/explicitgrid.hh
@@ -12,4 +12,7 @@ using LevelGridView = Grid::LevelGridView;
using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
using DeformedGrid = DeformedGridComplex::DeformedGridType;
+using DefLeafGridView = DeformedGrid::LeafGridView;
+using DefLevelGridView = DeformedGrid::LevelGridView;
+
#endif
diff --git a/src/explicitvectors.hh b/src/explicitvectors.hh
index 5f2d887fe6fd41720409d537b9960d122d9fb468..548c681faa9e6e867641e591170ee219101fcf0b 100644
--- a/src/explicitvectors.hh
+++ b/src/explicitvectors.hh
@@ -3,6 +3,7 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
@@ -12,4 +13,6 @@ using Vector = Dune::BlockVector<LocalVector>;
using Matrix = Dune::BCRSMatrix<LocalMatrix>;
using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+using BitVector = Dune::BitSetVector<MY_DIM>;
+
#endif
diff --git a/src/factories/stackedblocksfactory.cc b/src/factories/stackedblocksfactory.cc
index 0891faeed91b05c974393734c48b20442ca01a2c..11a3f98c65428e55638aa5eca9f3f0da26c38f4b 100644
--- a/src/factories/stackedblocksfactory.cc
+++ b/src/factories/stackedblocksfactory.cc
@@ -163,27 +163,40 @@ void StackedBlocksFactory<GridType, dims>::setBoundaryConditions() {
std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->leafView()), neumannFunction, "neumann");
body->addBoundaryCondition(neumannBoundary);
- // Dirichlet Boundary
+ // right Dirichlet Boundary
// identify Dirichlet nodes on leaf level
- std::shared_ptr<Dune::BitSetVector<dims>> dirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(leafVertexCount);
+ std::shared_ptr<Dune::BitSetVector<dims>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(leafVertexCount);
for (int j=0; j<leafVertexCount; j++) {
if (leafFaces_[i]->right.containsVertex(j))
- (*dirichletNodes)[j][0] = true;
+ (*velocityDirichletNodes)[j][0] = true;
- if (leafFaces_[i]->lower.containsVertex(j))
- (*dirichletNodes)[j][0] = true;
-
- #if MY_DIM == 3
+ #if MY_DIM == 3 //TODO: wrong, needs revision
if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
- dirichletNodes->at(j)[2] = true;
+ zeroDirichletNodes->at(j)[2] = true;
#endif
}
- std::shared_ptr<LeafBoundaryCondition> dirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
- dirichletBoundary->setBoundaryPatch(body->leafView(), dirichletNodes);
- dirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
- body->addBoundaryCondition(dirichletBoundary);
+ std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ velocityDirichletBoundary->setBoundaryPatch(body->leafView(), velocityDirichletNodes);
+ velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+ body->addBoundaryCondition(velocityDirichletBoundary);
+ }
+
+ // 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);
+ for (int j=0; j<firstLeafVertexCount; j++) {
+ if (leafFaces_[0]->lower.containsVertex(j))
+ (*zeroDirichletNodes)[j][1] = true;
}
+
+ std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ zeroDirichletBoundary->setBoundaryPatch(firstBody->leafView(), zeroDirichletNodes);
+
+ std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
+ zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
+ firstBody->addBoundaryCondition(zeroDirichletBoundary);
}
#include "stackedblocksfactory_tmpl.cc"
diff --git a/src/io/hdf5-bodywriter.hh b/src/io/hdf5-bodywriter.hh
new file mode 100644
index 0000000000000000000000000000000000000000..549062c706b1c617b779f0f6c884ace4aafb6bf7
--- /dev/null
+++ b/src/io/hdf5-bodywriter.hh
@@ -0,0 +1,79 @@
+#ifndef SRC_HDF5_BODYWRITER_HH
+#define SRC_HDF5_BODYWRITER_HH
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include "hdf5/frictionalboundary-writer.hh"
+#include "hdf5/patchinfo-writer.hh"
+
+template <class ProgramState, class VertexBasis, class GridView>
+class HDF5BodyWriter {
+private:
+ using VertexCoordinates = typename ProgramState::Vector;
+ using FrictionPatches = std::vector<const BoundaryPatch<GridView>* >;
+ using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
+
+ using LocalVector = typename VertexCoordinates::block_type;
+
+ friend class HDF5LevelWriter<ProgramState, VertexBasis, GridView>;
+
+public:
+ HDF5BodyWriter(
+ HDF5::Grouplike &file,
+ const VertexCoordinates& vertexCoordinates,
+ const VertexBasis& vertexBasis,
+ const FrictionPatches& frictionPatches,
+ const WeakPatches& weakPatches) :
+ patchCount_(frictionPatches.size()),
+ file_(file),
+#if MY_DIM == 3
+ patchInfoWriters_(patchCount_),
+#endif
+ frictionBoundaryWriters_(patchCount_) {
+
+ assert(patchCount_ == weakPatches.size());
+
+ for (size_t i=0; i<patchCount_; i++) {
+#if MY_DIM == 3
+ patchInfoWriters_[i] = new PatchInfoWriter<ProgramState, VertexBasis, GridView>(file_, vertexBasis, *frictionPatches[i], *weakPatches[i], i);
+#endif
+ frictionBoundaryWriters_[i] = new FrictionalBoundaryWriter<ProgramState, GridView>(file_, vertexCoordinates, *frictionPatches[i], i);
+ }
+ }
+
+ ~HDF5BodyWriter() {
+ for (size_t i=0; i<patchCount_; i++) {
+#if MY_DIM == 3
+ delete patchInfoWriters_[i];
+#endif
+ delete frictionBoundaryWriters_[i];
+ }
+ }
+
+ template <class Friction>
+ void reportSolution(ProgramState const &programState,
+ // for the friction coefficient
+ std::vector<std::shared_ptr<Friction>>& friction) {
+
+ for (size_t i=0; i<patchCount_; i++) {
+#if MY_DIM == 3
+ patchInfoWriters_[i]->write(programState);
+#endif
+ frictionBoundaryWriters_[i]->write(programState, *friction[i]);
+ }
+ }
+
+private:
+ const size_t patchCount_;
+
+ HDF5::Grouplike &file_;
+
+#if MY_DIM == 3
+ std::vector<PatchInfoWriter<ProgramState, VertexBasis, GridView>* > patchInfoWriters_;
+#endif
+
+ std::vector<FrictionalBoundaryWriter<ProgramState, GridView>* > frictionBoundaryWriters_;
+};
+#endif
diff --git a/src/io/hdf5-levelwriter.hh b/src/io/hdf5-levelwriter.hh
new file mode 100644
index 0000000000000000000000000000000000000000..66e3c048b7787a742fd3cd4797346b8cbbf3ea16
--- /dev/null
+++ b/src/io/hdf5-levelwriter.hh
@@ -0,0 +1,73 @@
+#ifndef SRC_HDF5_LEVELWRITER_HH
+#define SRC_HDF5_LEVELWRITER_HH
+
+#include <dune/fufem/hdf5/file.hh>
+
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/time-writer.hh"
+#include "hdf5-bodywriter.hh"
+
+template <class ProgramState, class VertexBasis, class GridView>
+class HDF5LevelWriter {
+private:
+ using HDF5BodyWriter = HDF5BodyWriter<ProgramState, VertexBasis, GridView>;
+ using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
+ using VertexBases = std::vector<VertexBasis>;
+ using FrictionPatches = std::vector<typename HDF5BodyWriter::FrictionPatches>;
+ using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
+
+ friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
+
+public:
+
+ HDF5LevelWriter(HDF5::Grouplike &file,
+ const VertexCoordinates& vertexCoordinates,
+ const VertexBases& vertexBases,
+ const FrictionPatches& frictionBoundaries,
+ const WeakPatches& weakPatches)
+ : bodyCount_(vertexCoordinates.size()),
+ file_(file),
+ iterationWriter_(file_),
+ timeWriter_(file_),
+ bodyWriters_(bodyCount_) {
+
+ for (size_t i=0; i<bodyCount_; i++) {
+ bodyWriters_[i] = new HDF5BodyWriter<ProgramState, VertexBasis, GridView>(file_, vertexCoordinates[i], vertexBases[i], frictionBoundaries[i], weakPatches[i]);
+ }
+ }
+
+ ~HDF5LevelWriter() {
+ for (size_t i=0; i<bodyCount_; i++) {
+ delete bodyWriters_[i];
+ }
+ }
+
+ template <class Friction>
+ void reportSolution(
+ ProgramState const &programState,
+ // for the friction coefficient
+ std::vector<std::shared_ptr<Friction>>& friction) {
+ timeWriter_.write(programState);
+
+ for (size_t i=0; i<bodyCount_; i++) {
+ bodyWriters_[i]->reportSolution(programState, *friction[i]); //TODO
+ }
+ }
+
+ void reportIterations(
+ ProgramState const &programState,
+ IterationRegister const &iterationCount) {
+ iterationWriter_.write(programState.timeStep, iterationCount);
+ }
+
+private:
+ const size_t bodyCount_;
+
+ HDF5::Grouplike &file_;
+
+ IterationWriter iterationWriter_;
+ TimeWriter<ProgramState> timeWriter_;
+
+ std::vector<HDF5BodyWriter* > bodyWriters_;
+};
+#endif
diff --git a/src/io/hdf5-writer.hh b/src/io/hdf5-writer.hh
deleted file mode 100644
index d0351be71cdc9e4e61c82fa32520ef81b5ba1707..0000000000000000000000000000000000000000
--- a/src/io/hdf5-writer.hh
+++ /dev/null
@@ -1,85 +0,0 @@
-#ifndef SRC_HDF5_WRITER_HH
-#define SRC_HDF5_WRITER_HH
-
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/geometry/convexpolyhedron.hh>
-#include <dune/fufem/hdf5/file.hh>
-
-#include "hdf5/frictionalboundary-writer.hh"
-#include "hdf5/iteration-writer.hh"
-#include "hdf5/patchinfo-writer.hh"
-#include "hdf5/surface-writer.hh"
-#include "hdf5/time-writer.hh"
-
-template <class ProgramState, class VertexBasis, class GridView>
-class HDF5Writer {
-private:
- using Vector = typename ProgramState::Vector;
- using Patch = BoundaryPatch<GridView>;
-
- using LocalVector = typename Vector::block_type;
-
-public:
- HDF5Writer(HDF5::Grouplike &file,
- const std::vector<Vector>& vertexCoordinates,
- const std::vector<const VertexBasis* >& vertexBases,
- const std::vector<Patch>& frictionalBoundaries,
- const std::vector<ConvexPolyhedron<LocalVector>>& weakPatches)
- : bodyCount_(vertexCoordinates.size()),
- file_(file),
- iterationWriter_(file_),
- timeWriter_(file_),
-#if MY_DIM == 3
- patchInfoWriters_(bodyCount_),
-#endif
- frictionalBoundaryWriters_(bodyCount_) {
-
- for (size_t i=0; i<bodyCount_; i++) {
-#if MY_DIM == 3
- patchInfoWriters_[i] = new PatchInfoWriter<ProgramState, VertexBasis, GridView>(file_, *vertexBases[i], frictionalBoundaries[i], weakPatches[i], i);
-#endif
- frictionalBoundaryWriters_[i] = new FrictionalBoundaryWriter<ProgramState, GridView>(file_, vertexCoordinates[i], frictionalBoundaries[i], i);
- }
- }
-
- ~HDF5Writer() {
- for (size_t i=0; i<bodyCount_; i++) {
-#if MY_DIM == 3
- delete patchInfoWriters_[i];
-#endif
- delete frictionalBoundaryWriters_[i];
- }
- }
-
- template <class Friction>
- void reportSolution(ProgramState const &programState,
- // for the friction coefficient
- std::vector<std::shared_ptr<Friction>>& friction) {
- timeWriter_.write(programState);
-
- for (size_t i=0; i<bodyCount_; i++) {
-#if MY_DIM == 3
- patchInfoWriters_[i]->write(programState);
-#endif
- frictionalBoundaryWriters_[i]->write(programState, *friction[i]);
- }
- }
-
- void reportIterations(ProgramState const &programState,
- IterationRegister const &iterationCount) {
- iterationWriter_.write(programState.timeStep, iterationCount);
- }
-
-private:
- const size_t bodyCount_;
-
- HDF5::Grouplike &file_;
-
- IterationWriter iterationWriter_;
- TimeWriter<ProgramState> timeWriter_;
-#if MY_DIM == 3
- std::vector<PatchInfoWriter<ProgramState, VertexBasis, GridView>* > patchInfoWriters_;
-#endif
- std::vector<FrictionalBoundaryWriter<ProgramState, GridView>* > frictionalBoundaryWriters_;
-};
-#endif
diff --git a/src/io/hdf5/frictionalboundary-writer.hh b/src/io/hdf5/frictionalboundary-writer.hh
index 4372cdad928bb60643b8996d7571dbc7d14f6cc7..0977401382ef827d78499316607077f3b722bf32 100644
--- a/src/io/hdf5/frictionalboundary-writer.hh
+++ b/src/io/hdf5/frictionalboundary-writer.hh
@@ -11,7 +11,7 @@ template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
using Patch = BoundaryPatch<GridView>;
public:
- FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
+ FrictionalBoundaryWriter(HDF5::Grouplike& file, Vector const &vertexCoordinates,
Patch const &frictionalBoundary, size_t const id);
template <class Friction>
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index de02c49787df22f05ef2cc28ec4b38749533d36d..4d6e7173513a627a72fe179388b6df072e3b9d7b 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -30,6 +30,7 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/formatstring.hh>
+#include <dune/fufem/hdf5/file.hh>
#include <dune/solvers/norms/energynorm.hh>
@@ -46,12 +47,13 @@
#include <dune/contact/projections/normalprojection.hh>
#include <dune/tectonic/geocoordinate.hh>
-#include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/globalfriction.hh>
-#include <dune/fufem/hdf5/file.hh>
+
#include "assemblers.hh"
#include "gridselector.hh"
+#include "explicitgrid.hh"
+#include "explicitvectors.hh"
#include "data-structures/enumparser.hh"
#include "data-structures/enums.hh"
@@ -61,7 +63,7 @@
#include "factories/stackedblocksfactory.hh"
-#include "io/hdf5-writer.hh"
+//#include "io/hdf5-levelwriter.hh"
#include "io/hdf5/restart-io.hh"
#include "io/vtk.hh"
@@ -69,9 +71,11 @@
#include "multi-body-problem-data/mybody.hh"
#include "multi-body-problem-data/grid/mygrids.hh"
-//#include "spatial-solving/solverfactory.hh"
+#include "spatial-solving/tnnmg/functional.hh"
+#include "spatial-solving/tnnmg/localbisectionsolver.hh"
+#include "spatial-solving/solverfactory.hh"
-//#include "time-stepping/adaptivetimestepper.hh"
+#include "time-stepping/adaptivetimestepper.hh"
#include "time-stepping/rate.hh"
#include "time-stepping/state.hh"
#include "time-stepping/updaters.hh"
@@ -82,7 +86,8 @@
// for getcwd
#include <unistd.h>
-#define USE_OLD_TNNMG
+//#include <tbb/tbb.h> //TODO multi threading preconditioner?
+//#include <pthread.h>
size_t const dims = MY_DIM;
@@ -115,18 +120,7 @@ int main(int argc, char *argv[]) {
auto const parset = getParameters(argc, argv);
- using Vector = Dune::BlockVector<Dune::FieldVector<double, dims>>;
- using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
- using DeformedGrid = DeformedGridComplex::DeformedGridType;
-
- using DefLeafGridView = DeformedGrid::LeafGridView;
- using DefLevelGridView = DeformedGrid::LevelGridView;
-
using Assembler = MyAssembler<DefLeafGridView, dims>;
- using Matrix = Assembler::Matrix;
- using LocalVector = Vector::block_type;
- using ScalarMatrix = Assembler::ScalarMatrix;
- using ScalarVector = Assembler::ScalarVector;
using field_type = Matrix::field_type;
@@ -225,18 +219,20 @@ int main(int argc, char *argv[]) {
typename LevelContactNetwork::BoundaryPatches frictionBoundaries;
levelContactNetwork.boundaryPatches("friction", frictionBoundaries);
+ /*
auto dataWriter =
writeData ? std::make_unique<
HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
*dataFile, vertexCoordinates, vertexBases,
- frictionBoundaries, weakPatches)
- : nullptr;
-/*
+ frictionBoundaries) //, weakPatches)
+ : nullptr;*/
+
const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "body");
IterationRegister iterationCount;
auto const report = [&](bool initial = false) {
+ /*
if (writeData) {
dataWriter->reportSolution(programState, levelContactNetwork.globalFriction());
if (!initial)
@@ -255,7 +251,7 @@ int main(int argc, char *argv[]) {
<< ", time = " << std::setw(12) << programState.relativeTime
<< ", tau = " << std::setw(12) << programState.relativeTau
<< std::endl;
-
+ */
if (parset.get<bool>("io.vtk.write")) {
std::vector<ScalarVector> stress(bodyCount);
@@ -264,6 +260,7 @@ int main(int argc, char *argv[]) {
body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
body->data()->getPoissonRatio(),
programState.u[i], stress[i]);
+
}
vtkWriter.write(programState.timeStep, programState.u, programState.v,
@@ -271,7 +268,7 @@ int main(int argc, char *argv[]) {
}
};
report(true);
-*/
+
// -------------------
// Set up TNNMG solver
// -------------------
@@ -279,8 +276,9 @@ int main(int argc, char *argv[]) {
Dune::BitSetVector<dims> totalDirichletNodes;
levelContactNetwork.totalNodes("dirichlet", totalDirichletNodes);
- using NonlinearFactory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
- NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, totalDirichletNodes);
+
+ using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
+ using NonlinearFactory = SolverFactory<Functional, BitVector>;
using BoundaryFunctions = typename LevelContactNetwork::BoundaryFunctions;
using BoundaryNodes = typename LevelContactNetwork::BoundaryNodes;
@@ -344,11 +342,11 @@ int main(int argc, char *argv[]) {
typename LevelContactNetwork::ExternalForces externalForces;
levelContactNetwork.externalForces(externalForces);
- AdaptiveTimeStepper<NonlinearFactory, Updaters, EnergyNorm<ScalarMatrix, ScalarVector>>
- adaptiveTimeStepper(factory, parset, globalFriction, current,
+ AdaptiveTimeStepper<NonlinearFactory, decltype(nBodyAssembler), Updaters, EnergyNorm<ScalarMatrix, ScalarVector>>
+ adaptiveTimeStepper(parset, nBodyAssembler, globalFriction, current,
programState.relativeTime, programState.relativeTau,
externalForces, stateEnergyNorms, mustRefine);
-/*
+
while (!adaptiveTimeStepper.reachedEnd()) {
programState.timeStep++;
@@ -374,7 +372,7 @@ int main(int argc, char *argv[]) {
break;
}
}
-*/
+
std::cout.rdbuf(coutbuf); //reset to standard output again
} catch (Dune::Exception &e) {
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index b850b5cc829734e50342d1e812254778f74ad97b..80028801f239be2aeb666ece146963b4d9807b61 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -44,12 +44,12 @@ b = 0.005 # [ ]
scheme = newmark
[u0.solver]
-maximumIterations = 100000
-verbosity = quiet
+maximumIterations = 10000
+verbosity = full
[a0.solver]
-maximumIterations = 100000
-verbosity = quiet
+maximumIterations = 10000
+verbosity = full
[v.solver]
maximumIterations = 100000
@@ -60,7 +60,7 @@ maximumIterations = 10000
lambda = 0.5
[solver.tnnmg.linear]
-maxiumumIterations = 100000
+maximumIterations = 100
pre = 3
cycle = 1 # 1 = V, 2 = W, etc.
post = 3
diff --git a/src/one-body-problem-data/mygrid.cc b/src/one-body-problem-data/mygrid.cc
index bc7c77daf011d97b2b96bfa17ec9e6279fdaa86b..7188f73a0c0ab5b295a5aca2174c3087144ee1cc 100644
--- a/src/one-body-problem-data/mygrid.cc
+++ b/src/one-body-problem-data/mygrid.cc
@@ -6,7 +6,7 @@
#include "mygrid.hh"
#include "midpoint.hh"
-#include "../diameter.hh"
+#include "../utils/diameter.hh"
#if MY_DIM == 3
SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
diff --git a/src/one-body-problem.cfg b/src/one-body-problem.cfg
index d970f5e93a5f6dc77c3be6e3f053cf4a85e73529..7336759ce77ba8ff4d3183a3f34f6bdbf5af4eb5 100644
--- a/src/one-body-problem.cfg
+++ b/src/one-body-problem.cfg
@@ -59,7 +59,7 @@ maximumIterations = 10000
lambda = 0.5
[solver.tnnmg.linear]
-maxiumumIterations = 100000
+maximumIterations = 100000
pre = 3
cycle = 1 # 1 = V, 2 = W, etc.
post = 3
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 6b776feb982debaca743bd3bb6598d4a8ebdd8ca..1dbaec370013f5b97931d5b894164f41a1b92e66 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -27,19 +27,21 @@
#include "fixedpointiterator.hh"
+#include "../utils/debugutils.hh"
+
void FixedPointIterationCounter::operator+=(
FixedPointIterationCounter const &other) {
iterations += other.iterations;
multigridIterations += other.multigridIterations;
}
-template <class Factory, class Updaters, class ErrorNorm>
-FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
- Factory &factory, Dune::ParameterTree const &parset,
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::FixedPointIterator(
+ Dune::ParameterTree const &parset,
+ NBodyAssembler& nBodyAssembler,
GlobalFrictionContainer& globalFriction, const std::vector<const ErrorNorm* >& errorNorms)
- : factory_(factory),
- step_(factory_.getStep()),
- parset_(parset),
+ : parset_(parset),
+ nBodyAssembler_(nBodyAssembler),
globalFriction_(globalFriction),
fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
@@ -49,23 +51,49 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
errorNorms_(errorNorms) {}
-template <class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
FixedPointIterationCounter
-FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
+FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
- auto contactAssembler = factory_.getNBodyAssembler();
-
- Matrix mat;
std::vector<const Matrix*> matrices_ptr(velocityMatrices.size());
for (size_t i=0; i<matrices_ptr.size(); i++) {
matrices_ptr[i] = &velocityMatrices[i];
}
- contactAssembler.assembleJacobian(matrices_ptr, mat);
- EnergyNorm<Matrix, Vector> energyNorm(mat);
- LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
+ // assemble full global contact problem
+ Matrix bilinearForm;
+ nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
+
+ Vector totalRhs;
+ nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
+
+ Vector totalX;
+ nBodyAssembler_.nodalToTransformed(velocityIterates, totalX);
+
+ // get lower and upper obstacles
+ const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
+ Vector lower(totalObstacles.size());
+ Vector upper(totalObstacles.size());
+
+ for (size_t j=0; j<totalObstacles.size(); ++j) {
+ const auto& totalObstaclesj = totalObstacles[j];
+ auto& lowerj = lower[j];
+ auto& upperj = upper[j];
+ for (size_t d=0; d<dims; ++d) {
+ lowerj[d] = totalObstaclesj[d][0];
+ upperj[d] = totalObstaclesj[d][1];
+ }
+ }
+
+ // set up functional and TNMMG solver
+ Functional J(bilinearForm, totalRhs, ..., lower, upper);
+ Factory solverFactory(parset.sub("solver.tnnmg"), J);
+ auto step = solverFactory.step();
+
+ EnergyNorm<Matrix, Vector> energyNorm(bilinearForm);
+ LoopSolver<Vector> velocityProblemSolver(step.get(), velocityMaxIterations_,
velocityTolerance_, &energyNorm,
verbosity_, false); // absolute error
@@ -101,20 +129,24 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
// assemble full global contact problem
Matrix bilinearForm;
- contactAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+ nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
Vector totalRhs;
- contactAssembler.assembleRightHandSide(rhs, totalRhs);
+ nBodyAssembler_.assembleRightHandSide(rhs, totalRhs);
Vector totalVelocityIterate;
- contactAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+ nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
+
+ print(bilinearForm, "matrix: ");
+ print(totalRhs, "totalRhs: ");
+ print(totalVelocityIterate, "iterate: ");
// solve a velocity problem
- step_->setProblem(bilinearForm, totalVelocityIterate, totalRhs);
+ step->setProblem(bilinearForm, totalVelocityIterate, totalRhs);
velocityProblemSolver.preprocess();
velocityProblemSolver.solve();
- contactAssembler.postprocess(totalVelocityIterate, velocityIterates);
+ nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
multigridIterations += velocityProblemSolver.getResult().iterations;
@@ -167,8 +199,8 @@ std::ostream &operator<<(std::ostream &stream,
<< ")";
}
-template <class Factory, class Updaters, class ErrorNorm>
-void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const {
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const {
// init result
v_rel.resize(v.size());
for (size_t i=0; i<v_rel.size(); i++) {
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index 688a80b516c3a8dfe4e34764a9a48688477118a7..02ad3089bfb3597267d46042dfa81985d1f0c57b 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -22,15 +22,18 @@ struct FixedPointIterationCounter {
std::ostream &operator<<(std::ostream &stream,
FixedPointIterationCounter const &fpic);
-template <class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
class FixedPointIterator {
+ using Functional = typename Factory::Functional;
using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
- using Nonlinearity = typename Factory::Nonlinearity;
+ using Nonlinearity = typename Factory::Functional::Nonlinearity;
+
+ const static int dims = Vector::block_type::dimension;
// using Nonlinearity = typename ConvexProblem::NonlinearityType;
- using DeformedGrid = typename Factory::DeformedGrid;
+ // using DeformedGrid = typename Factory::DeformedGrid;
public:
using GlobalFrictionContainer = GlobalFrictionContainer<Nonlinearity, 2>;
@@ -39,7 +42,8 @@ class FixedPointIterator {
void relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const;
public:
- FixedPointIterator(Factory& factory, const Dune::ParameterTree& parset,
+ FixedPointIterator(const Dune::ParameterTree& parset,
+ NBodyAssembler& nBodyAssembler,
GlobalFrictionContainer& globalFriction,
const std::vector<const ErrorNorm* >& errorNorms);
@@ -49,9 +53,8 @@ class FixedPointIterator {
std::vector<Vector>& velocityIterates);
private:
- Factory& factory_;
- std::shared_ptr<typename Factory::Step> step_;
Dune::ParameterTree const &parset_;
+ NBodyAssembler& nBodyAssembler_;
GlobalFrictionContainer& globalFriction_;
size_t fixedPointMaxIterations_;
diff --git a/src/spatial-solving/fixedpointiterator_tmpl.cc b/src/spatial-solving/fixedpointiterator_tmpl.cc
index eec3461083148f6a5b8b372dcd7abe3cd655e61a..20322c8cc84035fe8792b21fb2998fdd4d52435f 100644
--- a/src/spatial-solving/fixedpointiterator_tmpl.cc
+++ b/src/spatial-solving/fixedpointiterator_tmpl.cc
@@ -8,25 +8,37 @@
#include <dune/common/function.hh>
#include <dune/solvers/norms/energynorm.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
+
#include <dune/contact/common/deformedcontinuacomplex.hh>
#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
#include "../data-structures/levelcontactnetwork.hh"
#include "../time-stepping/rate/rateupdater.hh"
#include "../time-stepping/state/stateupdater.hh"
#include "../time-stepping/updaters.hh"
+
#include "solverfactory.hh"
+#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+//#include "../../dune/tectonic/globalfriction.hh"
+//#include "functional.hh"
+//#include "localbisectionsolver.hh"
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
+using MyLevelContactNetwork = LevelContactNetwork<Grid, MY_DIM>;
+using BoundaryNodes = typename MyLevelContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyLevelContactNetwork::BoundaryFunctions;
-using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
+using NonlinearLocalSolver = LocalBisectionSolver;
+using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<MyFunctional, NonlinearLocalSolver>;
+using Factory = SolverFactory<MyFunctional, BitVector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
@@ -34,4 +46,4 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
-template class FixedPointIterator<Factory, MyUpdaters, ErrorNorm>;
+template class FixedPointIterator<Factory, typename MyLevelContactNetwork::NBodyAssembler, MyUpdaters, ErrorNorm>;
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 0bba4b88001d23557b8608aaef3f7f40f02fd42f..25bc75755807abb46b4201ed8f25aa319732cc21 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -2,77 +2,41 @@
#include "config.h"
#endif
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#include <dune/common/bitsetvector.hh>
-
-#include <dune/fufem/assemblers/transferoperatorassembler.hh>
-#include <dune/solvers/solvers/solver.hh>
+#include <dune/solvers/common/wrapownshare.hh>
#include "solverfactory.hh"
-template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
-SolverFactory<DeformedGrid, Nonlinearity, MatrixType, VectorType>::SolverFactory(
- const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
- const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes)
- : nBodyAssembler_(nBodyAssembler),
- baseEnergyNorm_(baseSolverStep_),
- baseSolver_(&baseSolverStep_,
- parset.get<size_t>("linear.maxiumumIterations"),
- parset.get<double>("linear.tolerance"), &baseEnergyNorm_,
- Solver::QUIET),
- transferOperators_(nBodyAssembler.getGrids().at(0)->maxLevel()) {
-
- /* baseSolverStep_.obstacles_ = using HasObstacle = Dune::BitSetVector<BlockSize>;
- using Obstacle = BoxConstraint<Field, BlockSize>;
- const HasObstacle* hasObstacle_;
- const std::vector<Obstacle>* obstacles_;*/
-
- multigridStep_ = std::make_shared<Step>();
- // tnnmg iteration step
- multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
- multigridStep_->setIgnore(ignoreNodes);
- multigridStep_->setBaseSolver(baseSolver_);
- multigridStep_->setSmoother(&presmoother_, &postsmoother_);
- multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
- multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
- multigridStep_->setVerbosity(NumProc::QUIET);
+template <class Functional, class BitVector>
+SolverFactory<Functional, BitVector>::SolverFactory(
+ const Dune::ParameterTree& parset,
+ Functional& J) :
+ J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))),
+ //linear solver
+ //preconditioner_(),
+ linearSolverStep_(),
+ energyNorm_(linearSolverStep_)
+ {
+ nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
- // create the transfer operators
- const int nCouplings = nBodyAssembler_.nCouplings();
- const auto grids = nBodyAssembler_.getGrids();
+ // linearSolverStep_.setPreconditioner(preconditioner_);
+ linearSolver_ = std::make_shared<LinearSolver>(linearSolverStep_, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), energyNorm_, Solver::QUIET);
- for (size_t i=0; i<transferOperators_.size(); i++) {
- transferOperators_[i] = std::make_shared<TransferOperator>();
- }
-
- std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nCouplings);
- std::vector<const MatrixType*> mortarTransfer(nCouplings);
- std::vector<std::array<int,2> > gridIdx(nCouplings);
-
- for (int j=0; j<nCouplings; j++) {
- fineHasObstacle[j] = nBodyAssembler_.contactCoupling_[j]->nonmortarBoundary().getVertices();
- mortarTransfer[j] = &nBodyAssembler_.contactCoupling_[j]->mortarLagrangeMatrix();
- gridIdx[j] = nBodyAssembler_.coupling_[j].gridIdx_;
- }
-
- TransferOperator::setupHierarchy(transferOperators_,
- grids,
- mortarTransfer,
- nBodyAssembler_.localCoordSystems_,
- fineHasObstacle,
- gridIdx);
+ tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_, DefectProjection(), LineSearchSolver());
+ tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
+}
- multigridStep_->setTransferOperators(transferOperators_);
+template <class Functional, class BitVector>
+void SolverFactory<Functional, BitVector>::setProblem(Vector& x) {
+ nonlinearSmoother_->setProblem(x);
+ tnnmgStep_->setProblem(x);
}
-template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
-auto SolverFactory<DeformedGrid, Nonlinearity, MatrixType, VectorType>::getStep()
- -> std::shared_ptr<Step> {
- return multigridStep_;
+
+template <class Functional, class BitVector>
+auto SolverFactory<Functional, BitVector>::step()
+-> std::shared_ptr<Step> {
+ return tnnmgStep_;
}
#include "solverfactory_tmpl.cc"
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 5e87ea27d7920dbab7b044a485c5a7d4f8f159e7..906649b6337a9e8cc10ba68bb25144c8028641ea 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -1,66 +1,66 @@
#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
#define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
+#define NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
-#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
-//#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
+#include <dune/tnnmg/projections/obstacledefectprojection.hh>
-#include <dune/contact/assemblers/nbodyassembler.hh>
+#include "tnnmg/functional.hh"
+#include "tnnmg/linearization.hh"
+#include "tnnmg/linesearchsolver.hh"
+#include "tnnmg/localbisectionsolver.hh"
-//#include <dune/contact/estimators/hierarchiccontactestimator.hh>
-#include <dune/contact/solvers/nsnewtonmgstep.hh>
-#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
-#include <dune/contact/solvers/contactobsrestrict.hh>
-#include <dune/contact/solvers/contacttransferoperatorassembler.hh>
-#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
+template <class FunctionalTEMPLATE, class BitVectorType>
+class SolverFactory {
+public:
+ using Functional = FunctionalTEMPLATE;
+ using Matrix = typename Functional::Matrix;
+ using Vector = typename Functional::Vector;
+ using BitVector = BitVectorType;
-#define USE_OLD_TNNMG //needed for old bisection.hh in tnnmg
+ using LocalSolver = LocalBisectionSolver;
+ using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalSolver>;
-template <class DeformedGridTEMPLATE, class NonlinearityTEMPLATE, class MatrixType, class VectorType>
-class SolverFactory {
-//protected:
-// const size_t dim = DeformedGrid::dimension;
+ using Linearization = typename Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<Functional, BitVector>;
+ //using Preconditioner = Dune::Solvers::Preconditioner; //TODO Platzhalter für PatchPreconditioner
+ using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
+ using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
-public:
- using Matrix = MatrixType;
- using Vector = VectorType;
- using DeformedGrid = DeformedGridTEMPLATE;
- using Nonlinearity = NonlinearityTEMPLATE;
+ using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
+
+ using Step = typename Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
-public:
- using Step = Dune::Contact::NonSmoothNewtonMGStep<MatrixType, VectorType>;
- using TransferOperator = Dune::Contact::NonSmoothNewtonContactTransfer<VectorType>;
- SolverFactory(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
- const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes);
+ SolverFactory(Dune::ParameterTree const &,
+ Functional&);
- std::shared_ptr<Step> getStep();
+ void setProblem(Vector& x);
- const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& getNBodyAssembler() const {
- return nBodyAssembler_;
- }
+ auto step() -> std::shared_ptr<Step>;
private:
- const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler_;
+ Vector dummyIterate_;
+ std::shared_ptr<const Functional> J_;
- //ProjectedBlockGSStep<MatrixType, VectorType> baseSolverStep_;
- BlockGSStep<MatrixType, VectorType> baseSolverStep_;
- EnergyNorm<MatrixType, VectorType> baseEnergyNorm_;
- LoopSolver<VectorType> baseSolver_;
+ // nonlinear smoother
+ std::shared_ptr<NonlinearSmoother> nonlinearSmoother_;
- ProjectedBlockGSStep<MatrixType, VectorType> presmoother_;
- ProjectedBlockGSStep<MatrixType, VectorType> postsmoother_;
- std::shared_ptr<Step> multigridStep_;
+ // linear solver
+ //Preconditioner preconditioner_;
+ LinearSolverStep linearSolverStep_;
+ EnergyNorm<Matrix, Vector> energyNorm_;
+ std::shared_ptr<LinearSolver> linearSolver_;
- std::vector<std::shared_ptr<TransferOperator>> transferOperators_;
+ // TNNMGStep
+ std::shared_ptr<Step> tnnmgStep_;
};
#endif
diff --git a/src/spatial-solving/solverfactory_old.cc b/src/spatial-solving/solverfactory_old.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cbbe5ef1d7cbabc1c2ecd56677dad245192ed600
--- /dev/null
+++ b/src/spatial-solving/solverfactory_old.cc
@@ -0,0 +1,78 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+#include <dune/solvers/solvers/solver.hh>
+
+#include "solverfactory_old.hh"
+
+
+template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
+SolverFactoryOld<DeformedGrid, Nonlinearity, MatrixType, VectorType>::SolverFactoryOld(
+ const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
+ const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes)
+ : nBodyAssembler_(nBodyAssembler),
+ baseEnergyNorm_(baseSolverStep_),
+ baseSolver_(&baseSolverStep_,
+ parset.get<size_t>("linear.maximumIterations"),
+ parset.get<double>("linear.tolerance"), &baseEnergyNorm_,
+ Solver::QUIET),
+ transferOperators_(nBodyAssembler.getGrids().at(0)->maxLevel()) {
+
+ /* baseSolverStep_.obstacles_ = using HasObstacle = Dune::BitSetVector<BlockSize>;
+ using Obstacle = BoxConstraint<Field, BlockSize>;
+ const HasObstacle* hasObstacle_;
+ const std::vector<Obstacle>* obstacles_;*/
+
+ multigridStep_ = std::make_shared<Step>();
+ // tnnmg iteration step
+ multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
+ multigridStep_->setIgnore(ignoreNodes);
+ multigridStep_->setBaseSolver(baseSolver_);
+ multigridStep_->setSmoother(&presmoother_, &postsmoother_);
+ multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
+ multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
+ multigridStep_->setVerbosity(NumProc::QUIET);
+
+ // create the transfer operators
+ const int nCouplings = nBodyAssembler_.nCouplings();
+ const auto grids = nBodyAssembler_.getGrids();
+
+ for (size_t i=0; i<transferOperators_.size(); i++) {
+ transferOperators_[i] = std::make_shared<TransferOperator>();
+ }
+
+ std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nCouplings);
+ std::vector<const MatrixType*> mortarTransfer(nCouplings);
+ std::vector<std::array<int,2> > gridIdx(nCouplings);
+
+ for (int j=0; j<nCouplings; j++) {
+ fineHasObstacle[j] = nBodyAssembler_.contactCoupling_[j]->nonmortarBoundary().getVertices();
+ mortarTransfer[j] = &nBodyAssembler_.contactCoupling_[j]->mortarLagrangeMatrix();
+ gridIdx[j] = nBodyAssembler_.coupling_[j].gridIdx_;
+ }
+
+ TransferOperator::setupHierarchy(transferOperators_,
+ grids,
+ mortarTransfer,
+ nBodyAssembler_.localCoordSystems_,
+ fineHasObstacle,
+ gridIdx);
+
+ multigridStep_->setTransferOperators(transferOperators_);
+}
+
+template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
+auto SolverFactoryOld<DeformedGrid, Nonlinearity, MatrixType, VectorType>::getStep()
+ -> std::shared_ptr<Step> {
+ return multigridStep_;
+}
+
+#include "solverfactory_tmpl_old.cc"
diff --git a/src/spatial-solving/solverfactory_old.hh b/src/spatial-solving/solverfactory_old.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1c383419e1ac32eced29512d04260881f7157a82
--- /dev/null
+++ b/src/spatial-solving/solverfactory_old.hh
@@ -0,0 +1,66 @@
+#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORYOLD_HH
+#define SRC_SPATIAL_SOLVING_SOLVERFACTORYOLD_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/blockgsstep.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+//#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
+//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+//#include <dune/contact/estimators/hierarchiccontactestimator.hh>
+#include <dune/contact/solvers/nsnewtonmgstep.hh>
+#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
+#include <dune/contact/solvers/contactobsrestrict.hh>
+#include <dune/contact/solvers/contacttransferoperatorassembler.hh>
+#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
+
+#define USE_OLD_TNNMG //needed for old bisection.hh in tnnmg
+
+template <class DeformedGridTEMPLATE, class NonlinearityTEMPLATE, class MatrixType, class VectorType>
+class SolverFactoryOld {
+//protected:
+// const size_t dim = DeformedGrid::dimension;
+
+public:
+ using Matrix = MatrixType;
+ using Vector = VectorType;
+ using DeformedGrid = DeformedGridTEMPLATE;
+ using Nonlinearity = NonlinearityTEMPLATE;
+
+public:
+ using Step = Dune::Contact::NonSmoothNewtonMGStep<MatrixType, VectorType>;
+ using TransferOperator = Dune::Contact::NonSmoothNewtonContactTransfer<VectorType>;
+
+ SolverFactoryOld(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
+ const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes);
+
+ std::shared_ptr<Step> getStep();
+
+ const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& getNBodyAssembler() const {
+ return nBodyAssembler_;
+ }
+
+private:
+ const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler_;
+
+ //ProjectedBlockGSStep<MatrixType, VectorType> baseSolverStep_;
+ BlockGSStep<MatrixType, VectorType> baseSolverStep_;
+ EnergyNorm<MatrixType, VectorType> baseEnergyNorm_;
+ LoopSolver<VectorType> baseSolver_;
+
+ ProjectedBlockGSStep<MatrixType, VectorType> presmoother_;
+ ProjectedBlockGSStep<MatrixType, VectorType> postsmoother_;
+ std::shared_ptr<Step> multigridStep_;
+
+ std::vector<std::shared_ptr<TransferOperator>> transferOperators_;
+};
+#endif
diff --git a/src/spatial-solving/solverfactory_tmpl.cc b/src/spatial-solving/solverfactory_tmpl.cc
index 7e06390bb45510eb2a9b4b209a23db01d29a4b1c..60fc8802b7242b1c27eae50ee0bde47a5e715708 100644
--- a/src/spatial-solving/solverfactory_tmpl.cc
+++ b/src/spatial-solving/solverfactory_tmpl.cc
@@ -5,16 +5,10 @@
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
+//#include "../../dune/tectonic/globalratestatefriction.hh"
+#include "tnnmg/functional.hh"
+#include "tnnmg/zerononlinearity.hh"
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
+using FrictionlessFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity, Vector&, Vector&, double>;
+template class SolverFactory<FrictionlessFunctional, BitVector>;
-template class SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
-
-/*template class SolverFactory<
- MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
- ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
- Grid>;*/
diff --git a/src/spatial-solving/solverfactory_tmpl_old.cc b/src/spatial-solving/solverfactory_tmpl_old.cc
new file mode 100644
index 0000000000000000000000000000000000000000..600132d31b71a6fda13c7d955374bf98d01f9412
--- /dev/null
+++ b/src/spatial-solving/solverfactory_tmpl_old.cc
@@ -0,0 +1,16 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+
+#include <dune/tectonic/globalfriction.hh>
+
+template class SolverFactoryOld<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+
+/*template class SolverFactory<
+ MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
+ ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
+ Grid>;*/
diff --git a/src/spatial-solving/tnnmg/functional.hh b/src/spatial-solving/tnnmg/functional.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8b41cd3e7002eaf287dfbdc0c3c26fcb91ad5825
--- /dev/null
+++ b/src/spatial-solving/tnnmg/functional.hh
@@ -0,0 +1,420 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=8 sw=2 et sts=2:
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_FUNCTIONAL_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_FUNCTIONAL_HH
+
+#include <cstddef>
+#include <type_traits>
+
+#include <dune/solvers/common/wrapownshare.hh>
+#include <dune/solvers/common/copyorreference.hh>
+#include <dune/solvers/common/interval.hh>
+#include <dune/solvers/common/resize.hh>
+
+#include <dune/matrix-vector/algorithm.hh>
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+
+/** \brief Coordinate restriction of box constrained quadratic functional with nonlinearity;
+ * mainly used for presmoothing in TNNMG algorithm
+ *
+ * \tparam M global matrix type
+ * \tparam V global vector type
+ * \tparam V nonlinearity type
+ * \tparam R Range type
+ */
+template<class M, class V, class N, class R>
+class ShiftedFunctional : public Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunctional<M,V,R> {
+ using Base = Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunctional<M,V,R>;
+
+public:
+ using Matrix = M;
+ using Vector = V;
+ using Nonlinearity = N;
+ using LowerObstacle = Vector;
+ using UpperObstacle = Vector;
+ using Range = R;
+
+ ShiftedFunctional(const Matrix& quadraticPart, const Vector& linearPart, const Nonlinearity& phi, const LowerObstacle& lowerObstacle, const UpperObstacle& upperObstacle, const Vector& origin) :
+ Base(quadraticPart, linearPart, lowerObstacle, upperObstacle, origin),
+ phi_(phi)
+ {}
+
+ Range operator()(const Vector& v) const
+ {
+ auto temp = Base::Base::origin();
+ temp += v;
+ if (checkInsideBox(temp, this->originalLowerObstacle(), this->originalUpperObstacle()))
+ return Base::Base::operator()(v) + phi_.operator()(temp);
+ return std::numeric_limits<Range>::max();
+ }
+
+ const auto& phi() const {
+ return phi_;
+ }
+
+protected:
+ const Nonlinearity& phi_;
+};
+
+
+
+
+
+/** \brief Coordinate restriction of box constrained quadratic functional with nonlinearity;
+ * mainly used for line search step in TNNMG algorithm
+ *
+ * \tparam M Global matrix type
+ * \tparam V Global vector type
+ * \tparam N nonlinearity type
+ * \tparam L Global lower obstacle type
+ * \tparam U Global upper obstacle type
+ * \tparam R Range type
+ */
+template<class M, class V, class N, class L, class U, class R=double>
+class DirectionalRestriction :
+ public Dune::TNNMG::BoxConstrainedQuadraticDirectionalRestriction<M,V,L,U,R>
+{
+ using Base = Dune::TNNMG::BoxConstrainedQuadraticDirectionalRestriction<M,V,L,U,R>;
+ using Nonlinearity = N;
+ using Interval = typename Dune::Solvers::Interval<R>;
+
+public:
+ using GlobalMatrix = typename Base::GlobalMatrix;
+ using GlobalVector = typename Base::GlobalVector;
+ using GlobalLowerObstacle = typename Base::GlobalLowerObstacle;
+ using GlobalUpperObstacle = typename Base::GlobalUpperObstacle;
+
+ using Matrix = typename Base::Matrix;
+ using Vector = typename Base::Vector;
+
+ DirectionalRestriction(const GlobalMatrix& matrix, const GlobalVector& linearTerm, const Nonlinearity& phi, const GlobalLowerObstacle& lower, const GlobalLowerObstacle& upper, const GlobalVector& origin, const GlobalVector& direction) :
+ Base(matrix, linearTerm, lower, upper, origin, direction),
+ origin_(origin),
+ direction_(direction),
+ phi_(phi)
+ {
+ phi_.directionalDomain(origin_, direction_, domain_);
+
+ if (domain_[0] < this->defectLower_) {
+ domain_[0] = this->defectLower_;
+ }
+
+ if (domain_[1] > this->defectUpper_) {
+ domain_[1] = this->defectUpper_;
+ }
+ }
+
+ auto subDifferential(double x) const {
+ Interval D;
+ GlobalVector uxv = origin_;
+
+ Dune::MatrixVector::addProduct(uxv, x, direction_);
+ phi_.directionalSubDiff(uxv, direction_, D);
+ auto const Axmb = this->quadraticPart_ * x - this->linearPart_;
+ D[0] += Axmb;
+ D[1] += Axmb;
+
+ return D;
+ }
+
+ auto domain() const {
+ return domain_;
+ }
+
+ const GlobalVector& origin() const {
+ return origin_;
+ }
+
+ const GlobalVector& direction() const {
+ return direction_;
+ }
+
+protected:
+ const GlobalVector& origin_;
+ const GlobalVector& direction_;
+
+ const Nonlinearity& phi_;
+ Interval domain_;
+};
+
+/** \brief Box constrained quadratic functional with nonlinearity
+ * Note: call setIgnore() to set up functional in affine subspace with ignore information
+ *
+ * \tparam V Vector type
+ * \tparam N Nonlinearity type
+ * \tparam L Lower obstacle type
+ * \tparam U Upper obstacle type
+ * \tparam R Range type
+ */
+template<class V, class N, class L, class U, class O, class D, class R>
+class FirstOrderModelFunctional {
+ using Interval = typename Dune::Solvers::Interval<R>;
+
+public:
+ using Vector = std::decay_t<V>;
+ using Nonlinearity = std::decay_t<N>;
+ using LowerObstacle = std::decay_t<L>;
+ using UpperObstacle = std::decay_t<U>;
+ using Range = R;
+
+ using field_type = typename Vector::field_type;
+
+private:
+ void init() {
+ auto origin = origin_.get();
+ auto direction = direction_.get();
+
+ // set defect obstacles
+ defectLower_ = -std::numeric_limits<field_type>::max();
+ defectUpper_ = std::numeric_limits<field_type>::max();
+ Dune::TNNMG::directionalObstacles(origin, direction, lower_.get(), upper_.get(), defectLower_, defectUpper_);
+
+ // set domain
+ phi_.get().directionalDomain(origin, direction, domain_);
+
+ if (domain_[0] < defectLower_) {
+ domain_[0] = defectLower_;
+ }
+
+ if (domain_[1] > defectUpper_) {
+ domain_[1] = defectUpper_;
+ }
+
+ // set quadratic and linear parts of functional
+ quadraticPart_ = direction*direction;
+ quadraticPart_ *= maxEigenvalue_;
+
+ linearPart_ = linearTerm_.get()*direction - maxEigenvalue_*(direction*origin);
+ }
+
+public:
+ template <class MM, class VV1, class NN, class LL, class UU, class VV2, class VV3>
+ FirstOrderModelFunctional(
+ const MM& matrix,
+ VV1&& linearTerm,
+ NN&& phi,
+ LL&& lower,
+ UU&& upper,
+ VV2&& origin,
+ VV3&& direction) :
+ linearTerm_(std::forward<VV1>(linearTerm)),
+ lower_(std::forward<LL>(lower)),
+ upper_(std::forward<UU>(upper)),
+ origin_(std::forward<VV2>(origin)),
+ direction_(std::forward<VV3>(direction)),
+ phi_(std::forward<NN>(phi)) {
+
+ // set maxEigenvalue_ from matrix
+ Vector eigenvalues;
+ Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
+ maxEigenvalue_ = *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+
+ init();
+ }
+
+ template <class VV1, class NN, class LL, class UU, class VV2, class VV3>
+ FirstOrderModelFunctional(
+ const Range& maxEigenvalue,
+ VV1&& linearTerm,
+ NN&& phi,
+ LL&& lower,
+ UU&& upper,
+ VV2&& origin,
+ VV3&& direction) :
+ linearTerm_(std::forward<VV1>(linearTerm)),
+ lower_(std::forward<LL>(lower)),
+ upper_(std::forward<UU>(upper)),
+ origin_(std::forward<VV2>(origin)),
+ direction_(std::forward<VV3>(direction)),
+ phi_(std::forward<NN>(phi)),
+ maxEigenvalue_(maxEigenvalue) {
+
+ init();
+ }
+
+ template <class BitVector>
+ auto getIgnoreFunctional(const BitVector& ignore) const {
+ Vector direction = direction_.get();
+ Vector origin = origin_.get();
+ // Dune::Solvers::resizeInitializeZero(direction, linearPart());
+ for (size_t j = 0; j < ignore.size(); ++j)
+ if (ignore[j])
+ direction[j] = 0; // makes sure result remains in subspace after correction
+ else
+ origin[j] = 0; // shift affine offset
+
+ return FirstOrderModelFunctional<Vector, Nonlinearity&, LowerObstacle, UpperObstacle, Vector, Vector, Range>(maxEigenvalue_, linearTerm_, phi_, lower_, upper_, std::move(origin), std::move(direction));
+ }
+
+ Range operator()(const Vector& v) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderModelFunctional not implemented");
+ }
+
+ auto subDifferential(double x) const {
+ Interval Di;
+ Vector uxv = origin_.get();
+
+ Dune::MatrixVector::addProduct(uxv, x, direction_.get());
+ phi_.get().directionalSubDiff(uxv, direction_.get(), Di);
+
+ const auto Axmb = quadraticPart_ * x - linearPart_;
+ Di[0] += Axmb;
+ Di[1] += Axmb;
+
+ return Di;
+ }
+
+ const Interval& domain() const {
+ return domain_;
+ }
+
+ const Vector& origin() const {
+ return origin_.get();
+ }
+
+ const Vector& direction() const {
+ return direction_.get();
+ }
+
+ const field_type& defectLower() const {
+ return defectLower_;
+ }
+
+ const field_type& defectUpper() const {
+ return defectUpper_;
+ }
+
+private:
+ Dune::Solvers::ConstCopyOrReference<V> linearTerm_;
+ Dune::Solvers::ConstCopyOrReference<L> lower_;
+ Dune::Solvers::ConstCopyOrReference<U> upper_;
+
+ Dune::Solvers::ConstCopyOrReference<O> origin_;
+ Dune::Solvers::ConstCopyOrReference<D> direction_;
+
+ Dune::Solvers::ConstCopyOrReference<N> phi_;
+
+ Interval domain_;
+
+ Range maxEigenvalue_;
+
+ field_type quadraticPart_;
+ field_type linearPart_;
+
+ field_type defectLower_;
+ field_type defectUpper_;
+};
+
+// \ToDo This should be an inline friend of ShiftedBoxConstrainedQuadraticFunctional
+// but gcc-4.9.2 shipped with debian jessie does not accept
+// inline friends with auto return type due to bug-59766.
+// Notice, that this is fixed in gcc-4.9.3.
+template<class Index, class M, class V, class Nonlinearity, class R>
+auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, const Index& i)
+{
+ using Range = R;
+ using LocalMatrix = std::decay_t<decltype(f.quadraticPart()[i][i])>;
+ using LocalVector = std::decay_t<decltype(f.originalLinearPart()[i])>;
+ using LocalLowerObstacle = std::decay_t<decltype(f.originalLowerObstacle()[i])>;
+ using LocalUpperObstacle = std::decay_t<decltype(f.originalUpperObstacle()[i])>;
+
+ using namespace Dune::MatrixVector;
+ namespace H = Dune::Hybrid;
+
+ const LocalMatrix* Aii_p = nullptr;
+
+ LocalVector ri = f.originalLinearPart()[i];
+ const auto& Ai = f.quadraticPart()[i];
+ sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+ // TODO Here we must implement a wrapper to guarantee that this will work with proxy matrices!
+ H::ifElse(H::equals(j, i), [&](auto&& id){
+ Aii_p = id(&Aij);
+ });
+ Dune::TNNMG::Imp::mmv(Aij, f.origin()[j], ri, Dune::PriorityTag<1>());
+ });
+
+ LocalLowerObstacle dli = f.originalLowerObstacle()[i];
+ LocalUpperObstacle dui = f.originalUpperObstacle()[i];
+ dli -= f.origin()[i];
+ dui -= f.origin()[i];
+
+ auto&& phii = f.phi().restriction(i);
+
+ auto v = ri;
+ double const vnorm = v.two_norm();
+ if (vnorm > 1.0)
+ v /= vnorm;
+
+ return FirstOrderModelFunctional<LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, LocalVector, Range>(*Aii_p, std::move(ri), std::move(phii), std::move(dli), std::move(dui), std::move(f.origin()[i]), std::move(v));
+}
+
+
+/** \brief Box constrained quadratic functional with nonlinearity
+ *
+ * \tparam M Matrix type
+ * \tparam V Vector type
+ * \tparam L Lower obstacle type
+ * \tparam U Upper obstacle type
+ * \tparam R Range type
+ */
+template<class M, class V, class N, class L, class U, class R>
+class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L, U, R>
+{
+private:
+ using Base = Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L, U, R>;
+
+
+public:
+ using Nonlinearity = std::decay_t<N>;
+
+ using Matrix = typename Base::Matrix;
+ using Vector = typename Base::Vector;
+ using Range = typename Base::Range;
+ using LowerObstacle = typename Base::LowerObstacle;
+ using UpperObstacle = typename Base::UpperObstacle;
+
+private:
+ Dune::Solvers::ConstCopyOrReference<N> phi_;
+
+public:
+
+ template <class MM, class VV, class NN, class LL, class UU>
+ Functional(
+ MM&& matrix,
+ VV&& linearPart,
+ NN&& phi,
+ LL& lower,
+ UU& upper) :
+ Base(std::forward<MM>(matrix), std::forward<VV>(linearPart), std::forward<LL>(lower), std::forward<UU>(upper)),
+ phi_(std::forward<NN>(phi))
+ {}
+
+ const Nonlinearity& phi() const {
+ return phi_.get();
+ }
+
+ Range operator()(const Vector& v) const
+ {
+ if (Dune::TNNMG::checkInsideBox(v, this->lower_.get(), this->upper_.get()))
+ return Base::Base::operator()(v) + phi_.get().operator()(v);
+ return std::numeric_limits<Range>::max();
+ }
+
+ friend auto directionalRestriction(const Functional& f, const Vector& origin, const Vector& direction)
+ -> DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>
+ {
+ return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin, direction);
+ }
+
+ friend auto shift(const Functional& f, const Vector& origin)
+ -> ShiftedFunctional<Matrix, Vector, Nonlinearity, Range>
+ {
+ return ShiftedFunctional<Matrix, Vector, Nonlinearity, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
+ }
+};
+
+
+#endif
diff --git a/src/spatial-solving/tnnmg/levelglobalpreconditioner.hh b/src/spatial-solving/tnnmg/levelglobalpreconditioner.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2b7b4e2bb80f6ad5ec302097654dbb5bc320deb4
--- /dev/null
+++ b/src/spatial-solving/tnnmg/levelglobalpreconditioner.hh
@@ -0,0 +1,158 @@
+#ifndef LEVEL_GLOBAL_PRECONDITIONER_HH
+#define LEVEL_GLOBAL_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+
+#include <dune/istl/umfpack.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.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>
+class LevelGlobalPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+
+public:
+ enum BoundaryMode {homogeneous, fromIterate};
+
+private:
+ typedef typename BasisType::GridView GridView;
+ typedef typename GridView::Grid GridType;
+
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
+ const LocalAssembler& localAssembler_;
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
+
+ const GridType& grid_;
+ const int level_;
+ const BasisType basis_;
+
+ Dune::BitSetVector<1> boundaryDofs_;
+ MatrixType matrix_;
+ VectorType rhs_;
+
+ bool requireBuild_;
+
+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_);
+
+ GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
+ globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
+
+ typedef typename MatrixType::row_type RowType;
+ typedef typename RowType::ConstIterator ColumnIterator;
+
+ const GridView& gridView = levelInterfaceNetwork_.levelGridView();
+
+ // set boundary conditions
+ BoundaryPatch<GridView> boundaryPatch(gridView, true);
+ constructBoundaryDofs(boundaryPatch, basis_, boundaryDofs_);
+
+ for(size_t i=0; i<boundaryDofs_.size(); i++) {
+ if(!boundaryDofs_[i][0])
+ continue;
+
+ RowType& row = matrix_[i];
+
+ ColumnIterator cIt = row.begin();
+ ColumnIterator cEndIt = row.end();
+
+ for(; cIt!=cEndIt; ++cIt) {
+ row[cIt.index()] = 0;
+ }
+
+ row[i] = 1;
+ }
+
+ 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!");
+
+ Dune::Timer timer;
+ timer.reset();
+ timer.start();
+
+ this->x_->resize(matrix_.M());
+ *(this->x_) = 0;
+
+ #if HAVE_UMFPACK
+ Dune::InverseOperatorResult res;
+ VectorType rhsCopy(rhs_);
+
+ Dune::UMFPack<MatrixType> solver(matrix_);
+ solver.apply(*(this->x_), rhsCopy, res);
+
+ #else
+ #error No UMFPack!
+ #endif
+
+ //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/tnnmg/levelpatchpreconditioner.hh b/src/spatial-solving/tnnmg/levelpatchpreconditioner.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bbe464fdec9c3d9db5c218cecdbe4c448cd06fad
--- /dev/null
+++ b/src/spatial-solving/tnnmg/levelpatchpreconditioner.hh
@@ -0,0 +1,261 @@
+#ifndef LEVEL_PATCH_PRECONDITIONER_HH
+#define LEVEL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+
+#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
+#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
+#include <dune/faultnetworks/localproblem.hh>
+#include <dune/faultnetworks/levelinterfacenetwork.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+
+
+template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+
+public:
+ enum Mode {additive, multiplicative};
+ enum BoundaryMode {homogeneous, fromIterate};
+
+private:
+ typedef typename BasisType::GridView GridView;
+ typedef typename GridView::Grid GridType;
+
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
+ const BasisType& patchLevelBasis_;
+ const LocalAssembler& localAssembler_;
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
+ const Mode mode_;
+
+ const GridType& grid_;
+ const int level_;
+ const BasisType basis_;
+
+ size_t patchDepth_;
+ BoundaryMode boundaryMode_;
+
+ MatrixType matrix_;
+ std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
+
+public:
+
+ // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
+ LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+ const BasisType& patchLevelBasis,
+ const LocalAssembler& localAssembler,
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+ const Mode mode = LevelPatchPreconditioner::Mode::additive) :
+ levelInterfaceNetwork_(levelInterfaceNetwork),
+ patchLevelBasis_(patchLevelBasis),
+ localAssembler_(localAssembler),
+ localInterfaceAssemblers_(localInterfaceAssemblers),
+ mode_(mode),
+ grid_(levelInterfaceNetwork_.grid()),
+ level_(levelInterfaceNetwork_.level()),
+ basis_(levelInterfaceNetwork_)
+ {
+ setPatchDepth();
+ setBoundaryMode();
+
+ assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
+ }
+
+ void build() {
+ // assemble stiffness matrix for entire level including all faults
+ GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
+ globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
+
+ // set boundary conditions
+ Dune::BitSetVector<1> globalBoundaryDofs;
+ BoundaryPatch<GridView> boundaryPatch(levelInterfaceNetwork_.levelGridView(), true);
+ constructBoundaryDofs(boundaryPatch, basis_, globalBoundaryDofs);
+
+ typedef typename MatrixType::row_type RowType;
+ typedef typename RowType::ConstIterator ColumnIterator;
+
+ for(size_t i=0; i<globalBoundaryDofs.size(); i++) {
+ if(!globalBoundaryDofs[i][0])
+ continue;
+
+ RowType& row = matrix_[i];
+
+ ColumnIterator cIt = row.begin();
+ ColumnIterator cEndIt = row.end();
+
+ for(; cIt!=cEndIt; ++cIt) {
+ row[cIt.index()] = 0;
+ }
+ row[i] = 1;
+ }
+
+ // init vertexInElements
+ const int dim = GridType::dimension;
+ typedef typename GridType::template Codim<0>::Entity EntityType;
+ std::vector<std::vector<EntityType>> vertexInElements;
+
+ const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
+ vertexInElements.resize(patchLevelGridView.size(dim));
+ for (size_t i=0; i<vertexInElements.size(); i++) {
+ vertexInElements[i].resize(0);
+ }
+
+ typedef typename BasisType::LocalFiniteElement FE;
+ typedef typename GridView::template Codim <0>::Iterator ElementLevelIterator;
+ ElementLevelIterator endElemIt = patchLevelGridView.template end <0>();
+ for (ElementLevelIterator it = patchLevelGridView. template begin <0>(); it!=endElemIt; ++it) {
+
+ // compute coarseGrid vertexInElements
+ const FE& coarseFE = patchLevelBasis_.getLocalFiniteElement(*it);
+
+ for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
+ int dofIndex = patchLevelBasis_.indexInGridView(*it, i);
+ vertexInElements[dofIndex].push_back(*it);
+ }
+ }
+
+ localProblems_.resize(vertexInElements.size());
+
+ std::cout << std::endl;
+ std::cout << "---------------------------------------------" << std::endl;
+ std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
+
+ // init local fine level corrections
+ Dune::Timer timer;
+ timer.reset();
+ timer.start();
+
+ const int patchLevel = patchLevelBasis_.faultNetwork().level();
+ for (size_t i=0; i<vertexInElements.size(); i++) {
+ //std::cout << i << std::endl;
+ //std::cout << "---------------" << std::endl;
+
+ SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, basis_, patchLevel, level_, vertexInElements, i, patchDepth_);
+ std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
+ Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
+ 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);
+
+ /*
+ if ((i+1) % 10 == 0) {
+ std::cout << '\r' << std::floor((i+1.0)/patchLevelBasis_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*patchLevelBasis_.size() << " seconds." << std::flush;
+ }*/
+ }
+
+ std::cout << std::endl;
+ std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
+ std::cout << "---------------------------------------------" << std::endl;
+ timer.stop();
+ }
+
+ void setPatchDepth(const size_t patchDepth = 0) {
+ patchDepth_ = patchDepth;
+ }
+
+ void setBoundaryMode(const BoundaryMode boundaryMode = LevelPatchPreconditioner::BoundaryMode::homogeneous) {
+ boundaryMode_ = boundaryMode;
+ }
+
+ 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);
+
+ for (size_t i=0; i<localProblems_.size(); i++) {
+ if (boundaryMode_ == BoundaryMode::homogeneous)
+ localProblems_[i]->updateRhs(rhs);
+ else
+ localProblems_[i]->updateRhsAndBoundary(rhs, x);
+ }
+ }
+
+ virtual void iterate() {
+ if (mode_ == additive)
+ iterateAdd();
+ else
+ iterateMult();
+ }
+
+ void iterateAdd() {
+ *(this->x_) = 0;
+
+ VectorType it, x;
+ for (size_t i=0; i<localProblems_.size(); i++) {
+ localProblems_[i]->solve(it);
+ localProblems_[i]->prolong(it, x);
+
+ /*if (i==5) {
+ writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
+ }*/
+
+ *(this->x_) += x;
+ }
+ }
+
+ void iterateMult() {
+ *(this->x_) = 0;
+
+ VectorType it, x;
+ for (size_t i=0; i<localProblems_.size(); i++) {
+ VectorType updatedRhs(*(this->rhs_));
+ matrix_.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));
+
+ if (boundaryMode_ == BoundaryMode::homogeneous)
+ localProblems_[i]->updateRhs(updatedRhs);
+ else
+ localProblems_[i]->updateRhsAndBoundary(updatedRhs, *(this->x_));
+
+ localProblems_[i]->solve(it);
+ localProblems_[i]->prolong(it, x);
+
+
+ //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;
+ }
+ }
+
+ const BasisType& basis() const {
+ return basis_;
+ }
+
+ const GridView& gridView() const {
+ return basis_.getGridView();
+ }
+
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
+ return levelInterfaceNetwork_;
+ }
+
+ size_t size() const {
+ return localProblems_.size();
+ }
+};
+
+#endif
+
diff --git a/src/spatial-solving/tnnmg/linearization.hh b/src/spatial-solving/tnnmg/linearization.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4e7ce13238b727fec2ebdf91da04cdfd51aa9a0a
--- /dev/null
+++ b/src/spatial-solving/tnnmg/linearization.hh
@@ -0,0 +1,99 @@
+#ifndef DUNE_TECTONIC_LINEARIZATION_HH
+#define DUNE_TECTONIC_LINEARIZATION_HH
+
+#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
+#include <dune/istl/matrixindexset.hh>
+
+//#include <>
+
+/*#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/fufem/arithmetic.hh>
+#include <dune/solvers/common/interval.hh>
+#include <dune/solvers/computeenergy.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+#include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
+
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/minimisation.hh>
+#include <dune/tectonic/mydirectionalconvexfunction.hh>
+#include <dune/tectonic/quadraticenergy.hh>*/
+
+
+template<class F, class BV>
+class Linearization : public Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<F,BV> {
+private:
+ using Base = Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<F,BV>;
+ using Vector = typename Base::Vector;
+ using BitVector = typename Base::BitVector;
+
+ template <class T>
+ void determineRegularityTruncation(const Vector& x, BitVector&& truncationFlags, const T& truncationTolerance)
+ {
+ for (size_t i = 0; i < x.size(); ++i) {
+ if (this->f_.phi().regularity(i, x[i]) > truncationTolerance) {
+ truncationFlags[i] = true;
+ }
+ }
+ }
+
+public:
+ Linearization(const F& f, const BitVector& ignore) :
+ Base(f, ignore)
+ {}
+
+ void bind(const Vector& x) {
+ auto&& A = this->f_.quadraticPart();
+ auto&& phi = this->f_.phi();
+
+ // compute hessian
+ // ---------------
+
+ // construct sparsity pattern for linearization
+ Dune::MatrixIndexSet indices(A.N(), A.M());
+ indices.import(A);
+ phi.addHessianIndices(indices);
+
+ // construct matrix from pattern and initialize it
+ indices.exportIdx(this->hessian_);
+ this->hessian_ = 0.0;
+
+ // quadratic part
+ for (size_t i = 0; i < A.N(); ++i) {
+ auto const end = std::end(A[i]);
+ for (auto it = std::begin(A[i]); it != end; ++it)
+ this->hessian_[i][it.index()] += *it;
+ }
+
+ // nonlinearity part
+ phi.addHessian(x, this->hessian_);
+
+ // compute gradient
+ // ----------------
+
+ // quadratic part
+ this->negativeGradient_ = derivative(this->f_)(x); // A*x - b
+
+ // nonlinearity part
+ phi.addGradient(x, this->negativeGradient_);
+
+ // -grad is needed for Newton step
+ this->negativeGradient_ *= -1;
+
+
+ // determine which components to truncate
+ // ----------------
+ this->truncationFlags_ = this->ignore_;
+
+ determineTruncation(x, this->f_.lowerObstacle(), this->f_.upperObstacle(), this->truncationFlags_, this->truncationTolerance_); // obstacle truncation
+ determineRegularityTruncation(x, this->truncationFlags_, this->truncationTolerance_); // truncation due to regularity deficit of nonlinearity
+
+ // truncate gradient and hessian
+ truncateVector(this->negativeGradient_, this->truncationFlags_);
+ truncateMatrix(this->hessian_, this->truncationFlags_, this->truncationFlags_);
+ }
+};
+
+#endif
diff --git a/src/spatial-solving/tnnmg/linesearchsolver.hh b/src/spatial-solving/tnnmg/linesearchsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..834740a22386c08f671547c4894fb619e126f475
--- /dev/null
+++ b/src/spatial-solving/tnnmg/linesearchsolver.hh
@@ -0,0 +1,35 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_LINESEARCHSOLVER_HH
+#define DUNE_TECTONIC_LINESEARCHSOLVER_HH
+
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+
+/**
+ * \brief A local solver for scalar quadratic obstacle problems with nonlinearity
+ * using bisection
+ */
+class LineSearchSolver
+{
+public:
+ template<class Vector, class Functional, class BitVector>
+ void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
+ x = 0;
+
+ Dune::Solvers::Interval<double> D;
+ D = f.subDifferential(0);
+ if (D[1] > 0) // NOTE: Numerical instability can actually get us here
+ return;
+
+ if (not ignore) {
+ int bisectionsteps = 0;
+ const Bisection globalBisection;
+ x = globalBisection.minimize(f, 1.0, 0.0, bisectionsteps);
+ // x = globalBisection.minimize(f, vNorm, 0.0, bisectionsteps) / vNorm; // used rescaled v in f?
+ }
+ }
+};
+
+#endif
+
diff --git a/src/spatial-solving/tnnmg/localbisectionsolver.hh b/src/spatial-solving/tnnmg/localbisectionsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..50618dd6c69a1a603ae49f68d5f43e476371353b
--- /dev/null
+++ b/src/spatial-solving/tnnmg/localbisectionsolver.hh
@@ -0,0 +1,66 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
+#define DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
+
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+/**
+ * \brief A local solver for quadratic obstacle problems with nonlinearity
+ * using bisection
+ */
+class LocalBisectionSolver
+{
+public:
+ template<class Vector, class Functional, class BitVector>
+ void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
+ if (ignore.all()) {
+ return;
+ }
+
+ double safety = 1e-14;
+ if (f.defectUpper() - f.defectLower() < safety) {
+ x = f.origin();
+ return;
+ }
+
+ int bisectionsteps = 0;
+ const Bisection globalBisection;
+
+ if (ignore.any()) {
+ auto&& restrictedf = f.getIgnoreFunctional(ignore);
+ const auto alpha = globalBisection.minimize(restrictedf, 0.0, 0.0, bisectionsteps);
+
+#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+ x = restrictedf.origin();
+ auto correction = restrictedf.direction();
+ correction *= alpha;
+ x += correction;
+#else
+ auto x = restrictedf.direction();
+ x *= alpha;
+#endif
+
+
+ } else { // ignore.none()
+ const auto alpha = globalBisection.minimize(f, 0.0, 0.0, bisectionsteps);
+
+#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+ x = f.origin();
+ auto correction = f.direction();
+ correction *= alpha;
+ x += correction;
+#else
+ auto x = f.direction();
+ x *= alpha;
+#endif
+
+ }
+ }
+};
+
+#endif
+
+
+
+
diff --git a/src/spatial-solving/tnnmg/multilevelpatchpreconditioner.hh b/src/spatial-solving/tnnmg/multilevelpatchpreconditioner.hh
new file mode 100644
index 0000000000000000000000000000000000000000..584cca16aa9c02877efcc778c2b7675ea3d6f297
--- /dev/null
+++ b/src/spatial-solving/tnnmg/multilevelpatchpreconditioner.hh
@@ -0,0 +1,277 @@
+#ifndef MULTILEVEL_PATCH_PRECONDITIONER_HH
+#define MULTILEVEL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+
+#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 <dune/faultnetworks/dgmgtransfer.hh>
+
+#include <dune/faultnetworks/utils/debugutils.hh>
+
+template <class BasisType, class LocalOperatorAssembler, class LocalInterfaceAssembler, 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;
+
+
+ typedef typename Dune::BitSetVector<1> BitVector;
+
+ const InterfaceNetwork<GridType>& interfaceNetwork_;
+ const BitVector& activeLevels_;
+ const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers_;
+ const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers_;
+ const Mode mode_;
+
+ int itCounter_;
+
+ std::shared_ptr<LevelGlobalPreconditionerType> coarseGlobalPreconditioner_;
+ std::vector<std::shared_ptr<LevelPatchPreconditionerType>> levelPatchPreconditioners_;
+ std::vector<VectorType> levelX_;
+ std::vector<VectorType> levelRhs_;
+
+ std::shared_ptr<LevelInterfaceNetwork<GridView>> allFaultLevelInterfaceNetwork_;
+ std::shared_ptr<BasisType> allFaultBasis_;
+
+ std::vector<std::shared_ptr<DGMGTransfer<BasisType>>> mgTransfer_;
+
+public:
+ MultilevelPatchPreconditioner(const InterfaceNetwork<GridType>& interfaceNetwork,
+ const BitVector& activeLevels,
+ const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers,
+ const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers,
+ const Mode mode = MultilevelPatchPreconditioner::Mode::additive) :
+ interfaceNetwork_(interfaceNetwork),
+ 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!");
+
+ assert(activeLevels.size() == localAssemblers_.size() && activeLevels.size() == localInterfaceAssemblers_.size());
+
+ // init level fault preconditioners and multigrid transfer
+ levelPatchPreconditioners_.resize(0);
+ mgTransfer_.resize(0);
+
+ int maxLevel = -1;
+ 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]);
+
+ maxLevel = i;
+ break;
+ }
+ }
+
+ typename LevelPatchPreconditionerType::Mode levelMode = LevelPatchPreconditionerType::Mode::additive;
+ if (mode_ != additive){
+ levelMode = LevelPatchPreconditionerType::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));
+ }
+
+ maxLevel = i;
+ }
+ }
+
+ levelX_.resize(levelPatchPreconditioners_.size()+1);
+ levelRhs_.resize(levelPatchPreconditioners_.size()+1);
+
+ allFaultBasis_ = std::make_shared<BasisType>(interfaceNetwork_.levelInterfaceNetwork(maxLevel));
+
+ // 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<DGMGTransfer<BasisType>>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
+
+ itCounter_ = 0;
+ }
+
+
+ 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);
+
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ mgTransfer_[i]->restrict(x, levelX_[i]);
+ mgTransfer_[i]->restrict(rhs, levelRhs_[i]);
+
+ levelPatchPreconditioners_[i]->setProblem(mat, levelX_[i], levelRhs_[i]);
+ }
+
+ size_t j = levelPatchPreconditioners_.size();
+ mgTransfer_[j]->restrict(x, levelX_[j]);
+ mgTransfer_[j]->restrict(rhs, levelRhs_[j]);
+
+ coarseGlobalPreconditioner_->setProblem(mat, levelX_[j], levelRhs_[j]);
+ }
+
+
+ virtual void iterate() {
+ if (mode_ == additive)
+ iterateAdd();
+ else
+ iterateMult();
+ }
+
+ void iterateAdd() {
+ *(this->x_) = 0;
+
+ VectorType x;
+
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ levelPatchPreconditioners_[i]->iterate();
+ const VectorType& it = levelPatchPreconditioners_[i]->getSol();
+
+ mgTransfer_[i]->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& it = coarseGlobalPreconditioner_->getSol();
+
+ mgTransfer_[levelPatchPreconditioners_.size()]->prolong(it, x);
+ *(this->x_) += x;
+
+
+ /*
+ levelPatchPreconditioners_[itCounter_]->iterate();
+ const VectorType& it = levelPatchPreconditioners_[itCounter_]->getSol();
+
+ 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;
+
+ */
+ }
+
+
+ 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();
+ const VectorType& it = levelPatchPreconditioners_[i]->getSol();
+
+ mgTransfer_[i]->prolong(it, x);
+
+ //writeToVTK(*allFaultBasis_, x, "/storage/mi/podlesny/data/faultnetworks/fineIterates/multilevel", "exactvertexdata_step_"+std::to_string(itCounter_));
+
+
+ *(this->x_) += x;
+ }
+
+ VectorType updatedCoarseRhs(*(this->rhs_));
+ this->mat_->mmv(*(this->x_), updatedCoarseRhs);
+
+ size_t j = levelPatchPreconditioners_.size();
+ 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_++;
+ }
+
+ void build() {
+ coarseGlobalPreconditioner_->build();
+
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ levelPatchPreconditioners_[i]->build();
+ }
+ }
+
+ std::shared_ptr<LevelPatchPreconditionerType> levelPatchPreconditioner(const int level) {
+ return levelPatchPreconditioners_[level];
+ }
+
+ const BasisType& basis() const {
+ return *allFaultBasis_;
+ }
+
+ size_t size() const {
+ return levelPatchPreconditioners_.size()+1;
+ }
+};
+
+#endif
+
diff --git a/src/spatial-solving/tnnmg/supportpatchfactory.hh b/src/spatial-solving/tnnmg/supportpatchfactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..20754c0d5c371b5d9a715778a2b9d7e87cace00d
--- /dev/null
+++ b/src/spatial-solving/tnnmg/supportpatchfactory.hh
@@ -0,0 +1,253 @@
+#ifndef SUPPORT_PATCH_FACTORY_HH
+#define SUPPORT_PATCH_FACTORY_HH
+
+#include<queue>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/referenceelementhelper.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dune/faultnetworks/hierarchicleveliterator.hh>
+
+template <class BasisType>
+class SupportPatchFactory
+{
+ protected:
+ typedef typename BasisType::GridView GV;
+ typedef typename BasisType::LocalFiniteElement LFE;
+
+ typedef typename GV::Grid GridType;
+ typedef typename GridType::ctype ctype;
+ static const int dim = GridType::dimension;
+ typedef typename GridType::template Codim<0>::Entity Element;
+ typedef typename GridType::template Codim<dim>::Entity Vertex;
+
+ typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
+ typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType, Dune::MCMGElementLayout > ElementMapper;
+
+ private:
+ const BasisType& coarseBasis_;
+ const BasisType& fineBasis_;
+
+ const int coarseLevel_;
+ const int fineLevel_;
+
+ const std::vector< std::vector <Element>>& vertexInElements_;
+ const int centerVertexIdx_;
+ const int patchDepth_;
+
+
+ const GridType& grid;
+ const GV& coarseGridView;
+ const GV& fineGridView;
+
+ std::vector<int> localToGlobal;
+ Dune::BitSetVector<1> boundaryDofs;
+ std::vector<Element> fineRegionElements;
+
+ ElementMapper coarseMapper;
+
+ public:
+ //setup
+ SupportPatchFactory(const BasisType& coarseBasis, const BasisType& fineBasis, const int coarseLevel, const int fineLevel, const std::vector< std::vector <Element>>& vertexInElements, const int centerVertexIdx, const int patchDepth) :
+ coarseBasis_(coarseBasis),
+ fineBasis_(fineBasis),
+ coarseLevel_(coarseLevel),
+ fineLevel_(fineLevel),
+ vertexInElements_(vertexInElements),
+ centerVertexIdx_(centerVertexIdx),
+ patchDepth_(patchDepth),
+ grid(coarseBasis_.getGridView().grid()),
+ coarseGridView(coarseBasis_.getGridView()),
+ fineGridView(fineBasis_.getGridView()),
+ coarseMapper(grid, coarseLevel_)
+ {
+ fineRegionElements.resize(0);
+
+ std::vector<Element> coarseElements;
+ Dune::BitSetVector<1> visited(coarseMapper.size());
+ Dune::BitSetVector<1> vertexVisited(vertexInElements_.size());
+ visited.unsetAll();
+ vertexVisited.unsetAll();
+
+ std::set<int> localDofs;
+ std::set<int> localBoundaryDofs;
+
+ addVertexSupport(centerVertexIdx_, coarseElements, visited, vertexVisited);
+
+ // construct coarse patch
+ for (int depth=1; depth<patchDepth_; depth++) {
+ std::vector<Element> coarseElementNeighbors;
+
+ for (size_t i=0; i<coarseElements.size(); ++i) {
+ const Element& elem = coarseElements[i];
+ const LFE& coarseLFE = coarseBasis_.getLocalFiniteElement(elem);
+
+ for (size_t j=0; j<coarseLFE.localBasis().size(); ++j) {
+ int dofIndex = coarseBasis_.indexInGridView(elem, j);
+
+ if (!vertexVisited[dofIndex][0]) {
+ addVertexSupport(dofIndex, coarseElementNeighbors, visited, vertexVisited);
+ }
+ }
+ }
+
+ coarseElements.insert(coarseElements.begin(), coarseElementNeighbors.begin(), coarseElementNeighbors.end());
+ }
+
+
+ // construct fine patch
+ for (size_t i=0; i<coarseElements.size(); ++i) {
+ const Element& elem = coarseElements[i];
+ addFinePatchElements(elem, visited, localDofs, localBoundaryDofs);
+ }
+
+
+
+ localToGlobal.resize(localDofs.size());
+ boundaryDofs.resize(localDofs.size());
+ boundaryDofs.unsetAll();
+
+ std::set<int>::iterator dofIt = localDofs.begin();
+ std::set<int>::iterator dofEndIt = localDofs.end();
+ size_t i=0;
+ for (; dofIt != dofEndIt; ++dofIt) {
+ int localDof = *dofIt;
+ localToGlobal[i] = localDof;
+
+ if (localBoundaryDofs.count(localDof)) {
+ boundaryDofs[i][0] = true;
+ }
+ i++;
+ }
+ }
+
+ std::vector<int>& getLocalToGlobal() {
+ return localToGlobal;
+ }
+
+ std::vector<Element>& getRegionElements() {
+ return fineRegionElements;
+ }
+
+ Dune::BitSetVector<1>& getBoundaryDofs() {
+ return boundaryDofs;
+ }
+
+
+private:
+
+
+ void print(const Dune::BitSetVector<1>& x, std::string message){
+ std::cout << message << std::endl;
+ for (size_t i=0; i<x.size(); i++) {
+ std::cout << x[i][0] << " ";
+ }
+ std::cout << std::endl << std::endl;
+ }
+
+ void addVertexSupport(const int vertexIdx, std::vector<Element>& coarseElements, Dune::BitSetVector<1>& visited, Dune::BitSetVector<1>& vertexVisited) {
+ const std::vector<Element>& vertexInElement = vertexInElements_[vertexIdx];
+
+ for (size_t i=0; i<vertexInElement.size(); i++) {
+ const Element& elem = vertexInElement[i];
+ const int elemIdx = coarseMapper.index(elem);
+
+ if (!visited[elemIdx][0]) {
+ coarseElements.push_back(elem);
+ visited[elemIdx][0] = true;
+ }
+ }
+ vertexVisited[vertexIdx][0] = true;
+ }
+
+ bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) {
+ return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
+ }
+
+
+ void addFinePatchElements(const Element& coarseElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
+ HierarchicLevelIteratorType endHierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
+ for (HierarchicLevelIteratorType hierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
+
+ const Element& fineElem = *hierIt;
+ fineRegionElements.push_back(fineElem);
+
+ addLocalDofs(coarseElem, fineElem, inCoarsePatch, localDofs, localBoundaryDofs);
+ }
+ }
+
+ void addLocalDofs(const Element& coarseElem, const Element& fineElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
+ const LFE& fineLFE = fineBasis_.getLocalFiniteElement(fineElem);
+
+ /*
+ const auto& fineGeometry = fineElem.geometry();
+ const auto& coarseGeometry = coarseElem.geometry();
+
+ const auto& ref = Dune::ReferenceElements<ctype, dim>::general(fineElem.type()); */
+
+ // insert local dofs
+ for (size_t i=0; i<fineLFE.localBasis().size(); ++i) {
+ int dofIndex = fineBasis_.index(fineElem, i);
+ localDofs.insert(dofIndex);
+ }
+
+
+ // search for parts of the global and inner boundary
+ typename GridType::LevelIntersectionIterator neighborIt = fineGridView.ibegin(fineElem);
+ typename GridType::LevelIntersectionIterator neighborItEnd = fineGridView.iend(fineElem);
+ for (; neighborIt!=neighborItEnd; ++neighborIt) {
+ const typename GridType::LevelIntersection& intersection = *neighborIt;
+
+ bool isInnerBoundary = false;
+ if (intersection.neighbor()) {
+ Element outsideFather = intersection.outside();
+
+ while (outsideFather.level()>coarseLevel_) {
+ outsideFather = outsideFather.father();
+ }
+
+ if (!inCoarsePatch[coarseMapper.index(outsideFather)][0]) {
+ isInnerBoundary = true;
+ }
+ }
+
+ if (intersection.boundary() or isInnerBoundary) {
+ typedef typename LFE::Traits::LocalCoefficientsType LocalCoefficients;
+ const LocalCoefficients* localCoefficients = &fineBasis_.getLocalFiniteElement(intersection.inside()).localCoefficients();
+
+ for (size_t i=0; i<localCoefficients->size(); i++) {
+ unsigned int entity = localCoefficients->localKey(i).subEntity();
+ unsigned int codim = localCoefficients->localKey(i).codim();
+
+ if (containsInsideSubentity(fineElem, intersection, entity, codim))
+ localBoundaryDofs.insert(fineBasis_.index(intersection.inside(), i));
+ }
+ }
+ }
+
+ /* add interior coarse level vertices to boundary dofs
+ for(int i=0; i<coarseGeometry.corners(); ++i) {
+ const auto& local = fineGeometry.local(coarseGeometry.corner(i));
+
+ if (!ref.checkInside(local))
+ continue;
+
+ const int localBasisSize = fineLFE.localBasis().size();
+ std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+ fineLFE.localBasis().evaluateFunction(local, localBasisRep);
+
+ for(int k=0; k<localBasisSize; ++k) {
+ if (localBasisRep[k]==1) {
+ int dofIndex = fineBasis_.index(fineElem, k);
+ localBoundaryDofs.insert(dofIndex);
+ break;
+ }
+ }
+ } */
+
+ }
+};
+#endif
diff --git a/src/spatial-solving/tnnmg/zerononlinearity.hh b/src/spatial-solving/tnnmg/zerononlinearity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4885d72a9c70ebe7923ae6434934d1e3092a72c7
--- /dev/null
+++ b/src/spatial-solving/tnnmg/zerononlinearity.hh
@@ -0,0 +1,76 @@
+#ifndef DUNE_TECTONIC_ZERO_NONLINEARITY_HH
+#define DUNE_TECTONIC_ZERO_NONLINEARITY_HH
+
+#include <dune/solvers/common/interval.hh>
+
+/** \file
+ * \brief A dummy nonlinearity class representing the zero function
+ */
+
+class ZeroNonlinearity
+{
+public:
+ ZeroNonlinearity()
+ {}
+
+ const auto& restriction(size_t i) const {
+ return *this;
+ }
+
+ /** \brief Returns zero */
+ template <class VectorType>
+ double operator()(const VectorType& v) const
+ {
+ return 0.0;
+ }
+
+ template <class VectorType>
+ void addGradient(const VectorType& v, VectorType& gradient) const {}
+
+ template <class VectorType, class MatrixType>
+ void addHessian(const VectorType& v, MatrixType& hessian) const {}
+
+ template <class IndexSet>
+ void addHessianIndices(IndexSet& indices) const {}
+
+ /** \brief Returns the interval \f$ [0,0]\f$ */
+ void subDiff(int i, double x, Dune::Solvers::Interval<double>& D, int j) const
+ {
+ D[0] = 0.0;
+ D[1] = 0.0;
+ }
+
+ template <class VectorType>
+ void directionalSubDiff(const VectorType& u, const VectorType& v, Dune::Solvers::Interval<double>& subdifferential) const
+ {
+ subdifferential[0] = 0.0;
+ subdifferential[1] = 0.0;
+ }
+
+ /** \brief Returns 0 */
+ double regularity(int i, double x, int j) const
+ {
+ return 0.0;
+ }
+
+ void domain(int i, Dune::Solvers::Interval<double>& dom, int j) const
+ {
+ dom[0] = -std::numeric_limits<double>::max();
+ dom[1] = std::numeric_limits<double>::max();
+ return;
+ }
+
+ template <class VectorType>
+ void directionalDomain(const VectorType& u, const VectorType& v, Dune::Solvers::Interval<double>& dom) const
+ {
+ dom[0] = -std::numeric_limits<double>::max();
+ dom[1] = std::numeric_limits<double>::max();
+ return;
+ }
+
+ template <class BitVector>
+ void setIgnore(const BitVector& ignore) {}
+};
+
+#endif
+
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index cde9e81d81cb6d445e4acdc8aec166be90c56d45..d9618460fc4dd793e82fd4c08663eeb5b2ca480d 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -17,9 +17,10 @@ void IterationRegister::reset() {
finalCount = FixedPointIterationCounter();
}
-template <class Factory, class Updaters, class ErrorNorm>
-AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
- Factory &factory, Dune::ParameterTree const &parset,
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::AdaptiveTimeStepper(
+ Dune::ParameterTree const &parset,
+ NBodyAssembler nBodyAssembler,
GlobalFrictionContainer& globalFriction, Updaters ¤t,
double relativeTime, double relativeTau,
ExternalForces& externalForces,
@@ -28,8 +29,8 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
: relativeTime_(relativeTime),
relativeTau_(relativeTau),
finalTime_(parset.get<double>("problem.finalTime")),
- factory_(factory),
parset_(parset),
+ nBodyAssembler_(nBodyAssembler),
globalFriction_(globalFriction),
current_(current),
R1_(),
@@ -37,13 +38,13 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
mustRefine_(mustRefine),
errorNorms_(errorNorms) {}
-template <class Factory, class Updaters, class ErrorNorm>
-bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+bool AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::reachedEnd() {
return relativeTime_ >= 1.0;
}
-template <class Factory, class Updaters, class ErrorNorm>
-IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::advance() {
/*
| 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
@@ -93,13 +94,13 @@ IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
return iterationRegister_;
}
-template <class Factory, class Updaters, class ErrorNorm>
-typename AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::UpdatersWithCount
-AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+typename AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::UpdatersWithCount
+AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(
Updaters const &oldUpdaters, double rTime, double rTau) {
UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
newUpdatersAndCount.count =
- MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
+ MyCoupledTimeStepper(finalTime_, parset_, nBodyAssembler_, globalFriction_,
newUpdatersAndCount.updaters, errorNorms_,
externalForces_)
.step(rTime, rTau);
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index 7b251fade0bb1709a209bacae861c696593a6050..6e70c1692a6380905a78f102d55f2f7dfd0efacc 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -15,7 +15,7 @@ struct IterationRegister {
FixedPointIterationCounter finalCount;
};
-template <class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
class AdaptiveTimeStepper {
struct UpdatersWithCount {
Updaters updaters;
@@ -26,13 +26,15 @@ class AdaptiveTimeStepper {
//using ConvexProblem = typename Factory::ConvexProblem;
//using Nonlinearity = typename Factory::Nonlinearity;
- using MyCoupledTimeStepper = CoupledTimeStepper<Factory, Updaters, ErrorNorm>;
+ using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>;
using GlobalFrictionContainer = typename MyCoupledTimeStepper::GlobalFrictionContainer;
using ExternalForces = typename MyCoupledTimeStepper::ExternalForces;
public:
- AdaptiveTimeStepper(Factory &factory, Dune::ParameterTree const &parset,
+ AdaptiveTimeStepper(
+ Dune::ParameterTree const &parset,
+ NBodyAssembler& nBodyAssembler,
GlobalFrictionContainer& globalFriction,
Updaters ¤t, double relativeTime,
double relativeTau,
@@ -51,8 +53,8 @@ class AdaptiveTimeStepper {
double rTau);
double finalTime_;
- Factory &factory_;
Dune::ParameterTree const &parset_;
+ NBodyAssembler& nBodyAssembler_;
GlobalFrictionContainer& globalFriction_;
Updaters ¤t_;
UpdatersWithCount R1_;
diff --git a/src/time-stepping/adaptivetimestepper_tmpl.cc b/src/time-stepping/adaptivetimestepper_tmpl.cc
index 4f11f38799fb9caa92e811342a2ac15da900a5e5..a390e8244d4a45ac8f887bd0ecdfd0598e531324 100644
--- a/src/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/src/time-stepping/adaptivetimestepper_tmpl.cc
@@ -11,7 +11,6 @@
#include <dune/tnnmg/problem-classes/convexproblem.hh>
#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
#include "../data-structures/levelcontactnetwork.hh"
#include "../spatial-solving/solverfactory.hh"
@@ -19,8 +18,9 @@
#include "state/stateupdater.hh"
#include "updaters.hh"
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
+/*
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryNodes;
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryFunctions;
using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
@@ -30,3 +30,4 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
template class AdaptiveTimeStepper<Factory, MyUpdaters, ErrorNorm>;
+*/
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index e06ce73233e419f737b8455c354b3fe3ad5e0c61..b948957208f231eb8f0bc8e9eabb4d9ee68b8b95 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -4,23 +4,24 @@
#include "coupledtimestepper.hh"
-template <class Factory, class Updaters, class ErrorNorm>
-CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
- double finalTime, Factory &factory, Dune::ParameterTree const &parset,
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::CoupledTimeStepper(
+ double finalTime, Dune::ParameterTree const &parset,
+ NBodyAssembler& nBodyAssembler,
GlobalFrictionContainer& globalFriction, Updaters updaters,
const std::vector<const ErrorNorm* >& errorNorms,
ExternalForces& externalForces)
: finalTime_(finalTime),
- factory_(factory),
parset_(parset),
+ nBodyAssembler_(nBodyAssembler),
globalFriction_(globalFriction),
updaters_(updaters),
externalForces_(externalForces),
errorNorms_(errorNorms) {}
-template <class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
FixedPointIterationCounter
-CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
+CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(double relativeTime,
double relativeTau) {
updaters_.state_->nextTimeStep();
updaters_.rate_->nextTimeStep();
@@ -39,8 +40,8 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
updaters_.state_->setup(tau);
updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
- FixedPointIterator<Factory, Updaters, ErrorNorm> fixedPointIterator(
- factory_, parset_, globalFriction_, errorNorms_);
+ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm> fixedPointIterator(
+ parset_, nBodyAssembler_, globalFriction_, errorNorms_);
auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
velocityRHS, velocityIterate);
return iterations;
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 81a53c6fd5bd286dee6951a4eb6a64fbae9fb923..9d129d2d38aad716c8d47be6c0f7a12a94eada76 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -8,18 +8,19 @@
#include "../spatial-solving/fixedpointiterator.hh"
-template <class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
class CoupledTimeStepper {
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
//using Nonlinearity = typename Factory::Nonlinearity;
public:
- using GlobalFrictionContainer = typename FixedPointIterator<Factory, Updaters, ErrorNorm>::GlobalFrictionContainer;
+ using GlobalFrictionContainer = typename FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::GlobalFrictionContainer;
using ExternalForces = std::vector<const std::function<void(double, Vector &)>*>;
public:
- CoupledTimeStepper(double finalTime, Factory &factory,
+ CoupledTimeStepper(double finalTime,
Dune::ParameterTree const &parset,
+ NBodyAssembler& nBodyAssembler,
GlobalFrictionContainer& globalFriction,
Updaters updaters, const std::vector<const ErrorNorm* >& errorNorms,
ExternalForces& externalForces);
@@ -28,8 +29,8 @@ class CoupledTimeStepper {
private:
double finalTime_;
- Factory &factory_;
Dune::ParameterTree const &parset_;
+ NBodyAssembler& nBodyAssembler_;
GlobalFrictionContainer& globalFriction_;
Updaters updaters_;
ExternalForces& externalForces_;
diff --git a/src/time-stepping/coupledtimestepper_tmpl.cc b/src/time-stepping/coupledtimestepper_tmpl.cc
index e16e9df9c948a02153eb818a75334410cd1fd88b..d5f6d1d14d9d597ad84f283b96d503371ebecf1a 100644
--- a/src/time-stepping/coupledtimestepper_tmpl.cc
+++ b/src/time-stepping/coupledtimestepper_tmpl.cc
@@ -11,7 +11,6 @@
#include <dune/tnnmg/problem-classes/convexproblem.hh>
#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
#include "../data-structures/levelcontactnetwork.hh"
#include "../spatial-solving/solverfactory.hh"
@@ -19,8 +18,9 @@
#include "state/stateupdater.hh"
#include "updaters.hh"
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
+/*
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryNodes;
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryFunctions;
using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
@@ -31,3 +31,4 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
template class CoupledTimeStepper<Factory, MyUpdaters, ErrorNorm>;
+*/
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index 2961e7e8885a92681c35fe0b8f48aa7d76d2cee4..4e61938bc64482db31a5f5abc54f94a6773be134 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -21,7 +21,13 @@ void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
std::vector<Vector>& rhs, std::vector<Vector>& iterate,
std::vector<Matrix>& AM) {
- for (size_t i=0; i<this->u.size(); i++) {
+ const size_t bodyCount = this->u.size();
+
+ rhs.resize(bodyCount);
+ iterate.resize(bodyCount);
+ AM.resize(bodyCount);
+
+ for (size_t i=0; i<bodyCount; i++) {
this->tau = _tau;
/* We start out with the formulation
@@ -77,7 +83,7 @@ void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
iterate = this->v_o;
- for (size_t i=0; i<iterate.size(); i++) {
+ for (size_t i=0; i<bodyCount; i++) {
auto& bodyIterate = iterate[i];
const auto& bodyDirichletNodes = this->dirichletNodes[i];
@@ -92,7 +98,7 @@ void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
for (size_t j=0; j<dim; ++j) {
if (bcDirichletNodes[k][j]) {
- iterate[k][j] = dirichletValue;
+ bodyIterate[k][j] = dirichletValue;
}
}
}
diff --git a/src/time-stepping/rate/rateupdater_tmpl.cc b/src/time-stepping/rate/rateupdater_tmpl.cc
index 273c261198849795ad1fd0d9f39b316b80f256b2..aa8a7168591684a2aecd6d3fc81b974514e645e0 100644
--- a/src/time-stepping/rate/rateupdater_tmpl.cc
+++ b/src/time-stepping/rate/rateupdater_tmpl.cc
@@ -7,8 +7,8 @@
#include "../../data-structures/levelcontactnetwork.hh"
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryNodes;
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryFunctions;
template class RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
template class Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
diff --git a/src/time-stepping/rate_tmpl.cc b/src/time-stepping/rate_tmpl.cc
index 0d1271ed0e11cf146900a568c435c6716ce521da..9627a8ce754aa52a9530ee61641b14a950b624e1 100644
--- a/src/time-stepping/rate_tmpl.cc
+++ b/src/time-stepping/rate_tmpl.cc
@@ -1,8 +1,8 @@
#include "../data-structures/levelcontactnetwork_tmpl.cc"
#include "../explicitvectors.hh"
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryFunctions;
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryNodes;
template
auto initRateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
diff --git a/src/time-stepping/state_tmpl.cc b/src/time-stepping/state_tmpl.cc
index 5046cf258356874450b8fb4ad5de5d58ed162795..42fa0b333c77d85d96396c4f86ad81d0f1eb4a9d 100644
--- a/src/time-stepping/state_tmpl.cc
+++ b/src/time-stepping/state_tmpl.cc
@@ -7,7 +7,7 @@
#include "../data-structures/levelcontactnetwork.hh"
-using BoundaryPatchNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryPatchNodes;
+using BoundaryPatchNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryPatchNodes;
template
auto initStateUpdater<ScalarVector, Vector, BoundaryPatchNodes>(
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index 4671ff68de50d8c51cf674d89f968e9097548f56..150aad1cf2cf0b99175678de146e437dc7abd4f7 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -29,25 +29,29 @@ using namespace std;
template <int dim, typename ctype=double>
void print(const Dune::FieldVector<ctype, dim>& x, std::string message){
- std::cout << message << std::endl;
+ if (message != "")
+ std::cout << message << std::endl;
+
for (size_t i=0; i<x.size(); i++) {
- std::cout << x[i] << " ";
+ std::cout << x[i] << ";" << std::endl;;
}
- std::cout << std::endl;
+
}
template <int dim, typename ctype=double>
void print(const Dune::BlockVector<Dune::FieldVector<ctype, dim>>& x, std::string message){
std::cout << message << " " << "size " << x.size() << std::endl;
for (size_t i=0; i<x.size(); i++) {
- print(x[i], "block " + std::to_string(i) + ":");
+ print(x[i], "");
}
std::cout << std::endl << std::endl;
}
template <int dim, typename ctype=double>
- void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message){
- std::cout << message << std::endl;
+ void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, bool endOfLine = true){
+ if (message != "")
+ std::cout << message << std::endl;
+
for (size_t i=0; i<mat.N(); i++) {
for (size_t j=0; j<mat.M(); j++) {
if (mat.exists(i,j))
@@ -55,20 +59,44 @@ using namespace std;
else
std::cout << "0 ";
}
- std::cout << std::endl;
+
+ if (endOfLine)
+ std::cout << ";"<< std::endl;
}
- std::cout << std::endl;
}
- template <int dim, typename ctype=double>
- void print(const Dune::BCRSMatrix<Dune::FieldMatrix<ctype, dim, dim>>& mat, std::string message){
- std::cout << message << " " << "size (" << mat.N() << ", " << mat.M() << ")" <<std::endl;
- for (size_t i=0; i<mat.N(); i++) {
- for (size_t j=0; j<mat.M(); j++) {
- if (mat.exists(i,j))
- print(mat[i][j], "block (" + std::to_string(i) + ", " + std::to_string(j) + "):");
- }
- std::cout << std::endl;
+ template <int dim, typename ctype=double>
+ void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, const size_t line, bool endOfLine = true){
+ if (message != "")
+ std::cout << message << std::endl;
+
+ assert(line<mat.N());
+ for (size_t j=0; j<mat.M(); j++) {
+ if (mat.exists(line,j))
+ std::cout << mat[line][j] << " ";
+ else
+ std::cout << "0 ";
+ }
+
+ if (endOfLine)
+ std::cout << ";"<< std::endl;
+ }
+
+ template <int dim, typename ctype=double>
+ void print(const Dune::BCRSMatrix<Dune::FieldMatrix<ctype, dim, dim>>& mat, std::string message){
+ std::cout << message << " " << "size (" << mat.N() << ", " << mat.M() << ")" <<std::endl;
+
+ Dune::FieldMatrix<ctype, dim, dim> zeroEntry = 0;
+
+ for (size_t i=0; i<mat.N(); i++) {
+ for (size_t d=0; d<dim; d++) {
+ for (size_t j=0; j<mat.M(); j++) {
+ if (mat.exists(i,j))
+ print(mat[i][j], "", d, j==mat.M()-1);
+ else
+ print(zeroEntry, "", d, j==mat.M()-1);
+ }
+ }
}
std::cout << std::endl;
}