From 3fdd548cd92ca422be1d9fbf50a03195fb2108c5 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 12 Apr 2004 22:18:02 +0000 Subject: [PATCH] Fix some issues. Get closer to a working example. git-svn-id: https://svn.dealii.org/trunk@8997 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-15/step-15.cc | 201 +++++++++++++++++----------- 1 file changed, 120 insertions(+), 81 deletions(-) diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index cdc12ba6b9..2f4e11d66e 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -4,7 +4,7 @@ /* $Id$ */ /* Version: $Name$ */ /* */ -/* Copyright (C) 2002, 2003 by the deal.II authors */ +/* Copyright (C) 2002, 2003, 2004 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -18,7 +18,6 @@ #include #include // -#include #include #include #include @@ -65,13 +64,17 @@ class MinimizationProblem private: void setup_system (); - void assemble_newton_step (const bool p); - void solve (); + void assemble_step (const bool p); + double line_search (const Vector & update) const; + void do_step (); void initialize (); void refine_grid (); - double energy () const; void output_results (const unsigned int cycle) const; + static double energy (const DoFHandler &dof_handler, + const Vector &function); + + Triangulation triangulation; FE_Q fe; @@ -80,11 +83,10 @@ class MinimizationProblem ConstraintMatrix hanging_node_constraints; SparsityPattern sparsity_pattern; - SparseMatrix newton_matrix; + SparseMatrix matrix; Vector present_solution; - Vector newton_residual; - Vector newton_update; + Vector residual; }; @@ -103,7 +105,7 @@ class InitializationValues : public Function<1> double InitializationValues::value (const Point<1> &p, const unsigned int) const { - return std::pow(p(0), 1./3.); + return std::pow(p(0), 1.); } @@ -155,12 +157,11 @@ double gradient_power (const Tensor<1,dim> &v, template -void MinimizationProblem::assemble_newton_step (const bool p) +void MinimizationProblem::assemble_step (const bool p) { if (p) - newton_matrix.reinit (sparsity_pattern); - newton_update.reinit (dof_handler.n_dofs()); - newton_residual.reinit (dof_handler.n_dofs()); + matrix.reinit (sparsity_pattern); + residual.reinit (dof_handler.n_dofs()); QGauss3 quadrature_formula; @@ -209,36 +210,41 @@ void MinimizationProblem::assemble_newton_step (const bool p) if (p) for (unsigned int j=0; j::assemble_newton_step (const bool p) { if (p) for (unsigned int j=0; j boundary_values; VectorTools::interpolate_boundary_values (dof_handler, @@ -285,43 +291,69 @@ void MinimizationProblem::assemble_newton_step (const bool p) 1, ZeroFunction(), boundary_values); + Vector dummy (residual.size()); MatrixTools::apply_boundary_values (boundary_values, - newton_matrix, - newton_update, - newton_residual); -// std::cout << " res=" << newton_residual.l2_norm() << std::endl; + matrix, + dummy, + residual); } template -void MinimizationProblem::solve () +double +MinimizationProblem::line_search (const Vector &update) const { - SolverControl solver_control (1000, - 1e-3*newton_residual.l2_norm()); - PrimitiveVectorMemory<> vector_memory; - SolverGMRES<> gmres (solver_control, vector_memory); - - PreconditionJacobi<> preconditioner; - preconditioner.initialize(newton_matrix, 1.2); - -// std::cout << "------------------------------" << std::endl; -// newton_matrix.print_formatted(cout); -// std::cout << std::endl; -// for (unsigned int i=0; i tmp (present_solution.size()); - gmres.solve (newton_matrix, newton_update, newton_residual, - preconditioner); + for (double a=.01; a<=10; a*=1.5) + { + tmp = present_solution; + tmp.add (a, update); + const double e = energy(dof_handler, tmp); - present_solution += newton_update; - hanging_node_constraints.distribute (present_solution); + std::cout << "XX" << a << ' ' << e << std::endl; + if (e < optimal_energy) + { + optimal_energy = e; + alpha = a; + } + } + + std::cout << " Step length : " << alpha << ' ' << optimal_energy << std::endl; + + return alpha; } + +template +void MinimizationProblem::do_step () +{ + assemble_step (true); + + Vector update (present_solution.size()); + { + SolverControl solver_control (1000, + 1e-3*residual.l2_norm()); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> solver (solver_control, vector_memory); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(matrix, 1.2); + + solver.solve (matrix, update, residual, + preconditioner); + hanging_node_constraints.distribute (update); + } + + present_solution.add (line_search (update), update); +} + + template void MinimizationProblem::initialize () { @@ -360,17 +392,25 @@ void MinimizationProblem::refine_grid () Vector tmp (dof_handler.n_dofs()); solution_transfer.interpolate (present_solution, tmp); present_solution = tmp; + + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + hanging_node_constraints.distribute (present_solution); } template -double MinimizationProblem::energy () const +double +MinimizationProblem::energy (const DoFHandler &dof_handler, + const Vector &function) { double energy = 0.; QGauss3 quadrature_formula; - FEValues fe_values (fe, quadrature_formula, + FEValues fe_values (dof_handler.get_fe(), quadrature_formula, UpdateFlags(update_values | update_gradients | update_q_points | @@ -387,9 +427,9 @@ double MinimizationProblem::energy () const for (; cell!=endc; ++cell) { fe_values.reinit (cell); - fe_values.get_function_values (present_solution, + fe_values.get_function_values (function, local_solution_values); - fe_values.get_function_grads (present_solution, + fe_values.get_function_grads (function, local_solution_grads); for (unsigned int q_point=0; q_point void MinimizationProblem::run () { GridGenerator::hyper_cube (triangulation, 0., 1.); - triangulation.refine_global (1); + triangulation.refine_global (4); initialize (); - for (unsigned int refinement_cycle=0; refinement_cycle<50; + for (unsigned int refinement_cycle=0; refinement_cycle<5; ++refinement_cycle) { std::cout << "Cycle " << refinement_cycle << ':' << std::endl; @@ -455,12 +495,11 @@ void MinimizationProblem::run () setup_system (); unsigned int iteration=0; - for (; iteration<1000; ++iteration) + for (; iteration<5; ++iteration) { - assemble_newton_step (true); - solve (); + do_step (); - if (newton_residual.l2_norm() < 1.e-12) + if (residual.l2_norm() < 1.e-4) break; }; output_results (refinement_cycle); @@ -469,7 +508,7 @@ void MinimizationProblem::run () << iteration << std::endl; std::cout << " Energy : " - << energy () + << energy (dof_handler, present_solution) << std::endl; refine_grid (); -- 2.39.5