diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index 078ec7b3b55951970ecf625678973d6cde3044fc..ecd7b7b9b18906627d3e72f956a4bcd19809448f 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -23,10 +23,28 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
typedef typename MyConvexProblemType::LocalMatrixType LocalMatrixType;
int static const block_size = MyConvexProblemType::block_size;
+ int static const coarse_block_size = block_size;
/** \brief Solves one local system using a modified gradient method */
class IterateObject;
+ struct Linearization {
+ static const int block_size = coarse_block_size;
+
+ typedef typename MyBlockProblem<
+ MyConvexProblemTypeTEMPLATE>::LocalMatrixType LocalMatrixType;
+ typedef Dune::BCRSMatrix<typename Linearization::LocalMatrixType>
+ MatrixType;
+ typedef Dune::BlockVector<
+ Dune::FieldVector<double, Linearization::block_size>> VectorType;
+ typedef Dune::BitSetVector<Linearization::block_size> BitVectorType;
+ typename Linearization::MatrixType A;
+ typename Linearization::VectorType b;
+ typename Linearization::BitVectorType ignore;
+
+ Dune::BitSetVector<Linearization::block_size> truncation;
+ };
+
MyBlockProblem(Dune::ParameterTree const &parset,
MyConvexProblemType &problem)
: parset(parset), problem(problem) {
@@ -38,6 +56,110 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
0); // safety: acceptance factor for inexact minimization
}
+ std::string virtual getOutput(bool header = false) const {
+ // TODO: implement (not urgent)
+ return std::string();
+ }
+
+ void projectCoarseCorrection(VectorType const &u,
+ typename Linearization::VectorType const &v,
+ VectorType &projected_v,
+ Linearization const &linearization) const {
+ // TODO: implement (not urgent)
+ projected_v = v;
+ };
+
+ double computeDampingParameter(VectorType const &u,
+ VectorType const &projected_v) const {
+ // TODO: implement
+ return 1.0;
+ };
+
+ void assembleTruncate(VectorType const &u, Linearization &linearization,
+ Dune::BitSetVector<block_size> const &ignore) const {
+ // we can just copy the ignore information
+ linearization.ignore = ignore;
+
+ // determine truncation pattern
+ linearization.truncation.resize(u.size());
+ linearization.truncation.unsetAll();
+ Interval<double> ab;
+ for (int i = 0; i < u.size(); ++i) {
+ for (int j = 0; j < block_size; ++j) {
+ // problem.phi.smoothDomain(i, u[i][j], regularity_tol, ab, j);
+
+ // // if (phi is not regular at u) or (u is not in smooth domain) or
+ // (component should be ignored)
+ // // truncate this component
+ // if ((problem.phi.regularity(i, u[i][j], j) > regularity_tol)
+ // or (u[i][j] < ab[0])
+ // or (u[i][j] > ab[1])
+ // or ignore[i][j])
+ // linearization.truncation[i][j] = true;
+ }
+ }
+
+ // construct sparsity pattern for linearization
+ Dune::MatrixIndexSet indices(problem.A.N(), problem.A.M());
+ indices.import(problem.A);
+ // problem.phi.addHessianIndices(indices);
+
+ // construct matrix from pattern and initialize it
+ indices.exportIdx(linearization.A);
+ linearization.A = 0.0;
+
+ // compute quadratic part of hessian (linearization.A += problem.A)
+ for (int i = 0; i < problem.A.N(); ++i) {
+ typename MatrixType::row_type::ConstIterator it = problem.A[i].begin();
+ typename MatrixType::row_type::ConstIterator end = problem.A[i].end();
+ for (; it != end; ++it)
+ StaticMatrix::axpy(linearization.A[i][it.index()], 1.0, *it);
+ }
+
+ // compute nonlinearity part of hessian
+ problem.phi.addHessian(u, linearization.A);
+
+ // compute quadratic part of gradient
+ linearization.b.resize(u.size());
+ linearization.b = 0.0;
+ problem.A.template umv<VectorType, VectorType>(u, linearization.b);
+ linearization.b -= problem.f;
+
+ // compute nonlinearity part of gradient
+ problem.phi.addGradient(u, linearization.b);
+
+ // -grad is needed for Newton step
+ linearization.b *= -1.0;
+
+ // apply truncation to system
+ typename Linearization::MatrixType::row_type::Iterator it;
+ typename Linearization::MatrixType::row_type::Iterator end;
+ for (int row = 0; row < linearization.A.N(); ++row) {
+ it = linearization.A[row].begin();
+ end = linearization.A[row].end();
+ for (; it != end; ++it) {
+ int col = it.index();
+ for (int i = 0; i < (*it).N(); ++i) {
+ typename Linearization::MatrixType::block_type::row_type::Iterator
+ blockIt = (*it)[i].begin();
+ typename Linearization::MatrixType::block_type::row_type::Iterator
+ blockEnd = (*it)[i].end();
+ for (; blockIt != blockEnd; ++blockIt)
+ if (linearization.truncation[row][i] or linearization
+ .truncation[col][blockIt.index()])
+ (*blockIt) = 0.0;
+ }
+ }
+
+ for (int j = 0; j < block_size; ++j)
+ if (linearization.truncation[row][j])
+ linearization.b[row][j] = 0.0;
+ }
+
+ // for(int j=0; j<block_size; ++j)
+ // outStream << std::setw(9) << linearization.truncation.countmasked(j);
+ }
+
/** \brief Constructs and returns an iterate object */
IterateObject getIterateObject() {
return IterateObject(parset, bisection, problem);
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index efd7d49f29e3863357c5a26bb68e064d3924aa23..478d22586f566a952168e645626706e505ced5f9 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -49,6 +49,11 @@
#include <dune/tectonic/myconvexproblem.hh>
#include <dune/tectonic/myblockproblem.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
int const dim = 2;
template <class GridView>
@@ -242,17 +247,20 @@ int main(int argc, char *argv[]) {
u1 = 0;
VectorType u2 = u1;
VectorType u3 = u1;
+ VectorType u4 = u1;
VectorType u1_diff_new(grid.size(grid.maxLevel(), dim));
u1_diff_new = 0;
VectorType u2_diff_new = u1_diff_new;
VectorType u3_diff_new = u1_diff_new;
+ VectorType u4_diff_new = u1_diff_new;
CellVectorType vonMisesStress;
VectorType b1;
VectorType b2;
VectorType b3;
+ VectorType b4;
auto nodalIntegrals =
Dune::make_shared<Dune::BlockVector<Dune::FieldVector<double, 1>>>();
@@ -272,17 +280,19 @@ int main(int argc, char *argv[]) {
leafView, p1Basis, neumannNodes, b1, run);
b2 = b1;
b3 = b1;
+ b4 = b1;
// b -= linear update
stiffnessMatrix.mmv(u1, b1);
stiffnessMatrix.mmv(u2, b2);
stiffnessMatrix.mmv(u3, b3);
+ stiffnessMatrix.mmv(u4, b4);
if (parset.get<bool>("useNonlinearGS")) {
MyConvexProblemType myConvexProblem(
stiffnessMatrix, *myGlobalNonlinearity, b1, u1_diff_new);
- MyBlockProblemType myBlockProblem(parset, myConvexProblem);
- nonlinearGSStep.setProblem(u1_diff_new, myBlockProblem);
+ auto myBlockProblem = new MyBlockProblemType(parset, myConvexProblem);
+ nonlinearGSStep.setProblem(u1_diff_new, *myBlockProblem);
LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations,
solver_tolerance, &energyNorm, verbosity);
@@ -291,6 +301,72 @@ int main(int argc, char *argv[]) {
u1 += u1_diff_new;
+ if (parset.get<bool>("useTNNMG")) {
+ MyConvexProblemType myConvexProblem(
+ stiffnessMatrix, *myGlobalNonlinearity, b4, u4_diff_new);
+ auto myBlockProblem = new MyBlockProblemType(parset, myConvexProblem);
+
+ // {{{ Linear Solver;
+ TruncatedBlockGSStep<OperatorType, VectorType> *linearBaseSolverStep =
+ new TruncatedBlockGSStep<OperatorType, VectorType>;
+
+ EnergyNorm<OperatorType, VectorType> *baseEnergyNorm =
+ new EnergyNorm<OperatorType, VectorType>(*linearBaseSolverStep);
+ LoopSolver<VectorType> *linearBaseSolver = new LoopSolver<VectorType>(
+ linearBaseSolverStep, solver_maxIterations, solver_tolerance,
+ baseEnergyNorm, Solver::QUIET);
+ // }}}
+
+ // {{{ Smoothers
+ TruncatedBlockGSStep<OperatorType, VectorType> *linearPresmoother =
+ new TruncatedBlockGSStep<OperatorType, VectorType>;
+ TruncatedBlockGSStep<OperatorType, VectorType> *linearPostsmoother =
+ new TruncatedBlockGSStep<OperatorType, VectorType>;
+ // }}}
+
+ MultigridStep<OperatorType, VectorType> *linearIterationStep =
+ new MultigridStep<OperatorType, VectorType>;
+ linearIterationStep->setNumberOfLevels(levels);
+ linearIterationStep->setMGType(1, 3, 3); // 1/3/3 (mu/nu1/nu2)
+ linearIterationStep->basesolver_ = linearBaseSolver;
+ linearIterationStep->setSmoother(linearPresmoother, linearPostsmoother);
+
+ // {{{ Transfer operators
+ linearIterationStep->mgTransfer_.resize(levels - 1);
+ for (size_t i = 0; i < linearIterationStep->mgTransfer_.size(); i++) {
+ CompressedMultigridTransfer<VectorType> *newTransferOp =
+ new CompressedMultigridTransfer<VectorType>;
+ newTransferOp->setup(grid, i, i + 1);
+ linearIterationStep->mgTransfer_[i] = newTransferOp;
+ }
+ // }}}
+
+ // // Now the actual nonlinear solver
+ typedef MyBlockProblemType TNNMGProblemType;
+ typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType;
+ typedef TruncatedNonsmoothNewtonMultigrid<
+ TNNMGProblemType, NonlinearSmootherType> TNNMGStepType;
+
+ TNNMGProblemType *tnnmgProblem = myBlockProblem;
+ NonlinearSmootherType *nonlinearSmoother = new NonlinearSmootherType;
+
+ auto multigridStep =
+ new TNNMGStepType(*linearIterationStep, *nonlinearSmoother);
+ multigridStep->setProblem(u4_diff_new, *tnnmgProblem);
+ multigridStep->setSmoothingSteps(3, 1, 3); // 3/1/3 (pre/linear/post)
+ multigridStep->ignoreNodes_ = &ignoreNodes;
+
+ auto energyNorm =
+ new EnergyNorm<OperatorType, VectorType>(stiffnessMatrix);
+
+ LoopSolver<VectorType> overallSolver(
+ multigridStep, solver_maxIterations, solver_tolerance, energyNorm,
+ verbosity);
+ overallSolver.solve();
+ }
+
+ u4 += u4_diff_new;
+
{ // Compute von Mises stress and write everything to a file
auto *displacement =
new BasisGridFunction<P1Basis, VectorType>(p1Basis, u1);
@@ -351,12 +427,19 @@ int main(int argc, char *argv[]) {
std::cout << "sup |u1 - u3| = " << diff3.infinity_norm() << ", "
<< "|u1 - u3| = " << diff3.two_norm() << std::endl;
+ VectorType diff4 = u4;
+ diff4 -= u1;
+ std::cout << "sup |u1 - u4| = " << diff4.infinity_norm() << ", "
+ << "|u1 - u4| = " << diff4.two_norm() << std::endl;
+
// Print displacement on frictional boundary
for (size_t i = 0; i < frictionalNodes.size(); ++i)
if (frictionalNodes[i][0])
std::cout << boost::format("u1[%02d] = %+3e, %|40t|u2[%02d] = %+3e "
- "%|80t|u3[%02d] = %+3e") %
- i % u1[i] % i % u2[i] % i % u3[i] << std::endl;
+ "%|80t|u3[%02d] = %+3e %|120t|u4[%02d] = "
+ "%+3e") %
+ i % u1[i] % i % u2[i] % i % u3[i] % i %
+ u4[i] << std::endl;
}
catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index caef4f15a928c88678af871d6e3dc1989a520ee8..d49f051d8a2c6f3e46e59aef31a109ea2285973d 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -6,6 +6,7 @@ verbose = false
useNonlinearGS = true
useGS = false
useTruncatedGS = false
+useTNNMG = true
[grid]
refinements = 5