From 5feb394403ccb34570818e657bf3ca537054d52e Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Wed, 15 May 2019 21:00:20 +0200
Subject: [PATCH] .
---
dune/tectonic/frictionpotential.hh | 3 +
src/data-structures/levelcontactnetwork.cc | 5 +
src/data-structures/program_state.hh | 10 +-
src/factories/stackedblocksfactory.cc | 31 ++--
src/io/vtk.cc | 3 +
src/multi-body-problem-2D.cfg | 4 +-
src/multi-body-problem.cc | 44 +++++-
src/multi-body-problem.cfg | 7 +-
src/spatial-solving/fixedpointiterator.cc | 4 +-
src/spatial-solving/solverfactory.cc | 7 +-
src/spatial-solving/solverfactory.hh | 9 +-
src/spatial-solving/tnnmg/linesearchsolver.hh | 8 +-
src/spatial-solving/tnnmg/localproblem.hh | 141 ++++++++++++++++++
src/spatial-solving/tnnmg/tnnmgstep.hh | 87 +++++++++--
src/time-stepping/adaptivetimestepper.cc | 10 +-
src/time-stepping/rate/backward_euler.cc | 2 +
src/time-stepping/rate/newmark.cc | 5 +
.../state/ageinglawstateupdater.cc | 1 +
src/time-stepping/state/stateupdater.hh | 2 +-
src/utils/debugutils.hh | 46 +++++-
src/utils/tobool.hh | 2 +-
21 files changed, 378 insertions(+), 53 deletions(-)
create mode 100644 src/spatial-solving/tnnmg/localproblem.hh
diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index 326a5249..fcfcc5be 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -106,6 +106,9 @@ class RegularisedRateState : public FrictionPotential {
class ZeroFunction : public FrictionPotential {
public:
+ template <typename... Args>
+ ZeroFunction(Args... args) {}
+
double evaluate(double) const override { return 0; }
double coefficientOfFriction(double s) const override { return 0; }
diff --git a/src/data-structures/levelcontactnetwork.cc b/src/data-structures/levelcontactnetwork.cc
index aa336513..454c2f63 100644
--- a/src/data-structures/levelcontactnetwork.cc
+++ b/src/data-structures/levelcontactnetwork.cc
@@ -8,6 +8,7 @@
#include <dune/fufem/functions/constantfunction.hh>
#include <dune/tectonic/globalratestatefriction.hh>
+#include <dune/tectonic/frictionpotential.hh>
#include "levelcontactnetwork.hh"
@@ -105,6 +106,10 @@ void LevelContactNetwork<GridType, dimension>::assembleFriction(
}
}
+ /* globalFriction_ = std::make_shared<GlobalRateStateFriction<
+ Matrix, Vector, ZeroFunction, DeformedGridType>>(
+ nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);*/
+
switch (frictionModel) {
case Config::Truncated:
globalFriction_ = std::make_shared<GlobalRateStateFriction<
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 298f40f7..1eceb99a 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -215,9 +215,9 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
std::cout << "Energy norm: " << norm.operator ()(diff) << std::endl;*/
};
- timeStep = 0;
- relativeTime = 0.0;
- relativeTau = 1e-6;
+ timeStep = parset.get<size_t>("initialTime.timeStep");
+ relativeTime = parset.get<double>("initialTime.relativeTime");
+ relativeTau = parset.get<double>("initialTime.relativeTau");
std::vector<Vector> ell0(bodyCount_);
for (size_t i=0; i<bodyCount_; i++) {
@@ -251,6 +251,8 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
solveLinearProblem(dirichletNodes, levelContactNetwork.matrices().elasticity, ell0, u,
parset.sub("u0.solver"));
+ 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;
@@ -279,6 +281,8 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
Dune::BitSetVector<dims> noNodes(dirichletNodes.size(), false);
solveLinearProblem(noNodes, levelContactNetwork.matrices().mass, accelerationRHS, a,
parset.sub("a0.solver"));
+
+ print(a, "initial a:");
}
private:
diff --git a/src/factories/stackedblocksfactory.cc b/src/factories/stackedblocksfactory.cc
index 746691ab..320a7426 100644
--- a/src/factories/stackedblocksfactory.cc
+++ b/src/factories/stackedblocksfactory.cc
@@ -152,19 +152,27 @@ void StackedBlocksFactory<GridType, dims>::setBoundaryConditions() {
// identify Dirichlet nodes on leaf level
std::shared_ptr<Dune::BitSetVector<dims>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(leafVertexCount);
for (int j=0; j<leafVertexCount; j++) {
- if (leafFaces_[i]->right.containsVertex(j))
- (*velocityDirichletNodes)[j][0] = true;
+ if (leafFaces_[i]->upper.containsVertex(j))
+ (*velocityDirichletNodes)[j][0] = true;
- #if MY_DIM == 3 //TODO: wrong, needs revision
- if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
- zeroDirichletNodes->at(j)[2] = true;
- #endif
+ #if MY_DIM == 3 //TODO: wrong, needs revision
+ if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
+ zeroDirichletNodes->at(j)[2] = true;
+ #endif
}
+ if (i>0) {
std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+
velocityDirichletBoundary->setBoundaryPatch(body->leafView(), velocityDirichletNodes);
- velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
- body->addBoundaryCondition(velocityDirichletBoundary);
+
+ velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+ body->addBoundaryCondition(velocityDirichletBoundary);
+ } else {
+ //velocityDirichletBoundary->setBoundaryFunction(neumannFunction);
+ }
+
+
}
// lower Dirichlet Boundary
@@ -172,8 +180,11 @@ void StackedBlocksFactory<GridType, dims>::setBoundaryConditions() {
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;
+ if (leafFaces_[0]->lower.containsVertex(j)) {
+ for (size_t d=0; d<dims; d++) {
+ (*zeroDirichletNodes)[j][d] = true;
+ }
+ }
}
std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
diff --git a/src/io/vtk.cc b/src/io/vtk.cc
index df3c1033..7aa72c5b 100644
--- a/src/io/vtk.cc
+++ b/src/io/vtk.cc
@@ -8,6 +8,8 @@
#include "vtk.hh"
+#include "../utils/debugutils.hh"
+
template <class VertexBasis, class CellBasis>
BodyVTKWriter<VertexBasis, CellBasis>::BodyVTKWriter(
CellBasis const & cellBasis, VertexBasis const & vertexBasis,
@@ -27,6 +29,7 @@ void BodyVTKWriter<VertexBasis, CellBasis>::write(
vertexBasis_, u, "displacement");
writer.addVertexData(displacementPointer);
+ print(v, "VTK Writer v:");
auto const velocityPointer =
std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
vertexBasis_, v, "velocity");
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index 8212c581..a308c80b 100644
--- a/src/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem-2D.cfg
@@ -1,9 +1,9 @@
# -*- mode:conf -*-
[boundary.friction]
-smallestDiameter = 0.5 # 2e-3 [m]
+smallestDiameter = 0.01 # 2e-3 [m]
[timeSteps]
-refinementTolerance = 1e-5
+refinementTolerance = 1e-5 # 1e-5
[u0.solver]
tolerance = 1e-8
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 878d2cb2..2be31cd8 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -134,7 +134,8 @@ int main(int argc, char *argv[]) {
LevelContactNetwork& levelContactNetwork = stackedBlocksFactory.levelContactNetwork();
const size_t bodyCount = levelContactNetwork.nBodies();
- for (size_t i=0; i<levelContactNetwork.deformedGrids().size(); i++) {
+ for (size_t i=0; i<bodyCount; i++) {
+ printDofLocation(levelContactNetwork.leafView(i));
/* Vector def(levelContactNetwork.deformedGrids()[i]->size(dims));
def = 1;
deformedGridComplex.setDeformation(def, i);*/
@@ -153,6 +154,7 @@ int main(int argc, char *argv[]) {
// ----------------------------
levelContactNetwork.assemble();
+ printMortarBasis<Vector>(levelContactNetwork.nBodyAssembler());
// -----------------
// init input/output
@@ -188,6 +190,14 @@ int main(int argc, char *argv[]) {
else
programState.setupInitialConditions(parset, levelContactNetwork);
+
+ auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
+ for (size_t i=0; i<bodyCount; i++) {
+ levelContactNetwork.body(i)->setDeformation(programState.u[i]);
+ }
+ nBodyAssembler.assembleTransferOperator();
+ nBodyAssembler.assembleObstacle();
+
// ------------------------
// assemble global friction
// ------------------------
@@ -271,7 +281,6 @@ int main(int argc, char *argv[]) {
// -------------------
// Set up TNNMG solver
// -------------------
- auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
BitVector totalDirichletNodes;
levelContactNetwork.totalNodes("dirichlet", totalDirichletNodes);
@@ -288,7 +297,7 @@ int main(int argc, char *argv[]) {
}
}*/
- //print(totalDirichletNodes, "totalDirichletNodes:");
+ print(totalDirichletNodes, "totalDirichletNodes:");
using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
using NonlinearFactory = SolverFactory<Functional, BitVector>;
@@ -304,9 +313,19 @@ int main(int argc, char *argv[]) {
BoundaryNodes dirichletNodes;
levelContactNetwork.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));
+ }
+ }
+
std::vector<const Dune::BitSetVector<1>*> frictionNodes;
levelContactNetwork.frictionNodes(frictionNodes);
+ 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"),
@@ -335,13 +354,17 @@ int main(int argc, char *argv[]) {
std::cout << "Step 1" << std::endl;
std::vector<ScalarVector> coarseAlpha;
+ coarseAlpha.resize(bodyCount);
coarseUpdater.state_->extractAlpha(coarseAlpha);
- std::cout << "Step 2" << std::endl;
+ print(coarseAlpha, "coarseAlpha:");
std::vector<ScalarVector> fineAlpha;
+ fineAlpha.resize(bodyCount);
fineUpdater.state_->extractAlpha(fineAlpha);
+ print(fineAlpha, "fineAlpha:");
+
std::cout << "Step 3" << std::endl;
ScalarVector::field_type energyNorm = 0;
@@ -352,8 +375,11 @@ int main(int argc, char *argv[]) {
if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
continue;
+
energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
}
+ std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance << std::endl;
+ std::cout << "must refine: " << (energyNorm > refinementTolerance) << std::endl;
return energyNorm > refinementTolerance;
};
@@ -382,6 +408,11 @@ int main(int argc, char *argv[]) {
current.rate_->extractAcceleration(programState.a);
current.state_->extractAlpha(programState.alpha);
+ print(programState.u, "current u:");
+ print(programState.v, "current v:");
+ print(programState.a, "current a:");
+ print(programState.alpha, "current alpha:");
+
for (size_t i=0; i<bodyCount; i++) {
levelContactNetwork.body(i)->setDeformation(programState.u[i]);
}
@@ -390,6 +421,11 @@ int main(int argc, char *argv[]) {
report();
+ if (programState.timeStep==2000) {
+ std::cout << "limit of timeSteps reached!" << std::endl;
+ break; // TODO remove after debugging
+ }
+
if (terminationRequested) {
std::cerr << "Terminating prematurely" << std::endl;
break;
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index 8aa21764..5fa9847b 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -10,7 +10,7 @@ restarts.write = false #true
vtk.write = true
[problem]
-finalTime = 1000 # [s]
+finalTime = 100 # [s] #1000
bodyCount = 2
[body]
@@ -40,6 +40,11 @@ b = 0.017 # [ ]
a = 0.020 # [ ]
b = 0.005 # [ ]
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 1e-4 # 1e-6
+
[timeSteps]
scheme = newmark
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index db3dd3cb..f7f01caa 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -77,12 +77,12 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
Matrix bilinearForm;
nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
- print(bilinearForm, "bilinearForm:");
+ //print(bilinearForm, "bilinearForm:");
Vector totalRhs;
nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
- print(totalRhs, "totalRhs:");
+ //print(totalRhs, "totalRhs:");
Vector totalX;
nBodyAssembler_.nodalToTransformed(velocityIterates, totalX);
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 001a3731..b40d02d1 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -4,6 +4,7 @@
#include <dune/solvers/common/wrapownshare.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>
+#include <dune/solvers/solvers/umfpacksolver.hh>
#include "solverfactory.hh"
@@ -22,10 +23,10 @@ SolverFactory<Functional, BitVector>::SolverFactory(
nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
// linearSolverStep_.setPreconditioner(preconditioner_);
- linearSolverStep_ = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::gs(0.0, 0.0));
+ linearSolverStep_ = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::ldlt());
energyNorm_ = std::make_shared<EnergyNorm<Matrix, Vector>>(*linearSolverStep_);
- linearSolver_ = std::make_shared<LinearSolver>(*linearSolverStep_, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), *energyNorm_, Solver::QUIET);
- //linearSolver_ = std::make_shared<LinearSolver>();
+ //linearSolver_ = std::make_shared<LinearSolver>(*linearSolverStep_, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), *energyNorm_, Solver::QUIET);
+ linearSolver_ = std::make_shared<Dune::Solvers::UMFPackSolver<Matrix, Vector>>();
tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_, DefectProjection(), LineSearchSolver());
tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index e39eb945..b4a5efd5 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -11,15 +11,14 @@
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
//#include <dune/solvers/iterationsteps/cgstep.hh>
//#include <dune/solvers/iterationsteps/amgstep.hh>
-//#include <dune/solvers/solvers/umfpacksolver.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/linearsolvers/fastamgsolver.hh>
//#include "tnnmg/functional.hh"
-//#include "tnnmg/tnnmgstep.hh"
+#include "tnnmg/tnnmgstep.hh"
#include "tnnmg/linearization.hh"
#include "tnnmg/linesearchsolver.hh"
#include "tnnmg/localbisectionsolver.hh"
@@ -40,7 +39,9 @@ class SolverFactory {
//using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
//using LinearSolverStep = typename Dune::Solvers::AMGStep<Matrix, Vector>;
using LinearSolverStep = LinearIterationStep<Matrix, Vector>;
- using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ //using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ using LinearSolver = typename Dune::Solvers::LinearSolver<Matrix, Vector>;
+
//using LinearSolver = typename Dune::Solvers::UMFPackSolver<Matrix, Vector>;
//using LinearSolver = typename Dune::TNNMG::FastAMGSolver<Matrix, Vector>;
diff --git a/src/spatial-solving/tnnmg/linesearchsolver.hh b/src/spatial-solving/tnnmg/linesearchsolver.hh
index 352005b2..93471396 100644
--- a/src/spatial-solving/tnnmg/linesearchsolver.hh
+++ b/src/spatial-solving/tnnmg/linesearchsolver.hh
@@ -31,9 +31,13 @@ class LineSearchSolver
if (not ignore) {
int bisectionsteps = 0;
const Bisection globalBisection;
- std::cout << "LineSearchSolver: starting bisection" << std::endl;
+ // std::cout << "LineSearchSolver: starting bisection" << std::endl;
x = globalBisection.minimize(f, 1.0, 0.0, bisectionsteps);
- std::cout << "LineSearchSolver: bisection terminated" << std::endl;
+
+ x = f.domain().projectIn(x);
+
+
+ //std::cout << "LineSearchSolver: bisection terminated" << std::endl;
// x = globalBisection.minimize(f, vNorm, 0.0, bisectionsteps) / vNorm; // used rescaled v in f?
}
}
diff --git a/src/spatial-solving/tnnmg/localproblem.hh b/src/spatial-solving/tnnmg/localproblem.hh
new file mode 100644
index 00000000..8d202afc
--- /dev/null
+++ b/src/spatial-solving/tnnmg/localproblem.hh
@@ -0,0 +1,141 @@
+#ifndef OSC_LOCAL_PROBLEM_HH
+#define OSC_LOCAL_PROBLEM_HH
+
+#include <math.h>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/istl/matrixindexset.hh>
+//#include <dune/istl/superlu.hh>
+#include <dune/istl/umfpack.hh>
+
+#include <dune/fufem/assemblers/localoperatorassembler.hh>
+
+#include "../../utils/debugutils.hh"
+
+template <class MatrixType, class DomainType, class RangeType = DomainType>
+class LocalProblem {
+
+private:
+ typedef typename MatrixType::block_type block_type;
+ typedef typename MatrixType::field_type ctype;
+
+ const static size_t dim = DomainType::block_type::dimension;
+
+ using BitVector = Dune::BitSetVector<dim>;
+
+ MatrixType localMat;
+ const RangeType& rhs_;
+ const BitVector& ignoreNodes_;
+
+ /* size_t flatIndex(const size_t blockIdx, const size_t localBlockIdx) {
+ return blockIdx*dim + localBlockIdx;
+ }
+
+ bool isAllIgnored(const size_t blockIdx) {
+ bool res = true;
+ size_t flatIdx = blockIdx*dim;
+
+ for (size_t d=0; d<dim; d++) {
+ res = res && (globalToLocal_[flatIdx+d] == -1);
+ }
+
+ return res;
+ }*/
+ template <class LocalMat, class LocalBitVector>
+ void computeLocalMat(LocalMat& localMat, const LocalMat& refMat, const LocalBitVector& localIgnore, bool isDiagonal=false) {
+ if (isDiagonal) {
+ for (size_t i=0; i<refMat.N(); i++) {
+ bool isIgnored = localIgnore[i];
+ for (size_t j=0; j<refMat.M(); j++) {
+ if (isIgnored || refMat[i][j]==0)
+ localMat[i][j] = (double) (i==j);
+ else
+ localMat[i][j] = refMat[i][j];
+ }
+ }
+ } else {
+ for (size_t i=0; i<refMat.N(); i++) {
+ bool isIgnored = localIgnore[i];
+ for (size_t j=0; j<refMat.M(); j++) {
+ if (isIgnored)
+ localMat[i][j] = 0;
+ else
+ localMat[i][j] = refMat[i][j];
+ }
+ }
+ }
+ }
+
+public:
+ LocalProblem(const MatrixType& mat,
+ const RangeType& rhs,
+ const BitVector& ignoreNodes) :
+ rhs_(rhs),
+ ignoreNodes_(ignoreNodes)
+ {
+
+ // construct globalToLocal map
+ /* std::vector<int, int> localToGlobal;
+ std::vector<int, int> globalToLocal(ignoreNodes.size()*dim, -1);
+ int localIdx = 0;
+ for (size_t i=0; i<ignoreNodes.size(); ++i) {
+ const auto& ignoreNode = ignoreNodes[i];
+ for (size_t d=0; d<dim; d++) {
+ if (not toBool(ignoreNode[d])) {
+ size_t flatIdx = flatIndex(i, d);
+ localToGlobal[localIdx] = flatIdx;
+ globalToLocal[flatIdx] = localIdx;
+ localIdx++;
+ }
+ }
+ }*/
+
+ // build local stiffness matrix
+ localMat = mat;
+
+ for(size_t rowIdx = 0; rowIdx<localMat.N(); rowIdx++) {
+ const auto& row = mat[rowIdx];
+
+ auto colIt = row.begin();
+ const auto& colEndIt = row.end();
+ for(; colIt!=colEndIt; ++colIt) {
+ const auto colIdx = colIt.index();
+ computeLocalMat(localMat[rowIdx][colIdx], row[colIdx], ignoreNodes_[rowIdx], rowIdx==colIdx);
+ }
+ }
+ }
+
+ MatrixType& getMat() {
+ return localMat;
+ }
+
+ void getLocalRhs(const DomainType& iterate, RangeType& newRhs) {
+ newRhs = rhs_;
+
+ for (size_t i=0; i<newRhs.size(); i++) {
+ for (size_t d=0; d<dim; d++) {
+ if (ignoreNodes_[i][d]) {
+ newRhs[i][d] = iterate[i][d];
+ }
+ }
+ }
+ }
+
+ /* void solve(DomainType& x){
+ #if HAVE_SUPERLU
+ RangeType localRhsCopy(localRhs);
+ Dune::InverseOperatorResult res;
+
+ x.resize(localMat.M());
+
+ Dune::UMFPack<MatrixType> directSolver(localMat);
+ directSolver.apply(x, localRhsCopy, res);
+ #else
+ #error No SuperLU!
+ #endif
+ }*/
+};
+
+#endif
diff --git a/src/spatial-solving/tnnmg/tnnmgstep.hh b/src/spatial-solving/tnnmg/tnnmgstep.hh
index 40f77d06..8c1ce45c 100644
--- a/src/spatial-solving/tnnmg/tnnmgstep.hh
+++ b/src/spatial-solving/tnnmg/tnnmgstep.hh
@@ -13,9 +13,17 @@
#include "dune/solvers/iterationsteps/lineariterationstep.hh"
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/solvers/linearsolver.hh>
+#include <dune/solvers/common/canignore.hh>
+
+#include "localproblem.hh"
+
+
+#include <dune/tnnmg/iterationsteps/linearcorrection.hh>
#include <dune/tectonic/../../src/utils/debugutils.hh>
+namespace Dune {
+namespace TNNMG {
/**
* \brief One iteration of the TNNMG method
@@ -27,9 +35,9 @@ template<class F, class BV, class Linearization,
class DefectProjection,
class LineSearchSolver>
class TNNMGStep :
- public Dune::Solvers::IterationStep<typename F::Vector, BV>
+ public IterationStep<typename F::Vector, BV>
{
- using Base = Dune::Solvers::IterationStep<typename F::Vector, BV>;
+ using Base = IterationStep<typename F::Vector, BV>;
public:
@@ -39,8 +47,8 @@ class TNNMGStep :
using BitVector = typename Base::BitVector;
using ConstrainedBitVector = typename Linearization::ConstrainedBitVector;
using Functional = F;
- using IterativeSolver = Dune::Solvers::IterativeSolver< ConstrainedVector, Dune::Solvers::DefaultBitVector_t<ConstrainedVector> >;
- using LinearSolver = Dune::Solvers::LinearSolver< ConstrainedMatrix, ConstrainedVector >;
+ using IterativeSolver = Solvers::IterativeSolver< ConstrainedVector, Solvers::DefaultBitVector_t<ConstrainedVector> >;
+ using LinearSolver = 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
@@ -57,6 +65,7 @@ class TNNMGStep :
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
+ //linearCorrection_(makeLinearCorrection<ConstrainedMatrix>(iterativeSolver)),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -76,6 +85,7 @@ class TNNMGStep :
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
+ linearSolver_(linearSolver),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -88,14 +98,15 @@ class TNNMGStep :
*/
TNNMGStep(const Functional& f,
Vector& x,
- std::shared_ptr<Dune::Solvers::IterationStep<Vector,BitVector> > nonlinearSmoother,
- std::shared_ptr<Dune::Solvers::LinearIterationStep<ConstrainedMatrix,ConstrainedVector> > linearIterationStep,
+ std::shared_ptr<Solvers::IterationStep<Vector,BitVector> > nonlinearSmoother,
+ std::shared_ptr<Solvers::LinearIterationStep<ConstrainedMatrix,ConstrainedVector> > linearIterationStep,
unsigned int noOfLinearIterationSteps,
const DefectProjection& projection,
const LineSearchSolver& lineSolver)
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
+ //linearCorrection_(makeLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -123,15 +134,18 @@ class TNNMGStep :
const auto& ignore = (*this->ignoreNodes_);
auto& x = *getIterate();
- print(f.quadraticPart(), "f.quadraticPart():");
+ //std::cout << "TNNMGStep::iterate " << std::endl;
- print(f.linearPart(), "f.linearPart():");
+ //print(f.quadraticPart(), "f.quadraticPart():");
+
+ //print(f.linearPart(), "f.linearPart():");
// Nonlinear presmoothing
for (std::size_t i=0; i<preSmoothingSteps_; ++i)
nonlinearSmoother_->iterate();
- print(x, "TNNMG iterate after smoothing:");
+ //std::cout << "- nonlinear presmoothing: success" << std::endl;
+ //print(x, "TNNMG iterate after smoothing:");
// Compute constraint/truncated linearization
if (not linearization_)
@@ -142,28 +156,66 @@ class TNNMGStep :
auto&& A = linearization_->hessian();
auto&& r = linearization_->negativeGradient();
- print(A, "TNNMG Linear problem matrix:");
- print(r, "TNNMG Linear problem rhs:");
+ //print(A, "TNNMG Linear problem matrix:");
+ //print(r, "TNNMG Linear problem rhs:");
// Compute inexact solution of the linearized problem
Solvers::resizeInitializeZero(correction_, x);
Solvers::resizeInitializeZero(constrainedCorrection_, r);
- linearCorrection_(A, constrainedCorrection_, r);
+ // std::cout << "- computing linear correction..." << std::endl;
+
+ /* auto canIgnoreCast = std::dynamic_pointer_cast<CanIgnore<Dune::Solvers::DefaultBitVector_t<Vector>>>( linearSolver_ );
+ if (canIgnoreCast)
+ canIgnoreCast->setIgnore(*this->ignoreNodes_);*/
+
+ LocalProblem<ConstrainedMatrix, ConstrainedVector> localProblem(A, r, *this->ignoreNodes_);
+ Vector newR;
+ localProblem.getLocalRhs(constrainedCorrection_, newR);
+
+ /*print(*this->ignoreNodes_, "ignoreNodes:");
+ print(A, "A:");
+ print(localProblem.getMat(), "localMat:");*/
+
+ linearSolver_->setProblem(localProblem.getMat(), constrainedCorrection_, newR);
+ linearSolver_->preprocess();
+ linearSolver_->solve();
+
+ // linearCorrection_(A, constrainedCorrection_, r);
+
+ // std::cout << "- linear correction: success" << std::endl;
- print(constrainedCorrection_, "contrained correction:");
+ // print(constrainedCorrection_, "contrained correction:");
linearization_->extendCorrection(constrainedCorrection_, correction_);
- print(correction_, "correction:");
+ // std::cout << "- extended correction: success" << std::endl;
+ // print(correction_, "correction:");
// Project onto admissible set
projection_(f, x, correction_);
+ // std::cout << "- projection onto admissible set: success" << std::endl;
+
+ // print(correction_, "correction:");
+
+ double const correctionNorm = correction_.two_norm();
+ correction_ /= correctionNorm;
+
+
// Line search
auto fv = directionalRestriction(f, x, correction_);
+
+ /* std::cout << fv.quadraticPart() << " f.quadraticPart():" << std::endl;
+ std::cout << fv.linearPart() << " f.linearPart():" << std::endl;
+ std::cout << fv.subDifferential(0) << " f.subDifferential(0):" << std::endl;
+ std::cout << fv.domain()[0] << " " << fv.domain()[0] << " f.domain():" << std::endl;
+
+ std::cout << "- setup directional restriction: success" << std::endl; */
dampingFactor_ = 0;
lineSolver_(dampingFactor_, fv, false);
+
+ // std::cout << "- computing damping factor: success" << std::endl;
if (std::isnan(dampingFactor_))
dampingFactor_ = 0;
correction_ *= dampingFactor_;
@@ -196,6 +248,9 @@ class TNNMGStep :
std::shared_ptr<Linearization> linearization_;
+ //! \brief linear correction
+ std::shared_ptr<LinearSolver> linearSolver_;
+
typename Linearization::ConstrainedVector constrainedCorrection_;
Vector correction_;
DefectProjection projection_;
@@ -203,4 +258,8 @@ class TNNMGStep :
double dampingFactor_;
};
+
+} // end namespace TNNMG
+} // end namespace Dune
+
#endif
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index 74396864..e5a016d1 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -71,7 +71,7 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
bool didCoarsen = false;
iterationRegister_.reset();
UpdatersWithCount R2;
- UpdatersWithCount C;
+ UpdatersWithCount C; /*
while (relativeTime_ + relativeTau_ <= 1.0) {
R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
@@ -104,7 +104,7 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
R1_ = F1;
R2 = F2;
}
- }
+ } */
std::cout << "AdaptiveTimeStepper::advance() ...";
@@ -112,9 +112,9 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
relativeTime_ += relativeTau_;
current_ = R1_.updaters;
- //UpdatersWithCount emptyR1;
- //R1_ = emptyR1;
- R1_ = R2;
+ UpdatersWithCount emptyR1;
+ R1_ = emptyR1;
+ //R1_ = R2;
std::cout << " done" << std::endl;
diff --git a/src/time-stepping/rate/backward_euler.cc b/src/time-stepping/rate/backward_euler.cc
index b7a915f9..95946e60 100644
--- a/src/time-stepping/rate/backward_euler.cc
+++ b/src/time-stepping/rate/backward_euler.cc
@@ -3,6 +3,8 @@
#include "backward_euler.hh"
+#include "../../utils/debugutils.hh"
+
template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::BackwardEuler(
const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index 4e61938b..73bbe3a1 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -95,15 +95,20 @@ void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
double dirichletValue;
bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
+ std::cout << "dirichletValue: " << dirichletValue << std::endl;
+
+
for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
for (size_t j=0; j<dim; ++j) {
if (bcDirichletNodes[k][j]) {
+ std::cout << k << " " << j << std::endl;
bodyIterate[k][j] = dirichletValue;
}
}
}
}
}
+ print(iterate, "iterate:");
}
template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
diff --git a/src/time-stepping/state/ageinglawstateupdater.cc b/src/time-stepping/state/ageinglawstateupdater.cc
index a4c47377..796f8878 100644
--- a/src/time-stepping/state/ageinglawstateupdater.cc
+++ b/src/time-stepping/state/ageinglawstateupdater.cc
@@ -66,6 +66,7 @@ template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
ScalarVector& alpha) {
+ std::cout << "alpha size: " << alpha.size() << " nodes_.size() " << nodes_.size() << std::endl;
if (alpha.size() != nodes_.size()) {
alpha.resize(nodes_.size());
}
diff --git a/src/time-stepping/state/stateupdater.hh b/src/time-stepping/state/stateupdater.hh
index 458fc81a..ed7a5241 100644
--- a/src/time-stepping/state/stateupdater.hh
+++ b/src/time-stepping/state/stateupdater.hh
@@ -57,7 +57,7 @@ template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
}
}
- void extractAlpha(std::vector<ScalarVector>& alpha) {
+ void extractAlpha(std::vector<ScalarVector>& alpha) {
for (size_t i=0; i<localStateUpdaters_.size(); i++) {
auto& localStateUpdater = localStateUpdaters_[i];
localStateUpdater->extractAlpha(alpha[localStateUpdater->bodyIndex()]);
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index cb289f16..5190c427 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -118,7 +118,7 @@
void print(const std::vector<T>& x, std::string 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 << std::endl;
}
@@ -325,4 +325,48 @@
std::cout << std::endl;
}
+
+ template <class GridView>
+ void printDofLocation(const GridView& gridView) {
+
+ std::cout << "-- GridView vertices --" << std::endl;
+ std::cout << "Index Coords " << std::endl;
+ std::cout << "-----------------------" << std::endl;
+
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
+
+ for (auto&& it : vertices(gridView)) {
+ const auto i = vertexMapper.index(it);
+ const auto& geometry = it.geometry();
+ const auto& corner = geometry.corner(0);
+
+ std::cout << std::setw(5) << i << " ";
+ for(size_t i=0; i<corner.size(); ++i) {
+ std::cout << std::setw(5) << std::setprecision(3) << corner[i] << " ";
+ }
+ std::cout << std::endl;
+ }
+ std::cout << std::endl;
+ }
+
+ template <class Vector, class NBodyAssembler>
+ void printMortarBasis(const NBodyAssembler& nBodyAssembler) {
+ std::cout << "-- Mortar basis in canonical coords --" << std::endl;
+ std::cout << "--------------------------------------" << std::endl;
+
+ const auto& dim = nBodyAssembler.getTransformationMatrix().N();
+
+ Vector mortarBasisVector(dim);
+ std::vector<Vector> canonicalBasisVector(nBodyAssembler.nGrids());
+ for (size_t i=0; i<dim; i++) {
+ mortarBasisVector = 0;
+ mortarBasisVector[i] = 1;
+
+ nBodyAssembler.postprocess(mortarBasisVector, canonicalBasisVector);
+ print(canonicalBasisVector, "canonicalBasisVector " + std::to_string(i));
+ }
+ std::cout << std::endl;
+ }
+
#endif
diff --git a/src/utils/tobool.hh b/src/utils/tobool.hh
index 6f699601..bd8912f8 100644
--- a/src/utils/tobool.hh
+++ b/src/utils/tobool.hh
@@ -4,7 +4,7 @@
#include <dune/common/bitsetvector.hh>
template <class Alloc>
-bool toBool(Dune::BitSetVectorConstReference<1, Alloc> x) {
+bool toBool(const Dune::BitSetVectorConstReference<1, Alloc> x) {
return x[0];
}
#endif
--
GitLab