<commit_before>/* $Id$ */
/* Author: <NAME>, University of Heidelberg, 1999 */

/*    $Id$       */
/*                                                                */
/*    Copyright (C) 1999-2012 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyrightG and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

// @sect3{Include files}

// We are using the the same
// include files as in step-41:

#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.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/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.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/block_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>

#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_block_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/timer.h>
#include <fstream>
#include <iostream>
#include <list>
#include <time.h>


#include <deal.II/base/logstream.h>

namespace Step42
{
  using namespace dealii;


  // @sect3{The <code>PlasticityContactProblem</code> class template}

  // This class has only the purpose
  // to read in data from a picture file
  // that has to be stored as a pbm ascii
  // format. This data will be bilinear
  // interpolated and provides in this way
  // a function which describes an obstacle.

  template <int dim> class Input;

  // This class provides an interface
  // for a constitutive law. In this
  // example we are using an elastic
  // plastic material with linear,
  // isotropic hardening.

  template <int dim> class ConstitutiveLaw;

  // This class supplies all function
  // and variables needed to describe
  // the nonlinear contact problem. It is
  // close to step-41 but with some additional
  // features like: handling hanging nodes,
  // a newton method, using Trilinos and p4est
  // for parallel distributed computing.
  // To deal with hanging nodes makes
  // life a bit more complicated since
  // we need an other ConstraintMatrix now.
  // We create a newton method for the
  // active set method for the contact
  // situation and to handle the nonlinear
  // operator for the constitutive law.

  template <int dim>
  class PlasticityContactProblem
  {
  public:
    PlasticityContactProblem (int _n_refinements_global);
    void run ();

  private:
    void make_grid ();
    void setup_system();
    void assemble_nl_system (TrilinosWrappers::MPI::Vector &u);
    void residual_nl_system (TrilinosWrappers::MPI::Vector &u);
    void assemble_mass_matrix_diagonal (TrilinosWrappers::SparseMatrix &mass_matrix);
    void update_solution_and_constraints ();
    void dirichlet_constraints ();
    void solve ();
    void solve_newton ();
    void refine_grid ();
    void move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const;
    void output_results (const std::string &title) const;

    int                  n_refinements_global;

    MPI_Comm             mpi_communicator;

    parallel::distributed::Triangulation<dim>   triangulation;

    FESystem<dim>        fe;
    DoFHandler<dim>      dof_handler;

    IndexSet             locally_owned_dofs;
    IndexSet             locally_relevant_dofs;

    unsigned int         number_iterations;
    std::vector<double>  run_time;

    ConstraintMatrix     constraints;
    ConstraintMatrix     constraints_hanging_nodes;
    ConstraintMatrix     constraints_dirichlet_hanging_nodes;

    TrilinosWrappers::SparseMatrix system_matrix_newton;

    TrilinosWrappers::MPI::Vector       solution;
    TrilinosWrappers::MPI::Vector       old_solution;
    TrilinosWrappers::MPI::Vector       system_rhs_newton;
    TrilinosWrappers::MPI::Vector       resid_vector;
    TrilinosWrappers::MPI::Vector       diag_mass_matrix_vector;
    Vector<float>                       cell_constitution;
    IndexSet                            active_set;

    ConditionalOStream pcout;

    TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
    TrilinosWrappers::PreconditionAMG preconditioner_u;
    TrilinosWrappers::PreconditionAMG preconditioner_t;

    std::unique_ptr<Input<dim> >               input_obstacle;
    std::unique_ptr<ConstitutiveLaw<dim> >     plast_lin_hard;

    double sigma_0;    // Yield stress
    double gamma;      // Parameter for the linear isotropic hardening
    double e_modul;    // E-Modul
    double nu;         // Poisson ratio
  };

  // As explained above this class
  // allocates the obstacle which
  // will come into contact with our
  // deformable body.
  // The data which we read in by the
  // function read_obstacle () from the file
  // "obstacle_file.pbm" will be stored
  // in a double std::vector named
  // obstacle_data.
  // This vector composes the base
  // to calculate a piecewise bilinear
  // function as a polynomial interpolation.
  // This will be done by obstacle_function ().
  //
  // In the function run () of the class
  // PlasticityContactProblem we create
  // an object of the class Input which will
  // be used in the class Obstacle to
  // supply the obstacle function in
  // update_solution_and_constraints () of
  // the class PlasticityContactProblem.

  template <int dim>
  class Input
  {
   public:
    Input (const char* _name) :
      name (_name),
      mpi_communicator (MPI_COMM_WORLD),
      pcout (std::cout,
               (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
      obstacle_data (0),
      lx (1.0), // length of the cube in x direction
      ly (1.0), // length of the cube in y direction
      hx (0), //
      hy (0),
      nx (0),
      ny (0)
      {read_obstacle (name);}

      double hv (int i, int j);

      double obstacle_function (double x,double y);

      void read_obstacle (const char* name);

   private:
      const char*          name;
      MPI_Comm             mpi_communicator;
      ConditionalOStream   pcout;
      std::vector<double>  obstacle_data;
      double               lx, ly;
      double               hx, hy;
      int                  nx, ny;
  };

  template <int dim>
  double Input<dim>::hv (int i, int j)
  {
    assert(i>=0 && i<nx);
    assert(j>=0 && j<ny);
    return obstacle_data[nx*(ny-1-j)+i]; // i indiziert x-werte, j indiziert y-werte
  }

  template <int dim>
  double Input<dim>::obstacle_function (double x,double y)
  {
    int ix = (int)(x/hx);
    int iy = (int)(y/hy);

    if (ix<0)
      ix = 0;

    if (iy<0)
      iy = 0;

    if (ix>=nx-1)
      ix = nx-2;

    if (iy>=ny-1)
      iy = ny-2;

    double val = 0.0;
    {
      FullMatrix<double> H(4,4);
      Vector<double>  X(4);
      Vector<double>  b(4);

      double xx = 0.0;
      double yy = 0.0;

      xx = ix*hx;
      yy = iy*hy;
      H(0,0) = xx;
      H(0,1) = yy;
      H(0,2) = xx*yy;
      H(0,3) = 1.0;
      b(0)   = hv (ix, iy);

      xx = (ix + 1)*hx;
      yy = iy*hy;
      H(1,0) = xx;
      H(1,1) = yy;
      H(1,2) = xx*yy;
      H(1,3) = 1.0;
      b(1)   = hv (ix + 1, iy);

      xx = (ix + 1)*hx;
      yy = (iy + 1)*hy;
      H(2,0) = xx;
      H(2,1) = yy;
      H(2,2) = xx*yy;
      H(2,3) = 1.0;
      b(2)   = hv (ix + 1, iy + 1);

      xx = ix*hx;
      yy = (iy + 1)*hy;
      H(3,0) = xx;
      H(3,1) = yy;
      H(3,2) = xx*yy;
      H(3,3) = 1.0;
      b(3)   = hv (ix, iy + 1);

      H.gauss_jordan ();
      H.vmult (X, b);

      val = X(0)*x + X(1)*y + X(2)*x*y + X(3);
    }

    return val;
  }

  template <int dim>
  void Input<dim>::read_obstacle (const char* name)
  {
    int SZ = 100000;
    FILE* fp = fopen (name, "r");
    char* hlp_str = new char[SZ];
    double hlp;

    fscanf (fp, "%s", hlp_str);
    fscanf (fp, "%d", &nx);
    fscanf (fp, "%d", &ny);

    assert(nx>0 && ny>0);

    for (int k=0; k<nx*ny; k++)
      {
        fscanf (fp, "%lf", &hlp);
        obstacle_data.push_back (hlp);
      }

    fclose (fp);
    delete[] hlp_str;

    hx = 1.0/(nx - 1);
    hy = 1.0/(ny - 1);

    pcout << "Resolution of the scanned obstacle picture: " << nx << " x " << ny << std::endl;
  }

  template <int dim>
  class ConstitutiveLaw
  {
  public:
    ConstitutiveLaw (double _E,
                     double _nu,
                     double _sigma_0,
                     double _gamma,
                     MPI_Comm _mpi_communicator,
                     ConditionalOStream _pcout);

    void plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor,
                                 SymmetricTensor<2,dim> &strain_tensor,
                                 unsigned int            &elast_points,
                                 unsigned int            &plast_points,
                                 double                  &yield);
    void linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized,
                                            SymmetricTensor<4,dim> &stress_strain_tensor,
                                            SymmetricTensor<2,dim> &strain_tensor);
    inline SymmetricTensor<2,dim> get_strain (const FEValues<dim> &fe_values,
                                              const unsigned int  shape_func,
                                              const unsigned int  q_point) const;
    void set_sigma_0 (double sigma_hlp) {sigma_0 = sigma_hlp;}

  private:
    SymmetricTensor<4,dim>  stress_strain_tensor_mu;
    SymmetricTensor<4,dim>  stress_strain_tensor_kappa;
    double E;
    double nu;
    double sigma_0;
    double gamma;
    double mu;
    double kappa;
    MPI_Comm mpi_communicator;
    ConditionalOStream pcout;
  };

  template <int dim>
  ConstitutiveLaw<dim>::ConstitutiveLaw(double _E, double _nu, double _sigma_0, double _gamma, MPI_Comm _mpi_communicator, ConditionalOStream _pcout)
    :E (_E),
     nu (_nu),
     sigma_0 (_sigma_0),
     gamma (_gamma),
     mpi_communicator (_mpi_communicator),
     pcout (_pcout)
  {
    mu = E/(2*(1+nu));
    kappa = E/(3*(1-2*nu));
    stress_strain_tensor_kappa = kappa*outer_product(unit_symmetric_tensor<dim>(), unit_symmetric_tensor<dim>());
    stress_strain_tensor_mu = 2*mu*(identity_tensor<dim>() - outer_product(unit_symmetric_tensor<dim>(), unit_symmetric_tensor<dim>())/3.0);
  }

  template <int dim>
  inline
  SymmetricTensor<2,dim> ConstitutiveLaw<dim>::get_strain (const FEValues<dim> &fe_values,
                                                           const unsigned int   shape_func,
                                                           const unsigned int   q_point) const
  {
    const FEValuesExtractors::Vector displacement (0);
    SymmetricTensor<2,dim> tmp;

    tmp = fe_values[displacement].symmetric_gradient (shape_func,q_point);

    return tmp;
  }

  template <int dim>
  void ConstitutiveLaw<dim>::plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor,
                                                     SymmetricTensor<2,dim> &strain_tensor,
                                                     unsigned int            &elast_points,
                                                     unsigned int            &plast_points,
                                                     double                  &yield)
  {
    if (dim == 3)
      {
        SymmetricTensor<2,dim> stress_tensor;
        stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor;

        SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor);

        double deviator_stress_tensor_norm = deviator_stress_tensor.norm ();

        yield = 0;
        stress_strain_tensor = stress_strain_tensor_mu;
        double beta = 1.0;
        if (deviator_stress_tensor_norm >
            sigma_0)
          {
//            beta = (sigma_0 + gamma/(2*mu + gamma)*(deviator_stress_tensor_norm - sigma_0))/deviator_stress_tensor_norm;
            beta = sigma_0/deviator_stress_tensor_norm;
            stress_strain_tensor *= (gamma + (1 - gamma)*beta);
            yield = 1;
            plast_points += 1;
          }
        else
          elast_points += 1;

        stress_strain_tensor += stress_strain_tensor_kappa;
      }
  }

  template <int dim>
  void ConstitutiveLaw<dim>::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized,
      SymmetricTensor<4,dim> &stress_strain_tensor,
      SymmetricTensor<2,dim> &strain_tensor)
  {
    if (dim == 3)
      {
        SymmetricTensor<2,dim> stress_tensor;
        stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor;

        SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor);

        double deviator_stress_tensor_norm = deviator_stress_tensor.norm ();

        stress_strain_tensor = stress_strain_tensor_mu;
        stress_strain_tensor_linearized = stress_strain_tensor_mu;
        double beta = 1.0;
        if (deviator_stress_tensor_norm >
            sigma_0)
          {
//            beta = (sigma_0 + gamma/(2*mu + gamma)*(deviator_stress_tensor_norm - sigma_0))/deviator_stress_tensor_norm;
            beta = sigma_0/deviator_stress_tensor_norm;
            stress_strain_tensor *= (gamma + (1 - gamma)*beta);
            stress_strain_tensor_linearized *= (gamma + (1 - gamma)*beta);
            deviator_stress_tensor /= deviator_stress_tensor_norm;
            stress_strain_tensor_linearized -= (1 - gamma)*beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor);

//            std::cout << "Plastisch mit beta = " << beta
//                << ", sigma_0 = " << sigma_0
//                << std::endl;
          }

        stress_strain_tensor += stress_strain_tensor_kappa;
        stress_strain_tensor_linearized += stress_strain_tensor_kappa;
      }
  }

  namespace EquationData
  {
    template <int dim>
    class RightHandSide : public Function<dim>
    {
    public:
      RightHandSide () : Function<dim>(dim) {}

      virtual double value (const Point<dim>   &p,
                            const unsigned int  component = 0) const;

      virtual void vector_value (const Point<dim> &p,
                                 Vector<double>   &values) const;
    };

    template <int dim>
    double RightHandSide<dim>::value (const Point<dim> &p,
                                      const unsigned int component) const
    {
      double return_value = 0.0;

      if (component == 0)
        return_value = 0.0;
      if (component == 1)
        return_value = 0.0;
      if (component == 2)
        // if ((p(0)-0.5)*(p(0)-0.5)+(p(1)-0.5)*(p(1)-0.5) < 0.2)
        //  return_value = -5000;
        // else
        return_value = 0.0;
      // for (unsigned int i=0; i<dim; ++i)
      //   return_value += 4*std::pow(p(i), 4);

      return return_value;
    }

    template <int dim>
    void RightHandSide<dim>::vector_value (const Point<dim> &p,
                                           Vector<double>   &values) const
    {
      for (unsigned int c=0; c<this->n_components; ++c)
        values(c) = RightHandSide<dim>::value (p, c);
    }


    template <int dim>
    class BoundaryValues : public Function<dim>
    {
    public:
      BoundaryValues () : Function<dim>(dim) {};

      virtual double value (const Point<dim>   &p,
                            const unsigned int  component = 0) const;

      virtual void vector_value (const Point<dim> &p,
                                 Vector<double>   &values) const;
    };

    template <int dim>
    double BoundaryValues<dim>::value (const Point<dim> &p,
                                       const unsigned int component) const
    {
      double return_value = 0;

      if (component == 0)
        return_value = 0.0;
      if (component == 1)
        return_value = 0.0;
      if (component == 2)
        return_value = 0.0;

      return return_value;
    }

    template <int dim>
    void BoundaryValues<dim>::vector_value (const Point<dim> &p,
                                            Vector<double>   &values) const
    {
      for (unsigned int c=0; c<this->n_components; ++c)
        values(c) = BoundaryValues<dim>::value (p, c);
    }


    template <int dim>
    class Obstacle : public Function<dim>
    {
    public:
      Obstacle (std::unique_ptr<Input<dim> > const &_input) :
        Function<dim>(dim),
        input_obstacle_copy(std::move (_input)) {};

      virtual double value (const Point<dim>   &p,
                            const unsigned int  component = 0) const;

      virtual void vector_value (const Point<dim> &p,
                                 Vector<double>   &values) const;

    private:
      std::unique_ptr<Input<dim> >  const &input_obstacle_copy;
    };

    template <int dim>
    double Obstacle<dim>::value (const Point<dim> &p,
                                 const unsigned int component) const
    {
      double R = 0.03;
      double return_value = 100.0;
      if (component == 0)
        return_value = p(0);
      if (component == 1)
        return_value = p(1);
      if (component == 2)
        {
          return_value = -std::sqrt (0.36 - (p(0)-0.5)*(p(0)-0.5) - (p(1)-0.5)*(p(1)-0.5)) + 1.59;
//          return_value = 1.999 - input_obstacle_copy->obstacle_function (p(0), p(1));
        }
      return return_value;
    }

    template <int dim>
    void Obstacle<dim>::vector_value (const Point<dim> &p,
                                      Vector<double> &values) const
    {
      for (unsigned int c=0; c<this->n_components; ++c)
        values(c) = Obstacle<dim>::value (p, c);
    }
  }


  // @sect3{Implementation of the <code>PlasticityContactProblem</code> class}

  // Next for the implementation of the class
  // template that makes use of the functions
  // above. As before, we will write everything

  template <int dim>
  PlasticityContactProblem<dim>::PlasticityContactProblem (int _n_refinements_global)
    :
    n_refinements_global (_n_refinements_global),
    mpi_communicator (MPI_COMM_WORLD),
    triangulation (mpi_communicator),
    fe (FE_Q<dim>(1), dim),
    dof_handler (triangulation),
    pcout (std::cout,
           (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
    sigma_0 (400),
    gamma (0.01),
    e_modul (2.0e+5),
    nu (0.3)
  {
    // double _E, double _nu, double _sigma_0, double _gamma
    plast_lin_hard.reset (new ConstitutiveLaw<dim> (e_modul, nu, sigma_0, gamma, mpi_communicator, pcout));
  }

  template <int dim>
  void PlasticityContactProblem<dim>::make_grid ()
  {
    std::vector<unsigned int> repet(3);
    repet[0] = 1;//20;
    repet[1] = 1;
    repet[2] = 1;

    Point<dim> p1 (0,0,0);
    Point<dim> p2 (1.0, 1.0, 1.0);
    GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2);

    Triangulation<3>::active_cell_iterator
    cell = triangulation.begin_active(),
    endc = triangulation.end();

    /* boundary_indicators:
              _______
             /  9    /|
            /______ / |
          8|       | 8|
           |   8   | /
           |_______|/
               6
     */

    for (; cell!=endc; ++cell)
      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
        {
          if (cell->face (face)->center ()[2] == p2(2))
            cell->face (face)->set_boundary_indicator (9);
          if (cell->face (face)->center ()[0] == p1(0) ||
              cell->face (face)->center ()[0] == p2(0) ||
              cell->face (face)->center ()[1] == p1(1) ||
              cell->face (face)->center ()[1] == p2(1))
            cell->face (face)->set_boundary_indicator (8);
          if (cell->face (face)->center ()[2] == p1(2))
            cell->face (face)->set_boundary_indicator (6);
        }

    triangulation.refine_global (n_refinements_global);
  }

  template <int dim>
  void PlasticityContactProblem<dim>::setup_system ()
  {
    // setup dofs
    {
      dof_handler.distribute_dofs (fe);

      locally_owned_dofs = dof_handler.locally_owned_dofs ();
      locally_relevant_dofs.clear();
      DoFTools::extract_locally_relevant_dofs (dof_handler,
                                               locally_relevant_dofs);
    }

    // setup hanging nodes and dirichlet constraints
    {
      constraints_hanging_nodes.clear ();
      constraints_hanging_nodes.reinit (locally_relevant_dofs);
      DoFTools::make_hanging_node_constraints (dof_handler,
                                               constraints_hanging_nodes);
      constraints_hanging_nodes.close ();

      pcout << "   Number of active cells: "
            << triangulation.n_active_cells()
            << std::endl
            << "   Total number of cells: "
            << triangulation.n_cells()
            << std::endl
            << "   Number of degrees of freedom: "
            << dof_handler.n_dofs ()
            << std::endl;

      dirichlet_constraints ();
    }

    // Initialzation for matrices and vectors
    {
      solution.reinit (locally_relevant_dofs, mpi_communicator);
      system_rhs_newton.reinit (locally_owned_dofs, mpi_communicator);
      old_solution.reinit (system_rhs_newton);
      resid_vector.reinit (system_rhs_newton);
      diag_mass_matrix_vector.reinit (system_rhs_newton);
      cell_constitution.reinit (triangulation.n_active_cells ());
      active_set.clear ();
      active_set.set_size (locally_relevant_dofs.size ());
    }

    // setup sparsity pattern
    {
      TrilinosWrappers::SparsityPattern sp (locally_owned_dofs,
                                            mpi_communicator);

      DoFTools::make_sparsity_pattern (dof_handler, sp, constraints_dirichlet_hanging_nodes, false,
                                       Utilities::MPI::this_mpi_process(mpi_communicator));

      sp.compress();

      system_matrix_newton.reinit (sp);

      TrilinosWrappers::SparseMatrix mass_matrix;
      mass_matrix.reinit (sp);
      assemble_mass_matrix_diagonal (mass_matrix);
      const unsigned int
      start = (system_rhs_newton.local_range().first),
      end   = (system_rhs_newton.local_range().second);
      for (unsigned int j=start; j<end; j++)
        diag_mass_matrix_vector (j) = mass_matrix.diag_element (j);
      number_iterations = 0;

      diag_mass_matrix_vector.compress (VectorOperation::insert);
    }
  }

  template <int dim>
  void PlasticityContactProblem<dim>::assemble_nl_system (TrilinosWrappers::MPI::Vector &u)
  {
    QGauss<dim>  quadrature_formula(2);
    QGauss<dim-1>  face_quadrature_formula(2);

    FEValues<dim> fe_values (fe, quadrature_formula,
                             UpdateFlags(update_values    |
                                         update_gradients |
                                         update_q_points  |
                                         update_JxW_values));

    FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
                                      update_values   | update_quadrature_points |
                                      update_JxW_values);

    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
    const unsigned int   n_q_points    = quadrature_formula.size ();
    const unsigned int   n_face_q_points = face_quadrature_formula.size();

    const EquationData::RightHandSide<dim> right_hand_side;
    std::vector<Vector<double> > right_hand_side_values (n_q_points,
                                                         Vector<double>(dim));
    std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
                                                              Vector<double>(dim));

    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
    Vector<double>       cell_rhs (dofs_per_cell);

    std::vector<unsigned int> local_dof_indices (dofs_per_cell);

    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                                   endc = dof_handler.end();

    const FEValuesExtractors::Vector displacement (0);

    const double kappa = 1.0;
    for (; cell!=endc; ++cell)
      if (cell->is_locally_owned())
        {
          fe_values.reinit (cell);
          cell_matrix = 0;
          cell_rhs = 0;

          right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
                                             right_hand_side_values);

          std::vector<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
          fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);

          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
            {
              SymmetricTensor<4,dim> stress_strain_tensor_linearized;
              SymmetricTensor<4,dim> stress_strain_tensor;
              SymmetricTensor<2,dim> stress_tensor;

              plast_lin_hard->linearized_plast_linear_hardening (stress_strain_tensor_linearized,
                                                                 stress_strain_tensor,
                                                                 strain_tensor[q_point]);

              for (unsigned int i=0; i<dofs_per_cell; ++i)
                {
                  stress_tensor = stress_strain_tensor_linearized * plast_lin_hard->get_strain(fe_values, i, q_point);

                  for (unsigned int j=0; j<dofs_per_cell; ++j)
                    {
                      cell_matrix(i,j) += (stress_tensor *
                                           plast_lin_hard->get_strain(fe_values, j, q_point) *
                                           fe_values.JxW (q_point));
                    }

                  // the linearized part a(v^i;v^i,v) of the rhs
                  cell_rhs(i) += (stress_tensor *
                                  strain_tensor[q_point] *
                                  fe_values.JxW (q_point));

                  // the residual part a(v^i;v) of the rhs
                  cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor *
                                  plast_lin_hard->get_strain(fe_values, i, q_point) *
                                  fe_values.JxW (q_point));

                  // the residual part F(v) of the rhs
                  Tensor<1,dim> rhs_values;
                  rhs_values = 0;
                  cell_rhs(i) += (fe_values[displacement].value (i, q_point) *
                                  rhs_values *
                                  fe_values.JxW (q_point));
                }
            }

          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (cell->face (face)->at_boundary()
                  && cell->face (face)->boundary_indicator () == 9)
                {
                  fe_values_face.reinit (cell, face);

                  right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(),
                                                     right_hand_side_values_face);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      Tensor<1,dim> rhs_values;
                      rhs_values = 0;
                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
                                        rhs_values *
                                        fe_values_face.JxW (q_point));
                    }
                }
            }

          cell->get_dof_indices (local_dof_indices);
          constraints.distribute_local_to_global (cell_matrix, cell_rhs,
                                                  local_dof_indices,
                                                  system_matrix_newton, system_rhs_newton, true);
        };

    system_matrix_newton.compress ();
    system_rhs_newton.compress (VectorOperation::add);
  }

  template <int dim>
  void PlasticityContactProblem<dim>::residual_nl_system (TrilinosWrappers::MPI::Vector &u)
  {
    QGauss<dim>  quadrature_formula(2);
    QGauss<dim-1> face_quadrature_formula(2);

    FEValues<dim> fe_values (fe, quadrature_formula,
                             UpdateFlags(update_values    |
                                         update_gradients |
                                         update_q_points  |
                                         update_JxW_values));

    FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
                                      update_values   | update_quadrature_points |
                                      update_JxW_values);

    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
    const unsigned int   n_q_points    = quadrature_formula.size ();
    const unsigned int   n_face_q_points = face_quadrature_formula.size();

    const EquationData::RightHandSide<dim> right_hand_side;
    std::vector<Vector<double> > right_hand_side_values (n_q_points,
                                                         Vector<double>(dim));
    std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
                                                              Vector<double>(dim));

    Vector<double>       cell_rhs (dofs_per_cell);

    std::vector<unsigned int> local_dof_indices (dofs_per_cell);

    const FEValuesExtractors::Vector displacement (0);

    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                                   endc = dof_handler.end();

    unsigned int elast_points = 0;
    unsigned int plast_points = 0;
    double       yield = 0;
    unsigned int cell_number = 0;
    cell_constitution = 0;

    for (; cell!=endc; ++cell)
      if (cell->is_locally_owned())
        {
          fe_values.reinit (cell);
          cell_rhs = 0;

          right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
                                             right_hand_side_values);

          std::vector<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
          fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);

          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
            {
              SymmetricTensor<4,dim> stress_strain_tensor;
              SymmetricTensor<2,dim> stress_tensor;

              plast_lin_hard->plast_linear_hardening (stress_strain_tensor, strain_tensor[q_point],
                                                      elast_points, plast_points, yield);
              cell_constitution (cell_number) += yield;
              for (unsigned int i=0; i<dofs_per_cell; ++i)
                {
                  cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * //(stress_tensor) *
                                  plast_lin_hard->get_strain(fe_values, i, q_point) *
                                  fe_values.JxW (q_point));

                  Tensor<1,dim> rhs_values;
                  rhs_values = 0;
                  cell_rhs(i) += ((fe_values[displacement].value (i, q_point) *
                                   rhs_values) *
                                  fe_values.JxW (q_point));
                };
            };

          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (cell->face (face)->at_boundary()
                  && cell->face (face)->boundary_indicator () == 9)
                {
                  fe_values_face.reinit (cell, face);

                  right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(),
                                                     right_hand_side_values_face);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      Tensor<1,dim> rhs_values;
                      rhs_values = 0;
                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
                                        rhs_values *
                                        fe_values_face.JxW (q_point));
                    }
                }
            }

          cell->get_dof_indices (local_dof_indices);
          constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs,
              local_dof_indices,
              system_rhs_newton);

          cell_number += 1;
        };

    cell_constitution /= n_q_points;
    system_rhs_newton.compress (VectorOperation::add);

    unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, mpi_communicator);
    unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, mpi_communicator);
    pcout << "      Number of elastic quadrature points: " << sum_elast_points
          << " and plastic quadrature points: " << sum_plast_points <<std::endl;
  }

  template <int dim>
  void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal (TrilinosWrappers::SparseMatrix &mass_matrix)
  {
    QTrapez<dim-1>  face_quadrature_formula;

    FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
                                      update_values   |
                                      update_quadrature_points |
                                      update_JxW_values);

    const unsigned int   dofs_per_cell      = fe.dofs_per_cell;
    const unsigned int   n_face_q_points    = face_quadrature_formula.size();

    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);

    std::vector<unsigned int> local_dof_indices (dofs_per_cell);

    const FEValuesExtractors::Vector displacement (0);

    typename DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();
