From 9e463158989ba11a2bc4b9bfd1985751cbf4ee45 Mon Sep 17 00:00:00 2001 From: frohne Date: Thu, 23 Aug 2012 23:52:33 +0000 Subject: [PATCH] Implementing the KellyEstimator in refine_grid () git-svn-id: https://svn.dealii.org/trunk@26107 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 243 +++++++++++++++------------- 1 file changed, 128 insertions(+), 115 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 3a9a898ef8..3085e6a0ce 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -54,6 +54,7 @@ #include #include +#include #include #include #include @@ -111,9 +112,9 @@ namespace Step42 void dirichlet_constraints (); void solve (); void solve_newton (); - void output_results (const std::string& title) const; + void refine_grid (); void move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const; - void output_results (TrilinosWrappers::MPI::Vector vector, const std::string& title) const; + void output_results (const std::string& title) const; int n_refinements_global; int n_refinements_local; @@ -518,47 +519,47 @@ namespace Step42 triangulation.refine_global (n_refinements_global); // Lokale Verfeinerung des Gitters - for (int step=0; step::faces_per_cell; ++face) - { - // if (cell->face (face)->at_boundary() - // && cell->face (face)->boundary_indicator () == 9) - // { - // cell->set_refine_flag (); - // break; - // } - // else if (cell->level () == n_refinements + n_refinements_local - 1) - // { - // cell->set_refine_flag (); - // break; - // } - - // if (cell->face (face)->at_boundary() - // && cell->face (face)->boundary_indicator () == 9) - // { - // if (cell->face (face)->vertex (0)(0) <= 0.7 && - // cell->face (face)->vertex (1)(0) >= 0.3 && - // cell->face (face)->vertex (0)(1) <= 0.875 && - // cell->face (face)->vertex (2)(1) >= 0.125) - // { - // cell->set_refine_flag (); - // break; - // } - // } - - if (step == 0 && - cell->center ()(2) < n_refinements_local*9.0/64.0) - { - cell->set_refine_flag (); - break; - } - }; - triangulation.execute_coarsening_and_refinement (); - }; +// for (int step=0; step::faces_per_cell; ++face) +// { +// // if (cell->face (face)->at_boundary() +// // && cell->face (face)->boundary_indicator () == 9) +// // { +// // cell->set_refine_flag (); +// // break; +// // } +// // else if (cell->level () == n_refinements + n_refinements_local - 1) +// // { +// // cell->set_refine_flag (); +// // break; +// // } +// +// // if (cell->face (face)->at_boundary() +// // && cell->face (face)->boundary_indicator () == 9) +// // { +// // if (cell->face (face)->vertex (0)(0) <= 0.7 && +// // cell->face (face)->vertex (1)(0) >= 0.3 && +// // cell->face (face)->vertex (0)(1) <= 0.875 && +// // cell->face (face)->vertex (2)(1) >= 0.125) +// // { +// // cell->set_refine_flag (); +// // break; +// // } +// // } +// +// if (step == 0 && +// cell->center ()(2) < n_refinements_local*9.0/64.0) +// { +// cell->set_refine_flag (); +// break; +// } +// }; +// triangulation.execute_coarsening_and_refinement (); +// }; } template @@ -576,6 +577,7 @@ namespace Step42 // setup hanging nodes and dirichlet constraints { + constraints_hanging_nodes.clear (); constraints_hanging_nodes.reinit (locally_relevant_dofs); DoFTools::make_hanging_node_constraints (dof_handler, constraints_hanging_nodes); @@ -601,6 +603,7 @@ namespace Step42 old_solution.reinit (system_rhs_newton); resid_vector.reinit (system_rhs_newton); diag_mass_matrix_vector.reinit (system_rhs_newton); + active_set.clear (); active_set.set_size (locally_relevant_dofs.size ()); } @@ -1230,20 +1233,6 @@ namespace Step42 active_set_old = active_set; } // End of active-set-loop - - pcout<< "Creating output." < vertex_touched (triangulation.n_vertices(), + false); + + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (); + cell != dof_handler.end(); ++cell) + if (cell->is_locally_owned()) + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + if (vertex_touched[cell->vertex_index(v)] == false) + { + vertex_touched[cell->vertex_index(v)] = true; + + Point vertex_displacement; + for (unsigned int d=0; dvertex_dof_index(v,d)) != 0) + vertex_displacement[d] + = _complete_displacement(cell->vertex_dof_index(v,d)); + } + + cell->vertex(v) += vertex_displacement; + } + } + } + template void PlasticityContactProblem::output_results (const std::string& title) const { @@ -1315,72 +1352,48 @@ namespace Step42 } template - void PlasticityContactProblem::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const + void PlasticityContactProblem::run () { - pcout<< "Moving mesh." < vertex_touched (triangulation.n_vertices(), - false); + pcout << "Solving problem in " << dim << " space dimensions." << std::endl; + + Timer t; + run_time.resize (8); - for (typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active (); - cell != dof_handler.end(); ++cell) - if (cell->is_locally_owned()) - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - if (vertex_touched[cell->vertex_index(v)] == false) - { - vertex_touched[cell->vertex_index(v)] = true; + const unsigned int n_cycles = 5; + for (unsigned int cycle=0; cycle vertex_displacement; - for (unsigned int d=0; dvertex_dof_index(v,d)) != 0) - vertex_displacement[d] - = _complete_displacement(cell->vertex_dof_index(v,d)); - } + MPI_Barrier (mpi_communicator); + t.restart(); - cell->vertex(v) += vertex_displacement; - } + if (cycle == 0) + { + make_grid(); } - } + else + refine_grid (); - template - void PlasticityContactProblem::output_results (TrilinosWrappers::MPI::Vector vector, - const std::string& title) const - { - DataOut data_out; + setup_system (); - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (vector, "vector_to_plot"); - - data_out.build_patches (); - - std::ofstream output_vtk (dim == 2 ? - (title + ".vtk").c_str () : - (title + ".vtk").c_str ()); - data_out.write_vtk (output_vtk); - } - - template - void PlasticityContactProblem::run () - { - pcout << "Solving problem in " << dim << " space dimensions." << std::endl; - - run_time.resize (8); - - clock_t start, end; - - start = clock(); - make_grid(); - // mesh_surface (); + MPI_Barrier (mpi_communicator); + t.stop(); + run_time[0] += t.wall_time();; - setup_system (); - - end = clock(); - run_time[0] = (double)(end-start)/CLOCKS_PER_SEC; + solve_newton (); - solve_newton (); + pcout<< "Creating output." <