From 86dae9a8e4b375e143a0feb69f47a1385031ea75 Mon Sep 17 00:00:00 2001 From: frohne Date: Tue, 21 Aug 2012 17:46:45 +0000 Subject: [PATCH] Cleaning the code; some documentation git-svn-id: https://svn.dealii.org/trunk@26059 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 283 ++++++++++------------------ 1 file changed, 97 insertions(+), 186 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 281ed94ac9..3a9a898ef8 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -12,11 +12,8 @@ // @sect3{Include files} - // The first few (many?) include - // files have already been used in - // the previous example, so we will - // not explain their meaning here - // again. + // We are using the the same + // include files as in step-42: #include #include @@ -72,8 +69,31 @@ namespace Step42 // @sect3{The PlasticityContactProblem class template} + // This class provides an interface + // for a constitutive law. In this + // example we are using an elastic + // plastic material with linear, + // isotropic hardening. + template class ConstitutiveLaw; + // @sect3{The PlasticityContactProblem class template} + + // This class supplies all function + // and variables needed to describe + // the nonlinear contact problem. It is + // close to step-41 but with some additional + // features like: handling hanging nodes, + // a newton method, using Trilinos and p4est + // for parallel distributed computing. + // To deal with hanging nodes makes + // life a bit more complicated since + // we need an other ConstraintMatrix now. + // We create a newton method for the + // active set method for the contact + // situation and to handle the nonlinear + // operator for the constitutive law. + template class PlasticityContactProblem { @@ -84,10 +104,9 @@ namespace Step42 private: void make_grid (); void setup_system(); - void assemble_mass_matrix (); void assemble_nl_system (TrilinosWrappers::MPI::Vector &u); - void residual_nl_system (TrilinosWrappers::MPI::Vector &u, - Vector &sigma_eff_vector); + void residual_nl_system (TrilinosWrappers::MPI::Vector &u); + void assemble_mass_matrix_diagonal (TrilinosWrappers::SparseMatrix &mass_matrix); void update_solution_and_constraints (); void dirichlet_constraints (); void solve (); @@ -95,7 +114,6 @@ namespace Step42 void output_results (const std::string& title) const; 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 (Vector vector, const std::string& title) const; int n_refinements_global; int n_refinements_local; @@ -118,7 +136,6 @@ namespace Step42 ConstraintMatrix constraints_dirichlet_hanging_nodes; TrilinosWrappers::SparseMatrix system_matrix_newton; - TrilinosWrappers::SparseMatrix mass_matrix; TrilinosWrappers::MPI::Vector solution; TrilinosWrappers::MPI::Vector old_solution; @@ -139,21 +156,23 @@ namespace Step42 double gamma; // Parameter for the linear isotropic hardening double e_modul; // E-Modul double nu; // Poisson ratio - - std_cxx1x::shared_ptr Mp_preconditioner; }; template class ConstitutiveLaw { public: - ConstitutiveLaw (double _E, double _nu, double _sigma_0, double _gamma, MPI_Comm _mpi_communicator, ConditionalOStream _pcout); - // ConstitutiveLaw (double mu, double kappa); + ConstitutiveLaw (double _E, + double _nu, + double _sigma_0, + double _gamma, + MPI_Comm _mpi_communicator, + ConditionalOStream _pcout); + void plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor, SymmetricTensor<2,dim> &strain_tensor, - unsigned int &elast_points, - unsigned int &plast_points, - double &sigma_eff, + unsigned int &elast_points, + unsigned int &plast_points, double &yield); void linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, SymmetricTensor<4,dim> &stress_strain_tensor, @@ -210,7 +229,6 @@ namespace Step42 SymmetricTensor<2,dim> &strain_tensor, unsigned int &elast_points, unsigned int &plast_points, - double &sigma_eff, double &yield) { if (dim == 3) @@ -236,10 +254,7 @@ namespace Step42 else elast_points += 1; - // std::cout<< beta < @@ -563,7 +576,6 @@ 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); @@ -604,65 +616,18 @@ namespace Step42 system_matrix_newton.reinit (sp); + TrilinosWrappers::SparseMatrix mass_matrix; mass_matrix.reinit (sp); + assemble_mass_matrix_diagonal (mass_matrix); + const unsigned int + start = (system_rhs_newton.local_range().first), + end = (system_rhs_newton.local_range().second); + for (unsigned int j=start; j - void PlasticityContactProblem::assemble_mass_matrix () - { - QTrapez 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; - const unsigned int dofs_per_face = fe.dofs_per_face; - const unsigned int n_face_q_points = face_quadrature_formula.size(); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - const FEValuesExtractors::Vector displacement (0); - - 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) - if (cell->face (face)->at_boundary() - && cell->face (face)->boundary_indicator () == 9) - { - 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); - } - - mass_matrix.compress (); } template @@ -792,8 +757,7 @@ namespace Step42 } template - void PlasticityContactProblem::residual_nl_system (TrilinosWrappers::MPI::Vector &u, - Vector &sigma_eff_vector) + void PlasticityContactProblem::residual_nl_system (TrilinosWrappers::MPI::Vector &u) { QGauss quadrature_formula(2); QGauss face_quadrature_formula(2); @@ -819,7 +783,6 @@ namespace Step42 Vector(dim)); Vector cell_rhs (dofs_per_cell); - Vector cell_sigma_eff (dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); @@ -830,7 +793,6 @@ namespace Step42 unsigned int elast_points = 0; unsigned int plast_points = 0; - double sigma_eff = 0; double yield = 0; unsigned int cell_number = 0; for (; cell!=endc; ++cell) @@ -851,13 +813,8 @@ namespace Step42 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; + elast_points, plast_points, yield); - /* if (q_point == 0) - std::cout<< stress_strain_tensor < - void PlasticityContactProblem::output_results (Vector vector, const std::string& title) const - { - DataOut data_out; - - 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 () { @@ -1402,79 +1386,6 @@ namespace Step42 // @sect3{The main function} - // And this is the main function. It also - // looks mostly like in step-3, but if you - // look at the code below, note how we first - // create a variable of type - // PlasticityContactProblem@<2@> (forcing - // the compiler to compile the class template - // with dim replaced by - // 2) and run a 2d simulation, - // and then we do the whole thing over in 3d. - // - // In practice, this is probably not what you - // would do very frequently (you probably - // either want to solve a 2d problem, or one - // in 3d, but not both at the same - // time). However, it demonstrates the - // mechanism by which we can simply change - // which dimension we want in a single place, - // and thereby force the compiler to - // recompile the dimension independent class - // templates for the dimension we - // request. The emphasis here lies on the - // fact that we only need to change a single - // place. This makes it rather trivial to - // debug the program in 2d where computations - // are fast, and then switch a single place - // to a 3 to run the much more computing - // intensive program in 3d for `real' - // computations. - // - // Each of the two blocks is enclosed in - // braces to make sure that the - // laplace_problem_2d variable - // goes out of scope (and releases the memory - // it holds) before we move on to allocate - // memory for the 3d case. Without the - // additional braces, the - // laplace_problem_2d variable - // would only be destroyed at the end of the - // function, i.e. after running the 3d - // problem, and would needlessly hog memory - // while the 3d run could actually use it. - // - // Finally, the first line of the function is - // used to suppress some output. Remember - // that in the previous example, we had the - // output from the linear solvers about the - // starting residual and the number of the - // iteration where convergence was - // detected. This can be suppressed through - // the deallog.depth_console(0) - // call. - // - // The rationale here is the following: the - // deallog (i.e. deal-log, not de-allog) - // variable represents a stream to which some - // parts of the library write output. It - // redirects this output to the console and - // if required to a file. The output is - // nested in a way so that each function can - // use a prefix string (separated by colons) - // for each line of output; if it calls - // another function, that may also use its - // prefix which is then printed after the one - // of the calling function. Since output from - // functions which are nested deep below is - // usually not as important as top-level - // output, you can give the deallog variable - // a maximal depth of nested output for - // output to console and file. The depth zero - // which we gave here means that no output is - // written. By changing it you can get more - // information about the innards of the - // library. int main (int argc, char *argv[]) { using namespace dealii; -- 2.39.5