From: frohne Date: Fri, 18 Nov 2011 20:04:51 +0000 (+0000) Subject: changes of stop criterion and the active set condition X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2031f441ebcae4bf5582e280a39ec7c608b33da9;p=dealii-svn.git changes of stop criterion and the active set condition git-svn-id: https://svn.dealii.org/trunk@24754 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-41/step-41.cc b/deal.II/examples/step-41/step-41.cc index c44c5f68ba..131b883ae4 100644 --- a/deal.II/examples/step-41/step-41.cc +++ b/deal.II/examples/step-41/step-41.cc @@ -64,6 +64,7 @@ class Step41 void make_grid (); void setup_system(); void assemble_system (); + void assemble_mass_matrix (); void projection_active_set (); void solve (); void output_results (const std::string& title) const; @@ -78,6 +79,7 @@ class Step41 SparsityPattern sparsity_pattern; TrilinosWrappers::SparseMatrix system_matrix; TrilinosWrappers::SparseMatrix system_matrix_complete; + TrilinosWrappers::SparseMatrix mass_matrix; TrilinosWrappers::Vector solution; TrilinosWrappers::Vector tmp_solution; @@ -220,6 +222,7 @@ void Step41::setup_system () system_matrix.reinit (sparsity_pattern); system_matrix_complete.reinit (sparsity_pattern); + mass_matrix.reinit (sparsity_pattern); solution.reinit (dof_handler.n_dofs()); tmp_solution.reinit (dof_handler.n_dofs()); @@ -297,6 +300,51 @@ void Step41::assemble_system () } } +template +void Step41::assemble_mass_matrix () +{ + QTrapez quadrature_formula; + + FEValues fe_values (fe, quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_matrix = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + // This function apply the constraints + // to the system matrix and system rhs. + // The true parameter is set to make sure + // that the system rhs contains correct + // values in the rows with inhomogeneity + // constraints. + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + mass_matrix); + } +} + // @sect4{Step41::projection_active_set} // Updating of the active set which means to @@ -324,6 +372,7 @@ void Step41::projection_active_set () // to find and supply the constraints for the // obstacle condition active_set = 0.0; + const double c = 100.0; for (; cell!=endc; ++cell) for (unsigned int v=0; v::vertices_per_cell; ++v) { @@ -340,7 +389,8 @@ void Step41::projection_active_set () // the diag-entry of the mass-matrix. // TODO: I have to check the condition - if ((resid_vector (index_x)*std::pow (2, 2*n_refinements)*diag_mass_matrix_vector (index_x) >= solution_index_x - obstacle_value)) + if (resid_vector (index_x) + + diag_mass_matrix_vector (index_x)*c*(obstacle_value - solution_index_x) > 0) { constraints.add_line (index_x); constraints.set_inhomogeneity (index_x, obstacle_value); @@ -370,7 +420,7 @@ void Step41::projection_active_set () template void Step41::solve () { - ReductionControl reduction_control (100, 1e-12, 1e-2); + ReductionControl reduction_control (100, 1e-12, 1e-3); SolverCG solver (reduction_control); TrilinosWrappers::PreconditionAMG precondition; precondition.initialize (system_matrix); @@ -382,6 +432,8 @@ void Step41::solve () << " CG iterations needed to obtain convergence with an error: " << reduction_control.last_value() << std::endl; + + constraints.distribute (solution); } // @sect4{Step41::output_results} @@ -395,14 +447,14 @@ void Step41::output_results (const std::string& title) const DataOut data_out; data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (tmp_solution, "Displacement"); - data_out.add_data_vector (resid_vector, "Residual"); + // data_out.add_data_vector (tmp_solution, "Displacement"); + // data_out.add_data_vector (resid_vector, "Residual"); data_out.add_data_vector (active_set, "ActiveSet"); data_out.build_patches (); std::ofstream output_vtk ((title + ".vtk").c_str ()); - data_out.write_vtk (output_vtk); + data_out.write_gnuplot (output_vtk); } @@ -438,19 +490,24 @@ void Step41::run () // step of the active-set-iteration system_matrix_complete.copy_from (system_matrix); system_rhs_complete = system_rhs; - + + // to compute the factor which is used + // to scale the residual. You can consider + // this diagonal matrix as the discretization + // of a lagrange multiplier for the + // contact force + assemble_mass_matrix (); + for (unsigned int j=0; j::run () resid_vector = 0; resid_vector -= system_rhs_complete; system_matrix_complete.vmult_add (resid_vector, solution); - for (unsigned int k = 0; k 0) - resid_vector (k) = 0; std::cout<< "Update Active Set:"<