]> https://gitweb.dealii.org/ - dealii.git/commitdiff
change damped from int to bool; some documentation
authorJoerg Frohne <frohne@mathematik.uni-siegen.de>
Mon, 1 Apr 2013 22:45:30 +0000 (22:45 +0000)
committerJoerg Frohne <frohne@mathematik.uni-siegen.de>
Mon, 1 Apr 2013 22:45:30 +0000 (22:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@29134 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-42/step-42.cc

index d6c50ee10b25086ecc3e50b54e637e6191ce7687..cdd7085f1351f457d6aef147ba58efb935de6ed7 100644 (file)
@@ -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<double>  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 <int dim>
   double Input<dim>::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 <int dim>
   double Input<dim>::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 <int dim>
   void Input<dim>::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 <int dim>
   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 <int dim>
   ConstitutiveLaw<dim>::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<dim>() - outer_product(unit_symmetric_tensor<dim>(), unit_symmetric_tensor<dim>())/3.0);
   }
 
+  // Calculates the strain for the shape functions.
   template <int dim>
   inline
   SymmetricTensor<2,dim> ConstitutiveLaw<dim>::get_strain (const FEValues<dim> &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 <int dim>
   void ConstitutiveLaw<dim>::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 <int dim>
   void ConstitutiveLaw<dim>::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 <int dim>
     class RightHandSide : public Function<dim>
     {
@@ -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<dim; ++i)
-      //   return_value += 4*std::pow(p(i), 4);
 
       return return_value;
     }
@@ -417,7 +436,8 @@ namespace Step42
         values(c) = RightHandSide<dim>::value (p, c);
     }
 
-
+    // This function class is used to describe the prescribed displacements
+    // at the boundary. But again we set this to zero.
     template <int dim>
     class BoundaryValues : public Function<dim>
     {
@@ -455,7 +475,11 @@ namespace Step42
         values(c) = BoundaryValues<dim>::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 <int dim>
     class Obstacle : public Function<dim>
     {
@@ -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 <int dim>
   void PlasticityContactProblem<dim>::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<double>(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<resid_old)
-              damped=1;
+              damped=true;
 
             computing_timer.exit_section("Residual and lambda");
 
@@ -1363,7 +1398,7 @@ namespace Step42
         if (num_changed==0 && resid < 1e-8)
                    break;
         active_set_old = active_set;
-      } // End of active-set-loop
+      }
 
     pcout << "" << std::endl
           << "      Number of assembled systems = " << number_assemble_system

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.