ElastoplasticTorsion.h

Go to the documentation of this file.

 
 *   
 *   
 *   /*
 *    This piece of software solves the elliptic p-laplacian
 *    boundary-value problems:
 *   
 *      Min {∫ 1/2 W(|Du|²)+ 1/p |Du|^p -fu : u=g on ∂S }   (1)
 *       u
 *   
 *    for large values of p, which approximates (see Alvarez & Flores 2015)
 *   
 *      Min {∫ 1/2 W(|Du|²) -fu : |Du|<1 a.s. on S, u=g on ∂S }
 *       u
 *   
 *    By default W(t)=t and S=unit disk.
 *   
 *    Large portions of this code are borrowed from the deal.ii tutorials
 *   
 *         @ref step_15 "step-15" @ref step_29 "step-29".
 *   
 *    For further details see the technical report
 *    "Solving variational problems with uniform gradient bounds by p-Laplacian
 *     approximation: Elastoplastic torsion implementation using the deal.II
 *     library"
 *    available at the documentation and at http://www.dim.uchile.cl/~sflores.
 *   
 *   */
 *   
 * @endcode
 * 
 * Include files
 * 
 
 * 
 * 
 * @code
 *   #include <deal.II/base/quadrature_lib.h>
 *   #include <deal.II/base/function.h>
 *   #include <deal.II/base/logstream.h>
 *   #include <deal.II/base/utilities.h>
 *   #include <deal.II/base/convergence_table.h>
 *   #include <deal.II/base/smartpointer.h>
 *   #include <deal.II/base/parameter_handler.h>
 *   #include <deal.II/base/timer.h>
 *   
 *   #include <deal.II/lac/vector.h>
 *   #include <deal.II/lac/full_matrix.h>
 *   #include <deal.II/lac/sparse_matrix.h>
 *   #include <deal.II/lac/solver_cg.h>
 *   #include <deal.II/lac/precondition.h>
 *   #include <deal.II/lac/affine_constraints.h>
 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
 *   
 *   #include <deal.II/grid/tria.h>
 *   #include <deal.II/grid/grid_generator.h>
 *   #include <deal.II/grid/tria_accessor.h>
 *   #include <deal.II/grid/tria_iterator.h>
 *   #include <deal.II/grid/manifold_lib.h>
 *   #include <deal.II/grid/grid_refinement.h>
 *   #include <deal.II/grid/grid_in.h>
 *   
 *   #include <deal.II/dofs/dof_handler.h>
 *   #include <deal.II/dofs/dof_accessor.h>
 *   #include <deal.II/dofs/dof_tools.h>
 *   #include <deal.II/dofs/dof_renumbering.h>
 *   
 *   #include <deal.II/fe/fe_values.h>
 *   #include <deal.II/fe/fe_q.h>
 *   
 *   #include <deal.II/numerics/vector_tools.h>
 *   #include <deal.II/numerics/matrix_tools.h>
 *   #include <deal.II/numerics/data_out.h>
 *   #include <deal.II/numerics/error_estimator.h>
 *   #include <deal.II/numerics/solution_transfer.h>
 *   
 *   #include <typeinfo>
 *   #include <fstream>
 *   #include <iostream>
 *   
 *   #include <deal.II/numerics/solution_transfer.h>
 *   
 * @endcode
 * 
 * Open a namespace for this program and import everything from the
 * dealii namespace into it.
 * 
 * @code
 *   namespace nsp
 *   {
 *     using namespace dealii;
 *   
 * @endcode
 * 
 * ********************************************************
 * 
 * @code
 *     class ParameterReader : public Subscriptor
 *     {
 *     public:
 *       ParameterReader(ParameterHandler &);
 *       void read_parameters(const std::string);
 *     private:
 *       void declare_parameters();
 *       ParameterHandler &prm;
 *     };
 * @endcode
 * 
 * Constructor
 * 
 * @code
 *     ParameterReader::ParameterReader(ParameterHandler &paramhandler):
 *       prm(paramhandler)
 *     {}
 *   
 *     void ParameterReader::declare_parameters()
 *     {
 *   
 *       prm.enter_subsection ("Global Parameters");
 *       {
 *         prm.declare_entry("p", "100",Patterns::Double(2.1),
 *                           "Penalization parameter");
 *         prm.declare_entry("known_solution", "true",Patterns::Bool(),
 *                           "Whether the exact solution is known");
 *       }
 *       prm.leave_subsection ();
 *   
 *       prm.enter_subsection ("Mesh & Refinement Parameters");
 *       {
 *         prm.declare_entry("Code for the domain", "0",Patterns::Integer(0,2),
 *                           "Number identifying the domain in which we solve the problem");
 *         prm.declare_entry("No of initial refinements", "4",Patterns::Integer(0),
 *                           "Number of global mesh refinement steps applied to initial coarse grid");
 *         prm.declare_entry("No of adaptive refinements", "8",Patterns::Integer(0),
 *                           "Number of global adaptive mesh refinements");
 *         prm.declare_entry("top_fraction_of_cells", "0.25",Patterns::Double(0),
 *                           "refinement threshold");
 *         prm.declare_entry("bottom_fraction_of_cells", "0.05",Patterns::Double(0),
 *                           "coarsening threshold");
 *       }
 *       prm.leave_subsection ();
 *   
 *   
 *       prm.enter_subsection ("Algorithm Parameters");
 *       {
 *         prm.declare_entry("Descent_direction", "0",Patterns::Integer(0,1),
 *                           "0: Preconditioned descent, 1: Newton Method");
 *         prm.declare_entry("init_p", "10",Patterns::Double(2),
 *                           "Initial p");
 *         prm.declare_entry("delta_p", "50",Patterns::Double(0),
 *                           "increase of p");
 *         prm.declare_entry("Max_CG_it", "1500",Patterns::Integer(1),
 *                           "Maximum Number of CG iterations");
 *         prm.declare_entry("CG_tol", "1e-10",Patterns::Double(0),
 *                           "Tolerance for CG iterations");
 *         prm.declare_entry("max_LS_it", "45",Patterns::Integer(1),
 *                           "Maximum Number of LS iterations");
 *         prm.declare_entry("line_search_tolerence", "1e-6",Patterns::Double(0),
 *                           "line search tolerance constant (c1 in Nocedal-Wright)");
 *         prm.declare_entry("init_step_length", "1e-2",Patterns::Double(0),
 *                           "initial step length in line-search");
 *         prm.declare_entry("Max_inner", "800",Patterns::Integer(1),
 *                           "Maximum Number of inner iterations");
 *         prm.declare_entry("eps", "1.0e-8",Patterns::Double(0),
 *                           "Threshold on norm of the derivative to declare optimality achieved");
 *         prm.declare_entry("hi_eps", "1.0e-9",Patterns::Double(0),
 *                           "Threshold on norm of the derivative to declare optimality achieved in highly refined mesh");
 *         prm.declare_entry("hi_th", "8",Patterns::Integer(0),
 *                           "Number of adaptive refinement before change convergence threshold");
 *       }
 *       prm.leave_subsection ();
 *   
 *     }
 *     void ParameterReader::read_parameters (const std::string parameter_file)
 *     {
 *       declare_parameters();
 *       prm.parse_input (parameter_file);
 *     }
 *   
 * @endcode
 * 
 * ******************************************************************************************
 * The solution of the elastoplastic torsion problem on the unit disk with rhs=4.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     class Solution : public Function<dim>
 *     {
 *     public:
 *       Solution () : Function<dim>() {}
 *       virtual double value (const Point<dim> &pto, const unsigned int component = 0) const override;
 *       virtual Tensor<1,dim> gradient (const Point<dim> &pto, const unsigned int component = 0) const override;
 *     };
 *   
 *     template <int dim>
 *     double Solution<dim>::value (const Point<dim> &pto,const unsigned int) const
 *     {
 *       double r=sqrt(pto.square());
 *       if (r<0.5)
 *         return -1.0*std::pow(r,2.0)+0.75;
 *       else
 *         return 1.0-r;
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &pto,const unsigned int) const
 *     {
 *       double r=sqrt(pto.square());
 *       if (r<0.5)
 *         return -2.0*pto;
 *       else
 *         return  -1.0*pto/r;
 *     }
 *   
 *   
 *   
 *   
 * @endcode
 * 
 * ****************************************************************************************** 
 * 
 * @code
 *     /*                 Compute the Lagrange multiplier (as a derived quantity)                   */
 *   
 *   
 *     template <int dim>
 *     class ComputeMultiplier : public DataPostprocessor<dim>
 *     {
 *     private:
 *       double p;
 *     public:
 *       ComputeMultiplier (double pe);
 *   
 *       virtual
 *       void compute_derived_quantities_scalar (
 *         const std::vector< double > &,
 *         const std::vector< Tensor< 1, dim > > &,
 *         const std::vector< Tensor< 2, dim > > &,
 *         const std::vector< Point< dim > > &,
 *         const std::vector< Point< dim > > &,
 *         std::vector< Vector< double > > &
 *       ) const;
 *   
 *       virtual std::vector<std::string> get_names () const override;
 *   
 *       virtual
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *       get_data_component_interpretation () const override;
 *   
 *       virtual UpdateFlags get_needed_update_flags () const override;
 *     };
 *   
 *   
 *     template <int dim>
 *     ComputeMultiplier<dim>::ComputeMultiplier (double pe): p(pe)
 *     {}
 *   
 *   
 *     template <int dim>
 *     void ComputeMultiplier<dim>::compute_derived_quantities_scalar(
 *       const std::vector< double > &   /*uh*/,
 *       const std::vector< Tensor< 1, dim > >   &duh,
 *       const std::vector< Tensor< 2, dim > > &  /*dduh*/,
 *       const std::vector< Point< dim > > &  /* normals*/,
 *       const std::vector< Point< dim > > &  /*evaluation_points*/,
 *       std::vector< Vector< double > >    &computed_quantities )   const
 *     {
 *       const unsigned int n_quadrature_points = duh.size();
 *   
 *       for (unsigned int q=0; q<n_quadrature_points; ++q)
 *         {
 *           long  double sqrGrad=duh[q]* duh[q]; //squared norm of the gradient
 *           long double exponent=(p-2.0)/2*std::log(sqrGrad);
 *           computed_quantities[q](0) = std::sqrt(sqrGrad); // norm of the gradient
 *           computed_quantities[q](1)= std::exp(exponent); // multiplier
 *         }
 *     }
 *   
 *   
 *   
 *   
 *   
 *     template <int dim>
 *     std::vector<std::string>
 *     ComputeMultiplier<dim>::get_names() const
 *     {
 *       std::vector<std::string> solution_names;
 *       solution_names.push_back ("Gradient_norm");
 *       solution_names.push_back ("Lagrange_multiplier");
 *       return solution_names;
 *     }
 *   
 *   
 *     template <int dim>
 *     UpdateFlags
 *     ComputeMultiplier<dim>::get_needed_update_flags () const
 *     {
 *       return update_gradients;
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *     ComputeMultiplier<dim>:: get_data_component_interpretation () const
 *     {
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *       interpretation;
 * @endcode
 * 
 * norm of the gradient
 * 
 * @code
 *       interpretation.push_back (DataComponentInterpretation::component_is_scalar);
 * @endcode
 * 
 * Lagrange multiplier
 * 
 * @code
 *       interpretation.push_back (DataComponentInterpretation::component_is_scalar);
 *       return interpretation;
 *     }
 *   
 *   
 *   
 *   
 *   
 * @endcode
 * 
 * *************************************************************************************** 
 * 
 * @code
 *     template <int dim>
 *     class ElastoplasticTorsion
 *     {
 *     public:
 *       ElastoplasticTorsion (ParameterHandler &);
 *       ~ElastoplasticTorsion ();
 *       void run ();
 *   
 *     private:
 *       void setup_system (const bool initial_step);
 *       void assemble_system ();
 *       bool solve (const int inner_it);
 *       void init_mesh ();
 *       void refine_mesh ();
 *       void set_boundary_values ();
 *       double phi (const double alpha) const;
 *       bool checkWolfe(double &alpha, double &phi_alpha) const;
 *       bool determine_step_length (const int inner_it);
 *       void print_it_message (const int counter, bool ks);
 *       void output_results (unsigned int refinement) const;
 *       void format_convergence_tables();
 *       void process_solution (const unsigned int cycle);
 *       void process_multiplier (const unsigned int cycle,const int iter,double time);
 *       double dual_error () const;
 *       double dual_infty_error () const;
 *       double W (double Du2) const;
 *       double Wp (double Du2) const;
 *       double G (double Du2) const;
 *   
 *   
 *   
 *       ParameterHandler &prm;
 *       Triangulation<dim> triangulation;
 *       DoFHandler<dim> dof_handler;
 *       AffineConstraints<double> hanging_node_constraints;
 *       SparsityPattern sparsity_pattern;
 *       SparseMatrix<double> system_matrix;
 *       ConvergenceTable convergence_table;
 *       ConvergenceTable dual_convergence_table;
 *       Vector<double> present_solution;
 *       Vector<double> newton_update;
 *       Vector<double> system_rhs;
 *       Vector<double> grad_norm;
 *       Vector<double> lambda;
 *   
 *   
 *       double step_length,phi_zero,phi_alpha,phip,phip_zero;
 *       double old_step,old_phi_zero,old_phip;
 *       double L2_error;
 *       double H1_error;
 *       double Linfty_error;
 *       double dual_L1_error;
 *       double dual_L_infty_error;
 *       FE_Q<dim> fe;
 *       double p;
 *       double line_search_tolerence; // c_1 in Nocedal & Wright
 *       unsigned int dir_id;
 *       std::string elements;
 *       std::string Method;
 *   
 *     };
 *   
 *     /*******************************************************************************************/
 * @endcode
 * 
 * Boundary condition
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     class BoundaryValues : public Function<dim>
 *     {
 *     public:
 *       BoundaryValues () : Function<dim>() {}
 *   
 *       virtual double value (const Point<dim>   &p,
 *                             const unsigned int  component = 0) const override;
 *     };
 *   
 *   
 *     template <int dim>
 *     double BoundaryValues<dim>::value (const Point<dim> &/*pto*/,
 *                                        const unsigned int /*component*/) const
 *     {
 * @endcode
 * 
 * could be anything else (theory works provided |Dg|_infty < 1/2)
 * 
 * @code
 *       return 0.0;
 *   
 *       /* A challenging BC leading to overdetermined problems
 *          it is regulated by the parameter 0<eta<1.
 *          eta closer to 1 leads to more difficult problems.
 *   
 *       double pii=numbers::PI;
 *       double theta=std::atan2(p[1],p[0])+pii;
 *       double eta=0.9;
 *   
 *           if (theta <= 0.5)
 *           return eta*(theta*theta);
 *         else if ((theta >0.5) & (theta<= pii-0.5))
 *           return eta*(theta-0.25);
 *         else if ((theta>pii-0.5)&(theta<= pii+0.5))
 *           return eta*(pii-0.75-(theta-(pii-0.5))*(theta-(pii+0.5)));
 *         else if ((theta>pii+0.5)&(theta<= 2*pii-0.5))
 *           return eta*((2*pii-theta)-0.25);
 *         else
 *           return eta*((theta-2*pii)*(theta-2*pii) );*/
 *     }
 *   
 *   
 *   
 *     /******************************************************************************/
 * @endcode
 * 
 * Right-Hand Side
 * 
 * @code
 *     template <int dim>
 *     class RightHandSide : public Function<dim>
 *     {
 *     public:
 *       RightHandSide () : Function<dim>() {}
 *       virtual double value (const Point<dim> &p,
 *                             const unsigned int component = 0) const override;
 *     };
 *   
 *     template <int dim>
 *     double RightHandSide<dim>::value (const Point<dim> &/*p*/,
 *                                       const unsigned int /*component*/) const
 *     {
 * @endcode
 * 
 * set to constant = 4, for which explicit solution to compare exists
 * could be anything
 * 
 * @code
 *       double return_value = 4.0;
 *       return return_value;
 *     }
 *   
 *   
 *   
 *     /*******************************************************************/
 * @endcode
 * 
 * The ElastoplasticTorsion class implementation
 * 
 
 * 
 * Constructor of the class
 * 
 * @code
 *     template <int dim>
 *     ElastoplasticTorsion<dim>::ElastoplasticTorsion (ParameterHandler &param):
 *       prm(param),
 *       dof_handler (triangulation),
 *       L2_error(1.0),
 *       H1_error(1.0),
 *       Linfty_error(1.0),
 *       dual_L1_error(1.0),
 *       dual_L_infty_error(1.0),
 *       fe(2)
 *     {
 *       prm.enter_subsection ("Global Parameters");
 *       p=prm.get_double("p");
 *       prm.leave_subsection ();
 *       prm.enter_subsection ("Algorithm Parameters");
 *       line_search_tolerence=prm.get_double("line_search_tolerence");
 *       dir_id=prm.get_integer("Descent_direction");
 *       prm.leave_subsection ();
 *       if (fe.degree==1)
 *         elements="P1";
 *       else elements="P2";
 *   
 *       if (dir_id==0)
 *         Method="Precond";
 *       else
 *         Method="Newton";
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     ElastoplasticTorsion<dim>::~ElastoplasticTorsion ()
 *     {
 *       dof_handler.clear ();
 *     }
 *   
 *     /*****************************************************************************************/
 * @endcode
 * 
 * print iteration message
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::print_it_message (const int counter, bool ks)
 *     {
 *       if (ks)
 *         {
 *           process_solution (counter);
 *           std::cout << "iteration="<< counter+1 << "  J(u_h)= "<< phi_zero << ", H1 error: "
 *                     <<  H1_error  <<", W0-1,infty error: "<< Linfty_error<< " J'(u_h)(w)= "<< phip
 *                     << ", |J'(u_h)|= "<< system_rhs.l2_norm()<<std::endl;
 *         }
 *       else
 *         {
 *           std::cout << "iteration= " << counter+1 << " J(u_h)= "
 *                     << phi_alpha << " J'(u_h)= "<< phip<<std::endl;
 *         }
 *     }
 *   
 *   
 *     /*****************************************************************************************/
 * @endcode
 * 
 * Convergence Tables
 * 
 
 * 
 * 
 
 * 
 * 
 * @code
 *     /*************************************************************/
 * @endcode
 * 
 * formatting
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::format_convergence_tables()
 *     {
 *       convergence_table.set_precision("L2", 3);
 *       convergence_table.set_precision("H1", 3);
 *       convergence_table.set_precision("Linfty", 3);
 *       convergence_table.set_precision("function value", 3);
 *       convergence_table.set_precision("derivative", 3);
 *       dual_convergence_table.set_precision("dual_L1", 3);
 *       dual_convergence_table.set_precision("dual_Linfty", 3);
 *       dual_convergence_table.set_precision("L2", 3);
 *       dual_convergence_table.set_precision("H1", 3);
 *       dual_convergence_table.set_precision("Linfty", 3);
 *       convergence_table.set_scientific("L2", true);
 *       convergence_table.set_scientific("H1", true);
 *       convergence_table.set_scientific("Linfty", true);
 *       convergence_table.set_scientific("function value", true);
 *       convergence_table.set_scientific("derivative", true);
 *       dual_convergence_table.set_scientific("dual_L1", true);
 *       dual_convergence_table.set_scientific("dual_Linfty", true);
 *       dual_convergence_table.set_scientific("L2", true);
 *       dual_convergence_table.set_scientific("H1", true);
 *       dual_convergence_table.set_scientific("Linfty", true);
 *   
 *     }
 *   
 *     /****************************************/
 * @endcode
 * 
 * fill-in entry for the solution
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::process_solution (const unsigned int it)
 *     {
 *       Vector<float> difference_per_cell (triangulation.n_active_cells());
 *   
 * @endcode
 * 
 * compute L2 error (save to difference_per_cell)
 * 
 * @code
 *       VectorTools::integrate_difference (dof_handler,present_solution,
 *                                          Solution<dim>(),difference_per_cell,QGauss<dim>(3),VectorTools::L2_norm);
 *       L2_error = difference_per_cell.l2_norm();
 *   
 * @endcode
 * 
 * compute H1 error (save to difference_per_cell)
 * 
 * @code
 *       VectorTools::integrate_difference (dof_handler,present_solution,Solution<dim>(),
 *                                          difference_per_cell,QGauss<dim>(3),VectorTools::H1_seminorm);
 *       H1_error = difference_per_cell.l2_norm();
 *   
 * @endcode
 * 
 * compute W1infty error (save to difference_per_cell)
 * 
 * @code
 *       const QTrapezoid<1> q_trapez;
 *       const QIterated<dim> q_iterated (q_trapez, 5);
 *       VectorTools::integrate_difference (dof_handler,present_solution,Solution<dim>(),
 *                                          difference_per_cell,q_iterated,VectorTools::W1infty_seminorm);
 *       Linfty_error = difference_per_cell.linfty_norm();
 *   
 *   
 *       convergence_table.add_value("cycle", it);
 *       convergence_table.add_value("p", p);
 *       convergence_table.add_value("L2", L2_error);
 *       convergence_table.add_value("H1", H1_error);
 *       convergence_table.add_value("Linfty", Linfty_error);
 *       convergence_table.add_value("function value", phi_alpha);
 *       convergence_table.add_value("derivative", phip);
 *     }
 *   
 *   
 *     /***************************************/
 * @endcode
 * 
 * fill-in entry  for the multiplier
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::process_multiplier (const unsigned int cycle, const int iter,double time)
 *     {
 *       const unsigned int n_active_cells=triangulation.n_active_cells();
 *       const unsigned int n_dofs=dof_handler.n_dofs();
 *       dual_L1_error=dual_error();
 *       dual_L_infty_error=dual_infty_error();
 *   
 *   
 *       dual_convergence_table.add_value("cycle", cycle);
 *       dual_convergence_table.add_value("p", p);
 *       dual_convergence_table.add_value("iteration_number", iter);
 *       dual_convergence_table.add_value("cpu_time", time);
 *       dual_convergence_table.add_value("cells", n_active_cells);
 *       dual_convergence_table.add_value("dofs", n_dofs);
 *       dual_convergence_table.add_value("L2", L2_error);
 *       dual_convergence_table.add_value("H1", H1_error);
 *       dual_convergence_table.add_value("Linfty", Linfty_error);
 *       dual_convergence_table.add_value("dual_L1", dual_L1_error);
 *       dual_convergence_table.add_value("dual_Linfty", dual_L_infty_error);
 *   
 *     }
 *   
 *   
 *   
 *   
 *     /****************************************************************************************/
 * @endcode
 * 
 * ElastoplasticTorsion::setup_system
 * unchanged from @ref step_15 "step-15"
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::setup_system (const bool initial_step)
 *     {
 *       if (initial_step)
 *         {
 *           dof_handler.distribute_dofs (fe);
 *           present_solution.reinit (dof_handler.n_dofs());
 *           grad_norm.reinit (dof_handler.n_dofs());
 *           lambda.reinit (dof_handler.n_dofs());
 *   
 *           hanging_node_constraints.clear ();
 *           DoFTools::make_hanging_node_constraints (dof_handler,
 *                                                    hanging_node_constraints);
 *           hanging_node_constraints.close ();
 *         }
 *   
 *   
 * @endcode
 * 
 * The remaining parts of the function
 * 
 
 * 
 * 
 * @code
 *       newton_update.reinit (dof_handler.n_dofs());
 *       system_rhs.reinit (dof_handler.n_dofs());
 *       DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
 *       hanging_node_constraints.condense (c_sparsity);
 *       sparsity_pattern.copy_from(c_sparsity);
 *       system_matrix.reinit (sparsity_pattern);
 *     }
 *   
 *     /***************************************************************************************/
 *     /* the coefficients W, W' and G defining the problem.
 *   
 *        Min_u \int W(|Du|^2) dx
 *   
 *     They must be consistent as G(s)=W'(s)+2s W''(s) for any s>0.
 *     recall that they receive the SQUARED gradient. */
 *   
 *     template <int dim>
 *     double ElastoplasticTorsion<dim>::W (double Du2) const
 *     {
 *       return Du2;
 *     }
 *   
 *     template <int dim>
 *     double ElastoplasticTorsion<dim>::Wp (double /*Du2*/) const
 *     {
 *       return 1.0;
 *     }
 *   
 *     template <int dim>
 *     double ElastoplasticTorsion<dim>::G (double /*Du2*/) const
 *     {
 *       return 1.0;
 *     }
 *     /***************************************************************************************/
 *   
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::assemble_system ()
 *     {
 *   
 *   
 *       const QGauss<dim>  quadrature_formula(3);
 *       const RightHandSide<dim> right_hand_side;
 *       system_matrix = 0;
 *       system_rhs = 0;
 *   
 *       FEValues<dim> fe_values (fe, quadrature_formula,
 *                                update_gradients         |
 *                                update_values           |
 *                                update_quadrature_points |
 *                                update_JxW_values);
 *   
 *       const unsigned int           dofs_per_cell = fe.dofs_per_cell;
 *       const unsigned int           n_q_points    = quadrature_formula.size();
 *   
 *       FullMatrix<double>           cell_matrix (dofs_per_cell, dofs_per_cell);
 *       Vector<double>               cell_rhs (dofs_per_cell);
 *   
 *       std::vector<Tensor<1, dim> > old_solution_gradients(n_q_points);
 *       std::vector<types::global_dof_index>    local_dof_indices (dofs_per_cell);
 *   
 *   
 *       typename DoFHandler<dim>::active_cell_iterator
 *       cell = dof_handler.begin_active(),
 *       endc = dof_handler.end();
 *       for (; cell!=endc; ++cell)
 *         {
 *           cell_matrix = 0;
 *           cell_rhs = 0;
 *   
 *           fe_values.reinit (cell);
 *           fe_values.get_function_gradients(present_solution,
 *                                            old_solution_gradients);
 *   
 *           for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *             {
 *               long double coeff=0.0;
 *               long double a=old_solution_gradients[q_point] * old_solution_gradients[q_point];
 *               long double exponent=(p-2.0)/2*std::log(a);
 *               coeff= std::exp( exponent);
 *               for (unsigned int i=0; i<dofs_per_cell; ++i)
 *                 {
 *                   for (unsigned int j=0; j<dofs_per_cell; ++j)
 *                     {
 *                       if (dir_id==1)
 *                         {
 *                           cell_matrix(i, j) +=  fe_values.shape_grad(i, q_point) *  fe_values.shape_grad(j, q_point)
 *                                                 * (G(a)+(p-1.0)*coeff)    * fe_values.JxW(q_point);
 *                         }
 *                       else
 *                         {
 *                           cell_matrix(i, j) +=  fe_values.shape_grad(i, q_point) *  fe_values.shape_grad(j, q_point)
 *                                                 * (Wp(a)+coeff)
 *                                                 * fe_values.JxW(q_point);
 *                         }
 *                     }
 *   
 *                   cell_rhs(i) -= (  fe_values.shape_grad(i, q_point)
 *                                     * old_solution_gradients[q_point]
 *                                     * (Wp(a)+coeff)
 *                                     -right_hand_side.value(fe_values.quadrature_point(q_point))
 *                                     *fe_values.shape_value(i, q_point)
 *                                  )
 *                                  * fe_values.JxW(q_point);
 *                 }
 *             }
 *   
 *           cell->get_dof_indices (local_dof_indices);
 *           for (unsigned int i=0; i<dofs_per_cell; ++i)
 *             {
 *               for (unsigned int j=0; j<dofs_per_cell; ++j)
 *                 system_matrix.add (local_dof_indices[i],
 *                                    local_dof_indices[j],
 *                                    cell_matrix(i,j));
 *   
 *               system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *             }
 *         }
 *   
 *       hanging_node_constraints.condense (system_matrix);
 *       hanging_node_constraints.condense (system_rhs);
 *   
 *       std::map<types::global_dof_index,double> boundary_values;
 *       VectorTools::interpolate_boundary_values (dof_handler,
 *                                                 0,
 *                                                 Functions::ZeroFunction<dim>(),
 *                                                 boundary_values);
 *       MatrixTools::apply_boundary_values (boundary_values,
 *                                           system_matrix,
 *                                           newton_update,
 *                                           system_rhs);
 *     }
 *   
 *   
 *   
 *   
 *     /**********************************      Refine Mesh      ****************************************/
 * @endcode
 * 
 * unchanged from @ref step_15 "step-15"
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::refine_mesh ()
 *     {
 *       using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
 *   
 *       Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 *       KellyErrorEstimator<dim>::estimate (dof_handler,
 *                                           QGauss<dim-1>(3),
 *                                           FunctionMap(),
 *                                           present_solution,
 *                                           estimated_error_per_cell);
 *   
 *       prm.enter_subsection ("Mesh & Refinement Parameters");
 *       const double top_fraction=prm.get_double("top_fraction_of_cells");
 *       const double bottom_fraction=prm.get_double("bottom_fraction_of_cells");
 *       prm.leave_subsection ();
 *       GridRefinement::refine_and_coarsen_fixed_number (triangulation,
 *                                                        estimated_error_per_cell,
 *                                                        top_fraction, bottom_fraction);
 *   
 *       triangulation.prepare_coarsening_and_refinement ();
 *       SolutionTransfer<dim> solution_transfer(dof_handler);
 *       solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
 *       triangulation.execute_coarsening_and_refinement();
 *       dof_handler.distribute_dofs(fe);
 *   #  if DEAL_II_VERSION_GTE(9, 7, 0)
 *       present_solution.reinit(dof_handler.n_dofs());
 *       solution_transfer.interpolate(present_solution);
 *   #  else
 *       Vector<double> tmp(dof_handler.n_dofs());
 *       solution_transfer.interpolate(present_solution, tmp);
 *       present_solution = tmp;
 *   #  endif
 *       set_boundary_values ();
 *       hanging_node_constraints.clear();
 *   
 *       DoFTools::make_hanging_node_constraints(dof_handler,
 *                                               hanging_node_constraints);
 *       hanging_node_constraints.close();
 *       hanging_node_constraints.distribute (present_solution);
 *       setup_system (false);
 *     }
 *   
 *   
 *     /*******************************************************************************************/
 * @endcode
 * 
 * Dump the norm of the gradient and the lagrange multiplier in vtu format for visualization
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::output_results (unsigned int counter) const
 *     {
 * @endcode
 * 
 * multiplier object contains both |Du| and lambda.
 * 
 * @code
 *       ComputeMultiplier<dim> multiplier(p);
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler (dof_handler);
 *       data_out.add_data_vector (present_solution, "solution");
 *       data_out.add_data_vector (present_solution, multiplier);
 *       data_out.build_patches ();
 *       std::ostringstream p_str;
 *       p_str << p<<"-cycle-"<<counter;
 *       std::string str = p_str.str();
 *       const std::string filename = "solution-" + str+".vtu";
 *       std::ofstream output (filename.c_str());
 *       data_out.write_vtu (output);
 *     }
 *   
 *     /********************************************************************************************/
 * @endcode
 * 
 * unchanged from @ref step_15 "step-15"
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::set_boundary_values ()
 *     {
 *       std::map<types::global_dof_index, double> boundary_values;
 *       VectorTools::interpolate_boundary_values (dof_handler,
 *                                                 0,
 *                                                 BoundaryValues<dim>(),
 *                                                 boundary_values);
 *       for (std::map<types::global_dof_index, double>::const_iterator
 *            bp = boundary_values.begin();
 *            bp != boundary_values.end(); ++bp)
 *         present_solution(bp->first) = bp->second;
 *     }
 *   
 *   
 *     /****************************************************************************************/
 * @endcode
 * 
 * COMPUTE @f$\phi(\alpha)=J_p(u_h+\alpha w)@f$
 * 
 * @code
 *     template <int dim>
 *     double ElastoplasticTorsion<dim>::phi (const double alpha) const
 *     {
 *       double obj = 0.0;
 *       const RightHandSide<dim> right_hand_side;
 *       Vector<double> evaluation_point (dof_handler.n_dofs());
 *       evaluation_point = present_solution;  // copy of u_h
 *       evaluation_point.add (alpha, newton_update); // u_{n+1}=u_n+alpha w_n
 *   
 *       const QGauss<dim>  quadrature_formula(3);
 *       FEValues<dim> fe_values (fe, quadrature_formula,
 *                                update_gradients         |
 *                                update_values    |
 *                                update_quadrature_points |
 *                                update_JxW_values);
 *   
 *       const unsigned int           dofs_per_cell = fe.dofs_per_cell;
 *       const unsigned int           n_q_points    = quadrature_formula.size();
 *   
 *       Vector<double>               cell_residual (dofs_per_cell);
 *       std::vector<Tensor<1, dim> > gradients(n_q_points);
 *       std::vector<double> values(n_q_points);
 *   
 *   
 *       std::vector<types::global_dof_index>    local_dof_indices (dofs_per_cell);
 *   
 *       typename DoFHandler<dim>::active_cell_iterator
 *       cell = dof_handler.begin_active(),
 *       endc = dof_handler.end();
 *       for (; cell!=endc; ++cell)
 *         {
 *           cell_residual = 0;
 *           fe_values.reinit (cell);
 *           fe_values.get_function_gradients (evaluation_point, gradients);
 *           fe_values.get_function_values (evaluation_point, values);
 *   
 *   
 *           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *             {
 *               double Du2=gradients[q_point] *  gradients[q_point]; // Du2=|Du|^2
 *               double penalty;
 *               if (Du2<1.0e-10)
 *                 penalty=0.0;
 *               else
 *                 penalty=std::pow(Du2,p/2.0); // penalty=|Du|^p
 *   
 * @endcode
 * 
 * obj+= 1/2 W(|Du|^2)+1/p |Du|^p -fu (see (1))
 * 
 * @code
 *               obj+=(
 *                      (0.5*W(Du2)+penalty/p)- right_hand_side.value(fe_values.quadrature_point(q_point))*values[q_point]
 *                    ) * fe_values.JxW(q_point);
 *             }
 *   
 *         }
 *   
 *       return obj;
 *     }
 *   
 *   
 *     /***************************************************************************************************/
 * @endcode
 * 
 * Compute L^1 error norm of Lagrange Multiplier
 * with respect to exact solution (cf. Alvarez & Flores, 2015)
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     double ElastoplasticTorsion<dim>::dual_error () const
 *     {
 *       double obj = 0.0;
 *   
 *       const QGauss<dim>  quadrature_formula(3);
 *       FEValues<dim> fe_values (fe, quadrature_formula,
 *                                update_gradients         |
 *                                update_quadrature_points |
 *                                update_JxW_values);
 *   
 *       const unsigned int           dofs_per_cell = fe.dofs_per_cell;
 *       const unsigned int           n_q_points    = quadrature_formula.size();
 *   
 *       Vector<double>               cell_residual (dofs_per_cell);
 *       std::vector<Tensor<1, dim> > gradients(n_q_points);
 *   
 *       std::vector<types::global_dof_index>    local_dof_indices (dofs_per_cell);
 *   
 *       typename DoFHandler<dim>::active_cell_iterator
 *       cell = dof_handler.begin_active(),
 *       endc = dof_handler.end();
 *       for (; cell!=endc; ++cell)
 *         {
 *           cell_residual = 0;
 *           fe_values.reinit (cell);
 *           fe_values.get_function_gradients (present_solution, gradients);
 *   
 *           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *             {
 *               double coeff=gradients[q_point] *  gradients[q_point] ;
 *               if (coeff<1.0e-15)
 *                 coeff=0.0;
 *               else
 *                 coeff=std::pow(coeff,(p-2.0)/2.0); // |Du_p|^(p-2)
 *   
 *               double r=std::sqrt(fe_values.quadrature_point(q_point).square());
 *               double exact=0;
 *               if (r>0.5)
 *                 exact= 2*r-1;
 *   
 *               obj+=(     std::abs(coeff-exact) ) * fe_values.JxW(q_point);
 *             }
 *   
 *         }
 *   
 *       return obj;
 *     }
 *   
 *     /*******************************************************************************************/
 * @endcode
 * 
 * Compute L^infinity error norm of Lagrange Multiplier
 * with respect to exact solution (cf. Alvarez & Flores, 2015)
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     double ElastoplasticTorsion<dim>::dual_infty_error () const
 *     {
 *       double obj = 0.0;
 *       const QTrapezoid<1> q_trapez;
 *       const QIterated<dim> quadrature_formula (q_trapez, 10);
 *   
 *       FEValues<dim> fe_values (fe, quadrature_formula,
 *                                update_gradients         |
 *                                update_quadrature_points );
 *   
 *       const unsigned int           dofs_per_cell = fe.dofs_per_cell;
 *       const unsigned int           n_q_points    = quadrature_formula.size();
 *   
 *       Vector<double>               cell_residual (dofs_per_cell);
 *       std::vector<Tensor<1, dim> > gradients(n_q_points);
 *   
 *       std::vector<types::global_dof_index>    local_dof_indices (dofs_per_cell);
 *   
 *       typename DoFHandler<dim>::active_cell_iterator
 *       cell = dof_handler.begin_active(),
 *       endc = dof_handler.end();
 *       for (; cell!=endc; ++cell)
 *         {
 *           cell_residual = 0;
 *           fe_values.reinit (cell);
 *           fe_values.get_function_gradients (present_solution, gradients);
 *   
 *           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
 *             {
 *               long double sqdGrad=gradients[q_point] *  gradients[q_point] ;
 *               double r=std::sqrt(fe_values.quadrature_point(q_point).square());
 *               double exact=0;
 *               if (r>0.5)
 *                 exact= 2*r-1.0;
 * @endcode
 * 
 * compute |Du|^(p-2) as  exp(p-2/2*log(Du^2))
 * 
 * @code
 *               long double exponent=(p-2.0)/2*std::log(sqdGrad);
 *               long double coeff=std::exp(exponent);
 *   
 *               if (std::abs(coeff-exact)>obj )
 *                 obj=std::abs(coeff-exact);
 *             }
 *   
 *         }
 *   
 *       return obj;
 *     }
 *   
 *     /*****************************************************************************************/
 * @endcode
 * 
 * check whether putative step-length satisfies sufficient decrease conditions
 * 
 * @code
 *     template <int dim>
 *     bool ElastoplasticTorsion<dim>::checkWolfe(double &alpha, double &phi_alpha) const
 *     {
 *       if (phi_alpha< phi_zero+line_search_tolerence*phip*alpha )
 *         return true;
 *       else
 *         return false;
 *     }
 *   
 *   
 *     /*****************************************************************************************/
 * @endcode
 * 
 * Find a step-length satisfying sufficient decrease condition by line-search
 * uses quadratic interpolation
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     bool ElastoplasticTorsion<dim>::determine_step_length(const int inner_it)
 *     {
 *       unsigned int it=0;
 *       bool done;
 *       double alpha,nalpha;
 *       prm.enter_subsection ("Algorithm Parameters");
 *       const unsigned int max_LS_it=prm.get_integer("max_LS_it");
 *       double init_SL=prm.get_double("init_step_length");
 *       prm.leave_subsection ();
 *       if (inner_it==0)
 *         alpha=init_SL;
 *       else
 *         {
 *           alpha=std::min(1.45*old_step*old_phip/phip,1.0);
 *         }
 *       phi_alpha=phi(alpha);
 *       std::cerr << "Step length=" << alpha << ", Value= " << phi_alpha;
 * @endcode
 * 
 * check if step-size satisfies sufficient decrease condition
 * 
 * @code
 *       done=checkWolfe(alpha,phi_alpha);
 *       if (done)
 *         std::cerr << " accepted" << std::endl;
 *       else
 *         std::cerr << " rejected" ;
 *   
 *       while ((!done) & (it<max_LS_it))
 *         {
 * @endcode
 * 
 * new try obtained by quadratic interpolation
 * 
 * @code
 *           nalpha=-(phip*alpha*alpha)/(2*(phi_alpha-phi_zero-phip*alpha));
 *   
 *           if (nalpha<1e-3*alpha ||  std::abs(nalpha-alpha)/alpha<1e-8)
 *             nalpha=alpha/2;
 *           else if ( phi_alpha-phi_zero>1e3*std::abs(phi_zero) )
 *             nalpha=alpha/10;
 *           alpha=nalpha;
 *           phi_alpha=phi(alpha);
 *           done=checkWolfe(alpha,phi_alpha);
 *           if (done)
 *             std::cerr << ", finished with steplength= "<< alpha<< ", fcn value= "<< phi_alpha<<std::endl;
 *           it=it+1;
 *         }
 *       if (!done)
 *         {
 *           std::cerr << ", max. no. of iterations reached with steplength= "<< alpha
 *                     << ", fcn value= "<< phi_alpha<<std::endl;
 *           return false;
 *         }
 *       else
 *         {
 *           step_length=alpha;
 *           return true;
 *         }
 *   
 *     }
 *   
 *     /**************************************************************************************************/
 * @endcode
 * 
 * ElastoplasticTorsion::init_mesh()
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::init_mesh ()
 *     {
 * @endcode
 * 
 * get parameters
 * 
 * @code
 *       prm.enter_subsection ("Mesh & Refinement Parameters");
 *       const int domain_id=prm.get_integer("Code for the domain");
 *       const int init_ref=prm.get_integer("No of initial refinements");
 *       prm.leave_subsection ();
 *   
 *   
 *       if (domain_id==0)
 *         {
 * @endcode
 * 
 * For the unit disk around the origin
 * 
 * @code
 *           GridGenerator::hyper_ball (triangulation);
 *           static const SphericalManifold<dim> boundary;
 *           triangulation.set_manifold (0, boundary);
 *         }
 *       else if (domain_id==1)
 *         {
 * @endcode
 * 
 * For the unit square
 * 
 * @code
 *           GridGenerator::hyper_cube (triangulation, 0, 1);
 *         }
 *       else if (domain_id==2)
 *         {
 *           /* For Glowinski's domain
 *               ___    ___   __ 1
 *              |   |__|   |  __ .8
 *            |          |
 *            |          |
 *            |__________|  __ 0
 *   
 *              |   |  |   |
 *              0  .4 .6   1
 *   
 *           */
 *           Triangulation<dim> tria1;
 *           Triangulation<dim> tria2;
 *           Triangulation<dim> tria3;
 *           Triangulation<dim> tria4;
 *           Triangulation<dim> tria5;
 *           Triangulation<dim> tria6;
 *           GridGenerator::hyper_rectangle(tria1, Point<2>(0.0,0.0), Point<2>(0.4,0.8));
 *           GridGenerator::hyper_rectangle(tria2, Point<2>(0.0,0.8), Point<2>(0.4,1.0));
 *           GridGenerator::hyper_rectangle(tria3, Point<2>(0.4,0.0), Point<2>(0.6,0.8));
 *           GridGenerator::hyper_rectangle(tria4, Point<2>(0.6,0.0), Point<2>(1.0,0.8));
 *           GridGenerator::hyper_rectangle(tria5, Point<2>(0.6,0.8), Point<2>(1.0,1.0));
 *           GridGenerator::merge_triangulations (tria1, tria2, tria6);
 *           GridGenerator::merge_triangulations (tria6, tria3, tria6);
 *           GridGenerator::merge_triangulations (tria6, tria4, tria6);
 *           GridGenerator::merge_triangulations (tria6, tria5, triangulation);
 *         }
 * @endcode
 * 
 * perform initial refinements
 * 
 * @code
 *       triangulation.refine_global(init_ref);
 *     }
 *   
 *     /**************************************************************************************************/
 * @endcode
 * 
 * ElastoplasticTorsion::solve(inner_it)
 * Performs one inner iteration
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     bool ElastoplasticTorsion<dim>::solve (const int inner_it)
 *     {
 *       prm.enter_subsection ("Algorithm Parameters");
 *       const unsigned int max_CG_it=prm.get_integer("Max_CG_it");
 *       const double CG_tol=prm.get_double("CG_tol");
 *       prm.leave_subsection ();
 *   
 *       SolverControl solver_control (max_CG_it,CG_tol);
 *       SolverCG<>    solver (solver_control);
 *   
 *       PreconditionSSOR<> preconditioner;
 *       preconditioner.initialize(system_matrix,0.25);
 *   
 *       solver.solve (system_matrix, newton_update, system_rhs,
 *                     preconditioner);
 *       hanging_node_constraints.distribute (newton_update);
 *       /******  save current quantities for line-search  **** */
 * @endcode
 * 
 * Recall that phi(alpha)=J(u+alpha w)
 * 
 * @code
 *       old_step=step_length;
 *       old_phi_zero=phi_zero;
 *       phi_zero=phi(0); // phi(0)=J(u)
 *       old_phip=phip;
 *       phip=-1.0*(newton_update*system_rhs); //phi'(0)=J'(u) *w, rhs=-J'(u).
 *       if (inner_it==0)
 *         phip_zero=phip;
 *   
 *       if (phip>0)   // this should not happen, step back
 *         {
 *           std::cout << "Not a descent direction!" <<std::endl;
 *           present_solution.add (-1.0*step_length, newton_update);
 *           step_length=step_length/2;
 *           phip=old_phip;
 *           return false;
 *         }
 *       else
 *         {
 *           if (determine_step_length(inner_it))
 *             {
 * @endcode
 * 
 * update u_{n+1}=u_n+alpha w_n
 * 
 * @code
 *               present_solution.add (step_length, newton_update);
 *               return true;
 *             }
 *           else return false;
 *         }
 *     }
 *   
 *   
 *   
 *     /*************************************************************************************************************/
 * @endcode
 * 
 * ElastoplasticTorsion::run
 * 
 * @code
 *     template <int dim>
 *     void ElastoplasticTorsion<dim>::run ()
 *     {
 *   
 * @endcode
 * 
 * get parameters
 * 
 * @code
 *       prm.enter_subsection ("Mesh & Refinement Parameters");
 *       const int adapt_ref=prm.get_integer("No of adaptive refinements");
 *       prm.leave_subsection ();
 *       prm.enter_subsection ("Algorithm Parameters");
 *       const int max_inner=prm.get_integer("Max_inner");
 *       const double eps=prm.get_double("eps");
 *       const double hi_eps=prm.get_double("hi_eps");
 *       const int hi_th=prm.get_integer("hi_th");
 *       const double init_p=prm.get_double("init_p");
 *       const double delta_p=prm.get_double("delta_p");
 *       prm.leave_subsection ();
 *       prm.enter_subsection ("Global Parameters");
 *       bool known_solution=prm.get_bool("known_solution");
 *       double actual_p=prm.get_double("p");
 *       prm.leave_subsection ();
 *       /************************/
 *   
 * @endcode
 * 
 * init Timer
 * 
 * @code
 *       Timer timer;
 *       double ptime=0.0;
 *       timer.start ();
 *   
 * @endcode
 * 
 * initialize mesh for the selected domain
 * 
 * @code
 *       init_mesh();
 *   
 * @endcode
 * 
 * setup FE space
 * 
 * @code
 *       setup_system (true);
 *       set_boundary_values ();
 *   
 * @endcode
 * 
 * init counters
 * 
 * @code
 *       int global_it=0;    // Total inner iterations (counting both loops)
 *       int cycle=0;          // Total outer iterations (counting both loops)
 *       int refinement = 0;    //  Refinements  performed (adaptive) = outer iterations 2nd loop
 *   
 *   
 * @endcode
 * 
 * prepare to start first loop
 * 
 * @code
 *       p=init_p;
 *       bool well_solved=true;
 *   
 *       /*****************************          First loop      ***********************************/
 *       /****************** Prepare initial condition using increasing p  *************************/
 *       while (p<actual_p) // outer iteration, increasing p.
 *         {
 *           std::cout <<"--Preparing initial condition with p="<<p<<" iter.= " << global_it<<  "  .-- "<< std::endl;
 *           timer.restart();
 *           for (int inner_iteration=0; inner_iteration<max_inner; ++inner_iteration,++global_it)
 *             {
 *               assemble_system ();
 *               well_solved=solve (inner_iteration);
 *               print_it_message (global_it, known_solution);
 *               if (
 *                 ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) <1e-4) & (cycle<1)) |
 *                 ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) <1e-5) & (cycle>=1)) |
 *                 !well_solved
 *               )
 *                 break;
 *             }
 *           ptime=timer.cpu_time();
 *           if (well_solved)
 *             output_results (cycle);
 *   
 *           if (known_solution)
 *             {
 *               process_multiplier(cycle,global_it,ptime);
 * @endcode
 * 
 * dual_convergence_table.write_tex(dual_error_table_file);
 * 
 * @code
 *             }
 *           refine_mesh();
 *           cycle++;
 *           p+=delta_p;
 *         }
 *       /***************************    first loop finished        ********************/
 *   
 *   
 * @endcode
 * 
 * prepare for second loop
 * 
 * @code
 *       p=actual_p;
 *       well_solved=true;
 *   
 *   
 *       /*****************************        Second loop         *********************************/
 *       /**************************** Solve problem for target p  *********************************/
 *   
 *       std::cout << "============ Solving  problem with p="   <<p << " =================="  << std::endl;
 *       /*****    Outer iteration - refining mesh  ******************/
 *       while ((cycle<adapt_ref) & well_solved)
 *         {
 *           timer.restart();
 * @endcode
 * 
 * inner iteration
 * 
 * @code
 *           for (int inner_iteration=0; inner_iteration<max_inner; ++inner_iteration,++global_it)
 *             {
 *               assemble_system ();
 *               well_solved=solve (inner_iteration);
 *               print_it_message (global_it, known_solution);
 *   
 *               if (
 *                 ((system_rhs.l2_norm()/std::sqrt(system_rhs.size()) < eps) & (refinement<hi_th)) |
 *                 (( system_rhs.l2_norm()/ std::sqrt (system_rhs.size()) <hi_eps)  | (!well_solved))
 *               )
 *                 break;
 *             }
 * @endcode
 * 
 * inner iterations finished
 * 
 * @code
 *           ptime=timer.cpu_time();
 *           if (well_solved)
 *             output_results (cycle);
 *   
 * @endcode
 * 
 * compute and display error, if the explicit solution is known
 * 
 * @code
 *           if (known_solution)
 *             {
 *               process_multiplier(cycle,global_it,ptime);
 *               std::cout << "finished with H1 error: " <<  H1_error << ", dual error (L1): "
 *                         << dual_L1_error << "dual error (L infty): "<<dual_L_infty_error <<std::endl;
 *             }
 *   
 * @endcode
 * 
 * update counters
 * 
 * @code
 *           ++refinement;
 *           ++cycle;
 * @endcode
 * 
 * refine mesh
 * 
 * @code
 *           std::cout << "******** Refined mesh " << cycle    << " ********"  << std::endl;
 *           refine_mesh();
 *         }// second loop
 *   
 * @endcode
 * 
 * write convergence tables to file
 * 
 * @code
 *       if (known_solution)
 *         {
 *           format_convergence_tables();
 *           std::string error_filename = "error"+Method+elements+".tex";
 *           std::ofstream error_table_file(error_filename.c_str());
 *           std::string dual_error_filename = "dual_error"+Method+elements+".tex";
 *           std::ofstream dual_error_table_file(dual_error_filename.c_str());
 *           convergence_table.write_tex(error_table_file);
 *           dual_convergence_table.write_tex(dual_error_table_file);
 *         }
 *     }//run()
 *   
 *   }//namespace
 *   
 *   /**********************************************************************************************/
 * @endcode
 * 
 * The main function
 * 
 * @code
 *   int main ()
 *   {
 *     try
 *       {
 *         using namespace dealii;
 *         using namespace nsp;
 *         deallog.depth_console (0);
 *   
 *         ParameterHandler prm;
 *         ParameterReader param(prm);
 *         param.read_parameters("EPT.prm");
 *         ElastoplasticTorsion<2> ElastoplasticTorsionProblem(prm);
 *         ElastoplasticTorsionProblem .run ();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl << std::endl
 *                   << "----------------------------------------------------"    << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     return 0;
 *   }
 * @endcode
 
 
*/