From cfdb083aa5077f6a6980fd8ce7115a84de12e15a Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 14 Nov 2011 14:39:53 +0100
Subject: [PATCH] Bug fix: We're computation increments here!
Also, time-independent neumann conditions
---
src/one-body-sample.cc | 31 +++++++++++++++++++++++--------
1 file changed, 23 insertions(+), 8 deletions(-)
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 984eb1bc..4b04b4dd 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -96,8 +96,8 @@ void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
int run) { // constant sample function on neumann boundary
BoundaryPatch<GridView> neumannBoundary(gridView, neumannNodes);
LocalVectorType SampleVector(0);
- // FIXME: random values
- SampleVector[0] = 1 + run / 200;
+ // FIXME: random values (time-independent)
+ SampleVector[0] = 1;
SampleVector[1] = 0;
ConstantFunction<LocalVectorType, LocalVectorType> fNeumann(SampleVector);
NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
@@ -194,6 +194,15 @@ int main(int argc, char *argv[]) {
VectorType u1(grid.size(grid.maxLevel(), dim));
u1 = 0;
VectorType u3 = u1;
+
+ VectorType u1_diff_old(grid.size(grid.maxLevel(), dim));
+ u1_diff_old = 0;
+ VectorType u3_diff_old = u1_diff_old;
+
+ VectorType u1_diff_new(grid.size(grid.maxLevel(), dim));
+ u1_diff_new = 0;
+ VectorType u3_diff_new = u1_diff_new;
+
CellVectorType vonMisesStress;
VectorType b1;
@@ -213,8 +222,8 @@ int main(int argc, char *argv[]) {
b3 = b1;
// b += linear update
- stiffnessMatrix.umv(u1, b1);
- stiffnessMatrix.umv(u3, b3);
+ stiffnessMatrix.umv(u1_diff_old, b1);
+ stiffnessMatrix.umv(u3_diff_old, b3);
// {{{ Assemble terms for the nonlinearity
@@ -258,10 +267,10 @@ int main(int argc, char *argv[]) {
// }}}
{
- MyConvexProblemType myConvexProblem(stiffnessMatrix,
- *myGlobalNonlinearity, b1, u1);
+ MyConvexProblemType myConvexProblem(
+ stiffnessMatrix, *myGlobalNonlinearity, b1, u1_diff_new);
MyBlockProblemType myBlockProblem(parset, myConvexProblem);
- nonlinearGSStep.setProblem(u1, myBlockProblem);
+ nonlinearGSStep.setProblem(u1_diff_new, myBlockProblem);
LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations,
solver_tolerance, &energyNorm,
@@ -269,6 +278,9 @@ int main(int argc, char *argv[]) {
solver.solve();
}
+ u1 += u1_diff_new;
+ u1_diff_old = u1_diff_new;
+
auto *displacement =
new BasisGridFunction<P1Basis, VectorType>(p1Basis, u1);
VonMisesStressAssembler<GridType> localStressAssembler(E, nu,
@@ -294,7 +306,7 @@ int main(int argc, char *argv[]) {
// same results if phi vanishes (e.g. because the normalstress is zero)
{
TruncatedBlockGSStep<OperatorType, VectorType> blockGSStep(
- stiffnessMatrix, u3, b3);
+ stiffnessMatrix, u3_diff_new, b3);
blockGSStep.ignoreNodes_ = &ignoreNodes;
LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
@@ -302,6 +314,9 @@ int main(int argc, char *argv[]) {
Solver::QUIET);
solver.solve();
}
+
+ u3 += u3_diff_new;
+ u3_diff_old = u3_diff_new;
}
std::cout << std::endl;
--
GitLab