<commit_msg>fix the computation of cell_constitution for the parallel case<commit_after>/* $Id$ */
/* Author: <NAME>, University of Heidelberg, 1999 */

/*    $Id$       */
/*                                                                */
/*    Copyright (C) 1999-2012 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyrightG and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

// @sect3{Include files}

// We are using the the same
// include files as in step-41:

#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.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/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.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/block_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>

#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_block_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/timer.h>
#include <fstream>
#include <iostream>
#include <list>
#include <time.h>


#include <deal.II/base/logstream.h>

namespace Step42
{
  using namespace dealii;


  // @sect3{The <code>PlasticityContactProblem</code> class template}

  // This class has only the purpose
  // to read in data from a picture file
  // that has to be stored as a pbm ascii
  // format. This data will be bilinear
  // interpolated and provides in this way
  // a function which describes an obstacle.

  template <int dim> class Input;

  // This class provides an interface
  // for a constitutive law. In this
  // example we are using an elastic
  // plastic material with linear,
  // isotropic hardening.

  template <int dim> class ConstitutiveLaw;

  // This class supplies all function
  // and variables needed to describe
  // the nonlinear contact problem. It is
  // close to step-41 but with some additional
  // features like: handling hanging nodes,
  // a newton method, using Trilinos and p4est
  // for parallel distributed computing.
  // To deal with hanging nodes makes
  // life a bit more complicated since
  // we need an other ConstraintMatrix now.
  // We create a newton method for the
  // active set method for the contact
  // situation and to handle the nonlinear
  // operator for the constitutive law.

  template <int dim>
  class PlasticityContactProblem
  {
  public:
    PlasticityContactProblem (int _n_refinements_global);
    void run ();

  private:
    void make_grid ();
    void setup_system();
    void assemble_nl_system (TrilinosWrappers::MPI::Vector &u);
    void residual_nl_system (TrilinosWrappers::MPI::Vector &u);
    void assemble_mass_matrix_diagonal (TrilinosWrappers::SparseMatrix &mass_matrix);
    void update_solution_and_constraints ();
    void dirichlet_constraints ();
    void solve ();
    void solve_newton ();
    void refine_grid ();
    void move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const;
    void output_results (const std::string &title) const;

    int                  n_refinements_global;

    MPI_Comm             mpi_communicator;

    parallel::distributed::Triangulation<dim>   triangulation;

    FESystem<dim>        fe;
    DoFHandler<dim>      dof_handler;

    IndexSet             locally_owned_dofs;
    IndexSet             locally_relevant_dofs;

    unsigned int         number_iterations;
    std::vector<double>  run_time;

    ConstraintMatrix     constraints;
    ConstraintMatrix     constraints_hanging_nodes;
    ConstraintMatrix     constraints_dirichlet_hanging_nodes;

    TrilinosWrappers::SparseMatrix system_matrix_newton;

    TrilinosWrappers::MPI::Vector       solution;
    TrilinosWrappers::MPI::Vector       old_solution;
    TrilinosWrappers::MPI::Vector       system_rhs_newton;
    TrilinosWrappers::MPI::Vector       resid_vector;
    TrilinosWrappers::MPI::Vector       diag_mass_matrix_vector;
    Vector<float>                       cell_constitution;
    IndexSet                            active_set;

    ConditionalOStream pcout;

    TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
    TrilinosWrappers::PreconditionAMG preconditioner_u;
    TrilinosWrappers::PreconditionAMG preconditioner_t;

    std::unique_ptr<Input<dim> >               input_obstacle;
    std::unique_ptr<ConstitutiveLaw<dim> >     plast_lin_hard;

    double sigma_0;    // Yield stress
    double gamma;      // Parameter for the linear isotropic hardening
    double e_modul;    // E-Modul
    double nu;         // Poisson ratio
  };

  // As explained above this class
  // allocates the obstacle which
  // will come into contact with our
  // deformable body.
  // The data which we read in by the
  // function read_obstacle () from the file
  // "obstacle_file.pbm" will be stored
  // in a double std::vector named
  // obstacle_data.
  // This vector composes the base
  // to calculate a piecewise bilinear
  // function as a polynomial interpolation.
  // This will be done by obstacle_function ().
  //
  // In the function run () of the class
  // PlasticityContactProblem we create
  // an object of the class Input which will
  // be used in the class Obstacle to
  // supply the obstacle function in
  // update_solution_and_constraints () of
  // the class PlasticityContactProblem.

  template <int dim>
  class Input
  {
   public:
    Input (const char* _name) :
      name (_name),
      mpi_communicator (MPI_COMM_WORLD),
      pcout (std::cout,
               (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
      obstacle_data (0),
      lx (1.0), // length of the cube in x direction
      ly (1.0), // length of the cube in y direction
      hx (0), //
      hy (0),
      nx (0),
      ny (0)
      {read_obstacle (name);}

      double hv (int i, int j);

      double obstacle_function (double x,double y);

      void read_obstacle (const char* name);

   private:
      const char*          name;
      MPI_Comm             mpi_communicator;
      ConditionalOStream   pcout;
      std::vector<double>  obstacle_data;
      double               lx, ly;
      double               hx, hy;
      int                  nx, ny;
  };

  template <int dim>
  double Input<dim>::hv (int i, int j)
  {
    assert(i>=0 && i<nx);
    assert(j>=0 && j<ny);
    return obstacle_data[nx*(ny-1-j)+i]; // i indiziert x-werte, j indiziert y-werte
  }

  template <int dim>
  double Input<dim>::obstacle_function (double x,double y)
  {
    int ix = (int)(x/hx);
    int iy = (int)(y/hy);

    if (ix<0)
      ix = 0;

    if (iy<0)
      iy = 0;

    if (ix>=nx-1)
      ix = nx-2;

    if (iy>=ny-1)
      iy = ny-2;

    double val = 0.0;
    {
      FullMatrix<double> H(4,4);
      Vector<double>  X(4);
      Vector<double>  b(4);

      double xx = 0.0;
      double yy = 0.0;

      xx = ix*hx;
      yy = iy*hy;
      H(0,0) = xx;
      H(0,1) = yy;
      H(0,2) = xx*yy;
      H(0,3) = 1.0;
      b(0)   = hv (ix, iy);

      xx = (ix + 1)*hx;
      yy = iy*hy;
      H(1,0) = xx;
      H(1,1) = yy;
      H(1,2) = xx*yy;
      H(1,3) = 1.0;
      b(1)   = hv (ix + 1, iy);

      xx = (ix + 1)*hx;
      yy = (iy + 1)*hy;
      H(2,0) = xx;
      H(2,1) = yy;
      H(2,2) = xx*yy;
      H(2,3) = 1.0;
      b(2)   = hv (ix + 1, iy + 1);

      xx = ix*hx;
      yy = (iy + 1)*hy;
      H(3,0) = xx;
      H(3,1) = yy;
      H(3,2) = xx*yy;
      H(3,3) = 1.0;
      b(3)   = hv (ix, iy + 1);

      H.gauss_jordan ();
      H.vmult (X, b);

      val = X(0)*x + X(1)*y + X(2)*x*y + X(3);
    }

    return val;
  }

  template <int dim>
  void Input<dim>::read_obstacle (const char* name)
  {
    int SZ = 100000;
    FILE* fp = fopen (name, "r");
    char* hlp_str = new char[SZ];
    double hlp;

    fscanf (fp, "%s", hlp_str);
    fscanf (fp, "%d", &nx);
    fscanf (fp, "%d", &ny);

    assert(nx>0 && ny>0);

    for (int k=0; k<nx*ny; k++)
      {
        fscanf (fp, "%lf", &hlp);
        obstacle_data.push_back (hlp);
      }

    fclose (fp);
    delete[] hlp_str;

    hx = 1.0/(nx - 1);
    hy = 1.0/(ny - 1);

    pcout << "Resolution of the scanned obstacle picture: " << nx << " x " << ny << std::endl;
  }

  template <int dim>
  class ConstitutiveLaw
  {
  public:
    ConstitutiveLaw (double _E,
                     double _nu,
                     double _sigma_0,
                     double _gamma,
                     MPI_Comm _mpi_communicator,
                     ConditionalOStream _pcout);

    void plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor,
                                 SymmetricTensor<2,dim> &strain_tensor,
                                 unsigned int            &elast_points,
                                 unsigned int            &plast_points,
                                 double                  &yield);
    void linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized,
                                            SymmetricTensor<4,dim> &stress_strain_tensor,
                                            SymmetricTensor<2,dim> &strain_tensor);
    inline SymmetricTensor<2,dim> get_strain (const FEValues<dim> &fe_values,
                                              const unsigned int  shape_func,
                                              const unsigned int  q_point) const;
    void set_sigma_0 (double sigma_hlp) {sigma_0 = sigma_hlp;}

  private:
    SymmetricTensor<4,dim>  stress_strain_tensor_mu;
    SymmetricTensor<4,dim>  stress_strain_tensor_kappa;
    double E;
    double nu;
    double sigma_0;
    double gamma;
    double mu;
    double kappa;
    MPI_Comm mpi_communicator;
    ConditionalOStream pcout;
  };

  template <int dim>
  ConstitutiveLaw<dim>::ConstitutiveLaw(double _E, double _nu, double _sigma_0, double _gamma, MPI_Comm _mpi_communicator, ConditionalOStream _pcout)
    :E (_E),
     nu (_nu),
     sigma_0 (_sigma_0),
     gamma (_gamma),
     mpi_communicator (_mpi_communicator),
     pcout (_pcout)
  {
    mu = E/(2*(1+nu));
    kappa = E/(3*(1-2*nu));
    stress_strain_tensor_kappa = kappa*outer_product(unit_symmetric_tensor<dim>(), unit_symmetric_tensor<dim>());
    stress_strain_tensor_mu = 2*mu*(identity_tensor<dim>() - outer_product(unit_symmetric_tensor<dim>(), unit_symmetric_tensor<dim>())/3.0);
  }

  template <int dim>
  inline
  SymmetricTensor<2,dim> ConstitutiveLaw<dim>::get_strain (const FEValues<dim> &fe_values,
                                                           const unsigned int   shape_func,
                                                           const unsigned int   q_point) const
  {
    const FEValuesExtractors::Vector displacement (0);
    SymmetricTensor<2,dim> tmp;

    tmp = fe_values[displacement].symmetric_gradient (shape_func,q_point);

    return tmp;
  }

  template <int dim>
  void ConstitutiveLaw<dim>::plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor,
                                                     SymmetricTensor<2,dim> &strain_tensor,
                                                     unsigned int            &elast_points,
                                                     unsigned int            &plast_points,
                                                     double                  &yield)
  {
    if (dim == 3)
      {
        SymmetricTensor<2,dim> stress_tensor;
        stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor;

        SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor);

        double deviator_stress_tensor_norm = deviator_stress_tensor.norm ();

        yield = 0;
        stress_strain_tensor = stress_strain_tensor_mu;
        double beta = 1.0;
        if (deviator_stress_tensor_norm >
            sigma_0)
          {
//            beta = (sigma_0 + gamma/(2*mu + gamma)*(deviator_stress_tensor_norm - sigma_0))/deviator_stress_tensor_norm;
            beta = sigma_0/deviator_stress_tensor_norm;
            stress_strain_tensor *= (gamma + (1 - gamma)*beta);
            yield = 1;
            plast_points += 1;
          }
        else
          elast_points += 1;

        stress_strain_tensor += stress_strain_tensor_kappa;
      }
  }

  template <int dim>
  void ConstitutiveLaw<dim>::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized,
      SymmetricTensor<4,dim> &stress_strain_tensor,
      SymmetricTensor<2,dim> &strain_tensor)
  {
    if (dim == 3)
      {
        SymmetricTensor<2,dim> stress_tensor;
        stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor;

        SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor);

        double deviator_stress_tensor_norm = deviator_stress_tensor.norm ();

        stress_strain_tensor = stress_strain_tensor_mu;
        stress_strain_tensor_linearized = stress_strain_tensor_mu;
        double beta = 1.0;
        if (deviator_stress_tensor_norm >
            sigma_0)
          {
//            beta = (sigma_0 + gamma/(2*mu + gamma)*(deviator_stress_tensor_norm - sigma_0))/deviator_stress_tensor_norm;
            beta = sigma_0/deviator_stress_tensor_norm;
            stress_strain_tensor *= (gamma + (1 - gamma)*beta);
            stress_strain_tensor_linearized *= (gamma + (1 - gamma)*beta);
            deviator_stress_tensor /= deviator_stress_tensor_norm;
            stress_strain_tensor_linearized -= (1 - gamma)*beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor);

//            std::cout << "Plastisch mit beta = " << beta
//                << ", sigma_0 = " << sigma_0
//                << std::endl;
          }

        stress_strain_tensor += stress_strain_tensor_kappa;
        stress_strain_tensor_linearized += stress_strain_tensor_kappa;
      }
  }

  namespace EquationData
  {
    template <int dim>
    class RightHandSide : public Function<dim>
    {
    public:
      RightHandSide () : Function<dim>(dim) {}

      virtual double value (const Point<dim>   &p,
                            const unsigned int  component = 0) const;

      virtual void vector_value (const Point<dim> &p,
                                 Vector<double>   &values) const;
    };

    template <int dim>
    double RightHandSide<dim>::value (const Point<dim> &p,
                                      const unsigned int component) const
    {
      double return_value = 0.0;

      if (component == 0)
        return_value = 0.0;
      if (component == 1)
        return_value = 0.0;
      if (component == 2)
        // if ((p(0)-0.5)*(p(0)-0.5)+(p(1)-0.5)*(p(1)-0.5) < 0.2)
        //  return_value = -5000;
        // else
        return_value = 0.0;
      // for (unsigned int i=0; i<dim; ++i)
      //   return_value += 4*std::pow(p(i), 4);

      return return_value;
    }

    template <int dim>
    void RightHandSide<dim>::vector_value (const Point<dim> &p,
                                           Vector<double>   &values) const
    {
      for (unsigned int c=0; c<this->n_components; ++c)
        values(c) = RightHandSide<dim>::value (p, c);
    }


    template <int dim>
    class BoundaryValues : public Function<dim>
    {
    public:
      BoundaryValues () : Function<dim>(dim) {};

      virtual double value (const Point<dim>   &p,
                            const unsigned int  component = 0) const;

      virtual void vector_value (const Point<dim> &p,
                                 Vector<double>   &values) const;
    };

    template <int dim>
    double BoundaryValues<dim>::value (const Point<dim> &p,
                                       const unsigned int component) const
    {
      double return_value = 0;

      if (component == 0)
        return_value = 0.0;
      if (component == 1)
        return_value = 0.0;
      if (component == 2)
        return_value = 0.0;

      return return_value;
    }

    template <int dim>
    void BoundaryValues<dim>::vector_value (const Point<dim> &p,
                                            Vector<double>   &values) const
    {
      for (unsigned int c=0; c<this->n_components; ++c)
        values(c) = BoundaryValues<dim>::value (p, c);
    }


    template <int dim>
    class Obstacle : public Function<dim>
    {
    public:
      Obstacle (std::unique_ptr<Input<dim> > const &_input) :
        Function<dim>(dim),
        input_obstacle_copy(std::move (_input)) {};

      virtual double value (const Point<dim>   &p,
                            const unsigned int  component = 0) const;

      virtual void vector_value (const Point<dim> &p,
                                 Vector<double>   &values) const;

    private:
      std::unique_ptr<Input<dim> >  const &input_obstacle_copy;
    };

    template <int dim>
    double Obstacle<dim>::value (const Point<dim> &p,
                                 const unsigned int component) const
    {
      double R = 0.03;
      double return_value = 100.0;
      if (component == 0)
        return_value = p(0);
      if (component == 1)
        return_value = p(1);
      if (component == 2)
        {
          return_value = -std::sqrt (0.36 - (p(0)-0.5)*(p(0)-0.5) - (p(1)-0.5)*(p(1)-0.5)) + 1.59;
//          return_value = 1.999 - input_obstacle_copy->obstacle_function (p(0), p(1));
        }
      return return_value;
    }

    template <int dim>
    void Obstacle<dim>::vector_value (const Point<dim> &p,
                                      Vector<double> &values) const
    {
      for (unsigned int c=0; c<this->n_components; ++c)
        values(c) = Obstacle<dim>::value (p, c);
    }
  }


  // @sect3{Implementation of the <code>PlasticityContactProblem</code> class}

  // Next for the implementation of the class
  // template that makes use of the functions
  // above. As before, we will write everything

  template <int dim>
  PlasticityContactProblem<dim>::PlasticityContactProblem (int _n_refinements_global)
    :
    n_refinements_global (_n_refinements_global),
    mpi_communicator (MPI_COMM_WORLD),
    triangulation (mpi_communicator),
    fe (FE_Q<dim>(1), dim),
    dof_handler (triangulation),
    pcout (std::cout,
           (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
    sigma_0 (400),
    gamma (0.01),
    e_modul (2.0e+5),
    nu (0.3)
  {
    // double _E, double _nu, double _sigma_0, double _gamma
    plast_lin_hard.reset (new ConstitutiveLaw<dim> (e_modul, nu, sigma_0, gamma, mpi_communicator, pcout));
  }

  template <int dim>
  void PlasticityContactProblem<dim>::make_grid ()
  {
    std::vector<unsigned int> repet(3);
    repet[0] = 1;//20;
    repet[1] = 1;
    repet[2] = 1;

    Point<dim> p1 (0,0,0);
    Point<dim> p2 (1.0, 1.0, 1.0);
    GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2);

    Triangulation<3>::active_cell_iterator
    cell = triangulation.begin_active(),
    endc = triangulation.end();

    /* boundary_indicators:
              _______
             /  9    /|
            /______ / |
          8|       | 8|
           |   8   | /
           |_______|/
               6
     */

    for (; cell!=endc; ++cell)
      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
        {
          if (cell->face (face)->center ()[2] == p2(2))
            cell->face (face)->set_boundary_indicator (9);
          if (cell->face (face)->center ()[0] == p1(0) ||
              cell->face (face)->center ()[0] == p2(0) ||
              cell->face (face)->center ()[1] == p1(1) ||
              cell->face (face)->center ()[1] == p2(1))
            cell->face (face)->set_boundary_indicator (8);
          if (cell->face (face)->center ()[2] == p1(2))
            cell->face (face)->set_boundary_indicator (6);
        }

    triangulation.refine_global (n_refinements_global);
  }

  template <int dim>
  void PlasticityContactProblem<dim>::setup_system ()
  {
    // setup dofs
    {
      dof_handler.distribute_dofs (fe);

      locally_owned_dofs = dof_handler.locally_owned_dofs ();
      locally_relevant_dofs.clear();
      DoFTools::extract_locally_relevant_dofs (dof_handler,
                                               locally_relevant_dofs);
    }

    // setup hanging nodes and dirichlet constraints
    {
      constraints_hanging_nodes.clear ();
      constraints_hanging_nodes.reinit (locally_relevant_dofs);
      DoFTools::make_hanging_node_constraints (dof_handler,
                                               constraints_hanging_nodes);
      constraints_hanging_nodes.close ();

      pcout << "   Number of active cells: "
            << triangulation.n_active_cells()
            << std::endl
            << "   Total number of cells: "
            << triangulation.n_cells()
            << std::endl
            << "   Number of degrees of freedom: "
            << dof_handler.n_dofs ()
            << std::endl;

      dirichlet_constraints ();
    }

    // Initialzation for matrices and vectors
    {
      solution.reinit (locally_relevant_dofs, mpi_communicator);
      system_rhs_newton.reinit (locally_owned_dofs, mpi_communicator);
      old_solution.reinit (system_rhs_newton);
      resid_vector.reinit (system_rhs_newton);
      diag_mass_matrix_vector.reinit (system_rhs_newton);
      cell_constitution.reinit (triangulation.n_active_cells ());
      active_set.clear ();
      active_set.set_size (locally_relevant_dofs.size ());
    }

    // setup sparsity pattern
    {
      TrilinosWrappers::SparsityPattern sp (locally_owned_dofs,
                                            mpi_communicator);

      DoFTools::make_sparsity_pattern (dof_handler, sp, constraints_dirichlet_hanging_nodes, false,
                                       Utilities::MPI::this_mpi_process(mpi_communicator));

      sp.compress();

      system_matrix_newton.reinit (sp);

      TrilinosWrappers::SparseMatrix mass_matrix;
      mass_matrix.reinit (sp);
      assemble_mass_matrix_diagonal (mass_matrix);
      const unsigned int
      start = (system_rhs_newton.local_range().first),
      end   = (system_rhs_newton.local_range().second);
      for (unsigned int j=start; j<end; j++)
        diag_mass_matrix_vector (j) = mass_matrix.diag_element (j);
      number_iterations = 0;

      diag_mass_matrix_vector.compress (VectorOperation::insert);
    }
  }

  template <int dim>
  void PlasticityContactProblem<dim>::assemble_nl_system (TrilinosWrappers::MPI::Vector &u)
  {
    QGauss<dim>  quadrature_formula(2);
    QGauss<dim-1>  face_quadrature_formula(2);

    FEValues<dim> fe_values (fe, quadrature_formula,
                             UpdateFlags(update_values    |
                                         update_gradients |
                                         update_q_points  |
                                         update_JxW_values));

    FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
                                      update_values   | update_quadrature_points |
                                      update_JxW_values);

    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
    const unsigned int   n_q_points    = quadrature_formula.size ();
    const unsigned int   n_face_q_points = face_quadrature_formula.size();

    const EquationData::RightHandSide<dim> right_hand_side;
    std::vector<Vector<double> > right_hand_side_values (n_q_points,
                                                         Vector<double>(dim));
    std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
                                                              Vector<double>(dim));

    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
    Vector<double>       cell_rhs (dofs_per_cell);

    std::vector<unsigned int> local_dof_indices (dofs_per_cell);

    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                                   endc = dof_handler.end();

    const FEValuesExtractors::Vector displacement (0);

    const double kappa = 1.0;
    for (; cell!=endc; ++cell)
      if (cell->is_locally_owned())
        {
          fe_values.reinit (cell);
          cell_matrix = 0;
          cell_rhs = 0;

          right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
                                             right_hand_side_values);

          std::vector<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
          fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);

          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
            {
              SymmetricTensor<4,dim> stress_strain_tensor_linearized;
              SymmetricTensor<4,dim> stress_strain_tensor;
              SymmetricTensor<2,dim> stress_tensor;

              plast_lin_hard->linearized_plast_linear_hardening (stress_strain_tensor_linearized,
                                                                 stress_strain_tensor,
                                                                 strain_tensor[q_point]);

              for (unsigned int i=0; i<dofs_per_cell; ++i)
                {
                  stress_tensor = stress_strain_tensor_linearized * plast_lin_hard->get_strain(fe_values, i, q_point);

                  for (unsigned int j=0; j<dofs_per_cell; ++j)
                    {
                      cell_matrix(i,j) += (stress_tensor *
                                           plast_lin_hard->get_strain(fe_values, j, q_point) *
                                           fe_values.JxW (q_point));
                    }

                  // the linearized part a(v^i;v^i,v) of the rhs
                  cell_rhs(i) += (stress_tensor *
                                  strain_tensor[q_point] *
                                  fe_values.JxW (q_point));

                  // the residual part a(v^i;v) of the rhs
                  cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor *
                                  plast_lin_hard->get_strain(fe_values, i, q_point) *
                                  fe_values.JxW (q_point));

                  // the residual part F(v) of the rhs
                  Tensor<1,dim> rhs_values;
                  rhs_values = 0;
                  cell_rhs(i) += (fe_values[displacement].value (i, q_point) *
                                  rhs_values *
                                  fe_values.JxW (q_point));
                }
            }

          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (cell->face (face)->at_boundary()
                  && cell->face (face)->boundary_indicator () == 9)
                {
                  fe_values_face.reinit (cell, face);

                  right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(),
                                                     right_hand_side_values_face);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      Tensor<1,dim> rhs_values;
                      rhs_values = 0;
                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
                                        rhs_values *
                                        fe_values_face.JxW (q_point));
                    }
                }
            }

          cell->get_dof_indices (local_dof_indices);
          constraints.distribute_local_to_global (cell_matrix, cell_rhs,
                                                  local_dof_indices,
                                                  system_matrix_newton, system_rhs_newton, true);
        };

    system_matrix_newton.compress ();
    system_rhs_newton.compress (VectorOperation::add);
  }

  template <int dim>
  void PlasticityContactProblem<dim>::residual_nl_system (TrilinosWrappers::MPI::Vector &u)
  {
    QGauss<dim>  quadrature_formula(2);
    QGauss<dim-1> face_quadrature_formula(2);

    FEValues<dim> fe_values (fe, quadrature_formula,
                             UpdateFlags(update_values    |
                                         update_gradients |
                                         update_q_points  |
                                         update_JxW_values));

    FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
                                      update_values   | update_quadrature_points |
                                      update_JxW_values);

    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
    const unsigned int   n_q_points    = quadrature_formula.size ();
    const unsigned int   n_face_q_points = face_quadrature_formula.size();

    const EquationData::RightHandSide<dim> right_hand_side;
    std::vector<Vector<double> > right_hand_side_values (n_q_points,
                                                         Vector<double>(dim));
    std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
                                                              Vector<double>(dim));

    Vector<double>       cell_rhs (dofs_per_cell);

    std::vector<unsigned int> local_dof_indices (dofs_per_cell);

    const FEValuesExtractors::Vector displacement (0);

    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                                   endc = dof_handler.end();

    unsigned int elast_points = 0;
    unsigned int plast_points = 0;
    double       yield = 0;
    double       cell_number = 0;
    cell_constitution = 0;

    for (; cell!=endc; ++cell)
      if (cell->is_locally_owned())
        {
          fe_values.reinit (cell);
          cell_rhs = 0;

          right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
                                             right_hand_side_values);

          std::vector<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
          fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);

          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
            {
              SymmetricTensor<4,dim> stress_strain_tensor;
              SymmetricTensor<2,dim> stress_tensor;

              plast_lin_hard->plast_linear_hardening (stress_strain_tensor, strain_tensor[q_point],
                                                      elast_points, plast_points, yield);

              cell_constitution (cell_number) += yield;
              for (unsigned int i=0; i<dofs_per_cell; ++i)
                {
                  cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * //(stress_tensor) *
                                  plast_lin_hard->get_strain(fe_values, i, q_point) *
                                  fe_values.JxW (q_point));

                  Tensor<1,dim> rhs_values;
                  rhs_values = 0;
                  cell_rhs(i) += ((fe_values[displacement].value (i, q_point) *
                                   rhs_values) *
                                  fe_values.JxW (q_point));
                };
            };

          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (cell->face (face)->at_boundary()
                  && cell->face (face)->boundary_indicator () == 9)
                {
                  fe_values_face.reinit (cell, face);

                  right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(),
                                                     right_hand_side_values_face);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      Tensor<1,dim> rhs_values;
                      rhs_values = 0;
                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
                                        rhs_values *
                                        fe_values_face.JxW (q_point));
                    }
                }
            }

          cell->get_dof_indices (local_dof_indices);
          constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs,
              local_dof_indices,
              system_rhs_newton);

          cell_number += 1;
        }
      else
        {
          cell_constitution (cell_number) = 0;
          cell_number += 1;
        };

    cell_constitution /= n_q_points;
    cell_constitution.compress (VectorOperation::add);
    system_rhs_newton.compress (VectorOperation::add);

    unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, mpi_communicator);
    unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, mpi_communicator);
    pcout << "      Number of elastic quadrature points: " << sum_elast_points
          << " and plastic quadrature points: " << sum_plast_points <<std::endl;
  }

  template <int dim>
  void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal (TrilinosWrappers::SparseMatrix &mass_matrix)
  {
    QTrapez<dim-1>  face_quadrature_formula;

    FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
                                      update_values   |
                                      update_quadrature_points |
                                      update_JxW_values);

    const unsigned int   dofs_per_cell      = fe.dofs_per_cell;
    const unsigned int   n_face_q_points    = face_quadrature_formula.size();

    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);

    std::vector<unsigned int> local_dof_indices (dofs_per_cell);

    const FEValuesExtractors::Vector displacement (0);

    typename DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();

{'fraction_non_alphanumeric': 0.10453423856456191, 'fraction_numerical': 0.011068163340987862, 'mean_word_length': 1.9580626107501478, 'pattern_counts': {'":': 0, '<': 585, '<?xml version=': 0, '>': 539, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 226, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1062467', 'n_tokens_mistral': 22604, 'n_tokens_neox': 21170, 'n_words': 4365}