From: bangerth Date: Mon, 2 Apr 2012 15:23:16 +0000 (+0000) Subject: Fix which indices are allowed into the active set. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6027fed23a057e6ee7b14b80f02af3795b87803e;p=dealii-svn.git Fix which indices are allowed into the active set. git-svn-id: https://svn.dealii.org/trunk@25364 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 021669e43f..0b68743923 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -82,12 +82,12 @@ using namespace dealii; template class ConstitutiveLaw; template -class Step4 +class Step4 { public: Step4 (); void run (); - + private: void make_grid (); void setup_system(); @@ -104,48 +104,48 @@ private: void output_results (TrilinosWrappers::MPI::Vector vector, const std::string& title) const; void output_results (Vector vector, const std::string& title) const; - MPI_Comm mpi_communicator; + MPI_Comm mpi_communicator; parallel::distributed::Triangulation triangulation; FESystem fe; DoFHandler dof_handler; - + IndexSet locally_owned_dofs; IndexSet locally_relevant_dofs; - + int n_refinements; int n_refinements_local; unsigned int number_iterations; std::vector run_time; - + ConstraintMatrix constraints; ConstraintMatrix constraints_hanging_nodes; ConstraintMatrix constraints_dirichlet_hanging_nodes; - + TrilinosWrappers::SparseMatrix system_matrix_newton; TrilinosWrappers::SparseMatrix mass_matrix; - + TrilinosWrappers::MPI::Vector solution; TrilinosWrappers::MPI::Vector old_solution; TrilinosWrappers::MPI::Vector system_rhs_newton; TrilinosWrappers::MPI::Vector resid_vector; TrilinosWrappers::MPI::Vector diag_mass_matrix_vector; - IndexSet active_set; + IndexSet active_set; ConditionalOStream pcout; - + TrilinosWrappers::PreconditionAMG::AdditionalData additional_data; TrilinosWrappers::PreconditionAMG preconditioner_u; TrilinosWrappers::PreconditionAMG preconditioner_t; - + std::auto_ptr > plast_lin_hard; - + double sigma_0; // Yield stress double gamma; // Parameter for the linear isotropic hardening double e_modul; // E-Modul double nu; // Poisson ratio - + std_cxx1x::shared_ptr Mp_preconditioner; }; @@ -166,8 +166,8 @@ public: SymmetricTensor<2,dim> &strain_tensor); inline SymmetricTensor<2,dim> get_strain (const FEValues &fe_values, const unsigned int shape_func, - const unsigned int q_point) const; - + const unsigned int q_point) const; + private: SymmetricTensor<4,dim> stress_strain_tensor_mu; SymmetricTensor<4,dim> stress_strain_tensor_kappa; @@ -207,7 +207,7 @@ SymmetricTensor<2,dim> ConstitutiveLaw::get_strain (const FEValues &fe SymmetricTensor<2,dim> tmp; tmp = fe_values[displacement].symmetric_gradient (shape_func,q_point); - + return tmp; } @@ -226,13 +226,13 @@ void ConstitutiveLaw::plast_linear_hardening (SymmetricTensor<4,dim> &stre stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; double tmp = E/((1+nu)*(1-2*nu)); double stress_tensor_33 = 0.0;//tmp*(strain_tensor[0][0] + strain_tensor[1][1])*nu; - + SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); - + double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); deviator_stress_tensor_norm = std::sqrt (deviator_stress_tensor_norm*deviator_stress_tensor_norm + stress_tensor_33*stress_tensor_33); - + yield = 0; stress_strain_tensor = stress_strain_tensor_mu; double beta = 1.0; @@ -245,7 +245,7 @@ void ConstitutiveLaw::plast_linear_hardening (SymmetricTensor<4,dim> &stre } else elast_points += 1; - + // std::cout<< beta <::linearized_plast_linear_hardening (SymmetricTensor<4, stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; double tmp = E/((1+nu)*(1-2*nu)); double stress_tensor_33 = 0.0;//tmp*(strain_tensor[0][0] + strain_tensor[1][1])*nu; - + stress_strain_tensor = stress_strain_tensor_mu; stress_strain_tensor_linearized = stress_strain_tensor_mu; - + SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); - + double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); deviator_stress_tensor_norm = std::sqrt (deviator_stress_tensor_norm*deviator_stress_tensor_norm + stress_tensor_33*stress_tensor_33); double beta = 1.0; if (deviator_stress_tensor_norm >= sigma_0) - { + { beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; stress_strain_tensor *= beta; stress_strain_tensor_linearized *= beta; deviator_stress_tensor /= deviator_stress_tensor_norm; - stress_strain_tensor_linearized -= beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor); + stress_strain_tensor_linearized -= beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor); } - + stress_strain_tensor += stress_strain_tensor_kappa; stress_strain_tensor_linearized += stress_strain_tensor_kappa; } @@ -291,21 +291,21 @@ void ConstitutiveLaw::linearized_plast_linear_hardening (SymmetricTensor<4, namespace EquationData { template - class RightHandSide : public Function + class RightHandSide : public Function { public: RightHandSide () : Function(dim) {} - + virtual double value (const Point &p, const unsigned int component = 0) const; - + virtual void vector_value (const Point &p, Vector &values) const; }; template double RightHandSide::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { double return_value = 0.0; @@ -320,10 +320,10 @@ namespace EquationData return_value = 0.0; // for (unsigned int i=0; i void RightHandSide::vector_value (const Point &p, Vector &values) const @@ -334,21 +334,21 @@ namespace EquationData template - class BoundaryValues : public Function + class BoundaryValues : public Function { public: BoundaryValues () : Function(dim) {}; - + virtual double value (const Point &p, const unsigned int component = 0) const; - + virtual void vector_value (const Point &p, Vector &values) const; }; - + template double BoundaryValues::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { double return_value = 0; @@ -364,15 +364,15 @@ namespace EquationData template void BoundaryValues::vector_value (const Point &p, - Vector &values) const + Vector &values) const { for (unsigned int c=0; cn_components; ++c) values(c) = BoundaryValues::value (p, c); } - + template - class Obstacle : public Function + class Obstacle : public Function { public: Obstacle () : Function(dim) {}; @@ -386,7 +386,7 @@ namespace EquationData template double Obstacle::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { double R = 0.03; double return_value = 0.0; @@ -434,16 +434,16 @@ namespace EquationData return_value = 1e+5; } return return_value; - + // return 1e+10;//0.98; } - + template void Obstacle::vector_value (const Point &p, - Vector &values) const + Vector &values) const { for (unsigned int c=0; cn_components; ++c) - values(c) = Obstacle::value (p, c); + values(c) = Obstacle::value (p, c); } } @@ -453,7 +453,7 @@ namespace EquationData // Next for the implementation of the class // template that makes use of the functions // above. As before, we will write everything - + template Step4::Step4 () : @@ -483,11 +483,11 @@ void Step4::make_grid () Point p1 (0,0,0); Point p2 (1.0, 1.0, 1.0); GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2); - + Triangulation<3>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); - + /* boundary_indicators: _______ / 9 /| @@ -510,9 +510,9 @@ void Step4::make_grid () cell->face (face)->set_boundary_indicator (8); if (cell->face (face)->center ()[2] == p1(2)) cell->face (face)->set_boundary_indicator (6); - } - - n_refinements = 2; + } + + n_refinements = 3; n_refinements_local = 3; triangulation.refine_global (n_refinements); @@ -520,7 +520,7 @@ void Step4::make_grid () for (int step=0; step::setup_system () DoFTools::make_hanging_node_constraints (dof_handler, constraints_hanging_nodes); constraints_hanging_nodes.close (); - + pcout << "Number of active cells: " << triangulation.n_active_cells() << std::endl @@ -574,7 +574,7 @@ void Step4::setup_system () << "Number of degrees of freedom: " << dof_handler.n_dofs () << std::endl; - + dirichlet_constraints (); } @@ -595,12 +595,12 @@ void Step4::setup_system () DoFTools::make_sparsity_pattern (dof_handler, sp, constraints_dirichlet_hanging_nodes, false, Utilities::MPI::this_mpi_process(mpi_communicator)); - + sp.compress(); - + system_matrix_newton.reinit (sp); - - mass_matrix.reinit (sp); + + mass_matrix.reinit (sp); } assemble_mass_matrix (); @@ -613,11 +613,11 @@ void Step4::setup_system () } template -void Step4::assemble_mass_matrix () -{ +void Step4::assemble_mass_matrix () +{ QTrapez face_quadrature_formula; - FEFaceValues fe_values_face (fe, face_quadrature_formula, + FEFaceValues fe_values_face (fe, face_quadrature_formula, update_values | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; @@ -633,7 +633,7 @@ void Step4::assemble_mass_matrix () typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); - + for (; cell!=endc; ++cell) if (cell->is_locally_owned()) for (unsigned int face=0; face::faces_per_cell; ++face) @@ -642,15 +642,15 @@ void Step4::assemble_mass_matrix () { fe_values_face.reinit (cell, face); cell_matrix = 0; - + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - + constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_matrix, local_dof_indices, mass_matrix); @@ -660,18 +660,18 @@ void Step4::assemble_mass_matrix () } template -void Step4::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) +void Step4::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) { QGauss quadrature_formula(2); QGauss face_quadrature_formula(2); - FEValues fe_values (fe, quadrature_formula, + FEValues fe_values (fe, quadrature_formula, UpdateFlags(update_values | update_gradients | update_q_points | update_JxW_values)); - FEFaceValues fe_values_face (fe, face_quadrature_formula, + FEFaceValues fe_values_face (fe, face_quadrature_formula, update_values | update_quadrature_points | update_JxW_values); @@ -685,7 +685,7 @@ void Step4::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) std::vector > right_hand_side_values_face (n_face_q_points, Vector(dim)); - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); @@ -703,47 +703,47 @@ void Step4::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) fe_values.reinit (cell); cell_matrix = 0; cell_rhs = 0; - + right_hand_side.vector_value_list (fe_values.get_quadrature_points(), right_hand_side_values); - + std::vector > strain_tensor (n_q_points); fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor); - + for (unsigned int q_point=0; q_point stress_strain_tensor_linearized; SymmetricTensor<4,dim> stress_strain_tensor; SymmetricTensor<2,dim> stress_tensor; - + plast_lin_hard->linearized_plast_linear_hardening (stress_strain_tensor_linearized, stress_strain_tensor, strain_tensor[q_point]); - + // if (q_point == 0) // std::cout<< stress_strain_tensor_linearized <get_strain(fe_values, i, q_point); - + for (unsigned int j=0; jget_strain(fe_values, j, q_point) * fe_values.JxW (q_point)); } - + // the linearized part a(v^i;v^i,v) of the rhs cell_rhs(i) += (stress_tensor * strain_tensor[q_point] * fe_values.JxW (q_point)); - + // the residual part a(v^i;v) of the rhs cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * plast_lin_hard->get_strain(fe_values, i, q_point) * fe_values.JxW (q_point)); - + // the residual part F(v) of the rhs Tensor<1,dim> rhs_values; rhs_values = 0; @@ -752,21 +752,21 @@ void Step4::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) fe_values.JxW (q_point)); } } - + for (unsigned int face=0; face::faces_per_cell; ++face) { if (cell->face (face)->at_boundary() && cell->face (face)->boundary_indicator () == 9) { fe_values_face.reinit (cell, face); - + right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(), right_hand_side_values_face); - + for (unsigned int q_point=0; q_point rhs_values; - rhs_values = 0; + rhs_values = 0; for (unsigned int i=0; i::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) } } } - + cell->get_dof_indices (local_dof_indices); constraints.distribute_local_to_global (cell_matrix, cell_rhs, local_dof_indices, system_matrix_newton, system_rhs_newton, true); }; - + system_matrix_newton.compress (); system_rhs_newton.compress (); } @@ -792,13 +792,13 @@ void Step4::residual_nl_system (TrilinosWrappers::MPI::Vector &u, QGauss quadrature_formula(2); QGauss face_quadrature_formula(2); - FEValues fe_values (fe, quadrature_formula, + FEValues fe_values (fe, quadrature_formula, UpdateFlags(update_values | update_gradients | update_q_points | update_JxW_values)); - FEFaceValues fe_values_face (fe, face_quadrature_formula, + FEFaceValues fe_values_face (fe, face_quadrature_formula, update_values | update_quadrature_points | update_JxW_values); @@ -808,10 +808,10 @@ void Step4::residual_nl_system (TrilinosWrappers::MPI::Vector &u, const EquationData::RightHandSide right_hand_side; std::vector > right_hand_side_values (n_q_points, - Vector(dim)); + Vector(dim)); std::vector > right_hand_side_values_face (n_face_q_points, Vector(dim)); - + Vector cell_rhs (dofs_per_cell); Vector cell_sigma_eff (dofs_per_cell); @@ -832,24 +832,24 @@ void Step4::residual_nl_system (TrilinosWrappers::MPI::Vector &u, { fe_values.reinit (cell); cell_rhs = 0; - + right_hand_side.vector_value_list (fe_values.get_quadrature_points(), right_hand_side_values); - + std::vector > strain_tensor (n_q_points); fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor); - + for (unsigned int q_point=0; q_point stress_strain_tensor; SymmetricTensor<2,dim> stress_tensor; - + plast_lin_hard->plast_linear_hardening (stress_strain_tensor, strain_tensor[q_point], elast_points, plast_points, sigma_eff, yield); - + // sigma_eff_vector (cell_number) += sigma_eff; sigma_eff_vector (cell_number) += yield; - + /* if (q_point == 0) std::cout<< stress_strain_tensor <::residual_nl_system (TrilinosWrappers::MPI::Vector &u, cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * //(stress_tensor) * plast_lin_hard->get_strain(fe_values, i, q_point) * 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)); + fe_values.JxW (q_point)); }; }; - + for (unsigned int face=0; face::faces_per_cell; ++face) { if (cell->face (face)->at_boundary() && cell->face (face)->boundary_indicator () == 9) { fe_values_face.reinit (cell, face); - + right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(), right_hand_side_values_face); - + for (unsigned int q_point=0; q_point rhs_values; - rhs_values = 0; + rhs_values = 0; for (unsigned int i=0; i::residual_nl_system (TrilinosWrappers::MPI::Vector &u, } } } - + cell->get_dof_indices (local_dof_indices); constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs, local_dof_indices, @@ -918,7 +918,7 @@ void Step4::projection_active_set () typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); - + TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); distributed_solution = solution; TrilinosWrappers::MPI::Vector lambda (solution); @@ -937,30 +937,30 @@ void Step4::projection_active_set () for (unsigned int v=0; v::vertices_per_cell; ++v) { unsigned int index_z = cell->face (face)->vertex_dof_index (v,2); - + if (vertex_touched[cell->face (face)->vertex_index(v)] == false) vertex_touched[cell->face (face)->vertex_index(v)] = true; - else + else continue; - + // the local row where Point point (cell->face (face)->vertex (v)[0],/* + solution (index_x),*/ cell->face (face)->vertex (v)[1], cell->face (face)->vertex (v)[2]); - + double obstacle_value = obstacle.value (point, 2); double solution_index_z = solution (index_z); double gap = obstacle_value - point (2); - + if (lambda (index_z) + c*diag_mass_matrix_vector_relevant (index_z)*(solution_index_z - gap) > 0) { constraints.add_line (index_z); constraints.set_inhomogeneity (index_z, gap); - - distributed_solution (index_z) = gap; - - if (locally_owned_dofs.is_element (index_z)) + + distributed_solution (index_z) = gap; + + if (locally_relevant_dofs.is_element (index_z)) active_set.add_index (index_z); // std::cout<< index_z << ", " @@ -973,7 +973,7 @@ void Step4::projection_active_set () // <::projection_active_set () solution = distributed_solution; constraints.close (); - + const ConstraintMatrix::MergeConflictBehavior merge_conflict_behavior = ConstraintMatrix::left_object_wins; constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior); @@ -1013,7 +1013,7 @@ void Step4::dirichlet_constraints () EquationData::BoundaryValues(), constraints_dirichlet_hanging_nodes, component_mask); - + component_mask[0] = true; component_mask[1] = true; component_mask[2] = false; @@ -1026,7 +1026,7 @@ void Step4::dirichlet_constraints () } template -void Step4::solve () +void Step4::solve () { ReductionControl reduction_control (10000, 1e-15, 1e-4); @@ -1034,7 +1034,7 @@ void Step4::solve () distributed_solution = solution; constraints_hanging_nodes.set_zero (distributed_solution); - + // Solving iterative SolverCG solver (reduction_control, mpi_communicator); @@ -1096,17 +1096,17 @@ void Step4::solve_newton () assemble_nl_system (solution); //compute Newton-Matrix end = clock(); run_time[1] += (double)(end-start)/CLOCKS_PER_SEC; - + number_assemble_system += 1; - + start = clock(); solve (); end = clock(); run_time[2] += (double)(end-start)/CLOCKS_PER_SEC; - + TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); distributed_solution = solution; - + int damped = 0; tmp_vector = old_solution; double a = 0; @@ -1115,7 +1115,7 @@ void Step4::solve_newton () a=pow(0.5,i); old_solution = tmp_vector; old_solution.sadd(1-a,a, distributed_solution); - + start = clock(); system_rhs_newton = 0; sigma_eff_vector = 0; @@ -1136,10 +1136,10 @@ void Step4::solve_newton () // <::output_results (const std::string& title) const lambda = resid_vector; DataOut data_out; - + data_out.attach_dof_handler (dof_handler); - + data_out.add_data_vector (solution, "Displacement"); data_out.add_data_vector (lambda, "Residual"); data_out.add_data_vector (active_set, "ActiveSet"); - + Vector subdomain (triangulation.n_active_cells()); for (unsigned int i=0; i::output_results (const std::string& title) const filenames.push_back ("solution-" + Utilities::int_to_string (i, 4) + ".vtu"); - + std::ofstream master_output ((filename + ".pvtu").c_str()); data_out.write_pvtu_record (master_output, filenames); } - + TrilinosWrappers::MPI::Vector tmp (solution); tmp *= -1; move_mesh (tmp); @@ -1250,7 +1250,7 @@ void Step4::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displ if (vertex_touched[cell->vertex_index(v)] == false) { vertex_touched[cell->vertex_index(v)] = true; - + Point vertex_displacement; for (unsigned int d=0; d::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displ vertex_displacement[d] = _complete_displacement(cell->vertex_dof_index(v,d)); } - + cell->vertex(v) += vertex_displacement; } - } + } } template @@ -1298,7 +1298,7 @@ void Step4::output_results (Vector vector, const std::string& title } template -void Step4::run () +void Step4::run () { pcout << "Solving problem in " << dim << " space dimensions." << std::endl; @@ -1392,7 +1392,7 @@ void Step4::run () // written. By changing it you can get more // information about the innards of the // library. -int main (int argc, char *argv[]) +int main (int argc, char *argv[]) { deallog.depth_console (0); @@ -1408,6 +1408,6 @@ int main (int argc, char *argv[]) end = clock(); cout<< "%%%%%% Rechenzeit overall = " << (double)(end-start)/CLOCKS_PER_SEC <