From e9ac6528e05e886ffdd28f48d3d978eb7017d158 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 16 Nov 2011 00:03:48 +0100
Subject: [PATCH] Bring back BlockGSStep now that it works
Do not print vonMisesStress per node
That doesn't make sense -- we only have it for elements!
---
src/one-body-sample.cc | 37 +++++++++++++++++++++++++++++++------
1 file changed, 31 insertions(+), 6 deletions(-)
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 85d849f9..db4fe8ec 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -40,6 +40,7 @@
#include <dune/fufem/functionspacebases/p0basis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/solvers/common/numproc.hh> // Solver::FULL
+#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
@@ -199,19 +200,23 @@ int main(int argc, char *argv[]) {
VectorType u1(grid.size(grid.maxLevel(), dim));
u1 = 0;
+ VectorType u2 = u1;
VectorType u3 = u1;
VectorType u1_diff_old(grid.size(grid.maxLevel(), dim));
u1_diff_old = 0;
+ VectorType u2_diff_old = u1_diff_old;
VectorType u3_diff_old = u1_diff_old;
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;
CellVectorType vonMisesStress;
VectorType b1;
+ VectorType b2;
VectorType b3;
std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
@@ -225,10 +230,12 @@ int main(int argc, char *argv[]) {
// b = neumann
assemble_neumann<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
leafView, p1Basis, neumannNodes, b1, run);
+ b2 = b1;
b3 = b1;
// b += linear update
stiffnessMatrix.umv(u1_diff_old, b1);
+ stiffnessMatrix.umv(u2_diff_old, b2);
stiffnessMatrix.umv(u3_diff_old, b3);
// {{{ Assemble terms for the nonlinearity
@@ -309,6 +316,20 @@ int main(int argc, char *argv[]) {
writer.write(filename.c_str());
}
+ {
+ BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix,
+ u2_diff_new, b2);
+ blockGSStep.ignoreNodes_ = &ignoreNodes;
+
+ LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
+ solver_tolerance, &energyNorm,
+ Solver::QUIET); // Solver::QUIET);
+ solver.solve();
+ }
+
+ u2 += u2_diff_new;
+ u2_diff_old = u2_diff_new;
+
// Use a linear solver for comparison; should return roughly the
// same results if phi vanishes (e.g. because the normalstress is zero)
{
@@ -327,18 +348,22 @@ int main(int argc, char *argv[]) {
}
std::cout << std::endl;
+ VectorType diff2 = u2;
+ diff2 -= u1;
+ std::cout << "sup |u1 - u2| = " << diff2.infinity_norm() << ", "
+ << "|u1 - u2| = " << diff2.two_norm() << std::endl;
+
VectorType diff3 = u3;
diff3 -= u1;
- std::cout << "sup |u1 - u3| = " << diff3.infinity_norm() << std::endl;
- std::cout << "|u1 - u3| = " << diff3.two_norm() << std::endl;
+ std::cout << "sup |u1 - u3| = " << diff3.infinity_norm() << ", "
+ << "|u1 - u3| = " << diff3.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|u3[%02d] = %+3e "
- "%|80t|s[%02d] = %+3e") %
- i % u1[i] % i % u3[i] % i %
- vonMisesStress[i] << std::endl;
+ 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;
}
catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
--
GitLab