L-BFGS-phasefield-solver.h

Go to the documentation of this file.

 
 *   
 *   /* A monolithic scheme based on the L-BFGS method to solve the phase-field crack problem
 *    * 1. The phase-field formulation itself is based on "A phase field model for rate-independent
 *    *    crack propagation - Robust algorithmic implementation based on operator splits"
 *    *    by Christian Miehe , Martina Hofacker, Fabian Welschinger
 *    * 2. This code implements a monolithic approach. The phase-field irreversibility
 *    *    is enforced through the history field Phi_0^+ and the viscosity parameter.
 *    * 3. Using TBB for stiffness assembly and Gauss point calculation.
 *    * 4. Using adaptive mesh refinement.
 *    */
 *   
 *   #include <deal.II/grid/tria.h>
 *   #include <deal.II/grid/grid_generator.h>
 *   #include <deal.II/grid/grid_refinement.h>
 *   #include <deal.II/grid/grid_out.h>
 *   #include <deal.II/grid/grid_in.h>
 *   #include <deal.II/grid/manifold_lib.h>
 *   
 *   #include <deal.II/dofs/dof_handler.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_system.h>
 *   #include <deal.II/fe/fe_q.h>
 *   #include <deal.II/fe/fe_dgp_monomial.h>
 *   #include <deal.II/fe/mapping_q_eulerian.h>
 *   
 *   #include <deal.II/base/timer.h>
 *   #include <deal.II/base/quadrature_point_data.h>
 *   #include <deal.II/base/parameter_handler.h>
 *   
 *   #include <deal.II/lac/affine_constraints.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/dynamic_sparsity_pattern.h>
 *   #include <deal.II/lac/block_sparse_matrix.h>
 *   #include <deal.II/lac/block_vector.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/lac/solver_cg.h>
 *   #include <deal.II/lac/precondition.h>
 *   #include <deal.II/lac/linear_operator.h>
 *   #include <deal.II/lac/packaged_operation.h>
 *   #include <deal.II/lac/precondition_selector.h>
 *   #include <deal.II/lac/solver_selector.h>
 *   #include <deal.II/lac/sparse_direct.h>
 *   
 *   #include <deal.II/numerics/error_estimator.h>
 *   
 *   #include <deal.II/physics/elasticity/standard_tensors.h>
 *   
 *   #include <deal.II/base/quadrature_point_data.h>
 *   
 *   #include <deal.II/grid/grid_tools.h>
 *   
 *   #include <deal.II/base/work_stream.h>
 *   
 *   #include <deal.II/numerics/solution_transfer.h>
 *   
 *   #include <fstream>
 *   #include <iostream>
 *   #include <deal.II/base/logstream.h>
 *   
 *   #include "SpectrumDecomposition.h"
 *   #include "Utilities.h"
 *   
 *   namespace PhaseField
 *   {
 *     using namespace dealii;
 *   
 * @endcode
 * 
 * body force
 * 
 * @code
 *     template <int dim>
 *     void right_hand_side(const std::vector<Point<dim>> &points,
 *                          std::vector<Tensor<1, dim>> &  values,
 *                          const double fx,
 *                          const double fy,
 *                          const double fz)
 *     {
 *       Assert(values.size() == points.size(),
 *              ExcDimensionMismatch(values.size(), points.size()));
 *       Assert(dim >= 2, ExcNotImplemented());
 *   
 *       for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
 *         {
 *           if (dim == 2)
 *             {
 *               values[point_n][0] = fx;
 *               values[point_n][1] = fy;
 *             }
 *           else
 *             {
 *               values[point_n][0] = fx;
 *               values[point_n][1] = fy;
 *               values[point_n][2] = fz;
 *             }
 *         }
 *     }
 *   
 *     double degradation_function(const double d)
 *     {
 *       return (1.0 - d) * (1.0 - d);
 *     }
 *   
 *     double degradation_function_derivative(const double d)
 *     {
 *       return 2.0 * (d - 1.0);
 *     }
 *   
 *     double degradation_function_2nd_order_derivative(const double d)
 *     {
 *       (void) d;
 *       return 2.0;
 *     }
 *   
 *     namespace Parameters
 *     {
 *       struct Scenario
 *       {
 *         unsigned int m_scenario;
 *         std::string m_logfile_name;
 *         bool m_output_iteration_history;
 *         std::string m_type_nonlinear_solver;
 *         std::string m_type_line_search;
 *         std::string m_type_linear_solver;
 *         std::string m_refinement_strategy;
 *         unsigned int m_LBFGS_m;
 *         unsigned int m_global_refine_times;
 *         unsigned int m_local_prerefine_times;
 *         unsigned int m_max_adaptive_refine_times;
 *         int m_max_allowed_refinement_level;
 *         double m_phasefield_refine_threshold;
 *         double m_allowed_max_h_l_ratio;
 *         unsigned int m_total_material_regions;
 *         std::string m_material_file_name;
 *         int m_reaction_force_face_id;
 *   
 *         static void declare_parameters(ParameterHandler &prm);
 *         void parse_parameters(ParameterHandler &prm);
 *       };
 *   
 *       void Scenario::declare_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Scenario");
 *         {
 *           prm.declare_entry("Scenario number",
 *                             "1",
 *                             Patterns::Integer(0),
 *                             "Geometry, loading and boundary conditions scenario");
 *   
 *           prm.declare_entry("Log file name",
 *                             "Output.log",
 *                             Patterns::FileName(Patterns::FileName::input),
 *                             "Name of the file for log");
 *   
 *           prm.declare_entry("Output iteration history",
 *                             "yes",
 *                             Patterns::Selection("yes|no"),
 *                             "Shall we write iteration history to the log file?");
 *   
 *           prm.declare_entry("Nonlinear solver type",
 *                             "Newton",
 *                             Patterns::Selection("Newton|BFGS|LBFGS"),
 *                             "Type of solver used to solve the nonlinear system");
 *   
 *           prm.declare_entry("Line search type",
 *                             "GradientBased",
 *                             Patterns::Selection("GradientBased|StrongWolfe"),
 *                             "Type of line search method, the gradient-based method "
 *                             "should be preferred since it is generally faster");
 *   
 *           prm.declare_entry("Linear solver type",
 *                             "Direct",
 *                             Patterns::Selection("Direct|CG"),
 *                             "Type of solver used to solve the linear system B0");
 *   
 *           prm.declare_entry("Mesh refinement strategy",
 *                             "adaptive-refine",
 *                             Patterns::Selection("pre-refine|adaptive-refine"),
 *                             "Mesh refinement strategy: pre-refine or adaptive-refine");
 *   
 *           prm.declare_entry("LBFGS m",
 *                             "40",
 *                             Patterns::Integer(0),
 *                             "Number of vectors used for LBFGS");
 *   
 *           prm.declare_entry("Global refinement times",
 *                             "0",
 *                             Patterns::Integer(0),
 *                             "Global refinement times (across the entire domain)");
 *   
 *           prm.declare_entry("Local prerefinement times",
 *                             "0",
 *                             Patterns::Integer(0),
 *                             "Local pre-refinement times (assume crack path is known a priori), "
 *                             "only refine along the crack path.");
 *   
 *           prm.declare_entry("Max adaptive refinement times",
 *                             "100",
 *                             Patterns::Integer(0),
 *                             "Maximum number of adaptive refinement times allowed in each step");
 *   
 *           prm.declare_entry("Max allowed refinement level",
 *                             "100",
 *                             Patterns::Integer(0),
 *                             "Maximum allowed cell refinement level");
 *   
 *           prm.declare_entry("Phasefield refine threshold",
 *                             "0.8",
 *                             Patterns::Double(),
 *                             "Phasefield-based refinement threshold value");
 *   
 *           prm.declare_entry("Allowed max hl ratio",
 *                             "0.25",
 *                             Patterns::Double(),
 *                             "Allowed maximum ratio between mesh size h and length scale l");
 *   
 *           prm.declare_entry("Material regions",
 *                             "1",
 *                             Patterns::Integer(0),
 *                             "Number of material regions");
 *   
 *           prm.declare_entry("Material data file",
 *                             "1",
 *                             Patterns::FileName(Patterns::FileName::input),
 *                             "Material data file");
 *   
 *           prm.declare_entry("Reaction force face ID",
 *                             "1",
 *                             Patterns::Integer(),
 *                             "Face id where reaction forces should be calculated "
 *                             "(negative integer means not to calculate reaction force)");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       void Scenario::parse_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Scenario");
 *         {
 *           m_scenario = prm.get_integer("Scenario number");
 *           m_logfile_name = prm.get("Log file name");
 *           m_output_iteration_history = prm.get_bool("Output iteration history");
 *           m_type_nonlinear_solver = prm.get("Nonlinear solver type");
 *           m_type_line_search = prm.get("Line search type");
 *           m_type_linear_solver = prm.get("Linear solver type");
 *           m_refinement_strategy = prm.get("Mesh refinement strategy");
 *           m_LBFGS_m = prm.get_integer("LBFGS m");
 *           m_global_refine_times = prm.get_integer("Global refinement times");
 *           m_local_prerefine_times = prm.get_integer("Local prerefinement times");
 *           m_max_adaptive_refine_times = prm.get_integer("Max adaptive refinement times");
 *           m_max_allowed_refinement_level = prm.get_integer("Max allowed refinement level");
 *           m_phasefield_refine_threshold = prm.get_double("Phasefield refine threshold");
 *           m_allowed_max_h_l_ratio = prm.get_double("Allowed max hl ratio");
 *           m_total_material_regions = prm.get_integer("Material regions");
 *           m_material_file_name = prm.get("Material data file");
 *           m_reaction_force_face_id = prm.get_integer("Reaction force face ID");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       struct FESystem
 *       {
 *         unsigned int m_poly_degree;
 *         unsigned int m_quad_order;
 *   
 *         static void declare_parameters(ParameterHandler &prm);
 *   
 *         void parse_parameters(ParameterHandler &prm);
 *       };
 *   
 *   
 *       void FESystem::declare_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Finite element system");
 *         {
 *           prm.declare_entry("Polynomial degree",
 *                             "1",
 *                             Patterns::Integer(0),
 *                             "Phase field polynomial order");
 *   
 *           prm.declare_entry("Quadrature order",
 *                             "2",
 *                             Patterns::Integer(0),
 *                             "Gauss quadrature order");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       void FESystem::parse_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Finite element system");
 *         {
 *           m_poly_degree = prm.get_integer("Polynomial degree");
 *           m_quad_order  = prm.get_integer("Quadrature order");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 * @endcode
 * 
 * body force (N/m^3)
 * 
 * @code
 *       struct BodyForce
 *       {
 *         double m_x_component;
 *         double m_y_component;
 *         double m_z_component;
 *   
 *         static void declare_parameters(ParameterHandler &prm);
 *   
 *         void parse_parameters(ParameterHandler &prm);
 *       };
 *   
 *       void BodyForce::declare_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Body force");
 *         {
 *           prm.declare_entry("Body force x component",
 *                             "0.0",
 *                             Patterns::Double(),
 *                             "Body force x-component (N/m^3)");
 *   
 *           prm.declare_entry("Body force y component",
 *                             "0.0",
 *                             Patterns::Double(),
 *                             "Body force y-component (N/m^3)");
 *   
 *           prm.declare_entry("Body force z component",
 *                             "0.0",
 *                             Patterns::Double(),
 *                             "Body force z-component (N/m^3)");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       void BodyForce::parse_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Body force");
 *         {
 *           m_x_component = prm.get_double("Body force x component");
 *           m_y_component = prm.get_double("Body force y component");
 *           m_z_component = prm.get_double("Body force z component");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       struct NonlinearSolver
 *       {
 *         unsigned int m_max_iterations_NR;
 *         unsigned int m_max_iterations_BFGS;
 *         bool m_relative_residual;
 *   
 *         double       m_tol_u_residual;
 *         double       m_tol_d_residual;
 *         double       m_tol_u_incr;
 *         double       m_tol_d_incr;
 *   
 *         static void declare_parameters(ParameterHandler &prm);
 *   
 *         void parse_parameters(ParameterHandler &prm);
 *       };
 *   
 *       void NonlinearSolver::declare_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Nonlinear solver");
 *         {
 *           prm.declare_entry("Max iterations Newton-Raphson",
 *                             "10",
 *                             Patterns::Integer(0),
 *                             "Number of Newton-Raphson iterations allowed");
 *   
 *           prm.declare_entry("Max iterations BFGS",
 *                             "20",
 *                             Patterns::Integer(0),
 *                             "Number of BFGS iterations allowed");
 *   
 *           prm.declare_entry("Relative residual",
 *                             "yes",
 *                             Patterns::Selection("yes|no"),
 *                             "Shall we use relative residual for convergence?");
 *   
 *           prm.declare_entry("Tolerance displacement residual",
 *                             "1.0e-9",
 *                             Patterns::Double(0.0),
 *                             "Displacement residual tolerance");
 *   
 *           prm.declare_entry("Tolerance phasefield residual",
 *                             "1.0e-9",
 *                             Patterns::Double(0.0),
 *                             "Phasefield residual tolerance");
 *   
 *           prm.declare_entry("Tolerance displacement increment",
 *                             "1.0e-9",
 *                             Patterns::Double(0.0),
 *                             "Displacement increment tolerance");
 *   
 *           prm.declare_entry("Tolerance phasefield increment",
 *                             "1.0e-9",
 *                             Patterns::Double(0.0),
 *                             "Phasefield increment tolerance");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       void NonlinearSolver::parse_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Nonlinear solver");
 *         {
 *           m_max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
 *           m_max_iterations_BFGS = prm.get_integer("Max iterations BFGS");
 *           m_relative_residual = prm.get_bool("Relative residual");
 *   
 *           m_tol_u_residual           = prm.get_double("Tolerance displacement residual");
 *           m_tol_d_residual           = prm.get_double("Tolerance phasefield residual");
 *           m_tol_u_incr               = prm.get_double("Tolerance displacement increment");
 *           m_tol_d_incr               = prm.get_double("Tolerance phasefield increment");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       struct TimeInfo
 *       {
 *         double m_end_time;
 *         std::string m_time_file_name;
 *   
 *         static void declare_parameters(ParameterHandler &prm);
 *   
 *         void parse_parameters(ParameterHandler &prm);
 *       };
 *   
 *       void TimeInfo::declare_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Time");
 *         {
 *           prm.declare_entry("End time", "1", Patterns::Double(), "End time");
 *   
 *           prm.declare_entry("Time data file",
 *                             "1",
 *                             Patterns::FileName(Patterns::FileName::input),
 *                             "Time data file");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       void TimeInfo::parse_parameters(ParameterHandler &prm)
 *       {
 *         prm.enter_subsection("Time");
 *         {
 *           m_end_time = prm.get_double("End time");
 *           m_time_file_name = prm.get("Time data file");
 *         }
 *         prm.leave_subsection();
 *       }
 *   
 *       struct AllParameters : public Scenario,
 *                              public FESystem,
 *                              public BodyForce,
 *                              public NonlinearSolver,
 *                              public TimeInfo
 *       {
 *         AllParameters(const std::string &input_file);
 *   
 *         static void declare_parameters(ParameterHandler &prm);
 *   
 *         void parse_parameters(ParameterHandler &prm);
 *       };
 *   
 *       AllParameters::AllParameters(const std::string &input_file)
 *       {
 *         ParameterHandler prm;
 *         declare_parameters(prm);
 *         prm.parse_input(input_file);
 *         parse_parameters(prm);
 *       }
 *   
 *       void AllParameters::declare_parameters(ParameterHandler &prm)
 *       {
 *         Scenario::declare_parameters(prm);
 *         FESystem::declare_parameters(prm);
 *         BodyForce::declare_parameters(prm);
 *         NonlinearSolver::declare_parameters(prm);
 *         TimeInfo::declare_parameters(prm);
 *       }
 *   
 *       void AllParameters::parse_parameters(ParameterHandler &prm)
 *       {
 *         Scenario::parse_parameters(prm);
 *         FESystem::parse_parameters(prm);
 *         BodyForce::parse_parameters(prm);
 *         NonlinearSolver::parse_parameters(prm);
 *         TimeInfo::parse_parameters(prm);
 *       }
 *     } // namespace Parameters
 *   
 *     class Time
 *     {
 *     public:
 *       Time(const double time_end)
 *         : m_timestep(0)
 *         , m_time_current(0.0)
 *         , m_time_end(time_end)
 *         , m_delta_t(0.0)
 *         , m_magnitude(1.0)
 *       {}
 *   
 *       virtual ~Time() = default;
 *   
 *       double current() const
 *       {
 *         return m_time_current;
 *       }
 *       double end() const
 *       {
 *         return m_time_end;
 *       }
 *       double get_delta_t() const
 *       {
 *         return m_delta_t;
 *       }
 *       double get_magnitude() const
 *       {
 *         return m_magnitude;
 *       }
 *       unsigned int get_timestep() const
 *       {
 *         return m_timestep;
 *       }
 *       void increment(std::vector<std::array<double, 4>> time_table)
 *       {
 *         double t_1, t_delta, t_magnitude;
 *         for (auto & time_group : time_table)
 *           {
 *             t_1 = time_group[1];
 *             t_delta = time_group[2];
 *             t_magnitude = time_group[3];
 *   
 *             if (m_time_current < t_1 - 1.0e-6*t_delta)
 *               {
 *                 m_delta_t = t_delta;
 *                 m_magnitude = t_magnitude;
 *                 break;
 *               }
 *           }
 *   
 *         m_time_current += m_delta_t;
 *         ++m_timestep;
 *       }
 *   
 *     private:
 *       unsigned int m_timestep;
 *       double       m_time_current;
 *       const double m_time_end;
 *       double m_delta_t;
 *       double m_magnitude;
 *     };
 *   
 *     template <int dim>
 *     class LinearIsotropicElasticityAdditiveSplit
 *     {
 *     public:
 *       LinearIsotropicElasticityAdditiveSplit(const double lame_lambda,
 *                                              const double lame_mu,
 *                                              const double residual_k,
 *                                              const double length_scale,
 *                                              const double viscosity,
 *                                              const double gc)
 *         : m_lame_lambda(lame_lambda)
 *         , m_lame_mu(lame_mu)
 *         , m_residual_k(residual_k)
 *         , m_length_scale(length_scale)
 *         , m_eta(viscosity)
 *         , m_gc(gc)
 *         , m_phase_field_value(0.0)
 *         , m_grad_phasefield(Tensor<1, dim>())
 *         , m_strain(SymmetricTensor<2, dim>())
 *         , m_stress(SymmetricTensor<2, dim>())
 *         , m_stress_positive(SymmetricTensor<2, dim>())
 *         , m_mechanical_C(SymmetricTensor<4, dim>())
 *         , m_strain_energy_positive(0.0)
 *         , m_strain_energy_negative(0.0)
 *         , m_strain_energy_total(0.0)
 *         , m_crack_energy_dissipation(0.0)
 *       {
 *         Assert(  ( lame_lambda / (2*(lame_lambda + lame_mu)) <= 0.5)
 *                & ( lame_lambda / (2*(lame_lambda + lame_mu)) >=-1.0),
 *                ExcInternalError() );
 *       }
 *   
 *       const SymmetricTensor<4, dim> & get_mechanical_C() const
 *       {
 *         return m_mechanical_C;
 *       }
 *   
 *       const SymmetricTensor<2, dim> & get_cauchy_stress() const
 *       {
 *         return m_stress;
 *       }
 *   
 *       const SymmetricTensor<2, dim> & get_cauchy_stress_positive() const
 *       {
 *         return m_stress_positive;
 *       }
 *   
 *       double get_positive_strain_energy() const
 *       {
 *         return m_strain_energy_positive;
 *       }
 *   
 *       double get_negative_strain_energy() const
 *       {
 *         return m_strain_energy_negative;
 *       }
 *   
 *       double get_total_strain_energy() const
 *       {
 *         return m_strain_energy_total;
 *       }
 *   
 *       double get_crack_energy_dissipation() const
 *       {
 *         return m_crack_energy_dissipation;
 *       }
 *   
 *       double get_phase_field_value() const
 *       {
 *         return m_phase_field_value;
 *       }
 *   
 *       const Tensor<1, dim> get_phase_field_gradient() const
 *       {
 *         return m_grad_phasefield;
 *       }
 *   
 *       void update_material_data(const SymmetricTensor<2, dim> & strain,
 *                                 const double phase_field_value,
 *                                 const Tensor<1, dim> & grad_phasefield,
 *                                 const double phase_field_value_previous_step,
 *                                 const double delta_time)
 *       {
 *         m_strain = strain;
 *         m_phase_field_value = phase_field_value;
 *         m_grad_phasefield = grad_phasefield;
 *         Vector<double>              eigenvalues(dim);
 *         std::vector<Tensor<1, dim>> eigenvectors(dim);
 *         usr_spectrum_decomposition::spectrum_decomposition<dim>(m_strain,
 *                                                                     eigenvalues,
 *                                                                     eigenvectors);
 *   
 *         SymmetricTensor<2, dim> strain_positive, strain_negative;
 *         strain_positive = usr_spectrum_decomposition::positive_tensor(eigenvalues, eigenvectors);
 *         strain_negative = usr_spectrum_decomposition::negative_tensor(eigenvalues, eigenvectors);
 *   
 *         SymmetricTensor<4, dim> projector_positive, projector_negative;
 *         usr_spectrum_decomposition::positive_negative_projectors(eigenvalues,
 *                                                                      eigenvectors,
 *                                                                  projector_positive,
 *                                                                  projector_negative);
 *   
 *         SymmetricTensor<2, dim> stress_positive, stress_negative;
 *         const double degradation = degradation_function(m_phase_field_value) + m_residual_k;
 *         const double I_1 = trace(m_strain);
 *         stress_positive = m_lame_lambda * usr_spectrum_decomposition::positive_ramp_function(I_1)
 *                                         * Physics::Elasticity::StandardTensors<dim>::I
 *                         + 2 * m_lame_mu * strain_positive;
 *         stress_negative = m_lame_lambda * usr_spectrum_decomposition::negative_ramp_function(I_1)
 *                                         * Physics::Elasticity::StandardTensors<dim>::I
 *                               + 2 * m_lame_mu * strain_negative;
 *   
 *         m_stress = degradation * stress_positive + stress_negative;
 *         m_stress_positive = stress_positive;
 *   
 *         SymmetricTensor<4, dim> C_positive, C_negative;
 *         C_positive = m_lame_lambda * usr_spectrum_decomposition::heaviside_function(I_1)
 *                                    * Physics::Elasticity::StandardTensors<dim>::IxI
 *                    + 2 * m_lame_mu * projector_positive;
 *         C_negative = m_lame_lambda * usr_spectrum_decomposition::heaviside_function(-I_1)
 *                                    * Physics::Elasticity::StandardTensors<dim>::IxI
 *                          + 2 * m_lame_mu * projector_negative;
 *         m_mechanical_C = degradation * C_positive + C_negative;
 *   
 *         m_strain_energy_positive = 0.5 * m_lame_lambda * usr_spectrum_decomposition::positive_ramp_function(I_1)
 *                                                        * usr_spectrum_decomposition::positive_ramp_function(I_1)
 *                                  + m_lame_mu * strain_positive * strain_positive;
 *   
 *         m_strain_energy_negative = 0.5 * m_lame_lambda * usr_spectrum_decomposition::negative_ramp_function(I_1)
 *                                                        * usr_spectrum_decomposition::negative_ramp_function(I_1)
 *                                  + m_lame_mu * strain_negative * strain_negative;
 *   
 *         m_strain_energy_total = degradation * m_strain_energy_positive + m_strain_energy_negative;
 *   
 *         m_crack_energy_dissipation = m_gc * (  0.5 / m_length_scale * m_phase_field_value * m_phase_field_value
 *                                              + 0.5 * m_length_scale * m_grad_phasefield * m_grad_phasefield)
 * @endcode
 * 
 * the term due to viscosity regularization
 * 
 * @code
 *                                              + (m_phase_field_value - phase_field_value_previous_step)
 *                                              * (m_phase_field_value - phase_field_value_previous_step)
 *                                              * 0.5 * m_eta / delta_time;
 * @endcode
 * 
 * (void)delta_time;
 * (void)phase_field_value_previous_step;
 * 
 * @code
 *       }
 *   
 *     private:
 *       const double m_lame_lambda;
 *       const double m_lame_mu;
 *       const double m_residual_k;
 *       const double m_length_scale;
 *       const double m_eta;
 *       const double m_gc;
 *       double m_phase_field_value;
 *       Tensor<1, dim> m_grad_phasefield;
 *       SymmetricTensor<2, dim> m_strain;
 *       SymmetricTensor<2, dim> m_stress;
 *       SymmetricTensor<2, dim> m_stress_positive;
 *       SymmetricTensor<4, dim> m_mechanical_C;
 *       double m_strain_energy_positive;
 *       double m_strain_energy_negative;
 *       double m_strain_energy_total;
 *       double m_crack_energy_dissipation;
 *     };
 *   
 *   
 *     template <int dim>
 *     class PointHistory
 *     {
 *     public:
 *       PointHistory()
 *         : m_length_scale(0.0)
 *         , m_gc(0.0)
 *         , m_viscosity(0.0)
 *         , m_history_max_positive_strain_energy(0.0)
 *       {}
 *   
 *       virtual ~PointHistory() = default;
 *   
 *       void setup_lqp(const double lame_lambda,
 *                      const double lame_mu,
 *                      const double length_scale,
 *                      const double gc,
 *                      const double viscosity,
 *                      const double residual_k)
 *       {
 *         m_material =
 *                 std::make_shared<LinearIsotropicElasticityAdditiveSplit<dim>>(lame_lambda,
 *                                                                               lame_mu,
 *                                                                               residual_k,
 *                                                                               length_scale,
 *                                                                               viscosity,
 *                                                                               gc);
 *         m_history_max_positive_strain_energy = 0.0;
 *         m_length_scale = length_scale;
 *         m_gc = gc;
 *         m_viscosity = viscosity;
 *   
 *         update_field_values(SymmetricTensor<2, dim>(), 0.0, Tensor<1, dim>(), 0.0, 1.0);
 *       }
 *   
 *       void update_field_values(const SymmetricTensor<2, dim> & strain,
 *                                const double phase_field_value,
 *                                const Tensor<1, dim> & grad_phasefield,
 *                                const double phase_field_value_previous_step,
 *                                const double delta_time)
 *       {
 *         m_material->update_material_data(strain, phase_field_value, grad_phasefield,
 *                                          phase_field_value_previous_step, delta_time);
 *       }
 *   
 *       void update_history_variable()
 *       {
 *         double current_positive_strain_energy = m_material->get_positive_strain_energy();
 *         m_history_max_positive_strain_energy = std::fmax(m_history_max_positive_strain_energy,
 *                                                          current_positive_strain_energy);
 *       }
 *   
 *       double get_current_positive_strain_energy() const
 *       {
 *         return m_material->get_positive_strain_energy();
 *       }
 *   
 *       const SymmetricTensor<4, dim> & get_mechanical_C() const
 *       {
 *         return m_material->get_mechanical_C();
 *       }
 *   
 *       const SymmetricTensor<2, dim> & get_cauchy_stress() const
 *       {
 *         return m_material->get_cauchy_stress();
 *       }
 *   
 *       const SymmetricTensor<2, dim> & get_cauchy_stress_positive() const
 *       {
 *         return m_material->get_cauchy_stress_positive();
 *       }
 *   
 *       double get_total_strain_energy() const
 *       {
 *         return m_material->get_total_strain_energy();
 *       }
 *   
 *       double get_crack_energy_dissipation() const
 *       {
 *         return m_material->get_crack_energy_dissipation();
 *       }
 *   
 *       double get_phase_field_value() const
 *       {
 *         return m_material->get_phase_field_value();
 *       }
 *   
 *       const Tensor<1, dim> get_phase_field_gradient() const
 *       {
 *         return m_material->get_phase_field_gradient();
 *       }
 *   
 *       double get_history_max_positive_strain_energy() const
 *       {
 *         return m_history_max_positive_strain_energy;
 *       }
 *   
 *       double get_length_scale() const
 *       {
 *         return m_length_scale;
 *       }
 *   
 *       double get_critical_energy_release_rate() const
 *       {
 *         return m_gc;
 *       }
 *   
 *       double get_viscosity() const
 *       {
 *         return m_viscosity;
 *       }
 *     private:
 *       std::shared_ptr<LinearIsotropicElasticityAdditiveSplit<dim>> m_material;
 *       double m_length_scale;
 *       double m_gc;
 *       double m_viscosity;
 *       double m_history_max_positive_strain_energy;
 *     };
 *   
 *     template <int dim>
 *     class PhaseFieldMonolithicSolve
 *     {
 *     public:
 *       PhaseFieldMonolithicSolve(const std::string &input_file);
 *   
 *       virtual ~PhaseFieldMonolithicSolve() = default;
 *       void run();
 *   
 *     private:
 *       struct PerTaskData_ASM;
 *       struct ScratchData_ASM;
 *   
 *       struct PerTaskData_ASM_RHS_BFGS;
 *       struct ScratchData_ASM_RHS_BFGS;
 *   
 *       struct PerTaskData_UQPH;
 *       struct ScratchData_UQPH;
 *   
 *       Parameters::AllParameters m_parameters;
 *       Triangulation<dim> m_triangulation;
 *   
 *       CellDataStorage<typename Triangulation<dim>::cell_iterator,
 *                       PointHistory<dim>>
 *         m_quadrature_point_history;
 *   
 *       Time                m_time;
 *       std::ofstream m_logfile;
 *       mutable TimerOutput m_timer;
 *   
 *       DoFHandler<dim>                  m_dof_handler;
 *       FESystem<dim>                    m_fe;
 *       const unsigned int               m_dofs_per_cell;
 *       const FEValuesExtractors::Vector m_u_fe;
 *       const FEValuesExtractors::Scalar m_d_fe;
 *   
 *       static const unsigned int m_n_blocks          = 2;
 *       static const unsigned int m_n_components      = dim + 1;
 *       static const unsigned int m_first_u_component = 0;
 *       static const unsigned int m_d_component       = dim;
 *   
 *       enum
 *       {
 *         m_u_dof = 0,
 *         m_d_dof = 1
 *       };
 *   
 *       std::vector<types::global_dof_index> m_dofs_per_block;
 *   
 *       const QGauss<dim>     m_qf_cell;
 *       const QGauss<dim - 1> m_qf_face;
 *       const unsigned int    m_n_q_points;
 *   
 *       double m_vol_reference;
 *   
 *       AffineConstraints<double> m_constraints;
 *       BlockSparsityPattern      m_sparsity_pattern;
 *       BlockSparseMatrix<double> m_tangent_matrix;
 *       BlockVector<double>       m_system_rhs;
 *       BlockVector<double>       m_solution;
 *       SparseDirectUMFPACK       m_A_direct;
 *   
 *   
 *       std::map<unsigned int, std::vector<double>> m_material_data;
 *   
 *       std::vector<std::pair<double, std::vector<double>>> m_history_reaction_force;
 *       std::vector<std::pair<double, std::array<double, 3>>> m_history_energy;
 *   
 *   
 *       struct Errors
 *       {
 *         Errors()
 *           : m_norm(1.0)
 *           , m_u(1.0)
 *           , m_d(1.0)
 *         {}
 *   
 *         void reset()
 *         {
 *           m_norm = 1.0;
 *           m_u    = 1.0;
 *           m_d    = 1.0;
 *         }
 *   
 *         void normalize(const Errors &rhs)
 *         {
 *           if (rhs.m_norm != 0.0)
 *             m_norm /= rhs.m_norm;
 *           if (rhs.m_u != 0.0)
 *             m_u /= rhs.m_u;
 *           if (rhs.m_d != 0.0)
 *             m_d /= rhs.m_d;
 *         }
 *   
 *         double m_norm, m_u, m_d;
 *       };
 *   
 *       Errors m_error_residual, m_error_residual_0, m_error_residual_norm, m_error_update,
 *         m_error_update_0, m_error_update_norm;
 *   
 *       void get_error_residual(Errors &error_residual);
 *       void get_error_update(const BlockVector<double> &newton_update,
 *                             Errors & error_update);
 *   
 *       void make_grid();
 *       void make_grid_case_1();
 *       void make_grid_case_2();
 *       void make_grid_case_3();
 *       void make_grid_case_4();
 *       void make_grid_case_5();
 *       void make_grid_case_6();
 *       void make_grid_case_7();
 *       void make_grid_case_8();
 *       void make_grid_case_9();
 *       void make_grid_case_11();
 *   
 *       void setup_system();
 *   
 *       void determine_component_extractors();
 *   
 *       void make_constraints(const unsigned int it_nr);
 *   
 *       void assemble_system_newton(const BlockVector<double> & solution_old);
 *   
 *       void assemble_system_B0(const BlockVector<double> & solution_old);
 *   
 *       void assemble_system_newton_one_cell(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         ScratchData_ASM &                                     scratch,
 *         PerTaskData_ASM &                                     data) const;
 *   
 *       void assemble_system_B0_one_cell(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         ScratchData_ASM &                                     scratch,
 *         PerTaskData_ASM &                                     data) const;
 *   
 *       void assemble_system_rhs_BFGS_one_cell(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         ScratchData_ASM_RHS_BFGS &                           scratch,
 *         PerTaskData_ASM_RHS_BFGS &                           data) const;
 *   
 *       void assemble_system_rhs_BFGS(const BlockVector<double> & solution_old,
 *                                     BlockVector<double> & system_rhs);
 *   
 *       void assemble_system_rhs_BFGS_parallel(const BlockVector<double> & solution_old,
 *                                                  BlockVector<double> & system_rhs);
 *   
 *       bool solve_nonlinear_timestep_newton(BlockVector<double> &solution_delta);
 *   
 *       void solve_nonlinear_timestep_BFGS(BlockVector<double> &solution_delta);
 *   
 *       void solve_nonlinear_timestep_LBFGS(BlockVector<double> &solution_delta,
 *                                           BlockVector<double> & LBFGS_update_refine);
 *   
 *       double line_search_stepsize_strong_wolfe(const double phi_0,
 *                                                const double phi_0_prime,
 *                                                const BlockVector<double> & BFGS_p_vector,
 *                                                const BlockVector<double> & solution_delta);
 *   
 *       double line_search_stepsize_gradient_based(const BlockVector<double> & BFGS_p_vector,
 *                                                  const BlockVector<double> & solution_delta);
 *   
 *       double line_search_zoom_strong_wolfe(double phi_low, double phi_low_prime, double alpha_low,
 *                                            double phi_high, double phi_high_prime, double alpha_high,
 *                                            double phi_0, double phi_0_prime, const BlockVector<double> & BFGS_p_vector,
 *                                            double c1, double c2, unsigned int max_iter,
 *                                            const BlockVector<double> & solution_delta);
 *   
 *       double line_search_interpolation_cubic(const double alpha_0, const double phi_0, const double phi_0_prime,
 *                                              const double alpha_1, const double phi_1, const double phi_1_prime);
 *   
 *       std::pair<double, double> calculate_phi_and_phi_prime(const double alpha,
 *                                                             const BlockVector<double> & BFGS_p_vector,
 *                                                             const BlockVector<double> & solution_delta);
 *   
 *       std::vector<double> solve_linear_system(BlockVector<double> &newton_update);
 *   
 *       void LBFGS_B0(BlockVector<double> & LBFGS_r_vector,
 *                     BlockVector<double> & LBFGS_q_vector);
 *   
 *       void update_history_field_step();
 *   
 *       void output_results() const;
 *   
 *       void setup_qph();
 *   
 *       void update_qph_incremental(const BlockVector<double> &solution_delta,
 *                                   const BlockVector<double> &solution_old,
 *                                   const bool is_print);
 *   
 *       void update_qph_incremental_one_cell(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         ScratchData_UQPH &                                    scratch,
 *         PerTaskData_UQPH &                                    data);
 *   
 *       void copy_local_to_global_UQPH(const PerTaskData_UQPH & /*data*/)
 *       {}
 *   
 *       BlockVector<double>
 *       get_total_solution(const BlockVector<double> &solution_delta) const;
 *   
 * @endcode
 * 
 * Should not make this function const
 * 
 * @code
 *       void read_material_data(const std::string &data_file,
 *                               const unsigned int total_material_regions);
 *   
 *       void read_time_data(const std::string &data_file,
 *                               std::vector<std::array<double, 4>> & time_table);
 *   
 *       void print_conv_header_newton();
 *   
 *       void print_conv_header_BFGS();
 *   
 *       void print_conv_header_LBFGS();
 *   
 *       void print_parameter_information();
 *   
 *       void calculate_reaction_force(unsigned int face_ID);
 *   
 *       void write_history_data();
 *   
 *       double calculate_energy_functional() const;
 *   
 *       std::pair<double, double> calculate_total_strain_energy_and_crack_energy_dissipation() const;
 *   
 *       bool local_refine_and_solution_transfer(BlockVector<double> & solution_delta,
 *                                               BlockVector<double> & LBFGS_update_refine);
 *     }; // class PhaseFieldSplitSolve
 *   
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::get_error_residual(Errors &error_residual)
 *     {
 *       BlockVector<double> error_res(m_dofs_per_block);
 *   
 *       for (unsigned int i = 0; i < m_dof_handler.n_dofs(); ++i)
 *         if (!m_constraints.is_constrained(i))
 *           error_res(i) = m_system_rhs(i);
 *   
 *       error_residual.m_norm = error_res.l2_norm();
 *       error_residual.m_u    = error_res.block(m_u_dof).l2_norm();
 *       error_residual.m_d    = error_res.block(m_d_dof).l2_norm();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::get_error_update(const BlockVector<double> &newton_update,
 *                                                           Errors & error_update)
 *     {
 *       BlockVector<double> error_ud(m_dofs_per_block);
 *       for (unsigned int i = 0; i < m_dof_handler.n_dofs(); ++i)
 *         if (!m_constraints.is_constrained(i))
 *           error_ud(i) = newton_update(i);
 *   
 *       error_update.m_norm = error_ud.l2_norm();
 *       error_update.m_u    = error_ud.block(m_u_dof).l2_norm();
 *       error_update.m_d    = error_ud.block(m_d_dof).l2_norm();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::read_material_data(const std::string &data_file,
 *                                                        const unsigned int total_material_regions)
 *     {
 *       std::ifstream myfile (data_file);
 *   
 *       double lame_lambda, lame_mu, length_scale, gc, viscosity, residual_k;
 *       int material_region;
 *       double poisson_ratio;
 *       if (myfile.is_open())
 *         {
 *           m_logfile << "Reading material data file ..." << std::endl;
 *   
 *           while ( myfile >> material_region
 *                          >> lame_lambda
 *                          >> lame_mu
 *                          >> length_scale
 *                          >> gc
 *                          >> viscosity
 *                          >> residual_k)
 *             {
 *               m_material_data[material_region] = {lame_lambda,
 *                                                   lame_mu,
 *                                                   length_scale,
 *                                                   gc,
 *                                                   viscosity,
 *                                                   residual_k};
 *               poisson_ratio = lame_lambda / (2*(lame_lambda + lame_mu));
 *               Assert( (poisson_ratio <= 0.5)&(poisson_ratio >=-1.0) , ExcInternalError());
 *   
 *               m_logfile << "\tRegion " << material_region << " : " << std::endl;
 *               m_logfile << "\t\tLame lambda = " << lame_lambda << std::endl;
 *               m_logfile << "\t\tLame mu = "  << lame_mu << std::endl;
 *               m_logfile << "\t\tPoisson ratio = "  << poisson_ratio << std::endl;
 *               m_logfile << "\t\tPhase field length scale (l) = " << length_scale << std::endl;
 *               m_logfile << "\t\tCritical energy release rate (gc) = "  << gc << std::endl;
 *               m_logfile << "\t\tViscosity for regularization (eta) = "  << viscosity << std::endl;
 *               m_logfile << "\t\tResidual_k (k) = "  << residual_k << std::endl;
 *             }
 *   
 *           if (m_material_data.size() != total_material_regions)
 *             {
 *               m_logfile << "Material data file has " << m_material_data.size() << " rows. However, "
 *                         << "the mesh has " << total_material_regions << " material regions."
 *                         << std::endl;
 *               Assert(m_material_data.size() == total_material_regions,
 *                          ExcDimensionMismatch(m_material_data.size(), total_material_regions));
 *             }
 *           myfile.close();
 *         }
 *       else
 *         {
 *           m_logfile << "Material data file : " << data_file << " not exist!" << std::endl;
 *           Assert(false, ExcMessage("Failed to read material data file"));
 *         }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::read_time_data(const std::string &data_file,
 *                                                         std::vector<std::array<double, 4>> & time_table)
 *     {
 *       std::ifstream myfile (data_file);
 *   
 *       double t_0, t_1, delta_t, t_magnitude;
 *   
 *       if (myfile.is_open())
 *         {
 *           m_logfile << "Reading time data file ..." << std::endl;
 *   
 *           while ( myfile >> t_0
 *                          >> t_1
 *                          >> delta_t
 *                          >> t_magnitude)
 *             {
 *               Assert( t_0 < t_1,
 *                       ExcMessage("For each time pair, "
 *                                  "the start time should be smaller than the end time"));
 *               time_table.push_back({{t_0, t_1, delta_t, t_magnitude}});
 *             }
 *   
 *           Assert(std::fabs(t_1 - m_parameters.m_end_time) < 1.0e-9,
 *                  ExcMessage("End time in time table is inconsistent with input data in parameters.prm"));
 *   
 *           Assert(time_table.size() > 0,
 *                  ExcMessage("Time data file is empty."));
 *           myfile.close();
 *         }
 *       else
 *         {
 *           m_logfile << "Time data file : " << data_file << " not exist!" << std::endl;
 *           Assert(false, ExcMessage("Failed to read time data file"));
 *         }
 *   
 *       for (auto & time_group : time_table)
 *         {
 *           m_logfile << "\t\t"
 *                     << time_group[0] << ",\t"
 *                     << time_group[1] << ",\t"
 *                     << time_group[2] << ",\t"
 *                     << time_group[3] << std::endl;
 *         }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::setup_qph()
 *     {
 *       m_logfile << "\t\tSetting up quadrature point data ("
 *                 << m_n_q_points
 *                 << " points per cell)" << std::endl;
 *   
 *       m_quadrature_point_history.clear();
 *       for (auto const & cell : m_triangulation.active_cell_iterators())
 *         {
 *           m_quadrature_point_history.initialize(cell, m_n_q_points);
 *         }
 *   
 *       unsigned int material_id;
 *       double lame_lambda = 0.0;
 *       double lame_mu = 0.0;
 *       double length_scale = 0.0;
 *       double gc = 0.0;
 *       double viscosity = 0.0;
 *       double residual_k = 0.0;
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         {
 *           material_id = cell->material_id();
 *           if (m_material_data.find(material_id) != m_material_data.end())
 *             {
 *               lame_lambda                = m_material_data[material_id][0];
 *               lame_mu                    = m_material_data[material_id][1];
 *               length_scale               = m_material_data[material_id][2];
 *               gc                         = m_material_data[material_id][3];
 *               viscosity                  = m_material_data[material_id][4];
 *               residual_k                 = m_material_data[material_id][5];
 *             }
 *           else
 *             {
 *               m_logfile << "Could not find material data for material id: " << material_id << std::endl;
 *               AssertThrow(false, ExcMessage("Could not find material data for material id."));
 *             }
 *   
 *           const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
 *             m_quadrature_point_history.get_data(cell);
 *           Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *           for (unsigned int q_point = 0; q_point < m_n_q_points; ++q_point)
 *             lqph[q_point]->setup_lqp(lame_lambda, lame_mu, length_scale,
 *                                      gc, viscosity, residual_k);
 *         }
 *     }
 *   
 *     template <int dim>
 *     BlockVector<double> PhaseFieldMonolithicSolve<dim>::get_total_solution(
 *       const BlockVector<double> &solution_delta) const
 *     {
 *       BlockVector<double> solution_total(m_solution);
 *       solution_total += solution_delta;
 *       return solution_total;
 *     }
 *   
 *     template <int dim>
 *     void
 *     PhaseFieldMonolithicSolve<dim>::update_qph_incremental(const BlockVector<double> &solution_delta,
 *                                                            const BlockVector<double> &solution_old,
 *                                                            const bool is_print)
 *     {
 *       m_timer.enter_subsection("Update QPH data");
 *       if (is_print && m_parameters.m_output_iteration_history)
 *         m_logfile << " UQPH " << std::flush;
 *   
 *       const BlockVector<double> solution_total(get_total_solution(solution_delta));
 *   
 *       const UpdateFlags uf_UQPH(update_values | update_gradients);
 *       PerTaskData_UQPH  per_task_data_UQPH;
 *       ScratchData_UQPH  scratch_data_UQPH(m_fe,
 *                                           m_qf_cell,
 *                                           uf_UQPH,
 *                                           solution_total,
 *                                           solution_old,
 *                                           m_time.get_delta_t());
 *   
 *       auto worker = [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                            ScratchData_UQPH & scratch,
 *                            PerTaskData_UQPH & data)
 *         {
 *           this->update_qph_incremental_one_cell(cell, scratch, data);
 *         };
 *   
 *       auto copier = [this](const PerTaskData_UQPH &data)
 *         {
 *           this->copy_local_to_global_UQPH(data);
 *         };
 *   
 *       WorkStream::run(
 *           m_dof_handler.begin_active(),
 *           m_dof_handler.end(),
 *           worker,
 *           copier,
 *           scratch_data_UQPH,
 *           per_task_data_UQPH);
 *   
 *       m_timer.leave_subsection();
 *     }
 *   
 *     template <int dim>
 *     struct PhaseFieldMonolithicSolve<dim>::PerTaskData_UQPH
 *     {
 *       void reset()
 *       {}
 *     };
 *   
 *     template <int dim>
 *     struct PhaseFieldMonolithicSolve<dim>::ScratchData_UQPH
 *     {
 *       const BlockVector<double> & m_solution_UQPH;
 *   
 *       std::vector<SymmetricTensor<2, dim>> m_solution_symm_grads_u_cell;
 *       std::vector<double>         m_solution_values_phasefield_cell;
 *       std::vector<Tensor<1, dim>> m_solution_grad_phasefield_cell;
 *   
 *       FEValues<dim> m_fe_values;
 *   
 *       const BlockVector<double>&       m_solution_previous_step;
 *       std::vector<double>              m_phasefield_previous_step_cell;
 *   
 *       const double                     m_delta_time;
 *   
 *       ScratchData_UQPH(const FiniteElement<dim> & fe_cell,
 *                        const QGauss<dim> &        qf_cell,
 *                        const UpdateFlags          uf_cell,
 *                        const BlockVector<double> &solution_total,
 *                        const BlockVector<double> &solution_old,
 *                        const double delta_time)
 *         : m_solution_UQPH(solution_total)
 *         , m_solution_symm_grads_u_cell(qf_cell.size())
 *         , m_solution_values_phasefield_cell(qf_cell.size())
 *         , m_solution_grad_phasefield_cell(qf_cell.size())
 *         , m_fe_values(fe_cell, qf_cell, uf_cell)
 *         , m_solution_previous_step(solution_old)
 *         , m_phasefield_previous_step_cell(qf_cell.size())
 *         , m_delta_time(delta_time)
 *       {}
 *   
 *       ScratchData_UQPH(const ScratchData_UQPH &rhs)
 *         : m_solution_UQPH(rhs.m_solution_UQPH)
 *         , m_solution_symm_grads_u_cell(rhs.m_solution_symm_grads_u_cell)
 *         , m_solution_values_phasefield_cell(rhs.m_solution_values_phasefield_cell)
 *         , m_solution_grad_phasefield_cell(rhs.m_solution_grad_phasefield_cell)
 *         , m_fe_values(rhs.m_fe_values.get_fe(),
 *                       rhs.m_fe_values.get_quadrature(),
 *                       rhs.m_fe_values.get_update_flags())
 *         , m_solution_previous_step(rhs.m_solution_previous_step)
 *         , m_phasefield_previous_step_cell(rhs.m_phasefield_previous_step_cell)
 *         , m_delta_time(rhs.m_delta_time)
 *       {}
 *   
 *       void reset()
 *       {
 *         const unsigned int n_q_points = m_solution_symm_grads_u_cell.size();
 *         for (unsigned int q = 0; q < n_q_points; ++q)
 *           {
 *             m_solution_symm_grads_u_cell[q]  = 0.0;
 *             m_solution_values_phasefield_cell[q] = 0.0;
 *             m_solution_grad_phasefield_cell[q] = 0.0;
 *             m_phasefield_previous_step_cell[q] = 0.0;
 *           }
 *       }
 *     };
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::update_qph_incremental_one_cell(
 *       const typename DoFHandler<dim>::active_cell_iterator &cell,
 *       ScratchData_UQPH & scratch,
 *       PerTaskData_UQPH & /*data*/)
 *     {
 *       scratch.reset();
 *   
 *       scratch.m_fe_values.reinit(cell);
 *   
 *       const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
 *         m_quadrature_point_history.get_data(cell);
 *       Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *       const FEValuesExtractors::Vector displacement(0);
 *   
 *       scratch.m_fe_values[m_u_fe].get_function_symmetric_gradients(
 *         scratch.m_solution_UQPH, scratch.m_solution_symm_grads_u_cell);
 *       scratch.m_fe_values[m_d_fe].get_function_values(
 *         scratch.m_solution_UQPH, scratch.m_solution_values_phasefield_cell);
 *       scratch.m_fe_values[m_d_fe].get_function_gradients(
 *         scratch.m_solution_UQPH, scratch.m_solution_grad_phasefield_cell);
 *   
 *       scratch.m_fe_values[m_d_fe].get_function_values(
 *         scratch.m_solution_previous_step, scratch.m_phasefield_previous_step_cell);
 *   
 *       for (const unsigned int q_point :
 *            scratch.m_fe_values.quadrature_point_indices())
 *         lqph[q_point]->update_field_values(scratch.m_solution_symm_grads_u_cell[q_point],
 *                                            scratch.m_solution_values_phasefield_cell[q_point],
 *                                            scratch.m_solution_grad_phasefield_cell[q_point],
 *                                            scratch.m_phasefield_previous_step_cell[q_point],
 *                                            scratch.m_delta_time);
 *     }
 *   
 *     template <int dim>
 *     struct PhaseFieldMonolithicSolve<dim>::PerTaskData_ASM
 *     {
 *       FullMatrix<double>                   m_cell_matrix;
 *       Vector<double>                       m_cell_rhs;
 *       std::vector<types::global_dof_index> m_local_dof_indices;
 *   
 *       PerTaskData_ASM(const unsigned int dofs_per_cell)
 *         : m_cell_matrix(dofs_per_cell, dofs_per_cell)
 *         , m_cell_rhs(dofs_per_cell)
 *         , m_local_dof_indices(dofs_per_cell)
 *       {}
 *   
 *       void reset()
 *       {
 *         m_cell_matrix = 0.0;
 *         m_cell_rhs    = 0.0;
 *       }
 *     };
 *   
 *     template <int dim>
 *     struct PhaseFieldMonolithicSolve<dim>::PerTaskData_ASM_RHS_BFGS
 *     {
 *       Vector<double>                       m_cell_rhs;
 *       std::vector<types::global_dof_index> m_local_dof_indices;
 *   
 *       PerTaskData_ASM_RHS_BFGS(const unsigned int dofs_per_cell)
 *         : m_cell_rhs(dofs_per_cell)
 *         , m_local_dof_indices(dofs_per_cell)
 *       {}
 *   
 *       void reset()
 *       {
 *         m_cell_rhs    = 0.0;
 *       }
 *     };
 *   
 *     template <int dim>
 *     struct PhaseFieldMonolithicSolve<dim>::ScratchData_ASM
 *     {
 *       FEValues<dim>     m_fe_values;
 *       FEFaceValues<dim> m_fe_face_values;
 *   
 *       std::vector<std::vector<double>>                  m_Nx_phasefield;      // shape function values for phase-field
 *       std::vector<std::vector<Tensor<1, dim>>>          m_grad_Nx_phasefield; // gradient of shape function values for phase field
 *   
 *       std::vector<std::vector<Tensor<1, dim>>>          m_Nx_disp;       // shape function values for displacement
 *       std::vector<std::vector<Tensor<2, dim>>>          m_grad_Nx_disp;  // gradient of shape function values for displacement
 *       std::vector<std::vector<SymmetricTensor<2, dim>>> m_symm_grad_Nx_disp;  // symmetric gradient of shape function values for displacement
 *   
 *       const BlockVector<double>&       m_solution_previous_step;
 *       std::vector<double>              m_phasefield_previous_step_cell;
 *   
 *       ScratchData_ASM(const FiniteElement<dim> & fe_cell,
 *                       const QGauss<dim> &        qf_cell,
 *                       const UpdateFlags          uf_cell,
 *                       const QGauss<dim - 1> &    qf_face,
 *                       const UpdateFlags          uf_face,
 *                       const BlockVector<double>& solution_old)
 *         : m_fe_values(fe_cell, qf_cell, uf_cell)
 *         , m_fe_face_values(fe_cell, qf_face, uf_face)
 *         , m_Nx_phasefield(qf_cell.size(),
 *                           std::vector<double>(fe_cell.n_dofs_per_cell()))
 *         , m_grad_Nx_phasefield(qf_cell.size(),
 *                                std::vector<Tensor<1, dim>>(fe_cell.n_dofs_per_cell()))
 *         , m_Nx_disp(qf_cell.size(),
 *                     std::vector<Tensor<1, dim>>(fe_cell.n_dofs_per_cell()))
 *         , m_grad_Nx_disp(qf_cell.size(),
 *                          std::vector<Tensor<2, dim>>(fe_cell.n_dofs_per_cell()))
 *         , m_symm_grad_Nx_disp(qf_cell.size(),
 *                               std::vector<SymmetricTensor<2, dim>>(fe_cell.n_dofs_per_cell()))
 *         , m_solution_previous_step(solution_old)
 *         , m_phasefield_previous_step_cell(qf_cell.size())
 *       {}
 *   
 *       ScratchData_ASM(const ScratchData_ASM &rhs)
 *         : m_fe_values(rhs.m_fe_values.get_fe(),
 *                       rhs.m_fe_values.get_quadrature(),
 *                       rhs.m_fe_values.get_update_flags())
 *         , m_fe_face_values(rhs.m_fe_face_values.get_fe(),
 *                            rhs.m_fe_face_values.get_quadrature(),
 *                            rhs.m_fe_face_values.get_update_flags())
 *         , m_Nx_phasefield(rhs.m_Nx_phasefield)
 *         , m_grad_Nx_phasefield(rhs.m_grad_Nx_phasefield)
 *         , m_Nx_disp(rhs.m_Nx_disp)
 *         , m_grad_Nx_disp(rhs.m_grad_Nx_disp)
 *         , m_symm_grad_Nx_disp(rhs.m_symm_grad_Nx_disp)
 *         , m_solution_previous_step(rhs.m_solution_previous_step)
 *         , m_phasefield_previous_step_cell(rhs.m_phasefield_previous_step_cell)
 *       {}
 *   
 *       void reset()
 *       {
 *         const unsigned int n_q_points      = m_Nx_phasefield.size();
 *         const unsigned int n_dofs_per_cell = m_Nx_phasefield[0].size();
 *         for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *           {
 *             Assert(m_Nx_phasefield[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             Assert(m_grad_Nx_phasefield[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             Assert(m_Nx_disp[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             Assert(m_grad_Nx_disp[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             Assert(m_symm_grad_Nx_disp[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             m_phasefield_previous_step_cell[q_point] = 0.0;
 *             for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
 *               {
 *                 m_Nx_phasefield[q_point][k]           = 0.0;
 *                 m_grad_Nx_phasefield[q_point][k]      = 0.0;
 *                 m_Nx_disp[q_point][k]                 = 0.0;
 *                 m_grad_Nx_disp[q_point][k]            = 0.0;
 *                 m_symm_grad_Nx_disp[q_point][k]       = 0.0;
 *               }
 *           }
 *       }
 *     };
 *   
 *   
 *     template <int dim>
 *     struct PhaseFieldMonolithicSolve<dim>::ScratchData_ASM_RHS_BFGS
 *     {
 *       FEValues<dim>     m_fe_values;
 *       FEFaceValues<dim> m_fe_face_values;
 *   
 *       std::vector<std::vector<double>>                  m_Nx_phasefield;      // shape function values for phase-field
 *       std::vector<std::vector<Tensor<1, dim>>>          m_grad_Nx_phasefield; // gradient of shape function values for phase field
 *   
 *       std::vector<std::vector<Tensor<1, dim>>>          m_Nx_disp;       // shape function values for displacement
 *       std::vector<std::vector<Tensor<2, dim>>>          m_grad_Nx_disp;  // gradient of shape function values for displacement
 *       std::vector<std::vector<SymmetricTensor<2, dim>>> m_symm_grad_Nx_disp;  // symmetric gradient of shape function values for displacement
 *   
 *       const BlockVector<double>&       m_solution_previous_step;
 *       std::vector<double>              m_phasefield_previous_step_cell;
 *   
 *       ScratchData_ASM_RHS_BFGS(const FiniteElement<dim> & fe_cell,
 *                                const QGauss<dim> &        qf_cell,
 *                                const UpdateFlags          uf_cell,
 *                                const QGauss<dim - 1> &    qf_face,
 *                                const UpdateFlags          uf_face,
 *                                const BlockVector<double>& solution_old)
 *         : m_fe_values(fe_cell, qf_cell, uf_cell)
 *         , m_fe_face_values(fe_cell, qf_face, uf_face)
 *         , m_Nx_phasefield(qf_cell.size(),
 *                           std::vector<double>(fe_cell.n_dofs_per_cell()))
 *         , m_grad_Nx_phasefield(qf_cell.size(),
 *                                std::vector<Tensor<1, dim>>(fe_cell.n_dofs_per_cell()))
 *         , m_Nx_disp(qf_cell.size(),
 *                     std::vector<Tensor<1, dim>>(fe_cell.n_dofs_per_cell()))
 *         , m_grad_Nx_disp(qf_cell.size(),
 *                          std::vector<Tensor<2, dim>>(fe_cell.n_dofs_per_cell()))
 *         , m_symm_grad_Nx_disp(qf_cell.size(),
 *                               std::vector<SymmetricTensor<2, dim>>(fe_cell.n_dofs_per_cell()))
 *         , m_solution_previous_step(solution_old)
 *         , m_phasefield_previous_step_cell(qf_cell.size())
 *       {}
 *   
 *       ScratchData_ASM_RHS_BFGS(const ScratchData_ASM_RHS_BFGS &rhs)
 *         : m_fe_values(rhs.m_fe_values.get_fe(),
 *                       rhs.m_fe_values.get_quadrature(),
 *                       rhs.m_fe_values.get_update_flags())
 *         , m_fe_face_values(rhs.m_fe_face_values.get_fe(),
 *                            rhs.m_fe_face_values.get_quadrature(),
 *                            rhs.m_fe_face_values.get_update_flags())
 *         , m_Nx_phasefield(rhs.m_Nx_phasefield)
 *         , m_grad_Nx_phasefield(rhs.m_grad_Nx_phasefield)
 *         , m_Nx_disp(rhs.m_Nx_disp)
 *         , m_grad_Nx_disp(rhs.m_grad_Nx_disp)
 *         , m_symm_grad_Nx_disp(rhs.m_symm_grad_Nx_disp)
 *         , m_solution_previous_step(rhs.m_solution_previous_step)
 *         , m_phasefield_previous_step_cell(rhs.m_phasefield_previous_step_cell)
 *       {}
 *   
 *       void reset()
 *       {
 *         const unsigned int n_q_points      = m_Nx_phasefield.size();
 *         const unsigned int n_dofs_per_cell = m_Nx_phasefield[0].size();
 *         for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *           {
 *             Assert(m_Nx_phasefield[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             Assert(m_grad_Nx_phasefield[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             Assert(m_Nx_disp[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             Assert(m_grad_Nx_disp[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             Assert(m_symm_grad_Nx_disp[q_point].size() == n_dofs_per_cell,
 *                    ExcInternalError());
 *   
 *             m_phasefield_previous_step_cell[q_point] = 0.0;
 *             for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
 *               {
 *                 m_Nx_phasefield[q_point][k]           = 0.0;
 *                 m_grad_Nx_phasefield[q_point][k]      = 0.0;
 *                 m_Nx_disp[q_point][k]                 = 0.0;
 *                 m_grad_Nx_disp[q_point][k]            = 0.0;
 *                 m_symm_grad_Nx_disp[q_point][k]       = 0.0;
 *               }
 *           }
 *       }
 *     };
 *   
 * @endcode
 * 
 * constructor has no return type
 * 
 * @code
 *     template <int dim>
 *     PhaseFieldMonolithicSolve<dim>::PhaseFieldMonolithicSolve(const std::string &input_file)
 *       : m_parameters(input_file)
 *       , m_triangulation(Triangulation<dim>::maximum_smoothing)
 *       , m_time(m_parameters.m_end_time)
 *       , m_logfile(m_parameters.m_logfile_name)
 *       , m_timer(m_logfile, TimerOutput::summary, TimerOutput::wall_times)
 *       , m_dof_handler(m_triangulation)
 *       , m_fe(FE_Q<dim>(m_parameters.m_poly_degree),
 *              dim, // displacement
 *              FE_Q<dim>(m_parameters.m_poly_degree),
 *              1)   // phasefield
 *       , m_dofs_per_cell(m_fe.n_dofs_per_cell())
 *       , m_u_fe(m_first_u_component)
 *       , m_d_fe(m_d_component)
 *       , m_dofs_per_block(m_n_blocks)
 *       , m_qf_cell(m_parameters.m_quad_order)
 *       , m_qf_face(m_parameters.m_quad_order)
 *       , m_n_q_points(m_qf_cell.size())
 *       , m_vol_reference(0.0)
 *     {}
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid()
 *     {
 *       if (m_parameters.m_scenario == 1)
 *         make_grid_case_1();
 *       else if (m_parameters.m_scenario == 2)
 *         make_grid_case_2();
 *       else if (m_parameters.m_scenario == 3)
 *         make_grid_case_3();
 *       else if (m_parameters.m_scenario == 4)
 *         make_grid_case_4();
 *       else if (m_parameters.m_scenario == 5)
 *         make_grid_case_5();
 *       else if (m_parameters.m_scenario == 6)
 *         make_grid_case_6();
 *       else if (m_parameters.m_scenario == 7)
 *         make_grid_case_7();
 *       else if (m_parameters.m_scenario == 8)
 *         make_grid_case_8();
 *       else if (m_parameters.m_scenario == 9)
 *         make_grid_case_9();
 *       else if (m_parameters.m_scenario == 11)
 *         make_grid_case_11();
 *       else
 *         Assert(false, ExcMessage("The scenario has not been implemented!"));
 *   
 *       m_logfile << "\t\tTriangulation:"
 *                 << "\n\t\t\tNumber of active cells: "
 *                 << m_triangulation.n_active_cells()
 *                 << "\n\t\t\tNumber of used vertices: "
 *                 << m_triangulation.n_used_vertices()
 *                 << std::endl;
 *   
 *       std::ofstream out("original_mesh.vtu");
 *       GridOut       grid_out;
 *       grid_out.write_vtu(m_triangulation, out);
 *   
 *       m_vol_reference = GridTools::volume(m_triangulation);
 *       m_logfile << "\t\tGrid:\n\t\t\tReference volume: " << m_vol_reference << std::endl;
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_1()
 *     {
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\tSquare tension (unstructured)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 *   
 *       GridIn<dim> gridin;
 *       gridin.attach_triangulation(m_triangulation);
 *       std::ifstream f("square_tension_unstructured.msh");
 *       gridin.read_msh(f);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[1] + 0.5 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[1] - 0.5 ) < 1.0e-9)
 *                   face->set_boundary_id(1);
 *                 else
 *                   face->set_boundary_id(2);
 *               }
 *           }
 *   
 *       m_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       if (m_parameters.m_refinement_strategy == "pre-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
 *             {
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (   std::fabs(cell->center()[1]) < 0.01
 *                       && cell->center()[0] > 0.495)
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         cell->set_refine_flag();
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (   std::fabs(cell->center()[1] - 0.0) < 0.05
 *                       && std::fabs(cell->center()[0] - 0.5) < 0.05)
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_2()
 *     {
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\t\tSquare shear (unstructured)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 *   
 *       GridIn<dim> gridin;
 *       gridin.attach_triangulation(m_triangulation);
 *       std::ifstream f("square_shear_unstructured.msh");
 *       gridin.read_msh(f);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[1] + 0.5 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[1] - 0.5 ) < 1.0e-9)
 *                   face->set_boundary_id(1);
 *                 else if (   (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9)
 *                          || (std::fabs(face->center()[0] - 1.0 ) < 1.0e-9))
 *                   face->set_boundary_id(2);
 *                 else
 *                   face->set_boundary_id(3);
 *               }
 *           }
 *   
 *       m_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       if (m_parameters.m_refinement_strategy == "pre-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
 *             {
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    (cell->center()[0] > 0.45)
 *                        && (cell->center()[1] < 0.05) )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         cell->set_refine_flag();
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    std::fabs(cell->center()[0] - 0.5) < 0.025
 *                        && cell->center()[1] < 0.0 && cell->center()[1] > -0.025)
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_3()
 *     {
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\tSquare tension (structured)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 *   
 *       GridIn<dim> gridin;
 *       gridin.attach_triangulation(m_triangulation);
 *       std::ifstream f("square_tension_structured.msh");
 *       gridin.read_msh(f);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[1] - 1.0 ) < 1.0e-9)
 *                   face->set_boundary_id(1);
 *                 else
 *                   face->set_boundary_id(2);
 *               }
 *           }
 *   
 *       m_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       if (m_parameters.m_refinement_strategy == "pre-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
 *             {
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    (std::fabs(cell->center()[1] - 0.5) < 0.025)
 *                        && (cell->center()[0] > 0.475) )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         cell->set_refine_flag();
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    std::fabs(cell->center()[0] - 0.5) < 0.025
 *                        && std::fabs(cell->center()[1] - 0.5) < 0.025 )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_4()
 *     {
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\t\tSquare shear (structured)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 *   
 *       GridIn<dim> gridin;
 *       gridin.attach_triangulation(m_triangulation);
 *       std::ifstream f("square_shear_structured.msh");
 *       gridin.read_msh(f);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[1] - 1.0 ) < 1.0e-9)
 *                   face->set_boundary_id(1);
 *                 else if (   (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9)
 *                          || (std::fabs(face->center()[0] - 1.0 ) < 1.0e-9))
 *                   face->set_boundary_id(2);
 *                 else
 *                   face->set_boundary_id(3);
 *               }
 *           }
 *   
 *       m_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       if (m_parameters.m_refinement_strategy == "pre-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
 *             {
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    (cell->center()[0] > 0.475)
 *                        && (cell->center()[1] < 0.525) )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         cell->set_refine_flag();
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    std::fabs(cell->center()[0] - 0.5) < 0.025
 *                        && cell->center()[1] < 0.5 && cell->center()[1] > 0.475 )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_5()
 *     {
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\t\tThree-point bending (structured)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 *   
 *       GridIn<dim> gridin;
 *       gridin.attach_triangulation(m_triangulation);
 *       std::ifstream f("three_point_bending_structured.msh");
 *       gridin.read_msh(f);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[1] - 2.0 ) < 1.0e-9)
 *                   face->set_boundary_id(1);
 *                 else
 *                   face->set_boundary_id(2);
 *               }
 *           }
 *   
 *       m_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       if (m_parameters.m_refinement_strategy == "pre-refine")
 *         {
 *           for (const auto &cell : m_triangulation.active_cell_iterators())
 *             {
 *               if (    std::fabs(cell->center()[0] - 4.0) < 0.075
 *                    && cell->center()[1] < 1.6)
 *                 {
 *                   cell->set_refine_flag();
 *                 }
 *             }
 *           m_triangulation.execute_coarsening_and_refinement();
 *   
 *           unsigned int material_id;
 *           double length_scale;
 *           for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
 *             {
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    std::fabs(cell->center()[0] - 4.0) < 0.05
 *                        && cell->center()[1] < 1.6)
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         cell->set_refine_flag();
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    std::fabs(cell->center()[0] - 4.0) < 0.075
 *                        && std::fabs(cell->center()[1] - 0.4) < 0.075 )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_6()
 *     {
 *       AssertThrow(dim==3, ExcMessage("The dimension has to be 3D!"));
 *   
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\t\tSphere inclusion (3D structured)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       Triangulation<dim> tria_inner;
 *       GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.5);
 *   
 *       Triangulation<dim> tria_outer;
 *       GridGenerator::hyper_shell(
 *         tria_outer, Point<dim>(), 0.5, std::sqrt(dim), 2 * dim);
 *   
 *       Triangulation<dim> tmp_triangulation;
 *   
 *       GridGenerator::merge_triangulations(tria_inner, tria_outer, tmp_triangulation);
 *   
 *       tmp_triangulation.reset_all_manifolds();
 *       tmp_triangulation.set_all_manifold_ids(0);
 *   
 *       for (const auto &cell : tmp_triangulation.cell_iterators())
 *         {
 *           for (const auto &face : cell->face_iterators())
 *             {
 *               bool face_at_sphere_boundary = true;
 *               for (const auto v : face->vertex_indices())
 *                 {
 *                   if (std::abs(face->vertex(v).norm_square() - 0.25) > 1e-12)
 *                     face_at_sphere_boundary = false;
 *                 }
 *               if (face_at_sphere_boundary)
 *                 face->set_all_manifold_ids(1);
 *             }
 *           if (cell->center().norm_square() < 0.25)
 *             cell->set_material_id(1);
 *           else
 *             cell->set_material_id(0);
 *         }
 *   
 *       tmp_triangulation.set_manifold(1, SphericalManifold<dim>());
 *   
 *       TransfiniteInterpolationManifold<dim> transfinite_manifold;
 *       transfinite_manifold.initialize(tmp_triangulation);
 *       tmp_triangulation.set_manifold(0, transfinite_manifold);
 *   
 *       tmp_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       std::set<typename Triangulation< dim >::active_cell_iterator >
 *         cells_to_remove;
 *   
 *       for (const auto &cell : tmp_triangulation.active_cell_iterators())
 *         {
 *           if (   cell->center()[0] < 0.0
 *               || cell->center()[1] < 0.0
 *               || cell->center()[2] < 0.0)
 *             {
 *               cells_to_remove.insert(cell);
 *             }
 *         }
 *   
 *       GridGenerator::create_triangulation_with_removed_cells(tmp_triangulation,
 *                                                              cells_to_remove,
 *                                                              m_triangulation);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9)
 *                   face->set_boundary_id(1);
 *                 else if (std::fabs(face->center()[2] - 0.0 ) < 1.0e-9)
 *                   face->set_boundary_id(2);
 *                 else if (std::fabs(face->center()[2] - 1.0 ) < 1.0e-9)
 *                   face->set_boundary_id(3);
 *                 else
 *                   face->set_boundary_id(4);
 *               }
 *           }
 *   
 *       if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    cell->center()[2] > 0.525
 *                        && cell->center()[2] < 0.575
 *                        && cell->center()[0] < 0.05
 *                        && cell->center()[1] < 0.05 )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::cbrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_7()
 *     {
 *       AssertThrow(dim==3, ExcMessage("The dimension has to be 3D!"));
 *   
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\t\tSphere inclusion (3D structured version 2)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       Triangulation<dim> tria_inner;
 *       GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.49);
 *   
 *       Triangulation<dim> tria_outer;
 *       GridGenerator::hyper_shell(
 *         tria_outer, Point<dim>(), 0.49, std::sqrt(dim)*0.5, 2 * dim);
 *   
 *       Triangulation<dim> cube1;
 *       GridGenerator::hyper_rectangle(cube1, Point<dim>(0, 0, 0.5), Point<dim>(1, 1, 1.5));
 *       Triangulation<dim> cube2;
 *       GridGenerator::hyper_rectangle(cube2, Point<dim>(0, 0.5, -0.5), Point<dim>(1, 1.5, 0.5));
 *       Triangulation<dim> cube3;
 *       GridGenerator::hyper_rectangle(cube3, Point<dim>(0.5, -0.5, -0.5), Point<dim>(1.5, 0.5, 0.5));
 *   
 *       Triangulation<dim> tmp_triangulation;
 *       GridGenerator::merge_triangulations({&tria_inner, &tria_outer,
 *                                            &cube1, &cube2, &cube3}, tmp_triangulation);
 *   
 *       tmp_triangulation.reset_all_manifolds();
 *       tmp_triangulation.set_all_manifold_ids(0);
 *   
 *       for (const auto &cell : tmp_triangulation.cell_iterators())
 *         {
 *           for (const auto &face : cell->face_iterators())
 *             {
 *               bool face_at_sphere_boundary = true;
 *               for (const auto v : face->vertex_indices())
 *                 {
 *                   if (std::abs(face->vertex(v).norm_square() - 0.49 * 0.49) > 1e-12)
 *                     face_at_sphere_boundary = false;
 *                 }
 *               if (face_at_sphere_boundary)
 *                 face->set_all_manifold_ids(1);
 *             }
 *           if (cell->center().norm_square() < 0.1)
 *             cell->set_material_id(1);
 *           else
 *             cell->set_material_id(0);
 *         }
 *   
 *       tmp_triangulation.set_manifold(1, SphericalManifold<dim>());
 *   
 *       TransfiniteInterpolationManifold<dim> transfinite_manifold;
 *       transfinite_manifold.initialize(tmp_triangulation);
 *       tmp_triangulation.set_manifold(0, transfinite_manifold);
 *   
 *       tmp_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       std::set<typename Triangulation< dim >::active_cell_iterator >
 *         cells_to_remove;
 *   
 *       for (const auto &cell : tmp_triangulation.active_cell_iterators())
 *         {
 *           if (   cell->center()[0] < 0.0
 *               || cell->center()[1] < 0.0
 *               || cell->center()[2] < 0.0
 *               || cell->center()[0] > 1.0
 *               || cell->center()[1] > 1.0
 *               || cell->center()[2] > 1.0)
 *             {
 *               cells_to_remove.insert(cell);
 *             }
 *         }
 *   
 *       GridGenerator::create_triangulation_with_removed_cells(tmp_triangulation,
 *                                                              cells_to_remove,
 *                                                              m_triangulation);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9)
 *                   face->set_boundary_id(1);
 *                 else if (std::fabs(face->center()[2] - 0.0 ) < 1.0e-9)
 *                   face->set_boundary_id(2);
 *                 else if (std::fabs(face->center()[2] - 1.0 ) < 1.0e-9)
 *                   face->set_boundary_id(3);
 *                 else
 *                   face->set_boundary_id(4);
 *               }
 *           }
 *   
 *       if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    cell->center()[2] > 0.505
 *                        && cell->center()[2] < 0.575
 *                        && cell->center()[0] < 0.05
 *                        && cell->center()[1] < 0.05 )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::cbrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_8()
 *     {
 *       AssertThrow(dim==3, ExcMessage("The dimension has to be 3D!"));
 *   
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\t\tSphere inclusion (3D structured version 2 with barriers)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       Triangulation<dim> tria_inner;
 *       GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.49);
 *   
 *       Triangulation<dim> tria_outer;
 *       GridGenerator::hyper_shell(
 *         tria_outer, Point<dim>(), 0.49, std::sqrt(dim)*0.5, 2 * dim);
 *   
 *       Triangulation<dim> cube1;
 *       GridGenerator::hyper_rectangle(cube1, Point<dim>(0, 0, 0.5), Point<dim>(1, 1, 1.5));
 *       Triangulation<dim> cube2;
 *       GridGenerator::hyper_rectangle(cube2, Point<dim>(0, 0.5, -0.5), Point<dim>(1, 1.5, 0.5));
 *       Triangulation<dim> cube3;
 *       GridGenerator::hyper_rectangle(cube3, Point<dim>(0.5, -0.5, -0.5), Point<dim>(1.5, 0.5, 0.5));
 *   
 *       Triangulation<dim> tmp_triangulation;
 *       GridGenerator::merge_triangulations({&tria_inner, &tria_outer,
 *                                            &cube1, &cube2, &cube3}, tmp_triangulation);
 *   
 *       tmp_triangulation.reset_all_manifolds();
 *       tmp_triangulation.set_all_manifold_ids(0);
 *   
 *       for (const auto &cell : tmp_triangulation.cell_iterators())
 *         {
 *           for (const auto &face : cell->face_iterators())
 *             {
 *               bool face_at_sphere_boundary = true;
 *               for (const auto v : face->vertex_indices())
 *                 {
 *                   if (std::abs(face->vertex(v).norm_square() - 0.49 * 0.49) > 1e-12)
 *                     face_at_sphere_boundary = false;
 *                 }
 *               if (face_at_sphere_boundary)
 *                 face->set_all_manifold_ids(1);
 *             }
 *           if (cell->center().norm_square() < 0.1)
 *             cell->set_material_id(1);
 *           else
 *             cell->set_material_id(0);
 *         }
 *   
 *       tmp_triangulation.set_manifold(1, SphericalManifold<dim>());
 *   
 *       TransfiniteInterpolationManifold<dim> transfinite_manifold;
 *       transfinite_manifold.initialize(tmp_triangulation);
 *       tmp_triangulation.set_manifold(0, transfinite_manifold);
 *   
 *       tmp_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 * @endcode
 * 
 * some extra barriers
 * 
 * @code
 *       for (const auto &cell : tmp_triangulation.cell_iterators())
 *         {
 *           if (    std::fabs(cell->center()[1] - 0.75) < 0.05
 *                && std::fabs(cell->center()[2] - 0.5625) < 0.05
 *                && std::fabs(cell->center()[0] - 0.0) < 0.2)
 *             cell->set_material_id(1);
 *   
 *           if (    std::fabs(cell->center()[1] - 0.0) < 0.2
 *                && std::fabs(cell->center()[2] - 0.5) < 0.1
 *                && std::fabs(cell->center()[0] - 0.75) < 0.05)
 *             cell->set_material_id(1);
 *         }
 *   
 *       std::set<typename Triangulation< dim >::active_cell_iterator >
 *         cells_to_remove;
 *   
 *       for (const auto &cell : tmp_triangulation.active_cell_iterators())
 *         {
 *           if (   cell->center()[0] < 0.0
 *               || cell->center()[1] < 0.0
 *               || cell->center()[2] < 0.0
 *               || cell->center()[0] > 1.0
 *               || cell->center()[1] > 1.0
 *               || cell->center()[2] > 1.0)
 *             {
 *               cells_to_remove.insert(cell);
 *             }
 *         }
 *   
 *       GridGenerator::create_triangulation_with_removed_cells(tmp_triangulation,
 *                                                              cells_to_remove,
 *                                                              m_triangulation);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9)
 *                   face->set_boundary_id(1);
 *                 else if (std::fabs(face->center()[2] - 0.0 ) < 1.0e-9)
 *                   face->set_boundary_id(2);
 *                 else if (std::fabs(face->center()[2] - 1.0 ) < 1.0e-9)
 *                   face->set_boundary_id(3);
 *                 else
 *                   face->set_boundary_id(4);
 *               }
 *           }
 *   
 *       if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    cell->center()[2] > 0.505
 *                        && cell->center()[2] < 0.575
 *                        && cell->center()[0] < 0.05
 *                        && cell->center()[1] < 0.05 )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::cbrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_9()
 *     {
 *       AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 *   
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\t\tL-shape bending (2D structured)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       GridIn<dim> gridin;
 *       gridin.attach_triangulation(m_triangulation);
 *       std::ifstream f("L-Shape.msh");
 *       gridin.read_msh(f);
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9 )
 *                   face->set_boundary_id(0);
 *                 else
 *                   face->set_boundary_id(1);
 *               }
 *           }
 *   
 *       m_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       if (m_parameters.m_refinement_strategy == "pre-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
 *             {
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (    (cell->center()[1] > 242.0)
 *                        && (cell->center()[1] < 312.5)
 *                        && (cell->center()[0] < 258.0) )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         cell->set_refine_flag();
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (             (cell->center()[0] - 250) < 0.0
 *                        &&          (cell->center()[0] - 240) > 0.0
 *                        && std::fabs(cell->center()[1] - 250) < 10.0 )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::sqrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_grid_case_11()
 *     {
 *       AssertThrow(dim==3, ExcMessage("The dimension has to be 3D!"));
 *   
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *       m_logfile << "\t\t\t\tBrokenshire torsion (3D structured)" << std::endl;
 *       for (unsigned int i = 0; i < 80; ++i)
 *         m_logfile << "*";
 *       m_logfile << std::endl;
 *   
 *       Triangulation<2> triangulation_2d;
 *   
 *       double const length = 200.0;
 *       double const width = 50.0;
 *       double const height = 50.0;
 *       double const delta_L = 25.0;
 *       double const tan_theta = delta_L / (0.5*width);
 *   
 *       std::vector<unsigned int> repetitions(2, 1);
 *       repetitions[0] = 20;
 *       repetitions[1] = 5;
 *   
 *       Point<2> point1(0.0, 0.0);
 *       Point<2> point2(length, width);
 *   
 *       GridGenerator::subdivided_hyper_rectangle(triangulation_2d,
 *                                                 repetitions,
 *                                                 point1,
 *                                                 point2 );
 *   
 *       typename Triangulation<2>::vertex_iterator vertex_ptr;
 *       vertex_ptr = triangulation_2d.begin_active_vertex();
 *       while (vertex_ptr != triangulation_2d.end_vertex())
 *         {
 *           Point<2> & vertex_point = vertex_ptr->vertex();
 *   
 *           const double delta_x = (vertex_point(1) - 0.5*width) * tan_theta;
 *   
 *           if (std::fabs(vertex_point(0) - 0.5*length) < 1.0e-6)
 *             {
 *               vertex_point(0) += delta_x;
 *             }
 *           else if (std::fabs(vertex_point(0) + length/repetitions[0] - 0.5*length) < 1.0e-6)
 *             {
 *               vertex_point(0) += (delta_x + length/repetitions[0]*0.5);
 *             }
 *           else if (std::fabs(vertex_point(0) - length/repetitions[0] - 0.5*length) < 1.0e-6)
 *             {
 *               vertex_point(0) += (delta_x - length/repetitions[0]*0.5);
 *             }
 *           else if (vertex_point(0) < 0.5*length - length/repetitions[0] - 1.0e-6)
 *             {
 *               vertex_point(0) += (delta_x + length/repetitions[0]*0.5) * vertex_point(0)/(0.5*length - length/repetitions[0]);
 *             }
 *           else if (vertex_point(0) > 0.5*length + length/repetitions[0] + 1.0e-6)
 *             {
 *               vertex_point(0) += (delta_x - length/repetitions[0]*0.5) * (length - vertex_point(0))/(0.5*length - length/repetitions[0]);
 *             }
 *   
 *           ++vertex_ptr;
 *         }
 *   
 *       Triangulation<dim> tmp_triangulation;
 *       const unsigned int n_layer = repetitions[1] + 1;
 *       GridGenerator::extrude_triangulation(triangulation_2d, n_layer, height, tmp_triangulation);
 *   
 *       tmp_triangulation.refine_global(m_parameters.m_global_refine_times);
 *   
 *       std::set<typename Triangulation< dim >::active_cell_iterator >
 *         cells_to_remove;
 *   
 *       for (const auto &cell : tmp_triangulation.active_cell_iterators())
 *         {
 *           if (    (std::fabs(cell->center()[0] - (cell->center()[1] - 0.5*width)*tan_theta - 0.5*length) < 2.5)
 *                && cell->center()[2] > 0.5* height  )
 *             {
 *               cells_to_remove.insert(cell);
 *             }
 *         }
 *   
 *       GridGenerator::create_triangulation_with_removed_cells(tmp_triangulation,
 *                                                              cells_to_remove,
 *                                                              m_triangulation);
 *   
 *       if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *         {
 *           unsigned int material_id;
 *           double length_scale;
 *           bool initiation_point_refine_unfinished = true;
 *           while (initiation_point_refine_unfinished)
 *             {
 *               initiation_point_refine_unfinished = false;
 *               for (const auto &cell : m_triangulation.active_cell_iterators())
 *                 {
 *                   if (  (std::fabs(cell->center()[0] - (cell->center()[1] - 0.5*width)*tan_theta - 0.5*length) < 5.0)
 *                       && cell->center()[2] <= 0.5*height
 *                       && cell->center()[2] > 0.5*height - 5.0 )
 *                     {
 *                       material_id = cell->material_id();
 *                       length_scale = m_material_data[material_id][2];
 *                       if (  std::cbrt(cell->measure())
 *                           > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                         {
 *                           cell->set_refine_flag();
 *                           initiation_point_refine_unfinished = true;
 *                         }
 *                     }
 *                 }
 *               m_triangulation.execute_coarsening_and_refinement();
 *             }
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected mesh refinement strategy not implemented!"));
 *         }
 *   
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             if (face->at_boundary() == true)
 *               {
 *                 if (std::fabs(face->center()[0] - length) < 1.0e-6 )
 *                   face->set_boundary_id(0);
 *                 else if (std::fabs(face->center()[0] - 0.0) < 1.0e-6 )
 *                   face->set_boundary_id(1);
 *                 else
 *                   face->set_boundary_id(2);
 *               }
 *           }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::setup_system()
 *     {
 *       m_timer.enter_subsection("Setup system");
 *   
 *       std::vector<unsigned int> block_component(m_n_components,
 *                                                 m_u_dof); // displacement
 *       block_component[m_d_component] = m_d_dof;           // phasefield
 *   
 *       m_dof_handler.distribute_dofs(m_fe);
 *       DoFRenumbering::Cuthill_McKee(m_dof_handler);
 *       DoFRenumbering::component_wise(m_dof_handler, block_component);
 *   
 *       m_constraints.clear();
 *       DoFTools::make_hanging_node_constraints(m_dof_handler, m_constraints);
 *       m_constraints.close();
 *   
 *       m_dofs_per_block =
 *         DoFTools::count_dofs_per_fe_block(m_dof_handler, block_component);
 *   
 *       m_logfile << "\t\tTriangulation:"
 *                 << "\n\t\t\t Number of active cells: "
 *                 << m_triangulation.n_active_cells()
 *                 << "\n\t\t\t Number of used vertices: "
 *                 << m_triangulation.n_used_vertices()
 *                 << "\n\t\t\t Number of active edges: "
 *                 << m_triangulation.n_active_lines()
 *                 << "\n\t\t\t Number of active faces: "
 *                 << m_triangulation.n_active_faces()
 *                 << "\n\t\t\t Number of degrees of freedom (total): "
 *                 << m_dof_handler.n_dofs()
 *                 << "\n\t\t\t Number of degrees of freedom (disp): "
 *                 << m_dofs_per_block[m_u_dof]
 *                 << "\n\t\t\t Number of degrees of freedom (phasefield): "
 *                 << m_dofs_per_block[m_d_dof]
 *                 << std::endl;
 *   
 *       m_tangent_matrix.clear();
 *       {
 *         BlockDynamicSparsityPattern dsp(m_dofs_per_block, m_dofs_per_block);
 *   
 *         Table<2, DoFTools::Coupling> coupling(m_n_components, m_n_components);
 *         for (unsigned int ii = 0; ii < m_n_components; ++ii)
 *           for (unsigned int jj = 0; jj < m_n_components; ++jj)
 *             coupling[ii][jj] = DoFTools::always;
 *   
 *         DoFTools::make_sparsity_pattern(
 *           m_dof_handler, coupling, dsp, m_constraints, false);
 *         m_sparsity_pattern.copy_from(dsp);
 *       }
 *   
 *       m_tangent_matrix.reinit(m_sparsity_pattern);
 *   
 *       m_system_rhs.reinit(m_dofs_per_block);
 *       m_solution.reinit(m_dofs_per_block);
 *   
 *       setup_qph();
 *   
 *       m_timer.leave_subsection();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::make_constraints(const unsigned int it_nr)
 *     {
 *       const bool apply_dirichlet_bc = (it_nr == 0);
 *   
 *       if (it_nr > 1)
 *         {
 *           if (m_parameters.m_output_iteration_history)
 *             m_logfile << " --- " << std::flush;
 *           return;
 *         }
 *   
 *       if (m_parameters.m_output_iteration_history)
 *         m_logfile << " CST " << std::flush;
 *   
 *       if (apply_dirichlet_bc)
 *         {
 *           m_constraints.clear();
 *           DoFTools::make_hanging_node_constraints(m_dof_handler,
 *                                                   m_constraints);
 *   
 *           const FEValuesExtractors::Scalar x_displacement(0);
 *           const FEValuesExtractors::Scalar y_displacement(1);
 *           const FEValuesExtractors::Scalar z_displacement(2);
 *   
 *           const FEValuesExtractors::Vector displacements(0);
 *   
 *           if (   m_parameters.m_scenario == 1
 *               || m_parameters.m_scenario == 3)
 *             {
 * @endcode
 * 
 * Dirichlet B,C. bottom surface
 * 
 * @code
 *               const int boundary_id_bottom_surface = 0;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_bottom_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(y_displacement));
 *   
 *               typename Triangulation<dim>::active_vertex_iterator vertex_itr;
 *               vertex_itr = m_triangulation.begin_active_vertex();
 *               std::vector<types::global_dof_index> node_xy(m_fe.dofs_per_vertex);
 *   
 *               for (; vertex_itr != m_triangulation.end_vertex(); ++vertex_itr)
 *                 {
 *                   if (   (std::fabs(vertex_itr->vertex()[0] - 0.0) < 1.0e-9)
 *                       && (std::fabs(vertex_itr->vertex()[1] - 0.0) < 1.0e-9) )
 *                     {
 *                       node_xy = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
 *                     }
 *                 }
 *               m_constraints.add_line(node_xy[0]);
 *               m_constraints.set_inhomogeneity(node_xy[0], 0.0);
 *   
 *               m_constraints.add_line(node_xy[1]);
 *               m_constraints.set_inhomogeneity(node_xy[1], 0.0);
 *   
 *               const int boundary_id_top_surface = 1;
 *               /*
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_top_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(x_displacement));
 *               */
 *               const double time_inc = m_time.get_delta_t();
 *               double disp_magnitude = m_time.get_magnitude();
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_top_surface,
 *                                                        Functions::ConstantFunction<dim>(
 *                                                          disp_magnitude*time_inc, m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(y_displacement));
 *             }
 *           else if (   m_parameters.m_scenario == 2
 *                    || m_parameters.m_scenario == 4)
 *             {
 * @endcode
 * 
 * Dirichlet B,C. bottom surface
 * 
 * @code
 *               const int boundary_id_bottom_surface = 0;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_bottom_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(displacements));
 *   
 *               const int boundary_id_top_surface = 1;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_top_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(y_displacement));
 *   
 *               const double time_inc = m_time.get_delta_t();
 *               double disp_magnitude = m_time.get_magnitude();
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_top_surface,
 *                                                        Functions::ConstantFunction<dim>(
 *                                                          disp_magnitude*time_inc, m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(x_displacement));
 *   
 *               const int boundary_id_side_surfaces = 2;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_side_surfaces,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(y_displacement));
 *             }
 *           else if (m_parameters.m_scenario == 5)
 *             {
 *               typename Triangulation<dim>::active_vertex_iterator vertex_itr;
 *               vertex_itr = m_triangulation.begin_active_vertex();
 *               std::vector<types::global_dof_index> node_bottomleft(m_fe.dofs_per_vertex);
 *               std::vector<types::global_dof_index> node_bottomright(m_fe.dofs_per_vertex);
 *               std::vector<types::global_dof_index> node_topcenter(m_fe.dofs_per_vertex);
 *   
 *               for (; vertex_itr != m_triangulation.end_vertex(); ++vertex_itr)
 *                 {
 *                   if (   (std::fabs(vertex_itr->vertex()[0] - 0.0) < 1.0e-9)
 *                       && (std::fabs(vertex_itr->vertex()[1] - 0.0) < 1.0e-9) )
 *                     {
 *                       node_bottomleft = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
 *                     }
 *                   if (   (std::fabs(vertex_itr->vertex()[0] - 8.0) < 1.0e-9)
 *                       && (std::fabs(vertex_itr->vertex()[1] - 0.0) < 1.0e-9) )
 *                     {
 *                       node_bottomright = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
 *                     }
 *                   if (   (std::fabs(vertex_itr->vertex()[0] - 4.0) < 1.0e-9)
 *                       && (std::fabs(vertex_itr->vertex()[1] - 2.0) < 1.0e-9) )
 *                     {
 *                       node_topcenter = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
 *                     }
 *                 }
 * @endcode
 * 
 * bottom-left node fixed in both x- and y-directions
 * 
 * @code
 *               m_constraints.add_line(node_bottomleft[0]);
 *               m_constraints.set_inhomogeneity(node_bottomleft[0], 0.0);
 *   
 *               m_constraints.add_line(node_bottomleft[1]);
 *               m_constraints.set_inhomogeneity(node_bottomleft[1], 0.0);
 *   
 * @endcode
 * 
 * bottom-right node only fixed in y-direction
 * 
 * @code
 *               m_constraints.add_line(node_bottomright[1]);
 *               m_constraints.set_inhomogeneity(node_bottomright[1], 0.0);
 *   
 * @endcode
 * 
 * top-center node applied with y-displacement
 * 
 * @code
 *               const double time_inc = m_time.get_delta_t();
 *               double disp_magnitude = m_time.get_magnitude();
 *   
 *               m_constraints.add_line(node_topcenter[1]);
 *               m_constraints.set_inhomogeneity(node_topcenter[1], disp_magnitude*time_inc);
 *             }
 *           else if (   m_parameters.m_scenario == 6
 *                    || m_parameters.m_scenario == 7
 *                    || m_parameters.m_scenario == 8)
 *             {
 *               const int x0_surface = 0;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        x0_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(x_displacement));
 *               const int y0_surface = 1;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        y0_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(y_displacement));
 *               const int z0_surface = 2;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        z0_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(z_displacement));
 *   
 *               const int z1_surface = 3;
 *               const double time_inc = m_time.get_delta_t();
 *               double disp_magnitude = m_time.get_magnitude();
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        z1_surface,
 *                                                        Functions::ConstantFunction<dim>(
 *                                                          disp_magnitude*time_inc, m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(z_displacement));
 *             }
 *           else if (m_parameters.m_scenario == 9)
 *             {
 * @endcode
 * 
 * Dirichlet B,C. bottom surface
 * 
 * @code
 *               const int boundary_id_bottom_surface = 0;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_bottom_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(displacements));
 *   
 *               typename Triangulation<dim>::active_vertex_iterator vertex_itr;
 *               vertex_itr = m_triangulation.begin_active_vertex();
 *               std::vector<types::global_dof_index> node_disp_control(m_fe.dofs_per_vertex);
 *   
 *               for (; vertex_itr != m_triangulation.end_vertex(); ++vertex_itr)
 *                 {
 *                   if (   (std::fabs(vertex_itr->vertex()[0] - 470.0) < 1.0e-9)
 *                       && (std::fabs(vertex_itr->vertex()[1] - 250.0) < 1.0e-9) )
 *                     {
 *                       node_disp_control = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
 * @endcode
 * 
 * node applied with y-displacement
 * 
 * @code
 *                       const double time_inc = m_time.get_delta_t();
 *                       double disp_magnitude = m_time.get_magnitude();
 *   
 *                       m_constraints.add_line(node_disp_control[1]);
 *                       m_constraints.set_inhomogeneity(node_disp_control[1], disp_magnitude*time_inc);
 *                     }
 *                 }
 *             }
 *           else if (m_parameters.m_scenario == 11)
 *             {
 * @endcode
 * 
 * Dirichlet B,C. right surface
 * 
 * @code
 *               const int boundary_id_right_surface = 0;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_right_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(displacements));
 *   
 * @endcode
 * 
 * Dirichlet B,C. left surface
 * 
 * @code
 *               const int boundary_id_left_surface = 1;
 *               VectorTools::interpolate_boundary_values(m_dof_handler,
 *                                                        boundary_id_left_surface,
 *                                                        Functions::ZeroFunction<dim>(m_n_components),
 *                                                        m_constraints,
 *                                                        m_fe.component_mask(x_displacement));
 *   
 *               typename Triangulation<dim>::active_vertex_iterator vertex_itr;
 *               vertex_itr = m_triangulation.begin_active_vertex();
 *               std::vector<types::global_dof_index> node_rotate(m_fe.dofs_per_vertex);
 *               double node_dist = 0.0;
 *               double disp_mag = 0.0;
 *               double angle_theta = 0.0;
 *               double disp_y = 0;
 *               double disp_z = 0;
 *   
 *               for (; vertex_itr != m_triangulation.end_vertex(); ++vertex_itr)
 *                 {
 *                   if (std::fabs(vertex_itr->vertex()[0] - 0.0) < 1.0e-9)
 *                     {
 *                       node_rotate = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
 *                       node_dist = std::sqrt(  vertex_itr->vertex()[1] * vertex_itr->vertex()[1]
 *                                             + vertex_itr->vertex()[2] * vertex_itr->vertex()[2]);
 *   
 *                       angle_theta = m_time.get_delta_t() * m_time.get_magnitude();
 *                       disp_mag = node_dist * std::tan(angle_theta);
 *   
 *                       if (node_dist > 0)
 *                         {
 *                           disp_y = vertex_itr->vertex()[2]/node_dist * disp_mag;
 *                           disp_z = -vertex_itr->vertex()[1]/node_dist * disp_mag;
 *                         }
 *                       else
 *                         {
 *                           disp_y = 0.0;
 *                           disp_z = 0.0;
 *                         }
 *   
 *                       m_constraints.add_line(node_rotate[1]);
 *                       m_constraints.set_inhomogeneity(node_rotate[1], disp_y);
 *   
 *                       m_constraints.add_line(node_rotate[2]);
 *                       m_constraints.set_inhomogeneity(node_rotate[2], disp_z);
 *                     }
 *                 }
 *             }
 *           else
 *             Assert(false, ExcMessage("The scenario has not been implemented!"));
 *         }
 *       else  // inhomogeneous constraints
 *         {
 *           if (m_constraints.has_inhomogeneities())
 *             {
 *               AffineConstraints<double> homogeneous_constraints(m_constraints);
 *               for (unsigned int dof = 0; dof != m_dof_handler.n_dofs(); ++dof)
 *                 if (homogeneous_constraints.is_inhomogeneously_constrained(dof))
 *                   homogeneous_constraints.set_inhomogeneity(dof, 0.0);
 *               m_constraints.clear();
 *               m_constraints.copy_from(homogeneous_constraints);
 *             }
 *         }
 *       m_constraints.close();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::assemble_system_newton(const BlockVector<double> & solution_old)
 *     {
 *       m_timer.enter_subsection("Assemble system");
 *   
 *       if (m_parameters.m_output_iteration_history)
 *         m_logfile << " ASM_SYS " << std::flush;
 *   
 *       m_tangent_matrix = 0.0;
 *       m_system_rhs    = 0.0;
 *   
 *       const UpdateFlags uf_cell(update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       const UpdateFlags uf_face(update_values | update_normal_vectors |
 *                                 update_JxW_values);
 *   
 *       PerTaskData_ASM per_task_data(m_fe.n_dofs_per_cell());
 *       ScratchData_ASM scratch_data(m_fe, m_qf_cell, uf_cell, m_qf_face, uf_face, solution_old);
 *   
 *       auto worker =
 *         [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                ScratchData_ASM & scratch,
 *                PerTaskData_ASM & data)
 *         {
 *           this->assemble_system_newton_one_cell(cell, scratch, data);
 *         };
 *   
 *       auto copier = [this](const PerTaskData_ASM &data)
 *         {
 *           this->m_constraints.distribute_local_to_global(data.m_cell_matrix,
 *                                                          data.m_cell_rhs,
 *                                                          data.m_local_dof_indices,
 *                                                          m_tangent_matrix,
 *                                                          m_system_rhs);
 *         };
 *   
 *       WorkStream::run(
 *         m_dof_handler.active_cell_iterators(),
 *         worker,
 *         copier,
 *         scratch_data,
 *         per_task_data);
 *   
 *       m_timer.leave_subsection();
 *     }
 *   
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::assemble_system_B0(const BlockVector<double> & solution_old)
 *     {
 *       m_timer.enter_subsection("Assemble B0");
 *   
 *       m_tangent_matrix = 0.0;
 *   
 *       const UpdateFlags uf_cell(update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       const UpdateFlags uf_face(update_values | update_normal_vectors |
 *                                 update_JxW_values);
 *   
 *       PerTaskData_ASM per_task_data(m_fe.n_dofs_per_cell());
 *       ScratchData_ASM scratch_data(m_fe, m_qf_cell, uf_cell, m_qf_face, uf_face, solution_old);
 *   
 *       auto worker =
 *         [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                ScratchData_ASM & scratch,
 *                PerTaskData_ASM & data)
 *         {
 *           this->assemble_system_B0_one_cell(cell, scratch, data);
 *         };
 *   
 *       auto copier = [this](const PerTaskData_ASM &data)
 *         {
 *           this->m_constraints.distribute_local_to_global(data.m_cell_matrix,
 *                                                          data.m_local_dof_indices,
 *                                                          m_tangent_matrix);
 *         };
 *   
 *       WorkStream::run(
 *         m_dof_handler.active_cell_iterators(),
 *         worker,
 *         copier,
 *         scratch_data,
 *         per_task_data);
 *   
 *       m_timer.leave_subsection();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::assemble_system_rhs_BFGS_parallel(const BlockVector<double> & solution_old,
 *                                                                            BlockVector<double> & system_rhs)
 *     {
 *       m_timer.enter_subsection("Assemble RHS");
 *   
 * @endcode
 * 
 * m_logfile << " A_RHS " << std::flush;
 * 
 
 * 
 * 
 * @code
 *       system_rhs = 0.0;
 *   
 *       const UpdateFlags uf_cell(update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       const UpdateFlags uf_face(update_values | update_normal_vectors |
 *                                 update_JxW_values);
 *   
 *       PerTaskData_ASM_RHS_BFGS per_task_data(m_fe.n_dofs_per_cell());
 *       ScratchData_ASM_RHS_BFGS scratch_data(m_fe, m_qf_cell, uf_cell, m_qf_face, uf_face, solution_old);
 *   
 *       auto worker =
 *         [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                ScratchData_ASM_RHS_BFGS & scratch,
 *                PerTaskData_ASM_RHS_BFGS & data)
 *         {
 *           this->assemble_system_rhs_BFGS_one_cell(cell, scratch, data);
 *         };
 *   
 *       auto copier = [this, &system_rhs](const PerTaskData_ASM_RHS_BFGS &data)
 *         {
 *           this->m_constraints.distribute_local_to_global(data.m_cell_rhs,
 *                                                          data.m_local_dof_indices,
 *                                                          system_rhs);
 *         };
 *   
 *       WorkStream::run(
 *         m_dof_handler.active_cell_iterators(),
 *         worker,
 *         copier,
 *         scratch_data,
 *         per_task_data);
 *   
 *       m_timer.leave_subsection();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::assemble_system_rhs_BFGS_one_cell(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         ScratchData_ASM_RHS_BFGS & scratch,
 *         PerTaskData_ASM_RHS_BFGS & data) const
 *     {
 *       data.reset();
 *       scratch.reset();
 *       scratch.m_fe_values.reinit(cell);
 *       cell->get_dof_indices(data.m_local_dof_indices);
 *   
 *       scratch.m_fe_values[m_d_fe].get_function_values(
 *         scratch.m_solution_previous_step, scratch.m_phasefield_previous_step_cell);
 *   
 *       const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *         m_quadrature_point_history.get_data(cell);
 *       Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *       const double time_ramp = (m_time.current() / m_time.end());
 *       std::vector<Tensor<1, dim>> rhs_values(m_n_q_points);
 *   
 *       right_hand_side(scratch.m_fe_values.get_quadrature_points(),
 *                       rhs_values,
 *                       m_parameters.m_x_component*1.0,
 *                       m_parameters.m_y_component*1.0,
 *                       m_parameters.m_z_component*1.0);
 *   
 *       const double delta_time = m_time.get_delta_t();
 *   
 *       for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
 *         {
 *           for (const unsigned int k : scratch.m_fe_values.dof_indices())
 *             {
 *               const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 *   
 *               if (k_group == m_u_dof)
 *                 {
 *                   scratch.m_Nx_disp[q_point][k] =
 *                     scratch.m_fe_values[m_u_fe].value(k, q_point);
 *                   scratch.m_grad_Nx_disp[q_point][k] =
 *                     scratch.m_fe_values[m_u_fe].gradient(k, q_point);
 *                   scratch.m_symm_grad_Nx_disp[q_point][k] =
 *                     symmetrize(scratch.m_grad_Nx_disp[q_point][k]);
 *                 }
 *               else if (k_group == m_d_dof)
 *                 {
 *                   scratch.m_Nx_phasefield[q_point][k] =
 *                     scratch.m_fe_values[m_d_fe].value(k, q_point);
 *                   scratch.m_grad_Nx_phasefield[q_point][k] =
 *                     scratch.m_fe_values[m_d_fe].gradient(k, q_point);
 *                 }
 *               else
 *                 Assert(k_group <= m_d_dof, ExcInternalError());
 *             }
 *         }
 *   
 *       for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
 *         {
 *           const double length_scale            = lqph[q_point]->get_length_scale();
 *           const double gc                      = lqph[q_point]->get_critical_energy_release_rate();
 *           const double eta                     = lqph[q_point]->get_viscosity();
 *           const double history_strain_energy   = lqph[q_point]->get_history_max_positive_strain_energy();
 *           const double current_positive_strain_energy = lqph[q_point]->get_current_positive_strain_energy();
 *   
 *           double history_value = history_strain_energy;
 *           if (current_positive_strain_energy > history_strain_energy)
 *             history_value = current_positive_strain_energy;
 *   
 *           const double phasefield_value        = lqph[q_point]->get_phase_field_value();
 *           const Tensor<1, dim> phasefield_grad = lqph[q_point]->get_phase_field_gradient();
 *   
 *           const std::vector<double>         &      N_phasefield = scratch.m_Nx_phasefield[q_point];
 *           const std::vector<Tensor<1, dim>> & grad_N_phasefield = scratch.m_grad_Nx_phasefield[q_point];
 *           const double                old_phasefield = scratch.m_phasefield_previous_step_cell[q_point];
 *   
 *           const SymmetricTensor<2, dim> & cauchy_stress = lqph[q_point]->get_cauchy_stress();
 *   
 *           const std::vector<Tensor<1,dim>> & N_disp = scratch.m_Nx_disp[q_point];
 *           const std::vector<SymmetricTensor<2, dim>> & symm_grad_N_disp =
 *             scratch.m_symm_grad_Nx_disp[q_point];
 *           const double JxW = scratch.m_fe_values.JxW(q_point);
 *   
 *           SymmetricTensor<2, dim> symm_grad_Nx_i_x_C;
 *   
 *           for (const unsigned int i : scratch.m_fe_values.dof_indices())
 *             {
 *               const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *               if (i_group == m_u_dof)
 *                 {
 *                   data.m_cell_rhs(i) += (symm_grad_N_disp[i] * cauchy_stress) * JxW;
 *   
 * @endcode
 * 
 * contributions from the body force to right-hand side
 * 
 * @code
 *                   data.m_cell_rhs(i) -= N_disp[i] * rhs_values[q_point] * JxW;
 *                 }
 *               else if (i_group == m_d_dof)
 *                 {
 *                       data.m_cell_rhs(i) += (    gc * length_scale * grad_N_phasefield[i] * phasefield_grad
 *                                               +  (   gc / length_scale * phasefield_value
 *                                                + eta / delta_time  * (phasefield_value - old_phasefield)
 *                                                + degradation_function_derivative(phasefield_value) * history_value )
 *                                             * N_phasefield[i]
 *                                         ) * JxW;
 *                 }
 *               else
 *                 Assert(i_group <= m_d_dof, ExcInternalError());
 *             }  // i
 *         }  // q_point
 *   
 * @endcode
 * 
 * if there is surface pressure, this surface pressure always applied to the
 * reference configuration
 * 
 * @code
 *       const unsigned int face_pressure_id = 100;
 *       const double p0 = 0.0;
 *   
 *       for (const auto &face : cell->face_iterators())
 *         if (face->at_boundary() && face->boundary_id() == face_pressure_id)
 *           {
 *             scratch.m_fe_face_values.reinit(cell, face);
 *   
 *             for (const unsigned int f_q_point : scratch.m_fe_face_values.quadrature_point_indices())
 *               {
 *                 const Tensor<1, dim> &N = scratch.m_fe_face_values.normal_vector(f_q_point);
 *   
 *                 const double         pressure  = p0 * time_ramp;
 *                 const Tensor<1, dim> traction  = pressure * N;
 *   
 *                 for (const unsigned int i : scratch.m_fe_values.dof_indices())
 *                   {
 *                     const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *                     if (i_group == m_u_dof)
 *                       {
 *                             const unsigned int component_i = m_fe.system_to_component_index(i).first;
 *                             const double Ni = scratch.m_fe_face_values.shape_value(i, f_q_point);
 *                             const double JxW = scratch.m_fe_face_values.JxW(f_q_point);
 *                             data.m_cell_rhs(i) -= (Ni * traction[component_i]) * JxW;
 *                       }
 *                   }
 *               }
 *           }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::assemble_system_newton_one_cell(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         ScratchData_ASM & scratch,
 *         PerTaskData_ASM & data) const
 *     {
 *       data.reset();
 *       scratch.reset();
 *       scratch.m_fe_values.reinit(cell);
 *       cell->get_dof_indices(data.m_local_dof_indices);
 *   
 *       scratch.m_fe_values[m_d_fe].get_function_values(
 *         scratch.m_solution_previous_step, scratch.m_phasefield_previous_step_cell);
 *   
 *       const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *         m_quadrature_point_history.get_data(cell);
 *       Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *       const double time_ramp = (m_time.current() / m_time.end());
 *       std::vector<Tensor<1, dim>> rhs_values(m_n_q_points);
 *   
 *       right_hand_side(scratch.m_fe_values.get_quadrature_points(),
 *                       rhs_values,
 *                       m_parameters.m_x_component*1.0,
 *                       m_parameters.m_y_component*1.0,
 *                       m_parameters.m_z_component*1.0);
 *   
 *       const double delta_time = m_time.get_delta_t();
 *   
 *       for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
 *         {
 *           for (const unsigned int k : scratch.m_fe_values.dof_indices())
 *             {
 *               const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 *   
 *               if (k_group == m_u_dof)
 *                 {
 *                   scratch.m_Nx_disp[q_point][k] =
 *                     scratch.m_fe_values[m_u_fe].value(k, q_point);
 *                   scratch.m_grad_Nx_disp[q_point][k] =
 *                     scratch.m_fe_values[m_u_fe].gradient(k, q_point);
 *                   scratch.m_symm_grad_Nx_disp[q_point][k] =
 *                     symmetrize(scratch.m_grad_Nx_disp[q_point][k]);
 *                 }
 *               else if (k_group == m_d_dof)
 *                 {
 *                   scratch.m_Nx_phasefield[q_point][k] =
 *                     scratch.m_fe_values[m_d_fe].value(k, q_point);
 *                   scratch.m_grad_Nx_phasefield[q_point][k] =
 *                     scratch.m_fe_values[m_d_fe].gradient(k, q_point);
 *                 }
 *               else
 *                 Assert(k_group <= m_d_dof, ExcInternalError());
 *             }
 *         }
 *   
 *       for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
 *         {
 *           const double length_scale            = lqph[q_point]->get_length_scale();
 *           const double gc                      = lqph[q_point]->get_critical_energy_release_rate();
 *           const double eta                     = lqph[q_point]->get_viscosity();
 *           const double history_strain_energy   = lqph[q_point]->get_history_max_positive_strain_energy();
 *           const double current_positive_strain_energy = lqph[q_point]->get_current_positive_strain_energy();
 *   
 *           double history_value = history_strain_energy;
 *           if (current_positive_strain_energy > history_strain_energy)
 *             history_value = current_positive_strain_energy;
 *   
 *           const double phasefield_value        = lqph[q_point]->get_phase_field_value();
 *           const Tensor<1, dim> phasefield_grad = lqph[q_point]->get_phase_field_gradient();
 *   
 *           const std::vector<double>         &      N_phasefield = scratch.m_Nx_phasefield[q_point];
 *           const std::vector<Tensor<1, dim>> & grad_N_phasefield = scratch.m_grad_Nx_phasefield[q_point];
 *           const double                old_phasefield = scratch.m_phasefield_previous_step_cell[q_point];
 *   
 *           const SymmetricTensor<2, dim> & cauchy_stress = lqph[q_point]->get_cauchy_stress();
 *           const SymmetricTensor<2, dim> & cauchy_stress_positive = lqph[q_point]->get_cauchy_stress_positive();
 *           const SymmetricTensor<4, dim> & mechanical_C  = lqph[q_point]->get_mechanical_C();
 *   
 *           const std::vector<Tensor<1,dim>> & N_disp = scratch.m_Nx_disp[q_point];
 *           const std::vector<SymmetricTensor<2, dim>> & symm_grad_N_disp =
 *             scratch.m_symm_grad_Nx_disp[q_point];
 *           const double JxW = scratch.m_fe_values.JxW(q_point);
 *   
 *           SymmetricTensor<2, dim> symm_grad_Nx_i_x_C;
 *   
 *           for (const unsigned int i : scratch.m_fe_values.dof_indices())
 *             {
 *               const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *               if (i_group == m_u_dof)
 *                 {
 *                   data.m_cell_rhs(i) -= (symm_grad_N_disp[i] * cauchy_stress) * JxW;
 *   
 * @endcode
 * 
 * contributions from the body force to right-hand side
 * 
 * @code
 *                   data.m_cell_rhs(i) += N_disp[i] * rhs_values[q_point] * JxW;
 *                 }
 *               else if (i_group == m_d_dof)
 *                 {
 *                       data.m_cell_rhs(i) -= (    gc * length_scale * grad_N_phasefield[i] * phasefield_grad
 *                                               +  (   gc / length_scale * phasefield_value
 *                                                + eta / delta_time  * (phasefield_value - old_phasefield)
 *                                                + degradation_function_derivative(phasefield_value) * history_value )
 *                                             * N_phasefield[i]
 *                                         ) * JxW;
 *                 }
 *               else
 *                 Assert(i_group <= m_d_dof, ExcInternalError());
 *   
 *               if (i_group == m_u_dof)
 *                 {
 *                   symm_grad_Nx_i_x_C = symm_grad_N_disp[i] * mechanical_C;
 *                 }
 *   
 *               for (const unsigned int j : scratch.m_fe_values.dof_indices())
 *                 {
 *                   const unsigned int j_group = m_fe.system_to_base_index(j).first.first;
 *   
 *                   if ((i_group == j_group) && (i_group == m_u_dof))
 *                     {
 *                       data.m_cell_matrix(i, j) += symm_grad_Nx_i_x_C * symm_grad_N_disp[j] * JxW;
 *                     }
 *                   else if ((i_group == j_group) && (i_group == m_d_dof))
 *                     {
 *                       data.m_cell_matrix(i, j) += (  (  gc/length_scale + eta/delta_time +
 *                                                         degradation_function_2nd_order_derivative(phasefield_value)
 *                                                       * history_value  )
 *                                                     * N_phasefield[i] * N_phasefield[j]
 *                                                     + gc * length_scale * grad_N_phasefield[i] * grad_N_phasefield[j]
 *                                                   ) * JxW;
 *                     }
 *                   else if ((i_group == m_u_dof) && (j_group == m_d_dof))
 *                     {
 *                       data.m_cell_matrix(i, j) +=  symm_grad_N_disp[i] * cauchy_stress_positive
 *                                                  * degradation_function_derivative(phasefield_value)
 *                                                  * N_phasefield[j] * JxW;
 *                     }
 *                   else if ((i_group == m_d_dof) && (j_group == m_u_dof))
 *                     {
 *                       if (current_positive_strain_energy > history_strain_energy)
 *                         data.m_cell_matrix(i, j) +=  N_phasefield[i]
 *                                                     * degradation_function_derivative(phasefield_value)
 *                                                    * cauchy_stress_positive
 *                                                    * symm_grad_N_disp[j]
 *                                                    * JxW;
 *                       else
 *                         data.m_cell_matrix(i, j) += 0.0;
 *                     }
 *                   else
 *                     Assert((i_group <= m_d_dof) && (j_group <= m_d_dof),
 *                            ExcInternalError());
 *                 } // j
 *             }  // i
 *         }  // q_point
 *   
 * @endcode
 * 
 * if there is surface pressure, this surface pressure always applied to the
 * reference configuration
 * 
 * @code
 *       const unsigned int face_pressure_id = 100;
 *       const double p0 = 0.0;
 *   
 *       for (const auto &face : cell->face_iterators())
 *         if (face->at_boundary() && face->boundary_id() == face_pressure_id)
 *           {
 *             scratch.m_fe_face_values.reinit(cell, face);
 *   
 *             for (const unsigned int f_q_point : scratch.m_fe_face_values.quadrature_point_indices())
 *               {
 *                 const Tensor<1, dim> &N = scratch.m_fe_face_values.normal_vector(f_q_point);
 *   
 *                 const double         pressure  = p0 * time_ramp;
 *                 const Tensor<1, dim> traction  = pressure * N;
 *   
 *                 for (const unsigned int i : scratch.m_fe_values.dof_indices())
 *                   {
 *                     const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *                     if (i_group == m_u_dof)
 *                       {
 *                             const unsigned int component_i = m_fe.system_to_component_index(i).first;
 *                             const double Ni = scratch.m_fe_face_values.shape_value(i, f_q_point);
 *                             const double JxW = scratch.m_fe_face_values.JxW(f_q_point);
 *                             data.m_cell_rhs(i) += (Ni * traction[component_i]) * JxW;
 *                       }
 *                   }
 *               }
 *           }
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::assemble_system_B0_one_cell(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         ScratchData_ASM & scratch,
 *         PerTaskData_ASM & data) const
 *     {
 *       data.reset();
 *       scratch.reset();
 *       scratch.m_fe_values.reinit(cell);
 *       cell->get_dof_indices(data.m_local_dof_indices);
 *   
 *       scratch.m_fe_values[m_d_fe].get_function_values(
 *         scratch.m_solution_previous_step, scratch.m_phasefield_previous_step_cell);
 *   
 *       const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *         m_quadrature_point_history.get_data(cell);
 *       Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *       const double delta_time = m_time.get_delta_t();
 *   
 *       for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
 *         {
 *           for (const unsigned int k : scratch.m_fe_values.dof_indices())
 *             {
 *               const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 *   
 *               if (k_group == m_u_dof)
 *                 {
 *                   scratch.m_Nx_disp[q_point][k] =
 *                     scratch.m_fe_values[m_u_fe].value(k, q_point);
 *                   scratch.m_grad_Nx_disp[q_point][k] =
 *                     scratch.m_fe_values[m_u_fe].gradient(k, q_point);
 *                   scratch.m_symm_grad_Nx_disp[q_point][k] =
 *                     symmetrize(scratch.m_grad_Nx_disp[q_point][k]);
 *                 }
 *               else if (k_group == m_d_dof)
 *                 {
 *                   scratch.m_Nx_phasefield[q_point][k] =
 *                     scratch.m_fe_values[m_d_fe].value(k, q_point);
 *                   scratch.m_grad_Nx_phasefield[q_point][k] =
 *                     scratch.m_fe_values[m_d_fe].gradient(k, q_point);
 *                 }
 *               else
 *                 Assert(k_group <= m_d_dof, ExcInternalError());
 *             }
 *         }
 *   
 *       for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
 *         {
 *           const double length_scale            = lqph[q_point]->get_length_scale();
 *           const double gc                      = lqph[q_point]->get_critical_energy_release_rate();
 *           const double eta                     = lqph[q_point]->get_viscosity();
 *           const double history_strain_energy   = lqph[q_point]->get_history_max_positive_strain_energy();
 *           const double current_positive_strain_energy = lqph[q_point]->get_current_positive_strain_energy();
 *   
 *           double history_value = history_strain_energy;
 *           if (current_positive_strain_energy > history_strain_energy)
 *             history_value = current_positive_strain_energy;
 *   
 *           const double phasefield_value        = lqph[q_point]->get_phase_field_value();
 *   
 *           const std::vector<double>         &      N_phasefield = scratch.m_Nx_phasefield[q_point];
 *           const std::vector<Tensor<1, dim>> & grad_N_phasefield = scratch.m_grad_Nx_phasefield[q_point];
 *   
 * @endcode
 * 
 * const SymmetricTensor<2, dim> & cauchy_stress_positive = lqph[q_point]->get_cauchy_stress_positive();
 * 
 * @code
 *           const SymmetricTensor<4, dim> & mechanical_C  = lqph[q_point]->get_mechanical_C();
 *   
 *           const std::vector<SymmetricTensor<2, dim>> & symm_grad_N_disp =
 *             scratch.m_symm_grad_Nx_disp[q_point];
 *           const double JxW = scratch.m_fe_values.JxW(q_point);
 *   
 *           SymmetricTensor<2, dim> symm_grad_Nx_i_x_C;
 *   
 *           for (const unsigned int i : scratch.m_fe_values.dof_indices())
 *             {
 *               const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *               if (i_group == m_u_dof)
 *                 {
 *                   symm_grad_Nx_i_x_C = symm_grad_N_disp[i] * mechanical_C;
 *                 }
 *   
 *               for (const unsigned int j : scratch.m_fe_values.dof_indices())
 *                 {
 *                   const unsigned int j_group = m_fe.system_to_base_index(j).first.first;
 *   
 *                   if ((i_group == j_group) && (i_group == m_u_dof))
 *                     {
 *                       data.m_cell_matrix(i, j) += symm_grad_Nx_i_x_C * symm_grad_N_disp[j] * JxW;
 *                     }
 *                   else if ((i_group == j_group) && (i_group == m_d_dof))
 *                     {
 *                       data.m_cell_matrix(i, j) += (  (   gc/length_scale + eta/delta_time
 *                                                        + degradation_function_2nd_order_derivative(phasefield_value)
 *                                                        * history_value  )
 *                                                     * N_phasefield[i] * N_phasefield[j]
 *                                                     + gc * length_scale * grad_N_phasefield[i] * grad_N_phasefield[j]
 *                                                   ) * JxW;
 *                     }
 *                   else
 *                     Assert((i_group <= m_d_dof) && (j_group <= m_d_dof),
 *                            ExcInternalError());
 *                 } // j
 *             }  // i
 *         }  // q_point
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::assemble_system_rhs_BFGS(const BlockVector<double> & solution_old,
 *                                                                   BlockVector<double> & system_rhs)
 *     {
 *       m_timer.enter_subsection("Assemble RHS");
 *   
 * @endcode
 * 
 * m_logfile << " A_RHS " << std::flush;
 * 
 
 * 
 * 
 * @code
 *       system_rhs = 0.0;
 *   
 *       Vector<double> cell_rhs(m_dofs_per_cell);
 *       std::vector<types::global_dof_index> local_dof_indices(m_dofs_per_cell);
 *   
 *       const double time_ramp = (m_time.current() / m_time.end());
 *       const double delta_time = m_time.get_delta_t();
 *   
 *       std::vector<Tensor<1, dim>> rhs_values(m_n_q_points);
 *       const UpdateFlags uf_cell(update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       const UpdateFlags uf_face(update_values | update_normal_vectors |
 *                                 update_JxW_values);
 *   
 *       FEValues<dim> fe_values(m_fe, m_qf_cell, uf_cell);
 *       FEFaceValues<dim> fe_face_values(m_fe, m_qf_face, uf_face);
 *   
 * @endcode
 * 
 * shape function values for displacement field
 * 
 * @code
 *       std::vector<std::vector<Tensor<1, dim>>>
 *         Nx_disp(m_qf_cell.size(), std::vector<Tensor<1, dim>>(m_dofs_per_cell));
 *       std::vector<std::vector<Tensor<2, dim>>>
 *         grad_Nx_disp(m_qf_cell.size(), std::vector<Tensor<2, dim>>(m_dofs_per_cell));
 *       std::vector<std::vector<SymmetricTensor<2, dim>>>
 *         symm_grad_Nx_disp(m_qf_cell.size(), std::vector<SymmetricTensor<2, dim>>(m_dofs_per_cell));
 *   
 * @endcode
 * 
 * shape function values for phase field
 * 
 * @code
 *       std::vector<std::vector<double>>
 *         Nx_phasefield(m_qf_cell.size(), std::vector<double>(m_dofs_per_cell));
 *       std::vector<std::vector<Tensor<1, dim>>>
 *         grad_Nx_phasefield(m_qf_cell.size(), std::vector<Tensor<1, dim>>(m_dofs_per_cell));
 *   
 *       std::vector<double> phasefield_previous_step_cell(m_qf_cell.size());
 *   
 *       for (const auto &cell : m_dof_handler.active_cell_iterators())
 *         {
 *           const std::vector<std::shared_ptr< PointHistory<dim>>> lqph =
 *             m_quadrature_point_history.get_data(cell);
 *           Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *           cell_rhs = 0.0;
 *           fe_values.reinit(cell);
 *           right_hand_side(fe_values.get_quadrature_points(),
 *                           rhs_values,
 *                           m_parameters.m_x_component*time_ramp,
 *                           m_parameters.m_y_component*time_ramp,
 *                           m_parameters.m_z_component*time_ramp);
 *   
 *           fe_values[m_d_fe].get_function_values(
 *               solution_old, phasefield_previous_step_cell);
 *   
 *           for (const unsigned int q_point : fe_values.quadrature_point_indices())
 *             {
 *               for (const unsigned int k : fe_values.dof_indices())
 *                 {
 *                   const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 *   
 *                   if (k_group == m_u_dof)
 *                     {
 *                       Nx_disp[q_point][k] = fe_values[m_u_fe].value(k, q_point);
 *                       grad_Nx_disp[q_point][k] = fe_values[m_u_fe].gradient(k, q_point);
 *                       symm_grad_Nx_disp[q_point][k] = symmetrize(grad_Nx_disp[q_point][k]);
 *                     }
 *                   else if (k_group == m_d_dof)
 *                     {
 *                       Nx_phasefield[q_point][k] = fe_values[m_d_fe].value(k, q_point);
 *                       grad_Nx_phasefield[q_point][k] = fe_values[m_d_fe].gradient(k, q_point);
 *                     }
 *                   else
 *                     Assert(k_group <= m_d_dof, ExcInternalError());
 *                 }
 *             }
 *   
 *           for (const unsigned int q_point : fe_values.quadrature_point_indices())
 *             {
 *               const double length_scale            = lqph[q_point]->get_length_scale();
 *               const double gc                      = lqph[q_point]->get_critical_energy_release_rate();
 *               const double eta                     = lqph[q_point]->get_viscosity();
 *               const double history_strain_energy   = lqph[q_point]->get_history_max_positive_strain_energy();
 *               const double current_positive_strain_energy = lqph[q_point]->get_current_positive_strain_energy();
 *   
 *               double history_value = history_strain_energy;
 *               if (current_positive_strain_energy > history_strain_energy)
 *                 history_value = current_positive_strain_energy;
 *   
 *               const double phasefield_value        = lqph[q_point]->get_phase_field_value();
 *               const Tensor<1, dim> phasefield_grad = lqph[q_point]->get_phase_field_gradient();
 *   
 *               const std::vector<double>         &      N_phasefield = Nx_phasefield[q_point];
 *               const std::vector<Tensor<1, dim>> & grad_N_phasefield = grad_Nx_phasefield[q_point];
 *               const double                old_phasefield = phasefield_previous_step_cell[q_point];
 *   
 *               const SymmetricTensor<2, dim> & cauchy_stress = lqph[q_point]->get_cauchy_stress();
 *   
 *               const std::vector<Tensor<1,dim>> & N = Nx_disp[q_point];
 *               const std::vector<SymmetricTensor<2, dim>> & symm_grad_N = symm_grad_Nx_disp[q_point];
 *               const double JxW = fe_values.JxW(q_point);
 *   
 *               for (const unsigned int i : fe_values.dof_indices())
 *                 {
 *                   const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *                   if (i_group == m_u_dof)
 *                     {
 *                       cell_rhs(i) += (symm_grad_N[i] * cauchy_stress) * JxW;
 * @endcode
 * 
 * contributions from the body force to right-hand side
 * 
 * @code
 *                       cell_rhs(i) -= N[i] * rhs_values[q_point] * JxW;
 *                     }
 *                   else if (i_group == m_d_dof)
 *                     {
 *                       cell_rhs(i) += (    gc * length_scale * grad_N_phasefield[i] * phasefield_grad
 *                                            +  (   gc / length_scale * phasefield_value
 *                                             + eta / delta_time  * (phasefield_value - old_phasefield)
 *                                             + degradation_function_derivative(phasefield_value) * history_value )
 *                                        * N_phasefield[i]
 *                                      ) * JxW;
 *                     }
 *                   else
 *                     Assert(i_group <= m_d_dof, ExcInternalError());
 *                 }
 *             }
 *   
 * @endcode
 * 
 * if there is surface pressure, this surface pressure always applied to the
 * reference configuration
 * 
 * @code
 *           const unsigned int face_pressure_id = 100;
 *           const double p0 = 0.0;
 *   
 *           for (const auto &face : cell->face_iterators())
 *             {
 *               if (face->at_boundary() && face->boundary_id() == face_pressure_id)
 *                 {
 *                   fe_face_values.reinit(cell, face);
 *   
 *                   for (const unsigned int f_q_point : fe_face_values.quadrature_point_indices())
 *                     {
 *                       const Tensor<1, dim> &N = fe_face_values.normal_vector(f_q_point);
 *   
 *                       const double         pressure  = p0 * time_ramp;
 *                       const Tensor<1, dim> traction  = pressure * N;
 *   
 *                       for (const unsigned int i : fe_values.dof_indices())
 *                         {
 *                           const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *                           if (i_group == m_u_dof)
 *                             {
 *                               const unsigned int component_i = m_fe.system_to_component_index(i).first;
 *                               const double Ni = fe_face_values.shape_value(i, f_q_point);
 *                               const double JxW = fe_face_values.JxW(f_q_point);
 *                               cell_rhs(i) -= (Ni * traction[component_i]) * JxW;
 *                             }
 *                         }
 *                     }
 *                 }
 *             }
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *           for (const unsigned int i : fe_values.dof_indices())
 *             system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *         } // for (const auto &cell : m_dof_handler.active_cell_iterators())
 *   
 *       m_timer.leave_subsection();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::update_history_field_step()
 *     {
 *       m_logfile << "\t\tUpdate history variable" << std::endl;
 *   
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         {
 *           std::vector<std::shared_ptr< PointHistory<dim>>> lqph =
 *             m_quadrature_point_history.get_data(cell);
 *           Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *           for (unsigned int q_point = 0; q_point < m_n_q_points; ++q_point)
 *             {
 *               lqph[q_point]->update_history_variable();
 *             }
 *         }
 *     }
 *   
 *     template <int dim>
 *     double PhaseFieldMonolithicSolve<dim>::line_search_stepsize_gradient_based(const BlockVector<double> & BFGS_p_vector,
 *                                                                                const BlockVector<double> & solution_delta)
 *     {
 *       BlockVector<double> g_old(m_system_rhs);
 *   
 * @endcode
 * 
 * BFGS_p_vector is the search direction
 * 
 * @code
 *       BlockVector<double> solution_delta_trial(solution_delta);
 * @endcode
 * 
 * take a full step size 1.0
 * 
 * @code
 *       solution_delta_trial.add(1.0, BFGS_p_vector);
 *   
 *       update_qph_incremental(solution_delta_trial, m_solution, false);
 *   
 *       BlockVector<double> g_new(m_dofs_per_block);
 *       assemble_system_rhs_BFGS_parallel(m_solution, g_new);
 *   
 *       BlockVector<double> y_old(m_dofs_per_block);
 *   
 *       y_old = g_new - g_old;
 *   
 *       double alpha = 1.0;
 *   
 *       double alpha_old = 0.0;
 *   
 *       double delta_alpha_old = alpha - alpha_old;
 *   
 *       double delta_alpha_new;
 *   
 *       unsigned int ls_max = 10;
 *   
 *       for (unsigned int i = 1; i <= ls_max; ++i)
 *         {
 *           delta_alpha_new = -delta_alpha_old
 *                           * (g_new * BFGS_p_vector)/(y_old * BFGS_p_vector);
 *           alpha += delta_alpha_new;
 *   
 *           if (std::fabs(delta_alpha_new) < 1.0e-5)
 *             break;
 *   
 *           if (i == ls_max)
 *             {
 *               alpha = 1.0;
 *               break;
 *             }
 *   
 *           g_old = g_new;
 *   
 * @endcode
 * 
 * BFGS_p_vector is the search direction
 * 
 * @code
 *           solution_delta_trial = solution_delta;
 *           solution_delta_trial.add(alpha, BFGS_p_vector);
 *           update_qph_incremental(solution_delta_trial, m_solution, false);
 *           assemble_system_rhs_BFGS_parallel(m_solution, g_new);
 *   
 *           y_old = g_new - g_old;
 *   
 *           delta_alpha_old = delta_alpha_new;
 *         }
 *   
 *       if (alpha < 1.0e-3)
 *         alpha = 1.0;
 *   
 *       return alpha;
 *     }
 *   
 *     template <int dim>
 *     double PhaseFieldMonolithicSolve<dim>::line_search_stepsize_strong_wolfe(const double phi_0,
 *                                                                              const double phi_0_prime,
 *                                                                              const BlockVector<double> & BFGS_p_vector,
 *                                                                              const BlockVector<double> & solution_delta)
 *     {
 * @endcode
 * 
 * AssertThrow(phi_0_prime < 0,
 * ExcMessage("The derivative of phi at alpha = 0 should be negative!"));
 * 
 
 * 
 * Some line search parameters
 * 
 * @code
 *       const double c1 = 0.0001;
 *       const double c2 = 0.9;
 *       const double alpha_max = 100.0;
 *       const unsigned int max_iter = 20;
 *       double alpha = 1.0;
 *   
 *       double phi_old = phi_0;
 *       double phi_prime_old = phi_0_prime;
 *       double alpha_old = 0.0;
 *   
 *       double phi, phi_prime;
 *   
 *       std::pair<double, double> current_phi_phi_prime;
 *   
 *       unsigned int i = 0;
 *       for (; i < max_iter; ++i)
 *         {
 *           current_phi_phi_prime = calculate_phi_and_phi_prime(alpha, BFGS_p_vector, solution_delta);
 *           phi = current_phi_phi_prime.first;
 *           phi_prime = current_phi_phi_prime.second;
 *   
 *           if (   ( phi > (phi_0 + c1 * alpha * phi_0_prime) )
 *               || ( i > 0 && phi > phi_old ) )
 *             {
 *               return line_search_zoom_strong_wolfe(phi_old, phi_prime_old, alpha_old,
 *                                                    phi,     phi_prime,     alpha,
 *                                                    phi_0,   phi_0_prime,   BFGS_p_vector,
 *                                                    c1,      c2,            max_iter, solution_delta);
 *             }
 *   
 *           if (std::fabs(phi_prime) <= c2 * std::fabs(phi_0_prime))
 *             {
 *               return alpha;
 *             }
 *   
 *           if (phi_prime >= 0)
 *             {
 *               return line_search_zoom_strong_wolfe(phi,     phi_prime,     alpha,
 *                                                    phi_old, phi_prime_old, alpha_old,
 *                                                    phi_0,   phi_0_prime,   BFGS_p_vector,
 *                                                    c1,      c2,            max_iter, solution_delta);
 *             }
 *   
 *           phi_old = phi;
 *           phi_prime_old = phi_prime;
 *           alpha_old = alpha;
 *   
 *           alpha = std::min(2.0*alpha, alpha_max);
 *   
 * @endcode
 * 
 * AssertThrow(alpha < alpha_max,
 * ExcMessage("alpha is bigger than alpha_max, line search failed!"));
 * 
 * @code
 *         }
 *   
 * @endcode
 * 
 * AssertThrow(i < max_iter,
 * ExcMessage("max number attempts arrived, line search failed!"));
 * Instead of terminating the program, we can just take a full step.
 * 
 * @code
 *       if (i == max_iter)
 *         alpha = 1.0;
 *   
 *       return alpha;
 *     }
 *   
 *     template <int dim>
 *     double PhaseFieldMonolithicSolve<dim>::
 *       line_search_zoom_strong_wolfe(double phi_low, double phi_low_prime, double alpha_low,
 *                                     double phi_high, double phi_high_prime, double alpha_high,
 *                                     double phi_0, double phi_0_prime, const BlockVector<double> & BFGS_p_vector,
 *                                     double c1, double c2, unsigned int max_iter, const BlockVector<double> & solution_delta)
 *     {
 *       double alpha = 0;
 *       std::pair<double, double> current_phi_phi_prime;
 *       double phi, phi_prime;
 *   
 *       unsigned int i = 0;
 *       for (; i < max_iter; ++i)
 *         {
 * @endcode
 * 
 * a simple bisection is faster than cubic interpolation
 * 
 * @code
 *           alpha = 0.5 * (alpha_low + alpha_high);
 * @endcode
 * 
 * alpha = line_search_interpolation_cubic(alpha_low, phi_low, phi_low_prime,
 * alpha_high, phi_high, phi_high_prime);
 * 
 * @code
 *           current_phi_phi_prime = calculate_phi_and_phi_prime(alpha, BFGS_p_vector, solution_delta);
 *           phi = current_phi_phi_prime.first;
 *           phi_prime = current_phi_phi_prime.second;
 *   
 *           if (   (phi > phi_0 + c1 * alpha * phi_0_prime)
 *               || (phi > phi_low) )
 *             {
 *               alpha_high = alpha;
 *               phi_high = phi;
 *               phi_high_prime = phi_prime;
 *             }
 *           else
 *             {
 *               if (std::fabs(phi_prime) <= c2 * std::fabs(phi_0_prime))
 *                 {
 * @endcode
 * 
 * if (alpha < 1.0e-3)
 * alpha = 1.0e-3;
 * 
 * @code
 *                   return alpha;
 *                 }
 *   
 *               if (phi_prime * (alpha_high - alpha_low) >= 0.0)
 *                 {
 *                   alpha_high = alpha_low;
 *                   phi_high_prime = phi_low_prime;
 *                   phi_high = phi_low;
 *                 }
 *   
 *               alpha_low = alpha;
 *               phi_low_prime = phi_prime;
 *               phi_low = phi;
 *             }
 *         }
 *   
 *       if (alpha < 1.0e-3)
 *         alpha = 1.0;
 *   
 * @endcode
 * 
 * avoid unused variable warnings from compiler
 * 
 * @code
 *       (void)phi_high;
 *       (void)phi_high_prime;
 *       return alpha;
 *     }
 *   
 *     template <int dim>
 *     double PhaseFieldMonolithicSolve<dim>::
 *       line_search_interpolation_cubic(const double alpha_0, const double phi_0, const double phi_0_prime,
 *                                         const double alpha_1, const double phi_1, const double phi_1_prime)
 *     {
 *       const double d1 = phi_0_prime + phi_1_prime - 3.0 * (phi_0 - phi_1) / (alpha_0 - alpha_1);
 *   
 *       const double temp = d1 * d1 - phi_0_prime * phi_1_prime;
 *   
 *       if (temp < 0.0)
 *         return 0.5 * (alpha_0 + alpha_1);
 *   
 *       int sign;
 *       if (alpha_1 > alpha_0)
 *         sign = 1;
 *       else
 *         sign = -1;
 *   
 *       const double d2 = sign * std::sqrt(temp);
 *   
 *       const double alpha = alpha_1 - (alpha_1 - alpha_0)
 *                          * (phi_1_prime + d2 - d1) / (phi_1_prime - phi_0_prime + 2*d2);
 *   
 *       if (    (alpha_1 > alpha_0)
 *            && (alpha > alpha_1 || alpha < alpha_0))
 *         return 0.5 * (alpha_0 + alpha_1);
 *   
 *       if (    (alpha_0 > alpha_1)
 *            && (alpha > alpha_0 || alpha < alpha_1))
 *         return 0.5 * (alpha_0 + alpha_1);
 *   
 *       return alpha;
 *     }
 *   
 *     template <int dim>
 *     std::pair<double, double> PhaseFieldMonolithicSolve<dim>::
 *       calculate_phi_and_phi_prime(const double alpha,
 *                                   const BlockVector<double> & BFGS_p_vector,
 *                                   const BlockVector<double> & solution_delta)
 *     {
 * @endcode
 * 
 * the first component is phi(alpha), the second component is phi_prime(alpha),
 * 
 * @code
 *       std::pair<double, double> phi_values;
 *   
 *       BlockVector<double> solution_delta_trial(solution_delta);
 *       solution_delta_trial.add(alpha, BFGS_p_vector);
 *   
 *       update_qph_incremental(solution_delta_trial, m_solution, false);
 *   
 *       BlockVector<double> system_rhs(m_dofs_per_block);
 *       assemble_system_rhs_BFGS_parallel(m_solution, system_rhs);
 * @endcode
 * 
 * m_constraints.condense(system_rhs);
 * 
 
 * 
 * 
 * @code
 *       phi_values.first = calculate_energy_functional();
 *       phi_values.second = system_rhs * BFGS_p_vector;
 *       return phi_values;
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::LBFGS_B0(BlockVector<double> & LBFGS_r_vector,
 *                                                   BlockVector<double> & LBFGS_q_vector)
 *     {
 *       m_timer.enter_subsection("Solve B0");
 *   
 *       assemble_system_B0(m_solution);
 *   
 *       if (m_parameters.m_type_linear_solver == "Direct")
 *         {
 *           SparseDirectUMFPACK A_direct;
 *           A_direct.initialize(m_tangent_matrix);
 *           A_direct.vmult(LBFGS_r_vector,
 *                          LBFGS_q_vector);
 *         }
 *       else if (m_parameters.m_type_linear_solver == "CG")
 *         {
 *   /*
 *           SolverControl            solver_control(1e6, 1e-9);
 *           SolverCG<BlockVector<double>> cg(solver_control);
 *   
 *           PreconditionJacobi<BlockSparseMatrix<double>> preconditioner;
 *           preconditioner.initialize(m_tangent_matrix, 1.0);
 *   
 *           cg.solve(m_tangent_matrix,
 *                    LBFGS_r_vector,
 *                    LBFGS_q_vector,
 *                    preconditioner);
 *   */
 *           SolverControl            solver_control_uu(1e6, 1e-9);
 *           SolverCG<Vector<double>> cg_uu(solver_control_uu);
 *   
 *           PreconditionJacobi<SparseMatrix<double>> preconditioner_uu;
 *           preconditioner_uu.initialize(m_tangent_matrix.block(m_u_dof, m_u_dof), 1.0);
 *           cg_uu.solve(m_tangent_matrix.block(m_u_dof, m_u_dof),
 *                       LBFGS_r_vector.block(m_u_dof),
 *                       LBFGS_q_vector.block(m_u_dof),
 *                       preconditioner_uu);
 *   
 *           SolverControl            solver_control_dd(1e6, 1e-15);
 *           SolverCG<Vector<double>> cg_dd(solver_control_dd);
 *   
 *           PreconditionJacobi<SparseMatrix<double>> preconditioner_dd;
 *           preconditioner_dd.initialize(m_tangent_matrix.block(m_d_dof, m_d_dof), 1.0);
 *           cg_dd.solve(m_tangent_matrix.block(m_d_dof, m_d_dof),
 *                       LBFGS_r_vector.block(m_d_dof),
 *                       LBFGS_q_vector.block(m_d_dof),
 *                       preconditioner_dd);
 *         }
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("Selected linear solver not implemented!"));
 *         }
 *   
 *       m_timer.leave_subsection();
 *     }
 *   
 *     template <int dim>
 *     std::vector<double>
 *       PhaseFieldMonolithicSolve<dim>::solve_linear_system(BlockVector<double> & newton_update)
 *     {
 *       m_timer.enter_subsection("Solve coupled linear system");
 *   
 *       if (m_parameters.m_output_iteration_history)
 *         m_logfile << " SLV " << std::flush;
 *   
 *       std::vector<double> linear_solver_parameters(3);
 *   /*
 *       {
 *         SolverControl            solver_control(1e6, 1e-9);
 *         SolverCG<Vector<double>> cg(solver_control);
 *         cg.connect_condition_number_slot(
 *             [&] (double condition_number)
 *             {
 *               linear_solver_parameters[0] = condition_number;
 * @endcode
 * 
 * m_logfile << "   Estimated condition number = "<< condition_number << std::endl;
 * 
 * @code
 *             },
 *             false);
 *   
 *         PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *         preconditioner.initialize(m_system_matrix_displacement, 1.2);
 *   
 *         cg.solve(m_system_matrix_displacement,
 *                  newton_update,
 *                  m_system_rhs_displacement,
 *                  preconditioner);
 *   
 * @endcode
 * 
 * m_logfile << "   " << solver_control.last_step()
 * << " CG iterations needed to obtain convergence." << std::endl;
 * 
 * @code
 *         linear_solver_parameters[1] = solver_control.last_step();
 *         linear_solver_parameters[2] = solver_control.last_value();
 *       }
 *   */
 *       SparseDirectUMFPACK A_direct;
 *       A_direct.initialize(m_tangent_matrix);
 *       A_direct.vmult(newton_update,
 *                      m_system_rhs);
 *   
 *       m_constraints.distribute(newton_update);
 *   
 *       m_timer.leave_subsection();
 *       return linear_solver_parameters;
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::print_conv_header_newton()
 *     {
 *       static const unsigned int l_width = 135;
 *       m_logfile << '\t' << '\t';
 *       for (unsigned int i = 0; i < l_width; ++i)
 *         m_logfile << '_';
 *       m_logfile << std::endl;
 *   
 *       m_logfile << "                  SOLVER STEP (Newton)      "
 *                 << " |        Cond No.   Lin_Iter   Lin_Res    Res_Norm  "
 *                 << " Res_u      Res_d      Inc_Norm  "
 *                 << " Inc_u      Inc_d" << std::endl;
 *   
 *       m_logfile << '\t' << '\t';
 *       for (unsigned int i = 0; i < l_width; ++i)
 *         m_logfile << '_';
 *       m_logfile << std::endl;
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::print_conv_header_BFGS()
 *     {
 *       static const unsigned int l_width = 125;
 *       m_logfile << '\t' << '\t';
 *       for (unsigned int i = 0; i < l_width; ++i)
 *         m_logfile << '_';
 *       m_logfile << std::endl;
 *   
 *       m_logfile << "                  SOLVER STEP (BFGS)   "
 *                 << " |    Line Search alpha    Energy    Res_Norm    "
 *                 << " Res_u      Res_d    Inc_Norm    "
 *                 << " Inc_u      Inc_d" << std::endl;
 *   
 *       m_logfile << '\t' << '\t';
 *       for (unsigned int i = 0; i < l_width; ++i)
 *         m_logfile << '_';
 *       m_logfile << std::endl;
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::print_conv_header_LBFGS()
 *     {
 *       static const unsigned int l_width = 120;
 *       m_logfile << '\t' << '\t';
 *       for (unsigned int i = 0; i < l_width; ++i)
 *         m_logfile << '_';
 *       m_logfile << std::endl;
 *   
 *       m_logfile << "                  SOLVER STEP (LBFGS)  "
 *                 << " |  LS-alpha     Energy      Res_Norm    "
 *                 << " Res_u      Res_d    Inc_Norm   "
 *                 << " Inc_u      Inc_d" << std::endl;
 *   
 *       m_logfile << '\t' << '\t';
 *       for (unsigned int i = 0; i < l_width; ++i)
 *         m_logfile << '_';
 *       m_logfile << std::endl;
 *     }
 *   
 *     template <int dim>
 *     bool PhaseFieldMonolithicSolve<dim>::
 *       solve_nonlinear_timestep_newton(BlockVector<double> & solution_delta)
 *     {
 *       BlockVector<double> newton_update(m_dofs_per_block);
 *   
 *       m_error_residual.reset();
 *       m_error_residual_0.reset();
 *       m_error_residual_norm.reset();
 *       m_error_update.reset();
 *       m_error_update_0.reset();
 *       m_error_update_norm.reset();
 *   
 *       if (m_parameters.m_output_iteration_history)
 *         print_conv_header_newton();
 *   
 *       unsigned int newton_iteration = 0;
 *       for (; newton_iteration < m_parameters.m_max_iterations_NR; ++newton_iteration)
 *         {
 *           if (m_parameters.m_output_iteration_history)
 *             m_logfile << '\t' << '\t' << std::setw(2) << newton_iteration << ' '
 *                       << std::flush;
 *   
 *           make_constraints(newton_iteration);
 *           assemble_system_newton(m_solution);
 *   
 *           get_error_residual(m_error_residual);
 *           if (newton_iteration == 0)
 *             m_error_residual_0 = m_error_residual;
 *   
 *           m_error_residual_norm = m_error_residual;
 *           m_error_residual_norm.normalize(m_error_residual_0);
 *   
 *           if (newton_iteration > 0 && m_error_update_norm.m_u <= m_parameters.m_tol_u_incr
 *                                    && m_error_residual_norm.m_u <= m_parameters.m_tol_u_residual
 *                                    && m_error_update_norm.m_d <= m_parameters.m_tol_d_incr
 *                                    && m_error_residual_norm.m_d <= m_parameters.m_tol_d_residual)
 *             {
 *               if (m_parameters.m_output_iteration_history)
 *                 {
 *                   m_logfile << " CONVERGED!";
 *                   m_logfile << " | " << std::fixed << std::setprecision(3) << std::setw(7)
 *                             << std::scientific
 *                         << "  " << "  ----   "
 *                         << "  " << "  ----   "
 *                         << "  " << "  ----   "
 *                         << "  " << m_error_residual_norm.m_norm
 *                         << "  " << m_error_residual_norm.m_u
 *                         << "  " << m_error_residual_norm.m_d
 *                         << "  " << m_error_update_norm.m_norm
 *                         << "  " << m_error_update_norm.m_u
 *                         << "  " << m_error_update_norm.m_d
 *                         << "  " << std::endl;
 *   
 *                   m_logfile << '\t' << '\t';
 *                   for (unsigned int i = 0; i < 135; ++i)
 *                     m_logfile << '_';
 *                   m_logfile << std::endl;
 *                 }
 *   
 *               m_logfile << "\t\tConvergence is reached after "
 *                         << newton_iteration << " Newton iterations."<< std::endl;
 *   
 *               m_logfile << "\t\tResidual information of convergence:" << std::endl;
 *   
 *               m_logfile << "\t\t\tRelative residual of disp. equation: "
 *                         << m_error_residual_norm.m_u << std::endl;
 *   
 *               m_logfile << "\t\t\tAbsolute residual of disp. equation: "
 *                         << m_error_residual_norm.m_u * m_error_residual_0.m_u << std::endl;
 *   
 *               m_logfile << "\t\t\tRelative residual of phasefield equation: "
 *                         << m_error_residual_norm.m_d << std::endl;
 *   
 *               m_logfile << "\t\t\tAbsolute residual of phasefield equation: "
 *                         << m_error_residual_norm.m_d * m_error_residual_0.m_d << std::endl;
 *   
 *               m_logfile << "\t\t\tRelative increment of disp.: "
 *                         << m_error_update_norm.m_u << std::endl;
 *   
 *               m_logfile << "\t\t\tAbsolute increment of disp.: "
 *                         << m_error_update_norm.m_u * m_error_update_0.m_u << std::endl;
 *   
 *               m_logfile << "\t\t\tRelative increment of phasefield: "
 *                         << m_error_update_norm.m_d << std::endl;
 *   
 *               m_logfile << "\t\t\tAbsolute increment of phasefield: "
 *                         << m_error_update_norm.m_d * m_error_update_0.m_d << std::endl;
 *   
 * @endcode
 * 
 * break;
 * 
 * @code
 *               return true;
 *             }
 *   
 *           std::vector<double> linear_solver_parameters(3);
 *   
 *           linear_solver_parameters = solve_linear_system(newton_update);
 *   
 *           get_error_update(newton_update, m_error_update);
 *           if (newton_iteration == 0)
 *             m_error_update_0 = m_error_update;
 *   
 *           m_error_update_norm = m_error_update;
 *           m_error_update_norm.normalize(m_error_update_0);
 *   
 *           solution_delta += newton_update;
 *           update_qph_incremental(solution_delta, m_solution, true);
 *   
 *           if (m_parameters.m_output_iteration_history)
 *             {
 *               m_logfile << " | " << std::fixed << std::setprecision(3) << std::setw(7)
 *                         << std::scientific
 *                         << "  " << linear_solver_parameters[0]
 *                         << "  " << linear_solver_parameters[1]
 *                         << "  " << linear_solver_parameters[2]
 *                         << "  " << m_error_residual_norm.m_norm
 *                         << "  " << m_error_residual_norm.m_u
 *                         << "  " << m_error_residual_norm.m_d
 *                         << "  " << m_error_update_norm.m_norm
 *                         << "  " << m_error_update_norm.m_u
 *                         << "  " << m_error_update_norm.m_d
 *                         << "  " << std::endl;
 *             }
 *         }
 *   
 * @endcode
 * 
 * AssertThrow(newton_iteration < m_parameters.m_max_iterations_NR,
 * ExcMessage("No convergence in Newton-Raphson nonlinear solver!"));
 * 
 * @code
 *       return false;
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::
 *     solve_nonlinear_timestep_BFGS(BlockVector<double> & solution_delta)
 *     {
 *       AssertThrow(false,
 *                   ExcMessage("BFGS requires too much memory. Please use L-BFGS!"));
 *   
 *       BlockVector<double> BFGS_update(m_dofs_per_block);
 *   
 *       m_error_residual.reset();
 *       m_error_residual_0.reset();
 *       m_error_residual_norm.reset();
 *       m_error_update.reset();
 *       m_error_update_0.reset();
 *       m_error_update_norm.reset();
 *   
 *       print_conv_header_BFGS();
 *   
 *       unsigned int BFGS_iteration = 0;
 *   
 * @endcode
 * 
 * Initial guess B_0, which is a full matrix and takes a lot of memory
 * 
 * @code
 *       FullMatrix<double> BFGS_matrix = IdentityMatrix(m_dof_handler.n_dofs());
 *       Vector<double> BFGS_r_vector(m_dof_handler.n_dofs());
 *       Vector<double> BFGS_p_vector(m_dof_handler.n_dofs());
 *       Vector<double> BFGS_y_vector(m_dof_handler.n_dofs());
 *       Vector<double> BFGS_temp_vector(m_dof_handler.n_dofs());
 *   
 *       double line_search_parameter, rho;
 *   
 * @endcode
 * 
 * Most likely, we will not be able to create a second full matrix since
 * we will run out of memory on a laptop workstation
 * 
 * @code
 *       FullMatrix<double> temp_matrix_1(m_dof_handler.n_dofs());
 *       FullMatrix<double> temp_matrix_2(m_dof_handler.n_dofs());
 *   
 *       for (; BFGS_iteration < m_parameters.m_max_iterations_BFGS; ++BFGS_iteration)
 *         {
 *           m_logfile << '\t' << '\t' << std::setw(2) << BFGS_iteration << ' '
 *                     << std::flush;
 *   
 *           make_constraints(BFGS_iteration);
 *   
 * @endcode
 * 
 * At the first step, we simply distribute the inhomogeneous part of
 * the constraints
 * 
 * @code
 *           if (BFGS_iteration == 0)
 *             {
 *               m_constraints.distribute(BFGS_update);
 *               solution_delta += BFGS_update;
 *               m_logfile << " --- " << std::flush;
 *               m_logfile << " --- " << std::flush;
 *               update_qph_incremental(solution_delta, m_solution, false);
 *               m_logfile << " ---  |" << std::flush;
 *               m_logfile << std::endl;
 *               continue;
 *             }
 *           else if (BFGS_iteration == 1)
 *             {
 * @endcode
 * 
 * Calculate the residual vector r. NOTICE that in the context of
 * BFGS, this r is the gradient of the energy functional (objective function),
 * NOT the negative gradient of the energy functional
 * 
 * @code
 *               assemble_system_rhs_BFGS(m_solution, m_system_rhs);
 *   
 * @endcode
 * 
 * We cannot simply zero out the dofs that are constrained, since we might
 * have hanging node constraints. In this case, we need to modify the RHS
 * as C^T * b, which C contains entries of 0.5 (x_3 = 0.5*x_1 + 0.5*x_2)
 * for (unsigned int i = 0; i < m_dof_handler.n_dofs(); ++i)
 * if (m_constraints.is_constrained(i))
 * m_system_rhs(i) = 0.0;
 * 
 
 * 
 * if m_constraints has inhomogeneity, we cannot call m_constraints.condense(m_system_rhs),
 * since the m_system_matrix needs to be provided to modify the RHS properly. However, this
 * error will not be detected in the release mode and only will be detected on the debug mode
 * 
 * @code
 *               m_constraints.condense(m_system_rhs);
 *             }
 *   
 *           m_logfile << " --- " << std::flush;
 *           m_logfile << " --- " << std::flush;
 *           m_logfile << " --- " << std::flush;
 *   
 *           get_error_residual(m_error_residual);
 *           if (BFGS_iteration == 1)
 *             m_error_residual_0 = m_error_residual;
 *   
 *           m_error_residual_norm = m_error_residual;
 *           m_error_residual_norm.normalize(m_error_residual_0);
 *   
 *           if (BFGS_iteration > 1 && m_error_update_norm.m_u <= m_parameters.m_tol_u_incr
 *                                  && m_error_residual_norm.m_u <= m_parameters.m_tol_u_residual
 *                                  && m_error_update_norm.m_d <= m_parameters.m_tol_d_incr
 *                                  && m_error_residual_norm.m_d <= m_parameters.m_tol_d_residual)
 *             {
 *               m_logfile << " CONVERGED!";
 *               m_logfile << "| " << std::fixed << std::setprecision(3) << std::setw(7)
 *                         << std::scientific
 *                         << "  " << "  ----   "
 *                         << "  " << "  ----   "
 *                         << "  " << "  ----   "
 *                         << "  " << m_error_residual_norm.m_norm
 *                         << "  " << m_error_residual_norm.m_u
 *                         << "  " << m_error_residual_norm.m_d
 *                     << "  " << m_error_update_norm.m_norm
 *                     << "  " << m_error_update_norm.m_u
 *                     << "  " << m_error_update_norm.m_d
 *                     << "  " << std::endl;
 *   
 *               m_logfile << '\t' << '\t';
 *               for (unsigned int i = 0; i < 135; ++i)
 *                 m_logfile << '_';
 *               m_logfile << std::endl;
 *   
 *               break;
 *             }
 *   
 * @endcode
 * 
 * BFGS algorithm
 * 
 * @code
 *           BFGS_r_vector = m_system_rhs;
 *           BFGS_matrix.vmult(BFGS_p_vector, BFGS_r_vector);
 *           BFGS_p_vector *= -1.0;
 *           m_constraints.distribute(BFGS_p_vector);
 *   
 * @endcode
 * 
 * We need a line search algorithm to decide line_search_parameter
 * 
 * @code
 *           const double phi_0 = calculate_energy_functional();
 *           const double phi_0_prime = BFGS_r_vector * BFGS_p_vector;
 *   
 *           BlockVector<double> BFGS_p_vector_block(m_dofs_per_block);
 *           BFGS_p_vector_block = BFGS_p_vector;
 *           line_search_parameter = line_search_stepsize_strong_wolfe(phi_0,
 *                                                                     phi_0_prime,
 *                                                                     BFGS_p_vector_block,
 *                                                                     solution_delta);
 *   
 *           BFGS_p_vector *= line_search_parameter;
 *           BFGS_update = BFGS_p_vector;
 *   
 *           get_error_update(BFGS_update, m_error_update);
 *           if (BFGS_iteration == 1)
 *             m_error_update_0 = m_error_update;
 *   
 *           m_error_update_norm = m_error_update;
 *           m_error_update_norm.normalize(m_error_update_0);
 *   
 *           solution_delta += BFGS_update;
 *           update_qph_incremental(solution_delta, m_solution, false);
 *   
 *           BFGS_y_vector = m_system_rhs;
 *           BFGS_y_vector *= -1.0;
 *           assemble_system_rhs_BFGS(m_solution, m_system_rhs);
 *           m_constraints.condense(m_system_rhs);
 *           BFGS_temp_vector = m_system_rhs;
 *           BFGS_y_vector += BFGS_temp_vector;
 *   
 * @endcode
 * 
 * rho should be positive with the proper line search
 * 
 * @code
 *           rho = BFGS_y_vector * BFGS_p_vector;
 *           rho = 1.0/rho;
 *   
 *           if (rho < 0)
 *             m_logfile << "Rho is negative!" << std::endl;
 *   
 * @endcode
 * 
 * In the first step, we scale the identity matrix as
 * the BFGS matrix
 * 
 * @code
 *           if (BFGS_iteration == 1)
 *             {
 *               double scale_parameter = (BFGS_y_vector * BFGS_p_vector) / (BFGS_y_vector.norm_sqr());
 *               BFGS_matrix *= scale_parameter;
 *             }
 *   
 *           temp_matrix_1.outer_product(BFGS_p_vector, BFGS_y_vector);
 *           temp_matrix_2 = IdentityMatrix(m_dof_handler.n_dofs());
 *           temp_matrix_2.add(-rho, temp_matrix_1);
 *   
 *           temp_matrix_2.mmult(temp_matrix_1, BFGS_matrix);
 *           temp_matrix_1.mTmult(BFGS_matrix, temp_matrix_2);
 *   
 *           temp_matrix_1.outer_product(BFGS_p_vector, BFGS_p_vector);
 *   
 *           BFGS_matrix.add(rho, temp_matrix_1);
 *   
 *           const double energy_functional = calculate_energy_functional();
 *   
 *           m_logfile << " | " << std::fixed << std::setprecision(3) << std::setw(7)
 *                     << std::scientific
 *                     << "  " << line_search_parameter
 *                     << "       " << energy_functional
 *                     << "  " << m_error_residual_norm.m_norm
 *                     << "  " << m_error_residual_norm.m_u
 *                     << "  " << m_error_residual_norm.m_d
 *                     << "  " << m_error_update_norm.m_norm
 *                     << "  " << m_error_update_norm.m_u
 *                     << "  " << m_error_update_norm.m_d
 *                     << "  " << std::endl;
 *         }
 *   
 *       AssertThrow(BFGS_iteration < m_parameters.m_max_iterations_BFGS,
 *                   ExcMessage("No convergence in BFGS nonlinear solver!"));
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::
 *     solve_nonlinear_timestep_LBFGS(BlockVector<double> & solution_delta,
 *                                    BlockVector<double> & LBFGS_update_refine)
 *     {
 *       BlockVector<double> LBFGS_update(m_dofs_per_block);
 *   
 *       LBFGS_update = 0.0;
 *   
 *       m_error_residual.reset();
 *       m_error_residual_0.reset();
 *       m_error_residual_norm.reset();
 *       m_error_update.reset();
 *       m_error_update_0.reset();
 *       m_error_update_norm.reset();
 *   
 *       if (m_parameters.m_output_iteration_history)
 *         print_conv_header_LBFGS();
 *   
 *       unsigned int LBFGS_iteration = 0;
 *   
 *       BlockVector<double> LBFGS_r_vector(m_dofs_per_block);
 *       BlockVector<double> LBFGS_y_vector(m_dofs_per_block);
 *       BlockVector<double> LBFGS_q_vector(m_dofs_per_block);
 *       BlockVector<double> LBFGS_s_vector(m_dofs_per_block);
 *       std::list<std::pair< std::pair<BlockVector<double>,
 *                                      BlockVector<double>>,
 *                            double>> LBFGS_vector_list;
 *   
 *       const unsigned int LBFGS_m = m_parameters.m_LBFGS_m;
 *       std::list<double> LBFGS_alpha_list;
 *   
 *       double line_search_parameter = 0.0;
 *       double LBFGS_beta = 0.0;
 *       double rho = 0.0;
 *   
 *       for (; LBFGS_iteration < m_parameters.m_max_iterations_BFGS; ++LBFGS_iteration)
 *         {
 *           if (m_parameters.m_output_iteration_history)
 *             m_logfile << '\t' << '\t' << std::setw(2) << LBFGS_iteration << ' '
 *                       << std::flush;
 *   
 *           make_constraints(LBFGS_iteration);
 *   
 * @endcode
 * 
 * At the first step, we simply distribute the inhomogeneous part of
 * the constraints
 * 
 * @code
 *           if (LBFGS_iteration == 0)
 *             {
 * @endcode
 * 
 * use the solution from the previous solve on the
 * refined mesh as initial guess
 * 
 * @code
 *               LBFGS_update = LBFGS_update_refine;
 *   
 *               m_constraints.distribute(LBFGS_update);
 *               solution_delta += LBFGS_update;
 *               if (m_parameters.m_output_iteration_history)
 *                 {
 *                   m_logfile << " --- " << std::flush;
 *                   m_logfile << " --- " << std::flush;
 *                 }
 *               update_qph_incremental(solution_delta, m_solution, false);
 *               if (m_parameters.m_output_iteration_history)
 *                 {
 *                   m_logfile << " ---  |" << std::flush;
 *                   m_logfile << std::endl;
 *                 }
 *               continue;
 *             }
 *           else if (LBFGS_iteration == 1)
 *             {
 * @endcode
 * 
 * Calculate the residual vector r. NOTICE that in the context of
 * BFGS, this r is the gradient of the energy functional (objective function),
 * NOT the negative gradient of the energy functional
 * 
 * @code
 *               assemble_system_rhs_BFGS_parallel(m_solution, m_system_rhs);
 *   
 * @endcode
 * 
 * We cannot simply zero out the dofs that are constrained, since we might
 * have hanging node constraints. In this case, we need to modify the RHS
 * as C^T * b, which C contains entries of 0.5 (x_3 = 0.5*x_1 + 0.5*x_2)
 * for (unsigned int i = 0; i < m_dof_handler.n_dofs(); ++i)
 * if (m_constraints.is_constrained(i))
 * m_system_rhs(i) = 0.0;
 * 
 
 * 
 * if m_constraints has inhomogeneity, we cannot call m_constraints.condense(m_system_rhs),
 * since the m_system_matrix needs to be provided to modify the RHS properly. However, this
 * error will not be detected in the release mode and only will be detected on the debug mode
 * if we use assemble_system_rhs_BFGS_parallel, then condense() is not necessary
 * m_constraints.condense(m_system_rhs);
 * 
 * @code
 *             }
 *           if (m_parameters.m_output_iteration_history)
 *             {
 *               m_logfile << " --- " << std::flush;
 *               m_logfile << " --- " << std::flush;
 *               m_logfile << " --- " << std::flush;
 *             }
 *   
 *           get_error_residual(m_error_residual);
 *           if (LBFGS_iteration == 1)
 *             m_error_residual_0 = m_error_residual;
 *   
 *           m_error_residual_norm = m_error_residual;
 * @endcode
 * 
 * For three-point bending problem and 3D problem, we use absolute residual
 * for convergence test
 * 
 * @code
 *           if (m_parameters.m_relative_residual)
 *             m_error_residual_norm.normalize(m_error_residual_0);
 *   
 *           if (LBFGS_iteration > 1 && m_error_update_norm.m_u <= m_parameters.m_tol_u_incr
 *                                   && m_error_residual_norm.m_u <= m_parameters.m_tol_u_residual
 *                                   && m_error_update_norm.m_d <= m_parameters.m_tol_d_incr
 *                                   && m_error_residual_norm.m_d <= m_parameters.m_tol_d_residual
 *                                   )
 *             {
 *               if (m_parameters.m_output_iteration_history)
 *                 {
 *                   m_logfile << " | ";
 *                   m_logfile << " CONVERGED! " << std::fixed << std::setprecision(3) << std::setw(7)
 *                             << std::scientific
 *                         << "    ----    "
 *                         << "  " << m_error_residual_norm.m_norm
 *                         << "  " << m_error_residual_norm.m_u
 *                         << "  " << m_error_residual_norm.m_d
 *                         << "  " << m_error_update_norm.m_norm
 *                         << "  " << m_error_update_norm.m_u
 *                         << "  " << m_error_update_norm.m_d
 *                         << "  " << std::endl;
 *   
 *                   m_logfile << '\t' << '\t';
 *                   for (unsigned int i = 0; i < 120; ++i)
 *                     m_logfile << '_';
 *                   m_logfile << std::endl;
 *                 }
 *   
 *               m_logfile << "\t\tConvergence is reached after "
 *                         << LBFGS_iteration << " L-BFGS iterations."<< std::endl;
 *   
 *               m_logfile << "\t\tResidual information of convergence:" << std::endl;
 *   
 *               if (m_parameters.m_relative_residual)
 *                 {
 *                   m_logfile << "\t\t\tRelative residual of disp. equation: "
 *                             << m_error_residual_norm.m_u << std::endl;
 *   
 *                   m_logfile << "\t\t\tAbsolute residual of disp. equation: "
 *                             << m_error_residual_norm.m_u * m_error_residual_0.m_u << std::endl;
 *   
 *                   m_logfile << "\t\t\tRelative residual of phasefield equation: "
 *                             << m_error_residual_norm.m_d << std::endl;
 *   
 *                   m_logfile << "\t\t\tAbsolute residual of phasefield equation: "
 *                             << m_error_residual_norm.m_d * m_error_residual_0.m_d << std::endl;
 *   
 *                   m_logfile << "\t\t\tRelative increment of disp.: "
 *                             << m_error_update_norm.m_u << std::endl;
 *   
 *                   m_logfile << "\t\t\tAbsolute increment of disp.: "
 *                             << m_error_update_norm.m_u * m_error_update_0.m_u << std::endl;
 *   
 *                   m_logfile << "\t\t\tRelative increment of phasefield: "
 *                             << m_error_update_norm.m_d << std::endl;
 *   
 *                   m_logfile << "\t\t\tAbsolute increment of phasefield: "
 *                             << m_error_update_norm.m_d * m_error_update_0.m_d << std::endl;
 *                 }
 *               else
 *                 {
 *                   m_logfile << "\t\t\tAbsolute residual of disp. equation: "
 *                             << m_error_residual_norm.m_u << std::endl;
 *   
 *                   m_logfile << "\t\t\tAbsolute residual of phasefield equation: "
 *                             << m_error_residual_norm.m_d << std::endl;
 *   
 *                   m_logfile << "\t\t\tAbsolute increment of disp.: "
 *                             << m_error_update_norm.m_u << std::endl;
 *   
 *                   m_logfile << "\t\t\tAbsolute increment of phasefield: "
 *                             << m_error_update_norm.m_d << std::endl;
 *                 }
 *   
 *               break;
 *             }
 *   
 * @endcode
 * 
 * LBFGS algorithm
 * 
 * @code
 *           LBFGS_q_vector = m_system_rhs;
 *   
 *           LBFGS_alpha_list.clear();
 *           for (auto itr = LBFGS_vector_list.begin(); itr != LBFGS_vector_list.end(); ++itr)
 *             {
 *               LBFGS_s_vector = (itr->first).first;
 *               LBFGS_y_vector = (itr->first).second;
 *               rho = itr->second;
 *   
 *               const double alpha = rho * (LBFGS_s_vector * LBFGS_q_vector);
 *               LBFGS_alpha_list.push_back(alpha);
 *   
 *               LBFGS_q_vector.add(-alpha, LBFGS_y_vector);
 *             }
 *   /*
 *           double scale_gamma = 0.0;
 *           if (LBFGS_iteration == 1)
 *             {
 *               scale_gamma = 1.0;
 *             }
 *           else
 *             {
 *               LBFGS_s_vector = LBFGS_vector_list.front().first.first;
 *               LBFGS_y_vector = LBFGS_vector_list.front().first.second;
 *               scale_gamma = (LBFGS_s_vector * LBFGS_y_vector)/(LBFGS_y_vector * LBFGS_y_vector);
 *             }
 *   
 *           LBFGS_q_vector *= scale_gamma;
 *           LBFGS_r_vector = LBFGS_q_vector;
 *   */
 *           LBFGS_B0(LBFGS_r_vector,
 *                    LBFGS_q_vector);
 *   
 *           for (auto itr = LBFGS_vector_list.rbegin(); itr != LBFGS_vector_list.rend(); ++itr)
 *             {
 *               LBFGS_s_vector = (itr->first).first;
 *               LBFGS_y_vector = (itr->first).second;
 *               rho = itr->second;
 *   
 *               LBFGS_beta = rho * (LBFGS_y_vector * LBFGS_r_vector);
 *   
 *               const double alpha = LBFGS_alpha_list.back();
 *               LBFGS_alpha_list.pop_back();
 *   
 *               LBFGS_r_vector.add(alpha - LBFGS_beta, LBFGS_s_vector);
 *             }
 *   
 *           LBFGS_r_vector *= -1.0; // this is the p_vector (search direction)
 *   
 *           m_constraints.distribute(LBFGS_r_vector);
 *   
 * @endcode
 * 
 * We need a line search algorithm to decide line_search_parameter
 * 
 * @code
 *           if(m_parameters.m_type_line_search == "StrongWolfe")
 *             {
 *               const double phi_0 = calculate_energy_functional();
 *               const double phi_0_prime = m_system_rhs * LBFGS_r_vector;
 *   
 *               line_search_parameter = line_search_stepsize_strong_wolfe(phi_0,
 *                                                                         phi_0_prime,
 *                                                                         LBFGS_r_vector,
 *                                                                         solution_delta);
 *             }
 *           else if(m_parameters.m_type_line_search == "GradientBased")
 *             {
 * @endcode
 * 
 * LBFGS_r_vector is the search direction
 * 
 * @code
 *               line_search_parameter = line_search_stepsize_gradient_based(LBFGS_r_vector,
 *                                                                           solution_delta);
 *             }
 *           else
 *             {
 *               Assert(false, ExcMessage("An unknown line search method is called!"));
 *             }
 *   
 *           LBFGS_r_vector *= line_search_parameter;
 *           LBFGS_update = LBFGS_r_vector;
 *   
 *           get_error_update(LBFGS_update, m_error_update);
 *           if (LBFGS_iteration == 1)
 *             m_error_update_0 = m_error_update;
 *   
 *           m_error_update_norm = m_error_update;
 * @endcode
 * 
 * For three-point bending problem and the sphere inclusion problem,
 * we use absolute residual for convergence test
 * 
 * @code
 *           if (m_parameters.m_relative_residual)
 *             m_error_update_norm.normalize(m_error_update_0);
 *   
 *           solution_delta += LBFGS_update;
 *           update_qph_incremental(solution_delta, m_solution, false);
 *   
 *           LBFGS_y_vector = m_system_rhs;
 *           LBFGS_y_vector *= -1.0;
 *           assemble_system_rhs_BFGS_parallel(m_solution, m_system_rhs);
 * @endcode
 * 
 * if we use assemble_system_rhs_BFGS_parallel, then condense() is not necessary
 * m_constraints.condense(m_system_rhs);
 * 
 * @code
 *           LBFGS_y_vector += m_system_rhs;
 *   
 *           LBFGS_s_vector = LBFGS_update;
 *   
 *           const double g_norm = m_system_rhs.l2_norm();
 *   
 *           const double yxs = LBFGS_y_vector * LBFGS_s_vector;
 *   
 *           const double sxs = LBFGS_s_vector * LBFGS_s_vector;
 *   
 *           if (yxs/sxs >= 1.0e-6 * g_norm)
 *             {
 *               if (LBFGS_iteration > LBFGS_m)
 *                 LBFGS_vector_list.pop_back();
 *   
 *               rho = 1.0 / yxs;
 *   
 *               LBFGS_vector_list.push_front(std::make_pair(std::make_pair(LBFGS_s_vector,
 *                                                                          LBFGS_y_vector),
 *                                                           rho));
 *             }
 *   
 *           if (m_parameters.m_output_iteration_history)
 *             {
 *               const double energy_functional = calculate_energy_functional();
 *   
 *               m_logfile << " | " << std::fixed << std::setprecision(3) << std::setw(1)
 *                         << std::scientific
 *                         << "" << line_search_parameter
 *                         << std::fixed << std::setprecision(6) << std::setw(1)
 *                                           << std::scientific
 *                         << "  " << energy_functional
 *                         << std::fixed << std::setprecision(3) << std::setw(1)
 *                                           << std::scientific
 *                         << "  " << m_error_residual_norm.m_norm
 *                         << "  " << m_error_residual_norm.m_u
 *                         << "  " << m_error_residual_norm.m_d
 *                         << "  " << m_error_update_norm.m_norm
 *                         << "  " << m_error_update_norm.m_u
 *                         << "  " << m_error_update_norm.m_d
 *                         << "  " << std::endl;
 *             }
 *         }
 *   
 *       AssertThrow(LBFGS_iteration < m_parameters.m_max_iterations_BFGS,
 *                   ExcMessage("No convergence in L-BFGS nonlinear solver!"));
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::output_results() const
 *     {
 *       m_timer.enter_subsection("Output results");
 *   
 *       DataOut<dim> data_out;
 *   
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *         data_component_interpretation(
 *           dim, DataComponentInterpretation::component_is_part_of_vector);
 *   
 *       data_component_interpretation.push_back(
 *         DataComponentInterpretation::component_is_scalar);
 *   
 *       std::vector<std::string> solution_name(dim, "displacement");
 *       solution_name.emplace_back("phasefield");
 *   
 *       data_out.attach_dof_handler(m_dof_handler);
 *       data_out.add_data_vector(m_solution,
 *                                solution_name,
 *                                DataOut<dim>::type_dof_data,
 *                                data_component_interpretation);
 *   
 *   
 *       Vector<double> cell_material_id(m_triangulation.n_active_cells());
 * @endcode
 * 
 * output material ID for each cell
 * 
 * @code
 *       for (const auto &cell : m_triangulation.active_cell_iterators())
 *         {
 *           cell_material_id(cell->active_cell_index()) = cell->material_id();
 *         }
 *       data_out.add_data_vector(cell_material_id, "materialID");
 *   
 * @endcode
 * 
 * Stress L2 projection
 * 
 * @code
 *       DoFHandler<dim> stresses_dof_handler_L2(m_triangulation);
 *       FE_Q<dim>     stresses_fe_L2(m_parameters.m_poly_degree); //FE_Q element is continuous
 *       stresses_dof_handler_L2.distribute_dofs(stresses_fe_L2);
 *       AffineConstraints<double> constraints;
 *       constraints.clear();
 *       DoFTools::make_hanging_node_constraints(stresses_dof_handler_L2, constraints);
 *       constraints.close();
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *             data_component_interpretation_stress(1,
 *                                                  DataComponentInterpretation::component_is_scalar);
 *   
 *       for (unsigned int i = 0; i < dim; ++i)
 *         for (unsigned int j = i; j < dim; ++j)
 *           {
 *             Vector<double> stress_field_L2;
 *             stress_field_L2.reinit(stresses_dof_handler_L2.n_dofs());
 *   
 *             MappingQ<dim> mapping(m_parameters.m_poly_degree + 1);
 *             VectorTools::project(mapping,
 *                                  stresses_dof_handler_L2,
 *                                  constraints,
 *                                  m_qf_cell,
 *                                  [&] (const typename DoFHandler<dim>::active_cell_iterator & cell,
 *                                       const unsigned int q) -> double
 *                                  {
 *                                    return m_quadrature_point_history.get_data(cell)[q]->get_cauchy_stress()[i][j];
 *                                  },
 *                                  stress_field_L2);
 *   
 *             std::string stress_name = "Cauchy_stress_" + std::to_string(i+1) + std::to_string(j+1)
 *                                     + "_L2";
 *   
 *             data_out.add_data_vector(stresses_dof_handler_L2,
 *                                      stress_field_L2,
 *                                      stress_name,
 *                                      data_component_interpretation_stress);
 *           }
 *   
 *       data_out.build_patches(m_parameters.m_poly_degree);
 *   
 *       std::ofstream output("Solution-" + std::to_string(dim) + "d-" +
 *                            Utilities::int_to_string(m_time.get_timestep(),4) + ".vtu");
 *   
 *       data_out.write_vtu(output);
 *       m_timer.leave_subsection();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::calculate_reaction_force(unsigned int face_ID)
 *     {
 *       m_timer.enter_subsection("Calculate reaction force");
 *   
 *       BlockVector<double>       system_rhs;
 *       system_rhs.reinit(m_dofs_per_block);
 *   
 *       Vector<double> cell_rhs(m_dofs_per_cell);
 *       std::vector<types::global_dof_index> local_dof_indices(m_dofs_per_cell);
 *   
 *       const double time_ramp = (m_time.current() / m_time.end());
 *       std::vector<Tensor<1, dim>> rhs_values(m_n_q_points);
 *       const UpdateFlags uf_cell(update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       const UpdateFlags uf_face(update_values | update_normal_vectors |
 *                                 update_JxW_values);
 *   
 *       FEValues<dim> fe_values(m_fe, m_qf_cell, uf_cell);
 *       FEFaceValues<dim> fe_face_values(m_fe, m_qf_face, uf_face);
 *   
 * @endcode
 * 
 * shape function values for displacement field
 * 
 * @code
 *       std::vector<std::vector<Tensor<1, dim>>>
 *         Nx(m_qf_cell.size(), std::vector<Tensor<1, dim>>(m_dofs_per_cell));
 *       std::vector<std::vector<Tensor<2, dim>>>
 *         grad_Nx(m_qf_cell.size(), std::vector<Tensor<2, dim>>(m_dofs_per_cell));
 *       std::vector<std::vector<SymmetricTensor<2, dim>>>
 *         symm_grad_Nx(m_qf_cell.size(), std::vector<SymmetricTensor<2, dim>>(m_dofs_per_cell));
 *   
 *       for (const auto &cell : m_dof_handler.active_cell_iterators())
 *         {
 * @endcode
 * 
 * if calculate_reaction_force() is defined as const, then
 * we also need to put a const in std::shared_ptr,
 * that is, std::shared_ptr<const PointHistory<dim>>
 * 
 * @code
 *           const std::vector<std::shared_ptr< PointHistory<dim>>> lqph =
 *             m_quadrature_point_history.get_data(cell);
 *           Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *           cell_rhs = 0.0;
 *           fe_values.reinit(cell);
 *           right_hand_side(fe_values.get_quadrature_points(),
 *                               rhs_values,
 *                               m_parameters.m_x_component*time_ramp,
 *                               m_parameters.m_y_component*time_ramp,
 *                               m_parameters.m_z_component*time_ramp);
 *   
 *           for (const unsigned int q_point : fe_values.quadrature_point_indices())
 *             {
 *               for (const unsigned int k : fe_values.dof_indices())
 *                 {
 *                   const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 *   
 *                   if (k_group == m_u_dof)
 *                     {
 *                           Nx[q_point][k] = fe_values[m_u_fe].value(k, q_point);
 *                           grad_Nx[q_point][k] = fe_values[m_u_fe].gradient(k, q_point);
 *                           symm_grad_Nx[q_point][k] = symmetrize(grad_Nx[q_point][k]);
 *                     }
 *                 }
 *             }
 *   
 *           for (const unsigned int q_point : fe_values.quadrature_point_indices())
 *             {
 *               const SymmetricTensor<2, dim> & cauchy_stress = lqph[q_point]->get_cauchy_stress();
 *   
 *               const std::vector<Tensor<1,dim>> & N = Nx[q_point];
 *               const std::vector<SymmetricTensor<2, dim>> & symm_grad_N = symm_grad_Nx[q_point];
 *               const double JxW = fe_values.JxW(q_point);
 *   
 *               for (const unsigned int i : fe_values.dof_indices())
 *                 {
 *                   const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *                   if (i_group == m_u_dof)
 *                     {
 *                       cell_rhs(i) -= (symm_grad_N[i] * cauchy_stress) * JxW;
 * @endcode
 * 
 * contributions from the body force to right-hand side
 * 
 * @code
 *                           cell_rhs(i) += N[i] * rhs_values[q_point] * JxW;
 *                     }
 *                 }
 *             }
 *   
 * @endcode
 * 
 * if there is surface pressure, this surface pressure always applied to the
 * reference configuration
 * 
 * @code
 *           const unsigned int face_pressure_id = 100;
 *           const double p0 = 0.0;
 *   
 *           for (const auto &face : cell->face_iterators())
 *             {
 *               if (face->at_boundary() && face->boundary_id() == face_pressure_id)
 *                 {
 *                   fe_face_values.reinit(cell, face);
 *   
 *                   for (const unsigned int f_q_point : fe_face_values.quadrature_point_indices())
 *                     {
 *                       const Tensor<1, dim> &N = fe_face_values.normal_vector(f_q_point);
 *   
 *                       const double         pressure  = p0 * time_ramp;
 *                       const Tensor<1, dim> traction  = pressure * N;
 *   
 *                       for (const unsigned int i : fe_values.dof_indices())
 *                         {
 *                           const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 *   
 *                           if (i_group == m_u_dof)
 *                             {
 *                               const unsigned int component_i = m_fe.system_to_component_index(i).first;
 *                               const double Ni = fe_face_values.shape_value(i, f_q_point);
 *                               const double JxW = fe_face_values.JxW(f_q_point);
 *                               cell_rhs(i) += (Ni * traction[component_i]) * JxW;
 *                             }
 *                         }
 *                     }
 *                 }
 *             }
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *           for (const unsigned int i : fe_values.dof_indices())
 *             system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *         } // for (const auto &cell : m_dof_handler.active_cell_iterators())
 *   
 * @endcode
 * 
 * The difference between the above assembled system_rhs and m_system_rhs
 * is that m_system_rhs is condensed by the m_constraints, which zero out
 * the rhs values associated with the constrained DOFs and modify the rhs
 * values associated with the unconstrained DOFs.
 * 
 
 * 
 * 
 * @code
 *       std::vector< types::global_dof_index > mapping;
 *       std::set<types::boundary_id> boundary_ids;
 *       boundary_ids.insert(face_ID);
 *       DoFTools::map_dof_to_boundary_indices(m_dof_handler,
 *                                             boundary_ids,
 *                                             mapping);
 *   
 *       std::vector<double> reaction_force(dim, 0.0);
 *   
 *       for (unsigned int i = 0; i < m_dofs_per_block[m_u_dof]; ++i)
 *         {
 *           if (mapping[i] != numbers::invalid_dof_index)
 *             {
 *               reaction_force[i % dim] += system_rhs.block(m_u_dof)(i);
 *             }
 *         }
 *   
 *       for (unsigned int i = 0; i < dim; i++)
 *         m_logfile << "\t\tReaction force in direction " << i << " on boundary ID " << face_ID
 *                   << " = "
 *                   << std::fixed << std::setprecision(3) << std::setw(1)
 *                   << std::scientific
 *                   << reaction_force[i] << std::endl;
 *   
 *       std::pair<double, std::vector<double>> time_force;
 *       time_force.first = m_time.current();
 *       time_force.second = reaction_force;
 *       m_history_reaction_force.push_back(time_force);
 *   
 *       m_timer.leave_subsection();
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::write_history_data()
 *     {
 *       m_logfile << "\t\tWrite history data ... \n"<<std::endl;
 *   
 *       std::ofstream myfile_reaction_force ("Reaction_force.hist");
 *       if (myfile_reaction_force.is_open())
 *       {
 *         myfile_reaction_force << 0.0 << "\t";
 *         if (dim == 2)
 *           myfile_reaction_force << 0.0 << "\t"
 *                  << 0.0 << std::endl;
 *         if (dim == 3)
 *           myfile_reaction_force << 0.0 << "\t"
 *                  << 0.0 << "\t"
 *                  << 0.0 << std::endl;
 *   
 *         for (const auto &time_force : m_history_reaction_force)
 *           {
 *             myfile_reaction_force << time_force.first << "\t";
 *             if (dim == 2)
 *               myfile_reaction_force << time_force.second[0] << "\t"
 *                      << time_force.second[1] << std::endl;
 *             if (dim == 3)
 *               myfile_reaction_force << time_force.second[0] << "\t"
 *                      << time_force.second[1] << "\t"
 *                      << time_force.second[2] << std::endl;
 *           }
 *         myfile_reaction_force.close();
 *       }
 *       else
 *         m_logfile << "Unable to open file";
 *   
 *       std::ofstream myfile_energy ("Energy.hist");
 *       if (myfile_energy.is_open())
 *       {
 *         myfile_energy << std::fixed << std::setprecision(10) << std::scientific
 *                       << 0.0 << "\t"
 *                       << 0.0 << "\t"
 *                       << 0.0 << "\t"
 *                       << 0.0 << std::endl;
 *   
 *         for (const auto &time_energy : m_history_energy)
 *           {
 *             myfile_energy << std::fixed << std::setprecision(10) << std::scientific
 *                           << time_energy.first     << "\t"
 *                           << time_energy.second[0] << "\t"
 *                           << time_energy.second[1] << "\t"
 *                           << time_energy.second[2] << std::endl;
 *           }
 *         myfile_energy.close();
 *       }
 *       else
 *         m_logfile << "Unable to open file";
 *     }
 *   
 *     template <int dim>
 *     double PhaseFieldMonolithicSolve<dim>::calculate_energy_functional() const
 *     {
 *       double energy_functional = 0.0;
 *   
 *       FEValues<dim> fe_values(m_fe, m_qf_cell, update_JxW_values);
 *   
 *       for (const auto &cell : m_dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *   
 *           const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *             m_quadrature_point_history.get_data(cell);
 *           Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *           for (unsigned int q_point = 0; q_point < m_n_q_points; ++q_point)
 *             {
 *               const double JxW = fe_values.JxW(q_point);
 *               energy_functional += lqph[q_point]->get_total_strain_energy() * JxW;
 *               energy_functional += lqph[q_point]->get_crack_energy_dissipation() * JxW;
 *             }
 *         }
 *   
 *       return energy_functional;
 *     }
 *   
 *     template <int dim>
 *     std::pair<double, double>
 *       PhaseFieldMonolithicSolve<dim>::calculate_total_strain_energy_and_crack_energy_dissipation() const
 *     {
 *       double total_strain_energy = 0.0;
 *       double crack_energy_dissipation = 0.0;
 *   
 *       FEValues<dim> fe_values(m_fe, m_qf_cell, update_JxW_values);
 *   
 *       for (const auto &cell : m_dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *   
 *           const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *             m_quadrature_point_history.get_data(cell);
 *           Assert(lqph.size() == m_n_q_points, ExcInternalError());
 *   
 *           for (unsigned int q_point = 0; q_point < m_n_q_points; ++q_point)
 *             {
 *               const double JxW = fe_values.JxW(q_point);
 *               total_strain_energy += lqph[q_point]->get_total_strain_energy() * JxW;
 *               crack_energy_dissipation += lqph[q_point]->get_crack_energy_dissipation() * JxW;
 *             }
 *         }
 *   
 *       return std::make_pair(total_strain_energy, crack_energy_dissipation);
 *     }
 *   
 *   
 *     template <int dim>
 *     bool PhaseFieldMonolithicSolve<dim>::local_refine_and_solution_transfer(BlockVector<double> & solution_delta,
 *                                                                             BlockVector<double> & LBFGS_update_refine)
 *     {
 * @endcode
 * 
 * This is the solution at (n+1) obtained from the old (coarse) mesh
 * 
 * @code
 *       BlockVector<double> solution_next_step(m_dofs_per_block);
 *       solution_next_step = m_solution + solution_delta;
 *       bool mesh_is_same = true;
 *       bool cell_refine_flag = true;
 *   
 *       unsigned int material_id;
 *       double length_scale;
 *       double cell_length;
 *       while(cell_refine_flag)
 *         {
 *           cell_refine_flag = false;
 *   
 *           std::vector<types::global_dof_index> local_dof_indices(m_fe.dofs_per_cell);
 *           for (const auto &cell : m_dof_handler.active_cell_iterators())
 *             {
 *               cell->get_dof_indices(local_dof_indices);
 *   
 *               for (unsigned int i = 0; i< m_fe.dofs_per_cell; ++i)
 *                 {
 *                   const unsigned int comp_i = m_fe.system_to_component_index(i).first;
 *                   if (comp_i == m_d_component) //phasefield component
 *                     {
 *                       if (  solution_next_step(local_dof_indices[i])
 *                           > m_parameters.m_phasefield_refine_threshold )
 *                         {
 *                           material_id = cell->material_id();
 *                           length_scale = m_material_data[material_id][2];
 *                           if (dim == 2)
 *                             cell_length = std::sqrt(cell->measure());
 *                           else
 *                             cell_length = std::cbrt(cell->measure());
 *                           if (  cell_length
 *                               > length_scale * m_parameters.m_allowed_max_h_l_ratio )
 *                             {
 *                               if (cell->level() < m_parameters.m_max_allowed_refinement_level)
 *                                 {
 *                                   cell->set_refine_flag();
 *                                   break;
 *                                 }
 *                             }
 *                         }
 *                     }
 *                 }
 *             }
 *   
 *           for (const auto &cell : m_dof_handler.active_cell_iterators())
 *             {
 *               if (cell->refine_flag_set())
 *                 {
 *                   cell_refine_flag = true;
 *                   break;
 *                 }
 *             }
 *   
 * @endcode
 * 
 * if any cell is refined, we need to project the solution
 * to the newly refined mesh
 * 
 * @code
 *           if (cell_refine_flag)
 *             {
 *               mesh_is_same = false;
 *   
 *               std::vector<BlockVector<double> > old_solutions(2);
 *               old_solutions[0] = solution_next_step;
 *               old_solutions[1] = m_solution;
 *   
 *               m_triangulation.prepare_coarsening_and_refinement();
 *               SolutionTransfer<dim, BlockVector<double>> solution_transfer(m_dof_handler);
 *               solution_transfer.prepare_for_coarsening_and_refinement(old_solutions);
 *               m_triangulation.execute_coarsening_and_refinement();
 *   
 *               setup_system();
 *   
 *               std::vector<BlockVector<double>> tmp_solutions(2);
 *               tmp_solutions[0].reinit(m_dofs_per_block);
 *               tmp_solutions[1].reinit(m_dofs_per_block);
 *               solution_transfer.interpolate(tmp_solutions);
 *               solution_next_step = tmp_solutions[0];
 *               m_solution = tmp_solutions[1];
 *   
 * @endcode
 * 
 * make sure the projected solutions still satisfy
 * hanging node constraints
 * 
 * @code
 *               m_constraints.distribute(solution_next_step);
 *               m_constraints.distribute(m_solution);
 *             } // if (cell_refine_flag)
 *         } // while(cell_refine_flag)
 *   
 * @endcode
 * 
 * calculate field variables for newly refined cells
 * 
 * @code
 *       if (!mesh_is_same)
 *         {
 *           BlockVector<double> temp_solution_delta(m_dofs_per_block);
 *           BlockVector<double> temp_previous_solution(m_dofs_per_block);
 *           temp_solution_delta = 0.0;
 *           temp_previous_solution = 0.0;
 *           update_qph_incremental(temp_solution_delta, temp_previous_solution, false);
 *           update_history_field_step();
 *   
 * @endcode
 * 
 * initial guess for the resolve on the refined mesh
 * 
 * @code
 *           LBFGS_update_refine = solution_next_step - m_solution;
 *         }
 *   
 *       return mesh_is_same;
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::print_parameter_information()
 *     {
 *       m_logfile << "Scenario number = " << m_parameters.m_scenario << std::endl;
 *       m_logfile << "Log file = " << m_parameters.m_logfile_name << std::endl;
 *       m_logfile << "Write iteration history to log file? = " << std::boolalpha
 *                 << m_parameters.m_output_iteration_history << std::endl;
 *       m_logfile << "Nonlinear solver type = " << m_parameters.m_type_nonlinear_solver << std::endl;
 *       m_logfile << "Line search type = " << m_parameters.m_type_line_search << std::endl;
 *       m_logfile << "Linear solver type = " << m_parameters.m_type_linear_solver << std::endl;
 *       m_logfile << "Mesh refinement strategy = " << m_parameters.m_refinement_strategy << std::endl;
 *       m_logfile << "L-BFGS_m = " << m_parameters.m_LBFGS_m << std::endl;
 *       m_logfile << "Global refinement times = " << m_parameters.m_global_refine_times << std::endl;
 *       m_logfile << "Local prerefinement times = " <<m_parameters. m_local_prerefine_times << std::endl;
 *       m_logfile << "Maximum adaptive refinement times allowed in each step = "
 *                 << m_parameters.m_max_adaptive_refine_times << std::endl;
 *       m_logfile << "Maximum allowed cell refinement level = "
 *                     << m_parameters.m_max_allowed_refinement_level << std::endl;
 *       m_logfile << "Phasefield-based refinement threshold value = "
 *                 << m_parameters.m_phasefield_refine_threshold << std::endl;
 *       m_logfile << "Allowed maximum h/l ratio = " << m_parameters.m_allowed_max_h_l_ratio << std::endl;
 *       m_logfile << "total number of material types = " << m_parameters.m_total_material_regions << std::endl;
 *       m_logfile << "material data file name = " << m_parameters.m_material_file_name << std::endl;
 *       if (m_parameters.m_reaction_force_face_id >= 0)
 *         m_logfile << "Calculate reaction forces on Face ID = " << m_parameters.m_reaction_force_face_id << std::endl;
 *       else
 *         m_logfile << "No need to calculate reaction forces." << std::endl;
 *   
 *       if (m_parameters.m_relative_residual)
 *         m_logfile << "Relative residual for convergence." << std::endl;
 *       else
 *         m_logfile << "Absolute residual for convergence." << std::endl;
 *   
 *       m_logfile << "Body force = (" << m_parameters.m_x_component << ", "
 *                                     << m_parameters.m_y_component << ", "
 *                                     << m_parameters.m_z_component << ") (N/m^3)"
 *                                     << std::endl;
 *   
 *       m_logfile << "End time = " << m_parameters.m_end_time << std::endl;
 *       m_logfile << "Time data file name = " << m_parameters.m_time_file_name << std::endl;
 *     }
 *   
 *     template <int dim>
 *     void PhaseFieldMonolithicSolve<dim>::run()
 *     {
 *       print_parameter_information();
 *   
 *       read_material_data(m_parameters.m_material_file_name,
 *                              m_parameters.m_total_material_regions);
 *   
 *       std::vector<std::array<double, 4>> time_table;
 *   
 *       read_time_data(m_parameters.m_time_file_name, time_table);
 *   
 *       make_grid();
 *       setup_system();
 *       output_results();
 *   
 *       m_time.increment(time_table);
 *   
 *       while(m_time.current() < m_time.end() + m_time.get_delta_t()*1.0e-6)
 *         {
 *           m_logfile << std::endl
 *                     << "Timestep " << m_time.get_timestep() << " @ " << m_time.current()
 *                     << 's' << std::endl;
 *   
 *           bool mesh_is_same = false;
 *   
 * @endcode
 * 
 * initial guess for the resolve on the refined mesh
 * 
 * @code
 *           BlockVector<double> LBFGS_update_refine(m_dofs_per_block);
 *           LBFGS_update_refine = 0.0;
 *   
 * @endcode
 * 
 * local adaptive mesh refinement loop
 * 
 * @code
 *           unsigned int adp_refine_iteration = 0;
 *           for (; adp_refine_iteration < m_parameters.m_max_adaptive_refine_times + 1; ++adp_refine_iteration)
 *             {
 *               if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *                 m_logfile << "\tAdaptive refinement-"<< adp_refine_iteration << ": " << std::endl;
 *   
 *               BlockVector<double> solution_delta(m_dofs_per_block);
 *               solution_delta = 0.0;
 *   
 *               if (m_parameters.m_type_nonlinear_solver == "Newton")
 *                 {
 *                   bool newton_success = false;
 *                   newton_success = solve_nonlinear_timestep_newton(solution_delta);
 *                   AssertThrow(newton_success,
 *                               ExcMessage("No convergence in Newton-Raphson nonlinear solver!"));
 *                   /*
 * @endcode
 * 
 * if Newton-Raphson failed, use LBFGS solver
 * 
 * @code
 *                   if (!newton_success)
 *                     {
 *                       solution_delta = 0.0;
 *                       solve_nonlinear_timestep_LBFGS(solution_delta, LBFGS_update_refine);
 *                     }
 *                   */
 *                 }
 *               else if (m_parameters.m_type_nonlinear_solver == "BFGS")
 *                 solve_nonlinear_timestep_BFGS(solution_delta);
 *               else if (m_parameters.m_type_nonlinear_solver == "LBFGS")
 *                 solve_nonlinear_timestep_LBFGS(solution_delta, LBFGS_update_refine);
 *               else
 *                 AssertThrow(false, ExcMessage("Nonlinear solver type not implemented"));
 *   
 *               if (m_parameters.m_refinement_strategy == "adaptive-refine")
 *                 {
 *   
 *                   if (adp_refine_iteration == m_parameters.m_max_adaptive_refine_times)
 *                     {
 *                       m_solution += solution_delta;
 *                       break;
 *                     }
 *   
 *                   mesh_is_same = local_refine_and_solution_transfer(solution_delta,
 *                                                                     LBFGS_update_refine);
 *   
 *                   if (mesh_is_same)
 *                     {
 *                       m_solution += solution_delta;
 *                       break;
 *                     }
 *                 }
 *               else if (m_parameters.m_refinement_strategy == "pre-refine")
 *                 {
 *                   m_solution += solution_delta;
 *                   break;
 *                 }
 *               else
 *                 {
 *                   AssertThrow(false,
 *                               ExcMessage("Selected mesh refinement strategy not implemented!"));
 *                 }
 *             } // for (; adp_refine_iteration < m_parameters.m_max_adaptive_refine_times; ++adp_refine_iteration)
 *   
 * @endcode
 * 
 * AssertThrow(adp_refine_iteration < m_parameters.m_max_adaptive_refine_times,
 * ExcMessage("Number of local adaptive mesh refinement exceeds allowed maximum times!"));
 * 
 
 * 
 * 
 * @code
 *           update_history_field_step();
 * @endcode
 * 
 * output vtk files every 10 steps if there are too
 * many time steps
 * if (m_time.get_timestep() % 10 == 0)
 * 
 * @code
 *           output_results();
 *   
 *           double energy_functional_current = calculate_energy_functional();
 *           m_logfile << "\t\tEnergy functional (J) = " << std::fixed << std::setprecision(10) << std::scientific
 *                     << energy_functional_current << std::endl;
 *   
 *           std::pair<double, double> energy_pair = calculate_total_strain_energy_and_crack_energy_dissipation();
 *           m_logfile << "\t\tTotal strain energy (J) = " << std::fixed << std::setprecision(10) << std::scientific
 *                     << energy_pair.first << std::endl;
 *           m_logfile << "\t\tCrack energy dissipation (J) = " << std::fixed << std::setprecision(10) << std::scientific
 *                     << energy_pair.second << std::endl;
 *   
 *           std::pair<double, std::array<double, 3>> time_energy;
 *           time_energy.first = m_time.current();
 *           time_energy.second[0] = energy_pair.first;
 *           time_energy.second[1] = energy_pair.second;
 *           time_energy.second[2] = energy_pair.first + energy_pair.second;
 *           m_history_energy.push_back(time_energy);
 *   
 *           int face_ID = m_parameters.m_reaction_force_face_id;
 *           if (face_ID >= 0)
 *             calculate_reaction_force(face_ID);
 *   
 *           write_history_data();
 *   
 *           m_time.increment(time_table);
 *         } // while(m_time.current() < m_time.end() + m_time.get_delta_t()*1.0e-6)
 *     }
 *   } // namespace PhaseField
 *   
 *   
 *   int main(int argc, char* argv[])
 *   {
 *   
 *     using namespace dealii;
 *   
 *     if (argc != 2)
 *       AssertThrow(false,
 *                       ExcMessage("The number of arguments provided to the program has to be 2!"));
 *   
 *     const unsigned int dim = std::stoi(argv[1]);
 *     if (dim == 2 )
 *       {
 *         PhaseField::PhaseFieldMonolithicSolve<2> FEQ1Full("parameters.prm");
 *         FEQ1Full.run();
 *       }
 *     else if (dim == 3)
 *       {
 *         PhaseField::PhaseFieldMonolithicSolve<3> SphereInclusion3D("parameters.prm");
 *         SphereInclusion3D.run();
 *       }
 *     else
 *       {
 *         AssertThrow(false,
 *                     ExcMessage("Dimension has to be either 2 or 3"));
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
 
 
*/