From 365f4abcf5d13b11e4f8171c607c2ecf6f1b314f Mon Sep 17 00:00:00 2001 From: frohne Date: Mon, 1 Apr 2013 22:45:30 +0000 Subject: [PATCH] change damped from int to bool; some documentation git-svn-id: https://svn.dealii.org/trunk@29134 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 79 +++++++++++++++++++++-------- 1 file changed, 57 insertions(+), 22 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index d6c50ee10b..cdd7085f13 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -73,7 +73,7 @@ namespace Step42 // This class has the the only purpose // to read in data from a picture file - // that has to be stored as a pbm ascii + // that has to be stored in pbm ascii // format. This data will be bilinear // interpolated and provides in this way // a function which describes an obstacle. @@ -106,9 +106,7 @@ namespace Step42 pcout (std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), obstacle_data (0), - lx (1.0), // length of the cube in x direction - ly (1.0), // length of the cube in y direction - hx (0), // + hx (0), hy (0), nx (0), ny (0) @@ -125,11 +123,12 @@ namespace Step42 MPI_Comm mpi_communicator; ConditionalOStream pcout; std::vector obstacle_data; - double lx, ly; double hx, hy; int nx, ny; }; + // This function is used in obstacle_function () + // to provide the proper value of the obstacle. template double Input::hv (int i, int j) { @@ -138,6 +137,8 @@ namespace Step42 return obstacle_data[nx*(ny-1-j)+i]; // i indiziert x-werte, j indiziert y-werte } + // obstacle_function () calculates the bilinear interpolated + // value in the point (x,y). template double Input::obstacle_function (double x,double y) { @@ -206,6 +207,9 @@ namespace Step42 return val; } + // As mentioned above this function reads in the + // obstacle datas and stores them in the std::vector + // obstacle_data. It will be used only in run (). template void Input::read_obstacle (const char* name) { @@ -239,10 +243,11 @@ namespace Step42 // This class provides an interface // for a constitutive law. In this - // example we are using an elastic - // plastic material with linear, + // example we are using an elasto + // plastic material behavior with linear, // isotropic hardening. - + // For gamma = 0 we obtain perfect elasto + // plasticity behavior. template class ConstitutiveLaw { @@ -280,6 +285,15 @@ namespace Step42 ConditionalOStream pcout; }; + // The constructor of the ConstitutiveLaw class sets the + // required material parameter for our deformable body: + // E -> elastic modulus + // nu -> Passion's number + // sigma_0 -> yield stress + // gamma -> hardening parameter. + // Also it supplies the stress strain tensor of forth order + // of the volumetric and deviator part. For further details + // see the documentation above. template ConstitutiveLaw::ConstitutiveLaw(double _E, double _nu, double _sigma_0, double _gamma, MPI_Comm _mpi_communicator, ConditionalOStream _pcout) :E (_E), @@ -295,6 +309,7 @@ namespace Step42 stress_strain_tensor_mu = 2*mu*(identity_tensor() - outer_product(unit_symmetric_tensor(), unit_symmetric_tensor())/3.0); } + // Calculates the strain for the shape functions. template inline SymmetricTensor<2,dim> ConstitutiveLaw::get_strain (const FEValues &fe_values, @@ -309,6 +324,11 @@ namespace Step42 return tmp; } + // This is the implemented constitutive law. It projects the + // deviator part of the stresses in a quadrature point back to + // the yield stress plus the linear isotropic hardening. + // Also we sum up the elastic and the plastic quadrature + // points. template void ConstitutiveLaw::plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor, const SymmetricTensor<2,dim> &strain_tensor, @@ -342,6 +362,8 @@ namespace Step42 } } + // This function returns the linearized stress strain tensor. + // It contains the derivative of the nonlinear constitutive law. template void ConstitutiveLaw::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, SymmetricTensor<4,dim> &stress_strain_tensor, @@ -375,6 +397,8 @@ namespace Step42 namespace EquationData { + // It possible to apply an additional body force + // but in here it is set to zero. template class RightHandSide : public Function { @@ -399,12 +423,7 @@ namespace Step42 if (component == 1) return_value = 0.0; if (component == 2) - // if ((p(0)-0.5)*(p(0)-0.5)+(p(1)-0.5)*(p(1)-0.5) < 0.2) - // return_value = -5000; - // else return_value = 0.0; - // for (unsigned int i=0; i::value (p, c); } - + // This function class is used to describe the prescribed displacements + // at the boundary. But again we set this to zero. template class BoundaryValues : public Function { @@ -455,7 +475,11 @@ namespace Step42 values(c) = BoundaryValues::value (p, c); } - + // This function is obviously implemented to + // define the obstacle that penetrates our deformable + // body. You can choose between two ways to define + // your obstacle: to read it from a file or to use + // a function (here a ball). template class Obstacle : public Function { @@ -1250,8 +1274,13 @@ namespace Step42 computing_timer.exit_section("Solve"); } + // @sect4{PlasticityContactProblem::solve_newton} - + // In this function the damped Newton method is implemented. + // That means two nested loops: the outer loop for the newton + // iteration and the inner loop for the damping steps which + // will be used only if necessary. To obtain a good and reasonable + // starting value we solve an elastic problem in very first step (j=1). template void PlasticityContactProblem::solve_newton () { @@ -1280,7 +1309,6 @@ namespace Step42 unsigned int number_assemble_system = 0; for (; j<=100; j++) { - // Solve an elastic problem to obtain a better start solution if (j == 1 && cycle == 0) plast_lin_hard->set_sigma_0 (1e+10); else if (j == 2 || cycle > 0) @@ -1305,16 +1333,23 @@ namespace Step42 TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); distributed_solution = solution; - int damped = 0; + + // We handle a highly nonlinear problem so we have to damp + // the Newtons method. We refer that we iterate the new solution + // in each Newton step and not only the solution update. + // Since the solution set is a convex set and not a space we + // compute for the damping a linear combination of the + // previous and the current solution to guarantee that the + // damped solution is in our solution set again. + // At most we apply 10 damping steps. + bool damped = false; tmp_vector = old_solution; -// constraints.distribute (old_solution); double a = 0; for (unsigned int i=0; (i<10)&&(!damped); i++) { a=std::pow(0.5, static_cast(i)); old_solution = tmp_vector; old_solution.sadd(1-a,a, distributed_solution); -// constraints.distribute (old_solution); old_solution.compress (VectorOperation::add); computing_timer.enter_section("Residual and lambda"); @@ -1337,7 +1372,7 @@ namespace Step42 resid = res.l2_norm (); if (resid