From 7bb4bca8ea1f083231b736ef0973138d185bb61a Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Wed, 24 Jul 2019 17:32:04 +0200
Subject: [PATCH] add solverfactorytest
---
src/CMakeLists.txt | 26 +
src/data-structures/program_state.hh | 2 +
src/factories/stackedblocksfactory.cc | 7 +-
src/multi-body-problem-2D.cfg | 2 +-
src/multi-body-problem-data/bc.hh | 3 -
src/obstaclesolver.hh | 48 ++
src/solverfactorytest.cc | 756 ++++++++++++++++++
src/spatial-solving/fixedpointiterator.cc | 3 +
src/spatial-solving/solverfactory.cc | 9 +-
src/spatial-solving/solverfactory.hh | 13 +-
src/spatial-solving/tnnmg/functional.hh | 17 +-
src/spatial-solving/tnnmg/linearcorrection.hh | 18 +-
src/spatial-solving/tnnmg/tnnmgstep.hh | 51 +-
13 files changed, 903 insertions(+), 52 deletions(-)
create mode 100644 src/obstaclesolver.hh
create mode 100644 src/solverfactorytest.cc
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index e7cf62d6..b9748ffc 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -64,17 +64,42 @@ set(UGW_SOURCE_FILES
one-body-problem-data/mygrid.cc
)
+set(SFT_SOURCE_FILES
+ assemblers.cc
+ data-structures/body.cc
+ data-structures/levelcontactnetwork.cc
+ data-structures/contactnetwork.cc
+ data-structures/enumparser.cc
+ factories/stackedblocksfactory.cc
+ io/vtk.cc
+ multi-body-problem-data/grid/cuboidgeometry.cc
+ multi-body-problem-data/grid/mygrids.cc
+ multi-body-problem-data/grid/simplexmanager.cc
+ #spatial-solving/solverfactory.cc
+ spatial-solving/fixedpointiterator.cc
+ #spatial-solving/solverfactory_old.cc
+ time-stepping/adaptivetimestepper.cc
+ time-stepping/coupledtimestepper.cc
+ time-stepping/rate.cc
+ time-stepping/rate/rateupdater.cc
+ time-stepping/state.cc
+ solverfactorytest.cc
+)
+
foreach(_dim 2 3)
set(_sw_target one-body-problem-${_dim}D)
set(_msw_target multi-body-problem-${_dim}D)
set(_ugw_target uniform-grid-writer-${_dim}D)
+ set(_sft_target solverfactorytest-${_dim}D)
add_executable(${_sw_target} ${SW_SOURCE_FILES})
add_executable(${_msw_target} ${MSW_SOURCE_FILES})
add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
+ add_executable(${_sft_target} ${SFT_SOURCE_FILES})
add_dune_ug_flags(${_sw_target})
add_dune_ug_flags(${_msw_target})
add_dune_ug_flags(${_ugw_target})
+ add_dune_ug_flags(${_sft_target})
add_dune_hdf5_flags(${_sw_target})
add_dune_hdf5_flags(${_msw_target})
@@ -82,4 +107,5 @@ foreach(_dim 2 3)
set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+ set_property(TARGET ${_sft_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
endforeach()
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 382c82b1..afe9075f 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -243,6 +243,8 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
solver.preprocess();
solver.solve();
+ std::cout << "ProgramState: Energy of TNNMG solution: " << J(tnnmgStep->getSol()) << std::endl;
+
nBodyAssembler.postprocess(tnnmgStep->getSol(), _x);
Vector res = totalRhs;
diff --git a/src/factories/stackedblocksfactory.cc b/src/factories/stackedblocksfactory.cc
index 9f7fc32a..648518f3 100644
--- a/src/factories/stackedblocksfactory.cc
+++ b/src/factories/stackedblocksfactory.cc
@@ -73,10 +73,13 @@ void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
for (size_t i=1; i<this->bodyCount_-1; i++) {
origins[i] = cuboidGeometries_[i-1]->upperLeft();
+ auto height_i = height/3.0;
+
GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
- GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];
+ GlobalCoords upperWeakOrigin_i = {upperWeakOrigin[0], height_i};
+ upperWeakOrigin_i += origins[i];
- cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height);
+ cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height_i);
cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index 9112827c..a4a77808 100644
--- a/src/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem-2D.cfg
@@ -1,6 +1,6 @@
# -*- mode:conf -*-
[boundary.friction]
-smallestDiameter = 0.1 # 2e-3 [m]
+smallestDiameter = 0.25 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/multi-body-problem-data/bc.hh b/src/multi-body-problem-data/bc.hh
index 1b574cd6..c7c919ae 100644
--- a/src/multi-body-problem-data/bc.hh
+++ b/src/multi-body-problem-data/bc.hh
@@ -7,7 +7,6 @@ class VelocityDirichletCondition
: public Dune::VirtualFunction<double, double> {
void evaluate(double const &relativeTime, double &y) const {
// Assumed to vanish at time zero
- /*
double const finalVelocity = -5e-5;
std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
@@ -20,8 +19,6 @@ class VelocityDirichletCondition
y = (relativeTime <= 0.1)
? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
: finalVelocity;
- */
- y = relativeTime;
}
};
diff --git a/src/obstaclesolver.hh b/src/obstaclesolver.hh
new file mode 100644
index 00000000..49b4d7cb
--- /dev/null
+++ b/src/obstaclesolver.hh
@@ -0,0 +1,48 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_OBSTACLESOLVER_HH
+#define DUNE_TECTONIC_OBSTACLESOLVER_HH
+
+#include "spatial-solving/tnnmg/localproblem.hh"
+
+
+/**
+ * \brief A local solver for quadratic obstacle problems
+ *
+ * \todo Add concept check for the function interface
+ */
+class ObstacleSolver
+{
+public:
+ template<class Vector, class Functional, class BitVector>
+ constexpr void operator()(Vector& x, const Functional& f, const BitVector& ignore) const
+ {
+ auto& A = f.quadraticPart();
+ auto& r = f.linearPart();
+
+ LocalProblem<std::decay_t<decltype(A)> , std::decay_t<decltype(r)>> localProblem(A, r, ignore);
+ Vector newR;
+ localProblem.getLocalRhs(x, newR);
+
+ auto directSolver = Dune::Solvers::UMFPackSolver<std::decay_t<decltype(A)>, std::decay_t<decltype(r)>>();
+
+ directSolver->setProblem(localProblem.getMat(), x, newR);
+ directSolver->preprocess();
+ directSolver->solve();
+
+ const auto& lower = f.lowerObstacle();
+ const auto& upper = f.upperObstacle();
+
+ for (size_t i=0; i<x.size(); i++) {
+ if (ignore[i][0])
+ continue;
+
+ if (x[i] < lower[i])
+ x[i] = lower[i];
+ if (x[i] > upper[i])
+ x[i] = upper[i];
+ }
+ }
+};
+
+#endif
diff --git a/src/solverfactorytest.cc b/src/solverfactorytest.cc
new file mode 100644
index 00000000..19b68965
--- /dev/null
+++ b/src/solverfactorytest.cc
@@ -0,0 +1,756 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#define MY_DIM 2
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <exception>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/stdstreams.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/fufem/formatstring.hh>
+
+#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
+#include <dune/tectonic/geocoordinate.hh>
+
+#include "assemblers.hh"
+#include "gridselector.hh"
+#include "explicitgrid.hh"
+#include "explicitvectors.hh"
+
+#include "data-structures/enumparser.hh"
+#include "data-structures/enums.hh"
+#include "data-structures/contactnetwork.hh"
+#include "data-structures/matrices.hh"
+#include "data-structures/program_state.hh"
+
+#include "io/vtk.hh"
+
+#include "spatial-solving/tnnmg/tnnmgstep.hh"
+#include "spatial-solving/tnnmg/linesearchsolver.hh"
+#include "spatial-solving/preconditioners/nbodycontacttransfer.hh"
+#include "obstaclesolver.hh"
+
+#include "factories/stackedblocksfactory.hh"
+
+#include "time-stepping/rate.hh"
+#include "time-stepping/state.hh"
+#include "time-stepping/updaters.hh"
+
+#include "utils/tobool.hh"
+#include "utils/debugutils.hh"
+
+
+const int dim = MY_DIM;
+
+const int timeSteps = 1;
+
+Dune::ParameterTree parset;
+size_t bodyCount;
+
+std::vector<BodyState<Vector, ScalarVector>* > bodyStates;
+std::vector<Vector> u;
+std::vector<Vector> v;
+std::vector<Vector> a;
+std::vector<ScalarVector> alpha;
+std::vector<ScalarVector> weightedNormalStress;
+double relativeTime;
+double relativeTau;
+size_t timeStep = 0;
+
+const std::string path = "";
+const std::string outputFile = "solverfactorytest.log";
+
+template <class ContactNetwork>
+void solveProblem(const ContactNetwork& contactNetwork,
+ const Matrix& mat, const Vector& rhs, Vector& x,
+ const BitVector& ignore, const Vector& lower, const Vector& upper) {
+
+ using Solver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ using Norm = EnergyNorm<Matrix, Vector>;
+
+ //using LocalSolver = LocalBisectionSolver;
+ //using Linearization = Linearization<Functional, BitVector>;
+
+ // set up reference solver
+ Vector refX = x;
+ using ContactFunctional = Dune::TNNMG::BoxConstrainedQuadraticFunctional<Matrix&, Vector&, Vector&, Vector&, double>;
+ auto I = ContactFunctional(mat, rhs, lower, upper);
+
+ using LocalSolver = Dune::TNNMG::ScalarObstacleSolver;
+ auto localSolver = Dune::TNNMG::gaussSeidelLocalSolver(LocalSolver());
+
+ using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<ContactFunctional, decltype(localSolver), BitVector>;
+ auto nonlinearSmoother = std::make_shared<NonlinearSmoother>(I, refX, localSolver);
+
+ using Linearization = Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<ContactFunctional, BitVector>;
+ using DefectProjection = Dune::TNNMG::ObstacleDefectProjection;
+ using Step = Dune::TNNMG::TNNMGStep<ContactFunctional, BitVector, Linearization, DefectProjection, LocalSolver>;
+
+ // set multigrid solver
+ auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+
+ using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+ using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+ TransferOperators transfer(contactNetwork.nLevels()-1);
+ for (size_t i=0; i<transfer.size(); i++) {
+ transfer[i] = std::make_shared<TransferOperator>();
+ transfer[i]->setup(contactNetwork, i, i+1);
+ }
+
+ // Remove any recompute filed so that initially the full transferoperator is assembled
+ for (size_t i=0; i<transfer.size(); i++)
+ std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+ auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+ linearMultigridStep->setMGType(1, 3, 3);
+ linearMultigridStep->setSmoother(&smoother);
+ linearMultigridStep->setTransferOperators(transfer);
+
+ int mu = parset.get<int>("solver.tnnmg.main.multi"); // #multigrid steps in Newton step
+ auto step = Step(I, refX, nonlinearSmoother, linearMultigridStep, mu, DefectProjection(), LocalSolver());
+
+ // compute reference solution with generic functional and solver
+ auto norm = Norm(mat);
+ auto refSolver = Solver(step, parset.get<size_t>("u0.solver.maximumIterations"),
+ parset.get<double>("u0.solver.tolerance"), norm, Solver::FULL);
+
+ step.setIgnore(ignore);
+ step.setPreSmoothingSteps(parset.get<int>("solver.tnnmg.main.pre"));
+
+ refSolver.addCriterion(
+ [&](){
+ return Dune::formatString(" % 12.5e", I(x));
+ },
+ " energy ");
+
+ double initialEnergy = I(x);
+ refSolver.addCriterion(
+ [&](){
+ static double oldEnergy=initialEnergy;
+ double currentEnergy = I(x);
+ double decrease = currentEnergy - oldEnergy;
+ oldEnergy = currentEnergy;
+ return Dune::formatString(" % 12.5e", decrease);
+ },
+ " decrease ");
+
+ refSolver.addCriterion(
+ [&](){
+ return Dune::formatString(" % 12.5e", step.lastDampingFactor());
+ },
+ " damping ");
+
+
+ refSolver.addCriterion(
+ [&](){
+ return Dune::formatString(" % 12d", step.linearization().truncated().count());
+ },
+ " truncated ");
+
+
+ std::vector<double> correctionNorms;
+ refSolver.addCriterion(Dune::Solvers::correctionNormCriterion(step, norm, 1e-10, correctionNorms));
+
+ refSolver.preprocess();
+ refSolver.solve();
+ std::cout << correctionNorms.size() << std::endl;
+
+
+ // set up solver factory solver
+
+ // set up functional
+ using MyFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+ MyFunctional J(mat, rhs, ZeroNonlinearity(), lower, upper);
+
+ // set up TNMMG solver
+ // dummy solver, uses direct solver for linear correction no matter what is set here
+ Norm mgNorm(*linearMultigridStep);
+ auto mgSolver = std::make_shared<Solver>(linearMultigridStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), mgNorm, Solver::QUIET);
+
+ using Factory = SolverFactory<MyFunctional, BitVector>;
+ Factory factory(parset.sub("solver.tnnmg"), J, *mgSolver, ignore);
+
+ /* std::vector<BitVector> bodyDirichletNodes;
+ nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
+ for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
+ print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
+ }*/
+
+ /* print(bilinearForm, "matrix: ");
+ print(totalX, "totalX: ");
+ print(totalRhs, "totalRhs: ");*/
+
+ auto tnnmgStep = factory.step();
+ factory.setProblem(x);
+
+ LoopSolver<Vector> solver(
+ tnnmgStep.get(), parset.get<size_t>("u0.solver.maximumIterations"),
+ parset.get<double>("u0.solver.tolerance"), &norm,
+ Solver::FULL, true, &refX); // absolute error
+
+ solver.preprocess();
+ solver.solve();
+
+ auto diff = x;
+ diff -= refX;
+ std::cout << "Solver error in energy norm: " << norm(diff) << std::endl;
+
+ x = refX;
+}
+
+
+template <class ContactNetwork>
+void setupInitialConditions(const ContactNetwork& contactNetwork) {
+ using Matrix = typename ContactNetwork::Matrix;
+ const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+
+ std::cout << std::endl << "-- setupInitialConditions --" << std::endl;
+ std::cout << "----------------------------" << std::endl;
+
+ // Solving a linear problem with a multigrid solver
+ auto const solveLinearProblem = [&](
+ const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+ const std::vector<Vector>& _rhs, std::vector<Vector>& _x) {
+
+ std::vector<const Matrix*> matrices_ptr(_matrices.size());
+ for (size_t i=0; i<matrices_ptr.size(); i++) {
+ matrices_ptr[i] = _matrices[i].get();
+ }
+
+ // assemble full global contact problem
+ Matrix bilinearForm;
+
+ nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+ Vector totalRhs, oldTotalRhs;
+ nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
+ oldTotalRhs = totalRhs;
+
+ 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<dim; ++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");*/
+
+ solveProblem(contactNetwork, bilinearForm, totalRhs, totalX, _dirichletNodes, lower, upper);
+
+ nBodyAssembler.postprocess(totalX, _x);
+ };
+
+ timeStep = parset.get<size_t>("initialTime.timeStep");
+ relativeTime = parset.get<double>("initialTime.relativeTime");
+ relativeTau = parset.get<double>("initialTime.relativeTau");
+
+ bodyStates.resize(bodyCount);
+ u.resize(bodyCount);
+ v.resize(bodyCount);
+ a.resize(bodyCount);
+ alpha.resize(bodyCount);
+ weightedNormalStress.resize(bodyCount);
+
+ for (size_t i=0; i<bodyCount; i++) {
+ size_t leafVertexCount = contactNetwork.body(i)->nVertices();
+
+ u[i].resize(leafVertexCount),
+ v[i].resize(leafVertexCount),
+ a[i].resize(leafVertexCount),
+ alpha[i].resize(leafVertexCount),
+ weightedNormalStress[i].resize(leafVertexCount),
+
+ bodyStates[i] = new BodyState<Vector, ScalarVector>(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
+ }
+
+ 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());
+ ell0[i] = 0.0;
+
+ contactNetwork.body(i)->externalForce()(relativeTime, ell0[i]);
+ }
+
+ // Initial displacement: Start from a situation of minimal stress,
+ // which is automatically attained in the case [v = 0 = a].
+ // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
+ BitVector dirichletNodes;
+ contactNetwork.totalNodes("dirichlet", dirichletNodes);
+ /*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];
+ }
+
+ dirichletNodes[i] = val;
+ for (size_t d=0; d<dims; d++) {
+ dirichletNodes[i][d] = val;
+ }
+ }*/
+
+ std::cout << "solving linear problem for u..." << std::endl;
+
+ solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, u);
+
+ //print(u, "initial u:");
+
+ // Initial acceleration: Computed in agreement with Ma = ell0 - Au
+ // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
+ std::vector<Vector> accelerationRHS = ell0;
+ for (size_t i=0; i<bodyCount; i++) {
+ // Initial state
+ alpha[i] = parset.get<double>("boundary.friction.initialAlpha");
+
+ // Initial normal stress
+
+ const auto& body = contactNetwork.body(i);
+ std::vector<std::shared_ptr<typename ContactNetwork::LeafBody::BoundaryCondition>> frictionBoundaryConditions;
+ body->boundaryConditions("friction", frictionBoundaryConditions);
+ for (size_t j=0; j<frictionBoundaryConditions.size(); j++) {
+ ScalarVector frictionBoundaryStress(weightedNormalStress[i].size());
+
+ body->assembler()->assembleWeightedNormalStress(
+ *frictionBoundaryConditions[j]->boundaryPatch(), frictionBoundaryStress, body->data()->getYoungModulus(),
+ body->data()->getPoissonRatio(), u[i]);
+
+ weightedNormalStress[i] += frictionBoundaryStress;
+ }
+
+ Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
+ }
+
+ std::cout << "solving linear problem for a..." << std::endl;
+
+ BitVector noNodes(dirichletNodes.size(), false);
+ solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a);
+
+ //print(a, "initial a:");
+}
+
+template <class ContactNetwork>
+void relativeVelocities(const ContactNetwork& contactNetwork, const Vector& v, std::vector<Vector>& v_rel) {
+
+ const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+ const size_t nBodies = nBodyAssembler.nGrids();
+ // const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+
+ std::vector<const Dune::BitSetVector<1>*> bodywiseNonmortarBoundaries;
+ contactNetwork.frictionNodes(bodywiseNonmortarBoundaries);
+
+ size_t globalIdx = 0;
+ v_rel.resize(nBodies);
+ for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+ const auto& nonmortarBoundary = *bodywiseNonmortarBoundaries[bodyIdx];
+ auto& v_rel_ = v_rel[bodyIdx];
+
+ v_rel_.resize(nonmortarBoundary.size());
+
+ for (size_t i=0; i<v_rel_.size(); i++) {
+ if (toBool(nonmortarBoundary[i])) {
+ v_rel_[i] = v[globalIdx];
+ }
+ globalIdx++;
+ }
+ }
+
+}
+
+template <class Updaters, class ContactNetwork>
+void run(Updaters& updaters, ContactNetwork& contactNetwork,
+ const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+ std::vector<Vector>& velocityIterates) {
+
+ const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+ const auto nBodies = nBodyAssembler.nGrids();
+
+ std::vector<const Matrix*> matrices_ptr(nBodies);
+ for (int i=0; i<nBodies; i++) {
+ matrices_ptr[i] = &velocityMatrices[i];
+ }
+
+ // assemble full global contact problem
+ Matrix bilinearForm;
+ nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+ Vector totalRhs;
+ nBodyAssembler.assembleRightHandSide(velocityRHSs, totalRhs);
+
+ // 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<dim; ++d) {
+ lowerj[d] = totalObstaclesj[d][0];
+ upperj[d] = totalObstaclesj[d][1];
+ }
+ }
+
+ // compute velocity obstacles
+ Vector vLower, vUpper;
+ std::vector<Vector> u0, v0;
+ updaters.rate_->extractOldVelocity(v0);
+ updaters.rate_->extractOldDisplacement(u0);
+
+ Vector totalu0, totalv0;
+ nBodyAssembler.concatenateVectors(u0, totalu0);
+ nBodyAssembler.concatenateVectors(v0, totalv0);
+
+ updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
+ updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);
+
+ const auto& errorNorms = contactNetwork.stateEnergyNorms();
+
+ auto& globalFriction = contactNetwork.globalFriction();
+
+ BitVector totalDirichletNodes;
+ contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+ size_t fixedPointMaxIterations = parset.get<size_t>("v.fpi.maximumIterations");
+ double fixedPointTolerance = parset.get<double>("v.fpi.tolerance");
+ double lambda = parset.get<double>("v.fpi.lambda");
+
+ size_t fixedPointIteration;
+ size_t multigridIterations = 0;
+ std::vector<ScalarVector> alpha(nBodies);
+ updaters.state_->extractAlpha(alpha);
+ for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations;
+ ++fixedPointIteration) {
+
+ // contribution from nonlinearity
+ globalFriction.updateAlpha(alpha);
+
+ Vector totalVelocityIterate;
+ nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+
+ solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, vLower, vUpper);
+
+ nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
+
+ std::vector<Vector> v_m;
+ updaters.rate_->extractOldVelocity(v_m);
+
+ for (size_t i=0; i<v_m.size(); i++) {
+ v_m[i] *= 1.0 - lambda;
+ Dune::MatrixVector::addProduct(v_m[i], lambda, velocityIterates[i]);
+ }
+
+ // extract relative velocities in mortar basis
+ std::vector<Vector> v_rel;
+ relativeVelocities(contactNetwork, totalVelocityIterate, v_rel);
+
+ std::cout << "- State problem set" << std::endl;
+
+ // solve a state problem
+ updaters.state_->solve(v_rel);
+
+ std::cout << "-- Solved" << std::endl;
+
+ std::vector<ScalarVector> newAlpha(nBodies);
+ updaters.state_->extractAlpha(newAlpha);
+
+ bool breakCriterion = true;
+ for (int i=0; i<nBodies; i++) {
+ if (alpha[i].size()==0 || newAlpha[i].size()==0)
+ continue;
+
+ print(alpha[i], "alpha i:");
+ print(newAlpha[i], "new alpha i:");
+ if (errorNorms[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance) {
+ breakCriterion = false;
+ std::cout << "fixedPoint error: " << errorNorms[i]->diff(alpha[i], newAlpha[i]) << std::endl;
+ break;
+ }
+ }
+
+ if (lambda < 1e-12 or breakCriterion) {
+ std::cout << "-FixedPointIteration finished! " << (lambda < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
+ fixedPointIteration++;
+ break;
+ }
+ alpha = newAlpha;
+ }
+
+ std::cout << "-FixedPointIteration finished! " << std::endl;
+
+ if (fixedPointIteration == fixedPointMaxIterations)
+ DUNE_THROW(Dune::Exception, "FPI failed to converge");
+
+ updaters.rate_->postProcess(velocityIterates);
+}
+
+
+template <class Updaters, class ContactNetwork>
+void step(Updaters& updaters, ContactNetwork& contactNetwork) {
+ updaters.state_->nextTimeStep();
+ updaters.rate_->nextTimeStep();
+
+ auto const newRelativeTime = relativeTime + relativeTau;
+ typename ContactNetwork::ExternalForces externalForces;
+ contactNetwork.externalForces(externalForces);
+ std::vector<Vector> ell(externalForces.size());
+ for (size_t i=0; i<externalForces.size(); i++) {
+ (*externalForces[i])(newRelativeTime, ell[i]);
+ }
+
+ std::vector<Matrix> velocityMatrix;
+ std::vector<Vector> velocityRHS;
+ std::vector<Vector> velocityIterate;
+
+ double finalTime = parset.get<double>("problem.finalTime");
+ auto const tau = relativeTau * finalTime;
+ updaters.state_->setup(tau);
+ updaters.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
+
+ run(updaters, contactNetwork, velocityMatrix, velocityRHS, velocityIterate);
+}
+
+template <class Updaters, class ContactNetwork>
+void advanceStep(Updaters& current, ContactNetwork& contactNetwork) {
+ step(current, contactNetwork);
+ relativeTime += relativeTau;
+}
+
+void getParameters(int argc, char *argv[]) {
+ Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem.cfg", parset);
+ Dune::ParameterTreeParser::readINITree(
+ Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem-%dD.cfg", dim), parset);
+ Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+}
+
+int main(int argc, char *argv[]) { try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ std::ofstream out(path + outputFile);
+ std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+ std::cout.rdbuf(out.rdbuf()); //redirect std::cout to outputFile
+
+ std::cout << "-------------------------" << std::endl;
+ std::cout << "-- SolverFactory Test: --" << std::endl;
+ std::cout << "-------------------------" << std::endl << std::endl;
+
+ getParameters(argc, argv);
+
+ // ----------------------
+ // set up contact network
+ // ----------------------
+ StackedBlocksFactory<Grid, Vector> stackedBlocksFactory(parset);
+ using ContactNetwork = typename StackedBlocksFactory<Grid, Vector>::ContactNetwork;
+ stackedBlocksFactory.build();
+
+ ContactNetwork& contactNetwork = stackedBlocksFactory.contactNetwork();
+ bodyCount = contactNetwork.nBodies();
+
+ for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+ const auto& level = *contactNetwork.level(i);
+
+ for (size_t j=0; j<level.nBodies(); j++) {
+ writeToVTK(level.body(j)->gridView(), "debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+ }
+ }
+
+ for (size_t i=0; i<bodyCount; i++) {
+ writeToVTK(contactNetwork.body(i)->gridView(), "debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
+ }
+
+ // ----------------------------
+ // assemble contactNetwork
+ // ----------------------------
+ contactNetwork.assemble();
+
+ //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+ // -----------------
+ // init input/output
+ // -----------------
+
+ setupInitialConditions(contactNetwork);
+
+ auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+ for (size_t i=0; i<bodyCount; i++) {
+ contactNetwork.body(i)->setDeformation(u[i]);
+ }
+ nBodyAssembler.assembleTransferOperator();
+ nBodyAssembler.assembleObstacle();
+
+ // ------------------------
+ // assemble global friction
+ // ------------------------
+ contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), weightedNormalStress);
+
+ auto& globalFriction = contactNetwork.globalFriction();
+ globalFriction.updateAlpha(alpha);
+
+ using Assembler = MyAssembler<DefLeafGridView, dim>;
+ using field_type = Matrix::field_type;
+ using MyVertexBasis = typename Assembler::VertexBasis;
+ using MyCellBasis = typename Assembler::CellBasis;
+ std::vector<Vector> vertexCoordinates(bodyCount);
+ std::vector<const MyVertexBasis* > vertexBases(bodyCount);
+ std::vector<const MyCellBasis* > cellBases(bodyCount);
+
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = contactNetwork.body(i);
+ vertexBases[i] = &(body->assembler()->vertexBasis);
+ cellBases[i] = &(body->assembler()->cellBasis);
+
+ auto& vertexCoords = vertexCoordinates[i];
+ vertexCoords.resize(body->nVertices());
+
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->gridView(), Dune::mcmgVertexLayout());
+ for (auto &&v : vertices(body->gridView()))
+ vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
+ }
+
+ const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
+
+ auto const report = [&](bool initial = false) {
+ if (parset.get<bool>("io.printProgress"))
+ std::cout << "timeStep = " << std::setw(6) << timeStep
+ << ", time = " << std::setw(12) << relativeTime
+ << ", tau = " << std::setw(12) << relativeTau
+ << std::endl;
+
+ if (parset.get<bool>("io.vtk.write")) {
+ std::vector<ScalarVector> stress(bodyCount);
+
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = contactNetwork.body(i);
+ body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
+ body->data()->getPoissonRatio(),
+ u[i], stress[i]);
+
+ }
+
+ vtkWriter.write(timeStep, u, v, alpha, stress);
+ }
+ };
+ report(true);
+
+ // -------------------
+ // Set up TNNMG solver
+ // -------------------
+
+ BitVector totalDirichletNodes;
+ contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+ //print(totalDirichletNodes, "totalDirichletNodes:");
+
+ //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
+ //using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
+ //using NonlinearFactory = SolverFactory<Functional, BitVector>;
+
+ using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+ using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+ using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+ StateUpdater<ScalarVector, Vector>>;
+
+ BoundaryFunctions velocityDirichletFunctions;
+ contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+ BoundaryNodes dirichletNodes;
+ contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+ /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+ for (size_t j=0; j<dirichletNodes[i].size(); j++) {
+ print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
+ }
+ }*/
+
+
+ /*for (size_t i=0; i<frictionNodes.size(); i++) {
+ print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
+ }*/
+
+ Updaters current(
+ initRateUpdater(
+ parset.get<Config::scheme>("timeSteps.scheme"),
+ velocityDirichletFunctions,
+ dirichletNodes,
+ contactNetwork.matrices(),
+ u,
+ v,
+ a),
+ initStateUpdater<ScalarVector, Vector>(
+ parset.get<Config::stateModel>("boundary.friction.stateModel"),
+ alpha,
+ nBodyAssembler.getContactCouplings(),
+ contactNetwork.couplings())
+ );
+
+
+
+ //const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
+
+ for (timeStep=1; timeStep<=timeSteps; timeStep++) {
+
+ advanceStep(current, contactNetwork);
+
+ relativeTime += relativeTau;
+
+ current.rate_->extractDisplacement(u);
+ current.rate_->extractVelocity(v);
+ current.rate_->extractAcceleration(a);
+ current.state_->extractAlpha(alpha);
+
+ contactNetwork.setDeformation(u);
+
+ report();
+ }
+
+ bool passed = true;
+
+ std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+ std::cout.rdbuf(coutbuf); //reset to standard output again
+ return passed ? 0 : 1;
+
+} catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+} catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+} // end try
+} // end main
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index d84f5f09..43df107c 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -210,6 +210,9 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
std::cout << "-- Solved" << std::endl;
const auto& tnnmgSol = step->getSol();
+
+ std::cout << "FixPointIterator: Energy of TNNMG solution: " << J(tnnmgSol) << std::endl;
+
nBodyAssembler.postprocess(tnnmgSol, velocityIterates);
//nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index cfe36305..5f02147e 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -19,7 +19,14 @@ SolverFactory<Functional, BitVector>::SolverFactory(
const BitVector& ignoreNodes) :
J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))) {
- nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
+ auto localSolver = Dune::TNNMG::gaussSeidelLocalSolver(LocalSolver());
+ //nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, localSolver);
+
+ //using MyNonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, decltype(localSolver), BitVector>;
+ nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, localSolver);
+
+ //nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
+
auto linearSolver_ptr = Dune::Solvers::wrap_own_share<std::decay_t<LinearSolver>>(std::forward<LinearSolver>(linearSolver));
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 5d4a5aa2..5fe7e62f 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -10,11 +10,12 @@
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/cgstep.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/tnnmg/localsolvers/scalarobstaclesolver.hh>
-#include "tnnmg/tnnmgstep.hh"
+//#include "tnnmg/tnnmgstep.hh"
#include "tnnmg/linearization.hh"
#include "tnnmg/linesearchsolver.hh"
#include "tnnmg/localbisectionsolver.hh"
@@ -27,13 +28,13 @@ class SolverFactory {
using Vector = typename Functional::Vector;
using BitVector = BitVectorType;
- using LocalSolver = LocalBisectionSolver;
- using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalSolver>;
+ using LocalSolver = Dune::TNNMG::ScalarObstacleSolver;//LocalBisectionSolver;
+ using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, Dune::TNNMG::GaussSeidelLocalSolver<LocalSolver>, BitVector>;
using Linearization = Linearization<Functional, BitVector>;
using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
- using Step = typename Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
-
+ //using Step = TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+ using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
template <class LinearSolver>
SolverFactory(const Dune::ParameterTree&,
diff --git a/src/spatial-solving/tnnmg/functional.hh b/src/spatial-solving/tnnmg/functional.hh
index 4bad3e3f..098c7208 100644
--- a/src/spatial-solving/tnnmg/functional.hh
+++ b/src/spatial-solving/tnnmg/functional.hh
@@ -16,6 +16,12 @@
#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+template<class M, class V, class N, class R>
+class ShiftedFunctional;
+
+template <class M, class V, class N, class L, class U, class R>
+class Functional;
+
/** \brief Coordinate restriction of box constrained quadratic functional with nonlinearity;
* mainly used for presmoothing in TNNMG algorithm
*
@@ -346,12 +352,13 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
auto&& phii = f.phi().restriction(i);
- auto v = ri;
+ /*auto v = ri;
double const vnorm = v.two_norm();
if (vnorm > 1.0)
- v /= vnorm;
+ 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));
+ //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));
+ return Functional<LocalMatrix&, LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, Range>(*Aii_p, std::move(ri), std::move(phii), std::move(dli), std::move(dui));
}
@@ -389,8 +396,8 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
MM&& matrix,
VV&& linearPart,
NN&& phi,
- LL& lower,
- UU& upper) :
+ 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))
{}
diff --git a/src/spatial-solving/tnnmg/linearcorrection.hh b/src/spatial-solving/tnnmg/linearcorrection.hh
index c26f13df..188d580f 100644
--- a/src/spatial-solving/tnnmg/linearcorrection.hh
+++ b/src/spatial-solving/tnnmg/linearcorrection.hh
@@ -1,5 +1,5 @@
-#ifndef DUNE_TNNMG_ITERATIONSTEPS_LINEARCORRECTION_HH
-#define DUNE_TNNMG_ITERATIONSTEPS_LINEARCORRECTION_HH 1
+#ifndef DUNE_TECTONIC_ITERATIONSTEPS_LINEARCORRECTION_HH
+#define DUNE_TECTONIC_ITERATIONSTEPS_LINEARCORRECTION_HH
#include <functional>
#include <memory>
@@ -38,7 +38,7 @@ emptyIgnore(const Vector& v)
template<typename Matrix, typename Vector>
LinearCorrection<Matrix, Vector>
-makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector> > linearSolver)
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector> > linearSolver)
{
return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
@@ -56,7 +56,7 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector
template<typename Matrix, typename Vector>
LinearCorrection<Matrix, Vector>
-makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > iterativeSolver)
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > iterativeSolver)
{
return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
@@ -76,16 +76,16 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > i
directSolver->preprocess();
directSolver->solve();
- x = directX;
- /*using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
+ //x = directX;
+ using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
auto linearIterationStep = dynamic_cast<LinearIterationStep*>(&iterativeSolver->getIterationStep());
if (not linearIterationStep)
DUNE_THROW(Dune::Exception, "iterative solver must use a linear iteration step");
auto empty = emptyIgnore(x);
-
linearIterationStep->setIgnore(empty);
+ //linearIterationStep->setIgnore(ignore);
linearIterationStep->setProblem(A, x, b);
iterativeSolver->preprocess();
iterativeSolver->solve();
@@ -94,13 +94,13 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > i
const auto& norm = iterativeSolver->getErrorNorm();
auto error = norm.diff(linearIterationStep->getSol(), directX);
- std::cout << "Linear solver iterations: " << iterativeSolver->getResult().iterations << " Error: " << error << std::endl;*/
+ std::cout << "Linear solver iterations: " << iterativeSolver->getResult().iterations << " Error: " << error << std::endl;
};
}
template<typename Matrix, typename Vector>
LinearCorrection<Matrix, Vector>
-makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearIterationStep<Matrix, Vector> > linearIterationStep, int nIterationSteps = 1)
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::LinearIterationStep<Matrix, Vector> > linearIterationStep, int nIterationSteps = 1)
{
return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
linearIterationStep->setIgnore(ignore);
diff --git a/src/spatial-solving/tnnmg/tnnmgstep.hh b/src/spatial-solving/tnnmg/tnnmgstep.hh
index 44979dcf..1d439664 100644
--- a/src/spatial-solving/tnnmg/tnnmgstep.hh
+++ b/src/spatial-solving/tnnmg/tnnmgstep.hh
@@ -1,5 +1,5 @@
-#ifndef DUNE_TNNMG_ITERATIONSTEPS_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
-#define DUNE_TNNMG_ITERATIONSTEPS_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
+#ifndef DUNE_TECTONIC_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
+#define DUNE_TECTONIC_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
#include <string>
#include <sstream>
@@ -24,8 +24,6 @@
#include <dune/tectonic/../../src/utils/debugutils.hh>
-namespace Dune {
-namespace TNNMG {
/**
* \brief One iteration of the TNNMG method
@@ -49,8 +47,8 @@ class TNNMGStep :
using BitVector = typename Base::BitVector;
using ConstrainedBitVector = typename Linearization::ConstrainedBitVector;
using Functional = F;
- using IterativeSolver = Solvers::IterativeSolver< ConstrainedVector, Solvers::DefaultBitVector_t<ConstrainedVector> >;
- using LinearSolver = Solvers::LinearSolver< ConstrainedMatrix, ConstrainedVector >;
+ using IterativeSolver = Dune::Solvers::IterativeSolver< ConstrainedVector, Dune::Solvers::DefaultBitVector_t<ConstrainedVector> >;
+ using LinearSolver = Dune::Solvers::LinearSolver< ConstrainedMatrix, ConstrainedVector >;
/** \brief Constructor with an iterative solver object for the linear correction
* \param iterativeSolver This is a callback used to solve the constrained linearized system
@@ -67,7 +65,7 @@ class TNNMGStep :
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
- linearCorrection_(makeLinearCorrection<ConstrainedMatrix>(iterativeSolver)),
+ linearCorrection_(buildLinearCorrection<ConstrainedMatrix>(iterativeSolver)),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -87,7 +85,7 @@ class TNNMGStep :
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
- linearCorrection_(makeLinearCorrection(linearSolver)),
+ linearCorrection_(buildLinearCorrection(linearSolver)),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -100,15 +98,15 @@ class TNNMGStep :
*/
TNNMGStep(const Functional& f,
Vector& x,
- std::shared_ptr<Solvers::IterationStep<Vector,BitVector> > nonlinearSmoother,
- std::shared_ptr<Solvers::LinearIterationStep<ConstrainedMatrix,ConstrainedVector> > linearIterationStep,
+ std::shared_ptr<Dune::Solvers::IterationStep<Vector,BitVector> > nonlinearSmoother,
+ std::shared_ptr<Dune::Solvers::LinearIterationStep<ConstrainedMatrix,ConstrainedVector> > linearIterationStep,
unsigned int noOfLinearIterationSteps,
const DefectProjection& projection,
const LineSearchSolver& lineSolver)
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
- linearCorrection_(makeLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
+ linearCorrection_(buildLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -145,24 +143,24 @@ class TNNMGStep :
const auto& ignore = (*this->ignoreNodes_);
auto& x = *getIterate();
- std::cout << "TNNMGStep::iterate " << std::endl;
+ //std::cout << "TNNMGStep::iterate " << std::endl;
//print(f.quadraticPart(), "f.quadraticPart():");
//print(f.linearPart(), "f.linearPart():");
- std::cout << "-- energy before smoothing: " << f(x) << std::endl;
+ //std::cout << "-- energy before smoothing: " << f(x) << std::endl;
- print(x, "TNNMG iterate after smoothing:");
+ //print(x, "TNNMG iterate after smoothing:");
// Nonlinear presmoothing
for (std::size_t i=0; i<preSmoothingSteps_; ++i)
nonlinearSmoother_->iterate();
//std::cout << "- nonlinear presmoothing: success" << std::endl;
- print(x, "TNNMG iterate after smoothing:");
+ //print(x, "TNNMG iterate after smoothing:");
- std::cout << "-- energy after presmoothing: " << f(x) << std::endl;
+ //std::cout << "-- energy after presmoothing: " << f(x) << std::endl;
// Compute constraint/truncated linearization
if (not linearization_)
@@ -173,20 +171,20 @@ class TNNMGStep :
auto&& A = linearization_->hessian();
auto&& r = linearization_->negativeGradient();
- print(A, "TNNMG Linear problem matrix:");
+ /*print(A, "TNNMG Linear problem matrix:");
print(r, "TNNMG Linear problem rhs:");
print(ignore, "ignore:");
- print(linearization_->truncated(), "truncation:");
+ print(linearization_->truncated(), "truncation:");*/
// Compute inexact solution of the linearized problem
- Solvers::resizeInitializeZero(correction_, x);
- Solvers::resizeInitializeZero(constrainedCorrection_, r);
+ Dune::Solvers::resizeInitializeZero(correction_, x);
+ Dune::Solvers::resizeInitializeZero(constrainedCorrection_, r);
//print(constrainedCorrection_compare, "direct solver solution:");
- std::cout << "- computing linear correction..." << std::endl;
+ //std::cout << "- computing linear correction..." << std::endl;
//linearCorrection_(A, constrainedCorrection_, r, ignore);
@@ -207,9 +205,13 @@ class TNNMGStep :
//std::cout << "- extended correction: success" << std::endl;
//print(correction_, "correction:");
+ //std::cout << "-- energy before projection: " << f(x) << std::endl;
+
// Project onto admissible set
projection_(f, x, correction_);
+ //std::cout << "-- energy after projection: " << f(x) << std::endl;
+
//std::cout << "- projection onto admissible set: success" << std::endl;
//print(correction_, "correction:");
@@ -231,6 +233,8 @@ class TNNMGStep :
dampingFactor_ = 0;
lineSolver_(dampingFactor_, fv, false);
+ //std::cout << "damping factor: " << dampingFactor_ << std::endl;
+
// std::cout << "- computing damping factor: success" << std::endl;
if (std::isnan(dampingFactor_))
dampingFactor_ = 0;
@@ -240,7 +244,7 @@ class TNNMGStep :
x += correction_;
- std::cout << "-- energy after correction: " << energy(x) << std::endl;
+ //std::cout << "-- energy after correction: " << energy(x) << std::endl;
}
/**
@@ -279,7 +283,4 @@ class TNNMGStep :
};
-} // end namespace TNNMG
-} // end namespace Dune
-
#endif
--
GitLab