From 42384154b6868676553c83fcc7feb74c85c0ba4b Mon Sep 17 00:00:00 2001 From: frohne Date: Mon, 2 Apr 2012 15:06:44 +0000 Subject: [PATCH] Clean up. git-svn-id: https://svn.dealii.org/trunk@25363 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 295 +++++++++------------------- 1 file changed, 96 insertions(+), 199 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index c4e9230ad1..021669e43f 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -426,7 +426,7 @@ namespace EquationData // return_value = 0.032 + data->dicke - input_copy->mikro_height (p(0) + shift_walze_x, p(1) + shift_walze_y, p(2)); // Ball with radius R - double R = 0.2; + double R = 0.5; if (std::pow ((p(0)-1.0/2.0), 2) + std::pow ((p(1)-1.0/2.0), 2) < R*R) return_value = 1.0 + R - 0.001 - sqrt (R*R - std::pow ((p(0)-1.0/2.0), 2) - std::pow ((p(1)-1.0/2.0), 2)); @@ -464,7 +464,7 @@ Step4::Step4 () pcout (std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), sigma_0 (400), - gamma (2000), + gamma (1.e-2), e_modul (2.e5), nu (0.3) { @@ -493,7 +493,7 @@ void Step4::make_grid () / 9 /| /______ / | 8| | 8| - | 7 | / + | 8 | / |_______|/ 6 */ @@ -504,11 +504,10 @@ void Step4::make_grid () if (cell->face (face)->center ()[2] == p2(2)) cell->face (face)->set_boundary_indicator (9); if (cell->face (face)->center ()[0] == p1(0) || - cell->face (face)->center ()[0] == p2(0)) - cell->face (face)->set_boundary_indicator (8); - if (cell->face (face)->center ()[1] == p1(1) || + cell->face (face)->center ()[0] == p2(0) || + cell->face (face)->center ()[1] == p1(1) || cell->face (face)->center ()[1] == p2(1)) - cell->face (face)->set_boundary_indicator (7); + cell->face (face)->set_boundary_indicator (8); if (cell->face (face)->center ()[2] == p1(2)) cell->face (face)->set_boundary_indicator (6); } @@ -531,20 +530,6 @@ void Step4::make_grid () { if (cell->face (face)->at_boundary() && cell->face (face)->boundary_indicator () == 9) - // Verfeinerung fuer Dortmund - // if(// Einschraenkung in x-Richtung linker Rand - // cell->face (face)->vertex (0)[0] <= 0.7 &&//p2(0)*0.5+0.4 - 0.4*hlp_refinement && - // // Einschraenkung in x-Richtung rechter Rand - // cell->face (face)->vertex (1)[0] >= 0.3 &&//p2(0)*0.5-0.3 + 0.2*hlp_refinement && - // // Einschraenkung in y-Richtung oberer Rand - // cell->face (face)->vertex (0)[1] <= 0.875 &&//p2(1)-0.6 + 0.4*hlp_refinement) - // // Einschraenkung in y-Richtung unterer Rand - // cell->face (face)->vertex (2)[1] >= 0.125) - // { - // cell->set_refine_flag (); - // break; - // } - // Verfeinerung TKSE { cell->set_refine_flag (); break; @@ -555,8 +540,6 @@ void Step4::make_grid () break; } }; - - // Markierte Zellen werden verfeinert triangulation.execute_coarsening_and_refinement (); }; } @@ -763,16 +746,10 @@ void Step4::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) // the residual part F(v) of the rhs Tensor<1,dim> rhs_values; - rhs_values = 0;//right_hand_side_values[q_point](0); - //rhs_values[1] = 0;//right_hand_side_values[q_point](1); + rhs_values = 0; cell_rhs(i) += (fe_values[displacement].value (i, q_point) * rhs_values * fe_values.JxW (q_point)); - - // if (i == 7) - // std::cout<< i << ". " << stress_tensor//->get_strain(fe_values, i, q_point) - // << ", " << stress_strain_tensor_linearized - // <::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) for (unsigned int q_point=0; q_point rhs_values; - rhs_values[0] = right_hand_side_values_face[q_point](0); - rhs_values[1] = right_hand_side_values_face[q_point](1); - rhs_values[2] = right_hand_side_values_face[q_point](2); - + rhs_values = 0; for (unsigned int i=0; i::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) system_matrix_newton.compress (); system_rhs_newton.compress (); - - // for (unsigned int i=0; i @@ -897,12 +858,11 @@ void Step4::residual_nl_system (TrilinosWrappers::MPI::Vector &u, plast_lin_hard->get_strain(fe_values, i, q_point) * fe_values.JxW (q_point)); - /* Tensor<1,dim> rhs_values; - rhs_values[0] = right_hand_side_values[q_point](0); - rhs_values[1] = right_hand_side_values[q_point](1); - cell_rhs(i) += ((fe_values[displacement].value (i, q_point) * - rhs_values) * - fe_values.JxW (q_point));*/ + Tensor<1,dim> rhs_values; + rhs_values = 0; + cell_rhs(i) += ((fe_values[displacement].value (i, q_point) * + rhs_values) * + fe_values.JxW (q_point)); }; }; @@ -919,10 +879,7 @@ void Step4::residual_nl_system (TrilinosWrappers::MPI::Vector &u, for (unsigned int q_point=0; q_point rhs_values; - rhs_values[0] = right_hand_side_values_face[q_point](0); - rhs_values[1] = right_hand_side_values_face[q_point](1); - rhs_values[2] = right_hand_side_values_face[q_point](2); - + rhs_values = 0; for (unsigned int i=0; i::residual_nl_system (TrilinosWrappers::MPI::Vector &u, constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs, local_dof_indices, system_rhs_newton); - // for (unsigned int i=0; i::projection_active_set () constraints.reinit(locally_relevant_dofs); active_set.clear (); const double c = 100.0*e_modul; + for (; cell!=endc; ++cell) if (cell->is_locally_owned()) for (unsigned int face=0; face::faces_per_cell; ++face) @@ -994,16 +950,7 @@ void Step4::projection_active_set () double obstacle_value = obstacle.value (point, 2); double solution_index_z = solution (index_z); - double gap = obstacle_value - point (2); - - // std::cout<< point << ", " - // << "Error: " << lambda (index_z) + - // diag_mass_matrix_vector (index_z)*c*(solution_index_z - gap) - // << ", " << lambda (index_z) - // << ", " << diag_mass_matrix_vector (index_z) - // << ", " << obstacle_value - // << ", " << solution_index_z - // < 0) @@ -1015,21 +962,15 @@ void Step4::projection_active_set () if (locally_owned_dofs.is_element (index_z)) active_set.add_index (index_z); - - // std::cout<< point[0] << " " - // << -lambda (index_z) << " " - // << diag_mass_matrix_vector (index_z) << " " - // << 1.0/pow (2.0,n_refinements + n_refinements_local) - // <::dirichlet_constraints () / 9 /| /______ / | 8| | 8| - | 7 | / + | 8 | / |_______|/ 6 */ - // constraints_dirichlet_hanging_nodes.clear (); constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs); constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes); std::vector component_mask (dim, true); - component_mask[0] = false; + component_mask[0] = true; component_mask[1] = true; component_mask[2] = true; VectorTools::interpolate_boundary_values (dof_handler, @@ -1074,23 +1014,14 @@ void Step4::dirichlet_constraints () constraints_dirichlet_hanging_nodes, component_mask); - component_mask[0] = false; - component_mask[1] = true; - component_mask[2] = false; - VectorTools::interpolate_boundary_values (dof_handler, - 7, - EquationData::BoundaryValues(), - constraints_dirichlet_hanging_nodes, - component_mask); - component_mask[0] = true; - component_mask[1] = false; + component_mask[1] = true; component_mask[2] = false; VectorTools::interpolate_boundary_values (dof_handler, 8, EquationData::BoundaryValues(), constraints_dirichlet_hanging_nodes, - component_mask); + component_mask); constraints_dirichlet_hanging_nodes.close (); } @@ -1102,20 +1033,13 @@ void Step4::solve () TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); distributed_solution = solution; - // constraints_hanging_nodes.set_zero (distributed_solution); - const unsigned int - start = (distributed_solution.local_range().first), - end = (distributed_solution.local_range().second); - for (unsigned int i=start; i solver (reduction_control, mpi_communicator); preconditioner_u.initialize (system_matrix_newton, additional_data); - // preconditioner_u.reinit (); solver.solve (system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u); pcout << "Initial error: " << reduction_control.initial_value() <::solve_newton () TrilinosWrappers::MPI::Vector tmp_vector (system_rhs_newton); clock_t start, end; - char* name = new char[30]; - sprintf (name,"Error-FPV.dat"); - FILE* fp = fopen(name,"w"); - delete[] name; - fclose(fp); - std::vector > constant_modes; std::vector components (dim,true); components[dim] = false; @@ -1158,10 +1076,7 @@ void Step4::solve_newton () additional_data.output_details = false; additional_data.smoother_sweeps = 2; additional_data.aggregation_threshold = 1e-2; - //additional_data.constant_modes; - // std::cout<< "Update Active Set in Dim = " << dim < sigma_eff_vector; sigma_eff_vector.reinit (triangulation.n_active_cells()); @@ -1170,103 +1085,85 @@ void Step4::solve_newton () for (; j<=100;j++) { pcout<< " " < " << resid <