diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 0dfd443e5f726fa2cb75a79b23fde91b1749101e..3bce9241bf16ba97d7d1b21f4c98b1954a3251be 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -83,11 +83,12 @@ int main() { for (size_t j = 0; j < dim; ++j) f[i][j] = i + j; - VectorType u(grid.size(grid.maxLevel(), dim)); + VectorType u1(grid.size(grid.maxLevel(), dim)); // Fill initial guess with other semi-random numbers for (size_t i = 0; i < f.size(); ++i) for (size_t j = 0; j < dim; ++j) - u[i][j] = 17 * (j + 1); + u1[i][j] = 17 * (j + 1); + VectorType u2 = u1; Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes( grid.size(grid.maxLevel(), dim), false); @@ -138,8 +139,6 @@ int main() { std::cout << "Number of extremal nodes: " << extremal_nodes << std::endl; } - VectorType copy_of_u = u; - { // experiment with convex problems and the like typedef ZeroNonlinearity<SmallVector, SmallMatrix> NonlinearityType; @@ -149,14 +148,14 @@ int main() { VectorType unused_vector(grid.size(grid.maxLevel(), dim), 0); ConvexProblemType myConvexProblem(1, stiffnessMatrix, 0, unused_vector, - phi, f, copy_of_u); + phi, f, u1); typedef MyBlockProblem<ConvexProblemType> MyBlockProblemType; MyBlockProblemType myBlockProblem(myConvexProblem); GenericNonlinearGS<MyBlockProblemType> nonlinearGSStep; nonlinearGSStep.ignoreNodes_ = &ignoreNodes; - nonlinearGSStep.setProblem(copy_of_u, myBlockProblem); + nonlinearGSStep.setProblem(u1, myBlockProblem); // FIXME: Does this make any sense? EnergyNorm<OperatorType, VectorType> energyNorm(stiffnessMatrix); @@ -168,7 +167,7 @@ int main() { } // TODO: Why does blockGSStep even provide a default constructor? - BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u, f); + BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u2, f); blockGSStep.ignoreNodes_ = &ignoreNodes; // FIXME: Does this make any sense? @@ -180,8 +179,8 @@ int main() { solver.solve(); - VectorType diff = u; - diff -= copy_of_u; + VectorType diff = u2; + diff -= u1; std::cout << diff << std::endl; std::cout << std::endl; std::cout << diff.infinity_norm() << std::endl;