Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_SOLVENONLINEARPROBLEM_HH
#define DUNE_TECTONIC_SPATIAL_SOLVING_SOLVENONLINEARPROBLEM_HH
#include <dune/common/exceptions.hh>
#include <dune/matrix-vector/axpy.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/contact/assemblers/nbodyassembler.hh>
#include <dune/contact/common/dualbasisadapter.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include "../data-structures/enums.hh"
#include "../data-structures/enumparser.hh"
#include "../utils/tobool.hh"
#include "../utils/debugutils.hh"
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include "../spatial-solving/contact/nbodycontacttransfer.hh"
#include "solverfactory.hh"
#include "solverfactory.cc"
#include <dune/tectonic/utils/reductionfactors.hh>
#include "fixedpointiterator.hh"
template <class Functional, class BitVector>
class NonlinearSolver {
protected:
using Matrix = typename Functional::Matrix;
using Vector = typename Functional::Vector;
using Factory = SolverFactory<Functional, BitVector>;
using Norm = EnergyNorm<Matrix, Vector>;
using SolverType = Dune::Solvers::LoopSolver<Vector>;
const static int dims = Vector::block_type::dimension;
public:
template <class LinearSolver>
NonlinearSolver(const Dune::ParameterTree& tnnmgParset, std::shared_ptr<LinearSolver> linearSolver, std::shared_ptr<Functional> functional, const BitVector& dirichletNodes) :
functional_(functional),
norm_(functional_->quadraticPart()) {
// set up TNMMG step
solverFactory_ = std::make_shared<Factory>(tnnmgParset, functional_, dirichletNodes);
solverFactory_->build(linearSolver);
}
auto solve(const Dune::ParameterTree& solverParset, Vector& x) {
auto tnnmgStep = solverFactory_->step();
SolverType solver(*tnnmgStep.get(), solverParset.get<size_t>("maximumIterations"),
solverParset.get<double>("tolerance"), norm_,
solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
const auto& lower = functional_->lowerObstacle();
const auto& upper = functional_->upperObstacle();
// project in onto admissible set
for (size_t i=0; i<x.size(); i++) {
for (size_t j=0; j<dims; j++) {
if (x[i][j] < lower[i][j]) {
x[i][j] = lower[i][j];
}
if (x[i][j] > upper[i][j]) {
x[i][j] = upper[i][j];
}
}
}
solverFactory_->setProblem(x);
solver.preprocess();
solver.solve();
return solver.getResult();
}
private:
std::shared_ptr<Functional> functional_;
std::shared_ptr<Factory> solverFactory_;
Norm norm_;
};
#endif