/* ---------------------------------------------------------------------
 *
 * Copyright (C) 2010 - 2015 by the deal.II authors
 *                and Salvador Flores.
 *
 *
 *
 * This is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * ---------------------------------------------------------------------
 *
 * Author: Salvador Flores,
 *         Center for Mathematical Modelling,
 *       Universidad de Chile, 2015.
 */


/*
 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

      step-15 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.

*/

// Include files

#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>

// Open a namespace for this program and import everything from the
// dealii namespace into it.
namespace nsp
{
  using namespace dealii;

// ********************************************************//
  class ParameterReader : public Subscriptor
  {
  public:
    ParameterReader(ParameterHandler &);
    void read_parameters(const std::string);
  private:
    void declare_parameters();
    ParameterHandler &prm;
  };
// Constructor
  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);
  }

// ******************************************************************************************//
// The solution of the elastoplastic torsion problem on the unit disk with rhs=4.

  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;
    virtual Tensor<1,dim> gradient (const Point<dim> &pto, const unsigned int component = 0) const;
  };

  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;
  }




// ****************************************************************************************** //
  /*                 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;

    virtual
    std::vector<DataComponentInterpretation::DataComponentInterpretation>
    get_data_component_interpretation () const;
    virtual UpdateFlags get_needed_update_flags () const;
  };


  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;
    // norm of the gradient
    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
    // Lagrange multiplier
    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
    return interpretation;
  }





// *************************************************************************************** //
  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;

  };

  /*******************************************************************************************/
//                              Boundary condition

  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;
  };


  template <int dim>
  double BoundaryValues<dim>::value (const Point<dim> &/*pto*/,
                                     const unsigned int /*component*/) const
  {
    // could be anything else (theory works provided |Dg|_infty < 1/2)
    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) );*/
  }



  /******************************************************************************/
//                       Right-Hand Side
  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;
  };

  template <int dim>
  double RightHandSide<dim>::value (const Point<dim> &/*p*/,
                                    const unsigned int /*component*/) const
  {
    // set to constant = 4, for which explicit solution to compare exists
    // could be anything
    double return_value = 4.0;
    return return_value;
  }



  /*******************************************************************/
// The ElastoplasticTorsion class implementation

// Constructor of the class
  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 ();
  }

  /*****************************************************************************************/
//  print iteration message

  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;
      }
  }


  /*****************************************************************************************/
//                            Convergence Tables


  /*************************************************************/
// formating

  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);

  }

  /****************************************/
// fill-in entry for the solution
  template <int dim>
  void ElastoplasticTorsion<dim>::process_solution (const unsigned int it)
  {
    Vector<float> difference_per_cell (triangulation.n_active_cells());

    // compute L2 error (save to difference_per_cell)
    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();

    // compute H1 error (save to difference_per_cell)
    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();

    // compute W1infty error (save to difference_per_cell)
    const QTrapez<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);
  }


  /***************************************/
// fill-in entry  for the multiplier
  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);

  }




  /****************************************************************************************/
// ElastoplasticTorsion::setup_system
// unchanged from step-15

  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 ();
      }


    // The remaining parts of the function

    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 coeffcients 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,
                                              ZeroFunction<dim>(),
                                              boundary_values);
    MatrixTools::apply_boundary_values (boundary_values,
                                        system_matrix,
                                        newton_update,
                                        system_rhs);
  }




  /**********************************      Refine Mesh      ****************************************/
// unchanged from step-15

  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);
    Vector<double> tmp(dof_handler.n_dofs());
    solution_transfer.interpolate(present_solution, tmp);
    present_solution = tmp;
    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);
  }


  /*******************************************************************************************/
// Dump the norm of the gradient and the lagrange multiplier in vtu format for visualization
  template <int dim>
  void ElastoplasticTorsion<dim>::output_results (unsigned int counter) const
  {
    // multiplier object contains both |Du| and lambda.
    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);
  }

  /********************************************************************************************/
// unchanged from step-15
  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;
  }


  /****************************************************************************************/
//  COMPUTE $\phi(\alpha)=J_p(u_h+\alpha w)$
  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

            // obj+= 1/2 W(|Du|^2)+1/p |Du|^p -fu (see (1))
            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;
  }


  /***************************************************************************************************/
// Compute L^1 error norm of Lagrange Multiplier
// with respect to exact solution (cf. Alvarez & Flores, 2015)

  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;
  }

  /*******************************************************************************************/
// Compute L^infinity error norm of Lagrange Multiplier
// with respect to exact solution (cf. Alvarez & Flores, 2015)

  template <int dim>
  double ElastoplasticTorsion<dim>::dual_infty_error () const
  {
    double obj = 0.0;
    const QTrapez<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;
            // compute |Du|^(p-2) as  exp(p-2/2*log(Du^2))
            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;
  }

  /*****************************************************************************************/
// check whether putative step-length satisfies sufficient decrease conditions
  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;
  }


  /*****************************************************************************************/
// Find a step-length satisfying sufficient decrease condition by line-search
// uses quadratic interpolation

  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;
    // check if step-size satisfies sufficient decrease condition
    done=checkWolfe(alpha,phi_alpha);
    if (done)
      std::cerr << " accepted" << std::endl;
    else
      std::cerr << " rejected" ;

    while ((!done) & (it<max_LS_it))
      {
        // new try obtained by quadratic interpolation
        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 wiht steplength= "<< alpha
                  << ", fcn value= "<< phi_alpha<<std::endl;
        return false;
      }
    else
      {
        step_length=alpha;
        return true;
      }

  }

  /**************************************************************************************************/
  // ElastoplasticTorsion::init_mesh()

  template <int dim>
  void ElastoplasticTorsion<dim>::init_mesh ()
  {
    // get parameters
    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)
      {
        // For the unit disk around the origin
        GridGenerator::hyper_ball (triangulation);
        static const SphericalManifold<dim> boundary;
        triangulation.set_manifold (0, boundary);
      }
    else if (domain_id==1)
      {
        // For the unit square
        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);
      }
    // perform initial refinements
    triangulation.refine_global(init_ref);
  }

  /**************************************************************************************************/
  // ElastoplasticTorsion::solve(inner_it)
  // Performs one inner iteration

  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  **** */
    // Recall that phi(alpha)=J(u+alpha w)
    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))
          {
            // update u_{n+1}=u_n+alpha w_n
            present_solution.add (step_length, newton_update);
            return true;
          }
        else return false;
      }
  }



  /*************************************************************************************************************/
// ElastoplasticTorsion::run
  template <int dim>
  void ElastoplasticTorsion<dim>::run ()
  {

    // get parameters
    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 ();
    /************************/

    // init Timer
    Timer timer;
    double ptime=0.0;
    timer.start ();

    // initalize mesh for the selected domain
    init_mesh();

    // setup FE space
    setup_system (true);
    set_boundary_values ();

    // init counters
    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


    // prepare to start first loop
    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);
            //dual_convergence_table.write_tex(dual_error_table_file);
          }
        refine_mesh();
        cycle++;
        p+=delta_p;
      }
    /***************************    first loop finished        ********************/


    // prepare for second loop
    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();
        // inner iteration
        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;
          }
        //inner iterations finished
        ptime=timer.cpu_time();
        if (well_solved)
          output_results (cycle);

        // compute and display error, if the explicit solution is known
        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;
          }

        // update counters
        ++refinement;
        ++cycle;
        // refine mesh
        std::cout << "******** Refined mesh " << cycle    << " ********"  << std::endl;
        refine_mesh();
      }// second loop

    // write convergence tables to file
    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

/**********************************************************************************************/
// The main function
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;
}