Quasi_static_Finite_strain_Compressible_Elasticity.h

Go to the documentation of this file.

 
 * 
 * /*
 *  * Authors: Jean-Paul Pelteret, University of Erlangen-Nuremberg,
 *  *          Andrew McBride, University of Cape Town, 2015, 2017
 *  */
 * 
 * 
 * @endcode
 * 
 * We start by including all the necessary deal.II header files and some C++
 * related ones. They have been discussed in detail in previous tutorial
 * programs, so you need only refer to past tutorials for details.
 * 
 * @code
 * #include <deal.II/base/function.h>
 * #include <deal.II/base/parameter_handler.h>
 * #include <deal.II/base/point.h>
 * #include <deal.II/base/quadrature_lib.h>
 * #include <deal.II/base/symmetric_tensor.h>
 * #include <deal.II/base/tensor.h>
 * #include <deal.II/base/timer.h>
 * #include <deal.II/base/work_stream.h>
 * #include <deal.II/base/quadrature_point_data.h>
 * 
 * #include <deal.II/dofs/dof_renumbering.h>
 * #include <deal.II/dofs/dof_tools.h>
 * 
 * #include <deal.II/grid/grid_generator.h>
 * #include <deal.II/grid/grid_tools.h>
 * #include <deal.II/grid/grid_in.h>
 * #include <deal.II/grid/tria.h>
 * 
 * #include <deal.II/fe/fe_dgp_monomial.h>
 * #include <deal.II/fe/fe_q.h>
 * #include <deal.II/fe/fe_system.h>
 * #include <deal.II/fe/fe_tools.h>
 * #include <deal.II/fe/fe_values.h>
 * #include <deal.II/fe/mapping_q_eulerian.h>
 * #include <deal.II/fe/mapping_q.h>
 * 
 * #include <deal.II/lac/affine_constraints.h>
 * #include <deal.II/lac/block_sparse_matrix.h>
 * #include <deal.II/lac/block_vector.h>
 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
 * #include <deal.II/lac/full_matrix.h>
 * #include <deal.II/lac/precondition_selector.h>
 * #include <deal.II/lac/solver_cg.h>
 * #include <deal.II/lac/sparse_direct.h>
 * 
 * #include <deal.II/numerics/data_out.h>
 * #include <deal.II/numerics/vector_tools.h>
 * 
 * #include <deal.II/base/config.h>
 * #if DEAL_II_VERSION_MAJOR >= 9 && defined(DEAL_II_WITH_TRILINOS)
 * #include <deal.II/differentiation/ad.h>
 * #define ENABLE_SACADO_FORMULATION
 * #endif
 * 
 * @endcode
 * 
 * These must be included below the AD headers so that
 * their math functions are available for use in the
 * definition of tensors and kinematic quantities
 * 
 * @code
 * #include <deal.II/physics/elasticity/kinematics.h>
 * #include <deal.II/physics/elasticity/standard_tensors.h>
 * 
 * #include <iostream>
 * #include <fstream>
 * #include <memory>
 * 
 * 
 * @endcode
 * 
 * We then stick everything that relates to this tutorial program into a
 * namespace of its own, and import all the deal.II function and class names
 * into it:
 * 
 * @code
 * namespace Cook_Membrane
 * {
 *   using namespace dealii;
 * 
 * @endcode
 * 
 * 
 * <a name="Runtimeparameters"></a> 
 * <h3>Run-time parameters</h3>
 * 
 
 * 
 * There are several parameters that can be set in the code so we set up a
 * ParameterHandler object to read in the choices at run-time.
 * 
 * @code
 *   namespace Parameters
 *   {
 * @endcode
 * 
 * 
 * <a name="Assemblymethod"></a> 
 * <h4>Assembly method</h4>
 * 
 
 * 
 * Here we specify whether automatic differentiation is to be used to assemble
 * the linear system, and if so then what order of differentiation is to be
 * employed.
 * 
 * @code
 *     struct AssemblyMethod
 *     {
 *       unsigned int automatic_differentiation_order;
 * 
 *       static void
 *       declare_parameters(ParameterHandler &prm);
 * 
 *       void
 *       parse_parameters(ParameterHandler &prm);
 *     };
 * 
 * 
 *     void AssemblyMethod::declare_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Assembly method");
 *       {
 *         prm.declare_entry("Automatic differentiation order", "0",
 *                           Patterns::Integer(0,2),
 *                           "The automatic differentiation order to be used in the assembly of the linear system.\n"
 *                           "# Order = 0: Both the residual and linearisation are computed manually.\n"
 *                           "# Order = 1: The residual is computed manually but the linearisation is performed using AD.\n"
 *                           "# Order = 2: Both the residual and linearisation are computed using AD.");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 *     void AssemblyMethod::parse_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Assembly method");
 *       {
 *         automatic_differentiation_order = prm.get_integer("Automatic differentiation order");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 * @endcode
 * 
 * 
 * <a name="FiniteElementsystem"></a> 
 * <h4>Finite Element system</h4>
 * 
 
 * 
 * Here we specify the polynomial order used to approximate the solution.
 * The quadrature order should be adjusted accordingly.
 * 
 * @code
 *     struct FESystem
 *     {
 *       unsigned int poly_degree;
 *       unsigned int 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", "2",
 *                           Patterns::Integer(0),
 *                           "Displacement system polynomial order");
 * 
 *         prm.declare_entry("Quadrature order", "3",
 *                           Patterns::Integer(0),
 *                           "Gauss quadrature order");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 *     void FESystem::parse_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Finite element system");
 *       {
 *         poly_degree = prm.get_integer("Polynomial degree");
 *         quad_order = prm.get_integer("Quadrature order");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 * @endcode
 * 
 * 
 * <a name="Geometry"></a> 
 * <h4>Geometry</h4>
 * 
 
 * 
 * Make adjustments to the problem geometry and its discretisation.
 * 
 * @code
 *     struct Geometry
 *     {
 *       unsigned int elements_per_edge;
 *       double       scale;
 * 
 *       static void
 *       declare_parameters(ParameterHandler &prm);
 * 
 *       void
 *       parse_parameters(ParameterHandler &prm);
 *     };
 * 
 *     void Geometry::declare_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Geometry");
 *       {
 *         prm.declare_entry("Elements per edge", "32",
 *                           Patterns::Integer(0),
 *                           "Number of elements per long edge of the beam");
 * 
 *         prm.declare_entry("Grid scale", "1e-3",
 *                           Patterns::Double(0.0),
 *                           "Global grid scaling factor");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 *     void Geometry::parse_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Geometry");
 *       {
 *         elements_per_edge = prm.get_integer("Elements per edge");
 *         scale = prm.get_double("Grid scale");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 * @endcode
 * 
 * 
 * <a name="Materials"></a> 
 * <h4>Materials</h4>
 * 
 
 * 
 * We also need the shear modulus @f$ \mu @f$ and Poisson ration @f$ \nu @f$ for the
 * neo-Hookean material.
 * 
 * @code
 *     struct Materials
 *     {
 *       double nu;
 *       double mu;
 * 
 *       static void
 *       declare_parameters(ParameterHandler &prm);
 * 
 *       void
 *       parse_parameters(ParameterHandler &prm);
 *     };
 * 
 *     void Materials::declare_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Material properties");
 *       {
 *         prm.declare_entry("Poisson's ratio", "0.3",
 *                           Patterns::Double(-1.0,0.5),
 *                           "Poisson's ratio");
 * 
 *         prm.declare_entry("Shear modulus", "0.4225e6",
 *                           Patterns::Double(),
 *                           "Shear modulus");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 *     void Materials::parse_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Material properties");
 *       {
 *         nu = prm.get_double("Poisson's ratio");
 *         mu = prm.get_double("Shear modulus");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 * @endcode
 * 
 * 
 * <a name="Linearsolver"></a> 
 * <h4>Linear solver</h4>
 * 
 
 * 
 * Next, we choose both solver and preconditioner settings.  The use of an
 * effective preconditioner is critical to ensure convergence when a large
 * nonlinear motion occurs within a Newton increment.
 * 
 * @code
 *     struct LinearSolver
 *     {
 *       std::string type_lin;
 *       double      tol_lin;
 *       double      max_iterations_lin;
 *       std::string preconditioner_type;
 *       double      preconditioner_relaxation;
 * 
 *       static void
 *       declare_parameters(ParameterHandler &prm);
 * 
 *       void
 *       parse_parameters(ParameterHandler &prm);
 *     };
 * 
 *     void LinearSolver::declare_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Linear solver");
 *       {
 *         prm.declare_entry("Solver type", "CG",
 *                           Patterns::Selection("CG|Direct"),
 *                           "Type of solver used to solve the linear system");
 * 
 *         prm.declare_entry("Residual", "1e-6",
 *                           Patterns::Double(0.0),
 *                           "Linear solver residual (scaled by residual norm)");
 * 
 *         prm.declare_entry("Max iteration multiplier", "1",
 *                           Patterns::Double(0.0),
 *                           "Linear solver iterations (multiples of the system matrix size)");
 * 
 *         prm.declare_entry("Preconditioner type", "ssor",
 *                           Patterns::Selection("jacobi|ssor"),
 *                           "Type of preconditioner");
 * 
 *         prm.declare_entry("Preconditioner relaxation", "0.65",
 *                           Patterns::Double(0.0),
 *                           "Preconditioner relaxation value");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 *     void LinearSolver::parse_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Linear solver");
 *       {
 *         type_lin = prm.get("Solver type");
 *         tol_lin = prm.get_double("Residual");
 *         max_iterations_lin = prm.get_double("Max iteration multiplier");
 *         preconditioner_type = prm.get("Preconditioner type");
 *         preconditioner_relaxation = prm.get_double("Preconditioner relaxation");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 * @endcode
 * 
 * 
 * <a name="Nonlinearsolver"></a> 
 * <h4>Nonlinear solver</h4>
 * 
 
 * 
 * A Newton-Raphson scheme is used to solve the nonlinear system of governing
 * equations.  We now define the tolerances and the maximum number of
 * iterations for the Newton-Raphson nonlinear solver.
 * 
 * @code
 *     struct NonlinearSolver
 *     {
 *       unsigned int max_iterations_NR;
 *       double       tol_f;
 *       double       tol_u;
 * 
 *       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("Tolerance force", "1.0e-9",
 *                           Patterns::Double(0.0),
 *                           "Force residual tolerance");
 * 
 *         prm.declare_entry("Tolerance displacement", "1.0e-6",
 *                           Patterns::Double(0.0),
 *                           "Displacement error tolerance");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 *     void NonlinearSolver::parse_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Nonlinear solver");
 *       {
 *         max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
 *         tol_f = prm.get_double("Tolerance force");
 *         tol_u = prm.get_double("Tolerance displacement");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 * @endcode
 * 
 * 
 * <a name="Time"></a> 
 * <h4>Time</h4>
 * 
 
 * 
 * Set the timestep size @f$ \varDelta t @f$ and the simulation end-time.
 * 
 * @code
 *     struct Time
 *     {
 *       double delta_t;
 *       double end_time;
 * 
 *       static void
 *       declare_parameters(ParameterHandler &prm);
 * 
 *       void
 *       parse_parameters(ParameterHandler &prm);
 *     };
 * 
 *     void Time::declare_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Time");
 *       {
 *         prm.declare_entry("End time", "1",
 *                           Patterns::Double(),
 *                           "End time");
 * 
 *         prm.declare_entry("Time step size", "0.1",
 *                           Patterns::Double(),
 *                           "Time step size");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 *     void Time::parse_parameters(ParameterHandler &prm)
 *     {
 *       prm.enter_subsection("Time");
 *       {
 *         end_time = prm.get_double("End time");
 *         delta_t = prm.get_double("Time step size");
 *       }
 *       prm.leave_subsection();
 *     }
 * 
 * @endcode
 * 
 * 
 * <a name="Allparameters"></a> 
 * <h4>All parameters</h4>
 * 
 
 * 
 * Finally we consolidate all of the above structures into a single container
 * that holds all of our run-time selections.
 * 
 * @code
 *     struct AllParameters :
 *       public AssemblyMethod,
 *       public FESystem,
 *       public Geometry,
 *       public Materials,
 *       public LinearSolver,
 *       public NonlinearSolver,
 *       public Time
 * 
 *     {
 *       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)
 *     {
 *       AssemblyMethod::declare_parameters(prm);
 *       FESystem::declare_parameters(prm);
 *       Geometry::declare_parameters(prm);
 *       Materials::declare_parameters(prm);
 *       LinearSolver::declare_parameters(prm);
 *       NonlinearSolver::declare_parameters(prm);
 *       Time::declare_parameters(prm);
 *     }
 * 
 *     void AllParameters::parse_parameters(ParameterHandler &prm)
 *     {
 *       AssemblyMethod::parse_parameters(prm);
 *       FESystem::parse_parameters(prm);
 *       Geometry::parse_parameters(prm);
 *       Materials::parse_parameters(prm);
 *       LinearSolver::parse_parameters(prm);
 *       NonlinearSolver::parse_parameters(prm);
 *       Time::parse_parameters(prm);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Timeclass"></a> 
 * <h3>Time class</h3>
 * 
 
 * 
 * A simple class to store time data. Its functioning is transparent so no
 * discussion is necessary. For simplicity we assume a constant time step
 * size.
 * 
 * @code
 *   class Time
 *   {
 *   public:
 *     Time (const double time_end,
 *           const double delta_t)
 *       :
 *       timestep(0),
 *       time_current(0.0),
 *       time_end(time_end),
 *       delta_t(delta_t)
 *     {}
 * 
 *     virtual ~Time()
 *     {}
 * 
 *     double current() const
 *     {
 *       return time_current;
 *     }
 *     double end() const
 *     {
 *       return time_end;
 *     }
 *     double get_delta_t() const
 *     {
 *       return delta_t;
 *     }
 *     unsigned int get_timestep() const
 *     {
 *       return timestep;
 *     }
 *     void increment()
 *     {
 *       time_current += delta_t;
 *       ++timestep;
 *     }
 * 
 *   private:
 *     unsigned int timestep;
 *     double       time_current;
 *     const double time_end;
 *     const double delta_t;
 *   };
 * 
 * @endcode
 * 
 * 
 * <a name="CompressibleneoHookeanmaterialwithinaonefieldformulation"></a> 
 * <h3>Compressible neo-Hookean material within a one-field formulation</h3>
 * 
 
 * 
 * As discussed in the literature and @ref step_44 "step-44", Neo-Hookean materials are a type
 * of hyperelastic materials.  The entire domain is assumed to be composed of a
 * compressible neo-Hookean material.  This class defines the behaviour of
 * this material within a one-field formulation.  Compressible neo-Hookean
 * materials can be described by a strain-energy function (SEF) @f$ \Psi =
 * \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(J)
 * @f$.
 * 
 
 * 
 * The isochoric response is given by @f$
 * \Psi_{\text{iso}}(\overline{\mathbf{b}}) = c_{1} [\overline{I}_{1} - 3] @f$
 * where @f$ c_{1} = \frac{\mu}{2} @f$ and @f$\overline{I}_{1}@f$ is the first
 * invariant of the left- or right-isochoric Cauchy-Green deformation tensors.
 * That is @f$\overline{I}_1 :=\textrm{tr}(\overline{\mathbf{b}})@f$.  In this
 * example the SEF that governs the volumetric response is defined as @f$
 * \Psi_{\text{vol}}(J) = \kappa \frac{1}{4} [ J^2 - 1
 * - 2\textrm{ln}\; J ]@f$,  where @f$\kappa:= \lambda + 2/3 \mu@f$ is
 * the <a href="http://en.wikipedia.org/wiki/Bulk_modulus">bulk modulus</a>
 * and @f$\lambda@f$ is <a
 * href="http://en.wikipedia.org/wiki/Lam%C3%A9_parameters">Lame's first
 * parameter</a>.
 * 
 
 * 
 * The following class will be used to characterize the material we work with,
 * and provides a central point that one would need to modify if one were to
 * implement a different material model. For it to work, we will store one
 * object of this type per quadrature point, and in each of these objects
 * store the current state (characterized by the values or measures  of the
 * displacement field) so that we can compute the elastic coefficients
 * linearized around the current state.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   class Material_Compressible_Neo_Hook_One_Field
 *   {
 *   public:
 *     Material_Compressible_Neo_Hook_One_Field(const double mu,
 *                                              const double nu)
 *       :
 *       kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
 *       c_1(mu / 2.0)
 *     {
 *       Assert(kappa > 0, ExcInternalError());
 *     }
 * 
 *     ~Material_Compressible_Neo_Hook_One_Field()
 *     {}
 * 
 * @endcode
 * 
 * The first function is the total energy
 * @f$\Psi = \Psi_{\textrm{iso}} + \Psi_{\textrm{vol}}@f$.
 * 
 * @code
 *     NumberType
 *     get_Psi(const NumberType                        &det_F,
 *             const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       return get_Psi_vol(det_F) + get_Psi_iso(b_bar);
 *     }
 * 
 * @endcode
 * 
 * The second function determines the Kirchhoff stress @f$\boldsymbol{\tau}
 * = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}@f$
 * 
 * @code
 *     SymmetricTensor<2,dim,NumberType>
 *     get_tau(const NumberType                        &det_F,
 *             const SymmetricTensor<2,dim,NumberType> &b_bar)
 *     {
 * @endcode
 * 
 * See Holzapfel p231 eq6.98 onwards
 * 
 * @code
 *       return get_tau_vol(det_F) + get_tau_iso(b_bar);
 *     }
 * 
 * @endcode
 * 
 * The fourth-order elasticity tensor in the spatial setting
 * @f$\mathfrak{c}@f$ is calculated from the SEF @f$\Psi@f$ as @f$ J
 * \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}@f$
 * where @f$ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial
 * \mathbf{C} \partial \mathbf{C}}@f$
 * 
 * @code
 *     SymmetricTensor<4,dim,NumberType>
 *     get_Jc(const NumberType                        &det_F,
 *            const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       return get_Jc_vol(det_F) + get_Jc_iso(b_bar);
 *     }
 * 
 *   private:
 * @endcode
 * 
 * Define constitutive model parameters @f$\kappa@f$ (bulk modulus) and the
 * neo-Hookean model parameter @f$c_1@f$:
 * 
 * @code
 *     const double kappa;
 *     const double c_1;
 * 
 * @endcode
 * 
 * Value of the volumetric free energy
 * 
 * @code
 *     NumberType
 *     get_Psi_vol(const NumberType &det_F) const
 *     {
 *       return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
 *     }
 * 
 * @endcode
 * 
 * Value of the isochoric free energy
 * 
 * @code
 *     NumberType
 *     get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       return c_1 * (trace(b_bar) - dim);
 *     }
 * 
 * @endcode
 * 
 * Derivative of the volumetric free energy with respect to
 * @f$J@f$ return @f$\frac{\partial
 * \Psi_{\text{vol}}(J)}{\partial J}@f$
 * 
 * @code
 *     NumberType
 *     get_dPsi_vol_dJ(const NumberType &det_F) const
 *     {
 *       return (kappa / 2.0) * (det_F - 1.0 / det_F);
 *     }
 * 
 * @endcode
 * 
 * The following functions are used internally in determining the result
 * of some of the public functions above. The first one determines the
 * volumetric Kirchhoff stress @f$\boldsymbol{\tau}_{\textrm{vol}}@f$.
 * Note the difference in its definition when compared to @ref step_44 "step-44".
 * 
 * @code
 *     SymmetricTensor<2,dim,NumberType>
 *     get_tau_vol(const NumberType &det_F) const
 *     {
 *       return NumberType(get_dPsi_vol_dJ(det_F) * det_F) * Physics::Elasticity::StandardTensors<dim>::I;
 *     }
 * 
 * @endcode
 * 
 * Next, determine the isochoric Kirchhoff stress
 * @f$\boldsymbol{\tau}_{\textrm{iso}} =
 * \mathcal{P}:\overline{\boldsymbol{\tau}}@f$:
 * 
 * @code
 *     SymmetricTensor<2,dim,NumberType>
 *     get_tau_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       return Physics::Elasticity::StandardTensors<dim>::dev_P * get_tau_bar(b_bar);
 *     }
 * 
 * @endcode
 * 
 * Then, determine the fictitious Kirchhoff stress
 * @f$\overline{\boldsymbol{\tau}}@f$:
 * 
 * @code
 *     SymmetricTensor<2,dim,NumberType>
 *     get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       return 2.0 * c_1 * b_bar;
 *     }
 * 
 * @endcode
 * 
 * Second derivative of the volumetric free energy wrt @f$J@f$. We
 * need the following computation explicitly in the tangent so we make it
 * public.  We calculate @f$\frac{\partial^2
 * \Psi_{\textrm{vol}}(J)}{\partial J \partial
 * J}@f$
 * 
 * @code
 *     NumberType
 *     get_d2Psi_vol_dJ2(const NumberType &det_F) const
 *     {
 *       return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
 *     }
 * 
 * @endcode
 * 
 * Calculate the volumetric part of the tangent @f$J
 * \mathfrak{c}_\textrm{vol}@f$. Again, note the difference in its
 * definition when compared to @ref step_44 "step-44". The extra terms result from two
 * quantities in @f$\boldsymbol{\tau}_{\textrm{vol}}@f$ being dependent on
 * @f$\boldsymbol{F}@f$.
 * 
 * @code
 *     SymmetricTensor<4,dim,NumberType>
 *     get_Jc_vol(const NumberType &det_F) const
 *     {
 * @endcode
 * 
 * See Holzapfel p265
 * 
 * @code
 *       return det_F
 *              * ( (get_dPsi_vol_dJ(det_F) + det_F * get_d2Psi_vol_dJ2(det_F))*Physics::Elasticity::StandardTensors<dim>::IxI
 *                  - (2.0 * get_dPsi_vol_dJ(det_F))*Physics::Elasticity::StandardTensors<dim>::S );
 *     }
 * 
 * @endcode
 * 
 * Calculate the isochoric part of the tangent @f$J
 * \mathfrak{c}_\textrm{iso}@f$:
 * 
 * @code
 *     SymmetricTensor<4,dim,NumberType>
 *     get_Jc_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       const SymmetricTensor<2, dim> tau_bar = get_tau_bar(b_bar);
 *       const SymmetricTensor<2, dim> tau_iso = get_tau_iso(b_bar);
 *       const SymmetricTensor<4, dim> tau_iso_x_I
 *         = outer_product(tau_iso,
 *                         Physics::Elasticity::StandardTensors<dim>::I);
 *       const SymmetricTensor<4, dim> I_x_tau_iso
 *         = outer_product(Physics::Elasticity::StandardTensors<dim>::I,
 *                         tau_iso);
 *       const SymmetricTensor<4, dim> c_bar = get_c_bar();
 * 
 *       return (2.0 / dim) * trace(tau_bar)
 *              * Physics::Elasticity::StandardTensors<dim>::dev_P
 *              - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso)
 *              + Physics::Elasticity::StandardTensors<dim>::dev_P * c_bar
 *              * Physics::Elasticity::StandardTensors<dim>::dev_P;
 *     }
 * 
 * @endcode
 * 
 * Calculate the fictitious elasticity tensor @f$\overline{\mathfrak{c}}@f$.
 * For the material model chosen this is simply zero:
 * 
 * @code
 *     SymmetricTensor<4,dim,double>
 *     get_c_bar() const
 *     {
 *       return SymmetricTensor<4, dim>();
 *     }
 *   };
 * 
 * @endcode
 * 
 * 
 * <a name="Quadraturepointhistory"></a> 
 * <h3>Quadrature point history</h3>
 * 
 
 * 
 * As seen in @ref step_18 "step-18", the <code> PointHistory </code> class offers a method
 * for storing data at the quadrature points.  Here each quadrature point
 * holds a pointer to a material description.  Thus, different material models
 * can be used in different regions of the domain.  Among other data, we
 * choose to store the Kirchhoff stress @f$\boldsymbol{\tau}@f$ and the tangent
 * @f$J\mathfrak{c}@f$ for the quadrature points.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   class PointHistory
 *   {
 *   public:
 *     PointHistory()
 *     {}
 * 
 *     virtual ~PointHistory()
 *     {}
 * 
 * @endcode
 * 
 * The first function is used to create a material object and to
 * initialize all tensors correctly: The second one updates the stored
 * values and stresses based on the current deformation measure
 * @f$\textrm{Grad}\mathbf{u}_{\textrm{n}}@f$.
 * 
 * @code
 *     void setup_lqp (const Parameters::AllParameters &parameters)
 *     {
 *       material.reset(new Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu,
 *                      parameters.nu));
 *     }
 * 
 * @endcode
 * 
 * We offer an interface to retrieve certain data.
 * This is the strain energy:
 * 
 * @code
 *     NumberType
 *     get_Psi(const NumberType                        &det_F,
 *             const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       return material->get_Psi(det_F,b_bar);
 *     }
 * 
 * @endcode
 * 
 * Here are the kinetic variables. These are used in the material and
 * global tangent matrix and residual assembly operations:
 * First is the Kirchhoff stress:
 * 
 * @code
 *     SymmetricTensor<2,dim,NumberType>
 *     get_tau(const NumberType                        &det_F,
 *             const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       return material->get_tau(det_F,b_bar);
 *     }
 * 
 * @endcode
 * 
 * And the tangent:
 * 
 * @code
 *     SymmetricTensor<4,dim,NumberType>
 *     get_Jc(const NumberType                        &det_F,
 *            const SymmetricTensor<2,dim,NumberType> &b_bar) const
 *     {
 *       return material->get_Jc(det_F,b_bar);
 *     }
 * 
 * @endcode
 * 
 * In terms of member functions, this class stores for the quadrature
 * point it represents a copy of a material type in case different
 * materials are used in different regions of the domain, as well as the
 * inverse of the deformation gradient...
 * 
 * @code
 *   private:
 *     std::shared_ptr< Material_Compressible_Neo_Hook_One_Field<dim,NumberType> > material;
 *   };
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Quasistaticcompressiblefinitestrainsolid"></a> 
 * <h3>Quasi-static compressible finite-strain solid</h3>
 * 
 
 * 
 * Forward declarations for classes that will
 * perform assembly of the linear system.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   struct Assembler_Base;
 *   template <int dim,typename NumberType>
 *   struct Assembler;
 * 
 * @endcode
 * 
 * The Solid class is the central class in that it represents the problem at
 * hand. It follows the usual scheme in that all it really has is a
 * constructor, destructor and a <code>run()</code> function that dispatches
 * all the work to private functions of this class:
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   class Solid
 *   {
 *   public:
 *     Solid(const Parameters::AllParameters &parameters);
 * 
 *     virtual
 *     ~Solid();
 * 
 *     void
 *     run();
 * 
 *   private:
 * 
 * @endcode
 * 
 * We start the collection of member functions with one that builds the
 * grid:
 * 
 * @code
 *     void
 *     make_grid();
 * 
 * @endcode
 * 
 * Set up the finite element system to be solved:
 * 
 * @code
 *     void
 *     system_setup();
 * 
 * @endcode
 * 
 * Several functions to assemble the system and right hand side matrices
 * using multithreading. Each of them comes as a wrapper function, one
 * that is executed to do the work in the WorkStream model on one cell,
 * and one that copies the work done on this one cell into the global
 * object that represents it:
 * 
 * @code
 *     void
 *     assemble_system(const BlockVector<double> &solution_delta);
 * 
 * @endcode
 * 
 * We use a separate data structure to perform the assembly. It needs access
 * to some low-level data, so we simply befriend the class instead of
 * creating a complex interface to provide access as necessary.
 * 
 * @code
 *     friend struct Assembler_Base<dim,NumberType>;
 *     friend struct Assembler<dim,NumberType>;
 * 
 * @endcode
 * 
 * Apply Dirichlet boundary conditions on the displacement field
 * 
 * @code
 *     void
 *     make_constraints(const int &it_nr);
 * 
 * @endcode
 * 
 * Create and update the quadrature points. Here, no data needs to be
 * copied into a global object, so the copy_local_to_global function is
 * empty:
 * 
 * @code
 *     void
 *     setup_qph();
 * 
 * @endcode
 * 
 * Solve for the displacement using a Newton-Raphson method. We break this
 * function into the nonlinear loop and the function that solves the
 * linearized Newton-Raphson step:
 * 
 * @code
 *     void
 *     solve_nonlinear_timestep(BlockVector<double> &solution_delta);
 * 
 *     std::pair<unsigned int, double>
 *     solve_linear_system(BlockVector<double> &newton_update);
 * 
 * @endcode
 * 
 * Solution retrieval as well as post-processing and writing data to file :
 * 
 * @code
 *     BlockVector<double>
 *     get_total_solution(const BlockVector<double> &solution_delta) const;
 * 
 *     void
 *     output_results() const;
 * 
 * @endcode
 * 
 * Finally, some member variables that describe the current state: A
 * collection of the parameters used to describe the problem setup...
 * 
 * @code
 *     const Parameters::AllParameters &parameters;
 * 
 * @endcode
 * 
 * ...the volume of the reference and current configurations...
 * 
 * @code
 *     double                           vol_reference;
 *     double                           vol_current;
 * 
 * @endcode
 * 
 * ...and description of the geometry on which the problem is solved:
 * 
 * @code
 *     Triangulation<dim>               triangulation;
 * 
 * @endcode
 * 
 * Also, keep track of the current time and the time spent evaluating
 * certain functions
 * 
 * @code
 *     Time                             time;
 *     TimerOutput                      timer;
 * 
 * @endcode
 * 
 * A storage object for quadrature point information. As opposed to
 * @ref step_18 "step-18", deal.II's native quadrature point data manager is employed here.
 * 
 * @code
 *     CellDataStorage<typename Triangulation<dim>::cell_iterator,
 *                     PointHistory<dim,NumberType> > quadrature_point_history;
 * 
 * @endcode
 * 
 * A description of the finite-element system including the displacement
 * polynomial degree, the degree-of-freedom handler, number of DoFs per
 * cell and the extractor objects used to retrieve information from the
 * solution vectors:
 * 
 * @code
 *     const unsigned int               degree;
 *     const FESystem<dim>              fe;
 *     DoFHandler<dim>                  dof_handler_ref;
 *     const unsigned int               dofs_per_cell;
 *     const FEValuesExtractors::Vector u_fe;
 * 
 * @endcode
 * 
 * Description of how the block-system is arranged. There is just 1 block,
 * that contains a vector DOF @f$\mathbf{u}@f$.
 * There are two reasons that we retain the block system in this problem.
 * The first is pure laziness to perform further modifications to the
 * code from which this work originated. The second is that a block system
 * would typically necessary when extending this code to multiphysics
 * problems.
 * 
 * @code
 *     static const unsigned int        n_blocks = 1;
 *     static const unsigned int        n_components = dim;
 *     static const unsigned int        first_u_component = 0;
 * 
 *     enum
 *     {
 *       u_dof = 0
 *     };
 * 
 *     std::vector<types::global_dof_index>  dofs_per_block;
 * 
 * @endcode
 * 
 * Rules for Gauss-quadrature on both the cell and faces. The number of
 * quadrature points on both cells and faces is recorded.
 * 
 * @code
 *     const QGauss<dim>                qf_cell;
 *     const QGauss<dim - 1>            qf_face;
 *     const unsigned int               n_q_points;
 *     const unsigned int               n_q_points_f;
 * 
 * @endcode
 * 
 * Objects that store the converged solution and right-hand side vectors,
 * as well as the tangent matrix. There is a AffineConstraints object used
 * to keep track of constraints.  We make use of a sparsity pattern
 * designed for a block system.
 * 
 * @code
 *     AffineConstraints<double>        constraints;
 *     BlockSparsityPattern             sparsity_pattern;
 *     BlockSparseMatrix<double>        tangent_matrix;
 *     BlockVector<double>              system_rhs;
 *     BlockVector<double>              solution_n;
 * 
 * @endcode
 * 
 * Then define a number of variables to store norms and update norms and
 * normalisation factors.
 * 
 * @code
 *     struct Errors
 *     {
 *       Errors()
 *         :
 *         norm(1.0), u(1.0)
 *       {}
 * 
 *       void reset()
 *       {
 *         norm = 1.0;
 *         u = 1.0;
 *       }
 *       void normalise(const Errors &rhs)
 *       {
 *         if (rhs.norm != 0.0)
 *           norm /= rhs.norm;
 *         if (rhs.u != 0.0)
 *           u /= rhs.u;
 *       }
 * 
 *       double norm, u;
 *     };
 * 
 *     Errors error_residual, error_residual_0, error_residual_norm, error_update,
 *            error_update_0, error_update_norm;
 * 
 * @endcode
 * 
 * Methods to calculate error measures
 * 
 * @code
 *     void
 *     get_error_residual(Errors &error_residual);
 * 
 *     void
 *     get_error_update(const BlockVector<double> &newton_update,
 *                      Errors &error_update);
 * 
 * @endcode
 * 
 * Print information to screen in a pleasing way...
 * 
 * @code
 *     static
 *     void
 *     print_conv_header();
 * 
 *     void
 *     print_conv_footer();
 * 
 *     void
 *     print_vertical_tip_displacement();
 *   };
 * 
 * @endcode
 * 
 * 
 * <a name="ImplementationofthecodeSolidcodeclass"></a> 
 * <h3>Implementation of the <code>Solid</code> class</h3>
 * 
 
 * 
 * 
 * <a name="Publicinterface"></a> 
 * <h4>Public interface</h4>
 * 
 
 * 
 * We initialise the Solid class using data extracted from the parameter file.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   Solid<dim,NumberType>::Solid(const Parameters::AllParameters &parameters)
 *     :
 *     parameters(parameters),
 *     vol_reference (0.0),
 *     vol_current (0.0),
 *     triangulation(Triangulation<dim>::maximum_smoothing),
 *     time(parameters.end_time, parameters.delta_t),
 *     timer(std::cout,
 *           TimerOutput::summary,
 *           TimerOutput::wall_times),
 *     degree(parameters.poly_degree),
 * @endcode
 * 
 * The Finite Element System is composed of dim continuous displacement
 * DOFs.
 * 
 * @code
 *     fe(FE_Q<dim>(parameters.poly_degree), dim), // displacement
 *     dof_handler_ref(triangulation),
 *     dofs_per_cell (fe.dofs_per_cell),
 *     u_fe(first_u_component),
 *     dofs_per_block(n_blocks),
 *     qf_cell(parameters.quad_order),
 *     qf_face(parameters.quad_order),
 *     n_q_points (qf_cell.size()),
 *     n_q_points_f (qf_face.size())
 *   {
 * 
 *   }
 * 
 * @endcode
 * 
 * The class destructor simply clears the data held by the DOFHandler
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   Solid<dim,NumberType>::~Solid()
 *   {
 *     dof_handler_ref.clear();
 *   }
 * 
 * 
 * @endcode
 * 
 * In solving the quasi-static problem, the time becomes a loading parameter,
 * i.e. we increasing the loading linearly with time, making the two concepts
 * interchangeable. We choose to increment time linearly using a constant time
 * step size.
 * 
 
 * 
 * We start the function with preprocessing, and then output the initial grid
 * before starting the simulation proper with the first time (and loading)
 * increment.
 * 
 
 * 
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::run()
 *   {
 *     make_grid();
 *     system_setup();
 *     output_results();
 *     time.increment();
 * 
 * @endcode
 * 
 * We then declare the incremental solution update @f$\varDelta
 * \mathbf{\Xi}:= \{\varDelta \mathbf{u}\}@f$ and start the loop over the
 * time domain.
 *     
 
 * 
 * At the beginning, we reset the solution update for this time step...
 * 
 * @code
 *     BlockVector<double> solution_delta(dofs_per_block);
 *     while (time.current() <= time.end())
 *       {
 *         solution_delta = 0.0;
 * 
 * @endcode
 * 
 * ...solve the current time step and update total solution vector
 * @f$\mathbf{\Xi}_{\textrm{n}} = \mathbf{\Xi}_{\textrm{n-1}} +
 * \varDelta \mathbf{\Xi}@f$...
 * 
 * @code
 *         solve_nonlinear_timestep(solution_delta);
 *         solution_n += solution_delta;
 * 
 * @endcode
 * 
 * ...and plot the results before moving on happily to the next time
 * step:
 * 
 * @code
 *         output_results();
 *         time.increment();
 *       }
 * 
 * @endcode
 * 
 * Lastly, we print the vertical tip displacement of the Cook cantilever
 * after the full load is applied
 * 
 * @code
 *     print_vertical_tip_displacement();
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Privateinterface"></a> 
 * <h3>Private interface</h3>
 * 
 
 * 
 * 
 * <a name="Solidmake_grid"></a> 
 * <h4>Solid::make_grid</h4>
 * 
 
 * 
 * On to the first of the private member functions. Here we create the
 * triangulation of the domain, for which we choose a scaled an anisotripically
 * discretised rectangle which is subsequently transformed into the correct
 * of the Cook cantilever. Each relevant boundary face is then given a boundary
 * ID number.
 * 
 
 * 
 * We then determine the volume of the reference configuration and print it
 * for comparison.
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   Point<dim> grid_y_transform (const Point<dim> &pt_in)
 *   {
 *     const double &x = pt_in[0];
 *     const double &y = pt_in[1];
 * 
 *     const double y_upper = 44.0 + (16.0/48.0)*x; // Line defining upper edge of beam
 *     const double y_lower =  0.0 + (44.0/48.0)*x; // Line defining lower edge of beam
 *     const double theta = y/44.0; // Fraction of height along left side of beam
 *     const double y_transform = (1-theta)*y_lower + theta*y_upper; // Final transformation
 * 
 *     Point<dim> pt_out = pt_in;
 *     pt_out[1] = y_transform;
 * 
 *     return pt_out;
 *   }
 * 
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::make_grid()
 *   {
 * @endcode
 * 
 * Divide the beam, but only along the x- and y-coordinate directions
 * 
 * @code
 *     std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge);
 * @endcode
 * 
 * Only allow one element through the thickness
 * (modelling a plane strain condition)
 * 
 * @code
 *     if (dim == 3)
 *       repetitions[dim-1] = 1;
 * 
 *     const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, -0.5) : Point<dim>(0.0, 0.0));
 *     const Point<dim> top_right = (dim == 3 ? Point<dim>(48.0, 44.0, 0.5) : Point<dim>(48.0, 44.0));
 * 
 *     GridGenerator::subdivided_hyper_rectangle(triangulation,
 *                                               repetitions,
 *                                               bottom_left,
 *                                               top_right);
 * 
 * @endcode
 * 
 * Since we wish to apply a Neumann BC to the right-hand surface, we
 * must find the cell faces in this part of the domain and mark them with
 * a distinct boundary ID number.  The faces we are looking for are on the
 * +x surface and will get boundary ID 11.
 * Dirichlet boundaries exist on the left-hand face of the beam (this fixed
 * boundary will get ID 1) and on the +Z and -Z faces (which correspond to
 * ID 2 and we will use to impose the plane strain condition)
 * 
 * @code
 *     const double tol_boundary = 1e-6;
 *     typename Triangulation<dim>::active_cell_iterator cell =
 *       triangulation.begin_active(), endc = triangulation.end();
 *     for (; cell != endc; ++cell)
 *       for (unsigned int face = 0;
 *            face < GeometryInfo<dim>::faces_per_cell; ++face)
 *         if (cell->face(face)->at_boundary() == true)
 *           {
 *             if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary)
 *               cell->face(face)->set_boundary_id(1); // -X faces
 *             else if (std::abs(cell->face(face)->center()[0] - 48.0) < tol_boundary)
 *               cell->face(face)->set_boundary_id(11); // +X faces
 *             else if (dim == 3 && std::abs(std::abs(cell->face(face)->center()[2]) - 0.5) < tol_boundary)
 *               cell->face(face)->set_boundary_id(2); // +Z and -Z faces
 *           }
 * 
 * @endcode
 * 
 * Transform the hyper-rectangle into the beam shape
 * 
 * @code
 *     GridTools::transform(&grid_y_transform<dim>, triangulation);
 * 
 *     GridTools::scale(parameters.scale, triangulation);
 * 
 *     vol_reference = GridTools::volume(triangulation);
 *     vol_current = vol_reference;
 *     std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidsystem_setup"></a> 
 * <h4>Solid::system_setup</h4>
 * 
 
 * 
 * Next we describe how the FE system is setup.  We first determine the number
 * of components per block. Since the displacement is a vector component, the
 * first dim components belong to it.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::system_setup()
 *   {
 *     timer.enter_subsection("Setup system");
 * 
 *     std::vector<unsigned int> block_component(n_components, u_dof); // Displacement
 * 
 * @endcode
 * 
 * The DOF handler is then initialised and we renumber the grid in an
 * efficient manner. We also record the number of DOFs per block.
 * 
 * @code
 *     dof_handler_ref.distribute_dofs(fe);
 *     DoFRenumbering::Cuthill_McKee(dof_handler_ref);
 *     DoFRenumbering::component_wise(dof_handler_ref, block_component);
 *     dofs_per_block = DoFTools::count_dofs_per_fe_block(dof_handler_ref, block_component);
 * 
 *     std::cout << "Triangulation:"
 *               << "\n\t Number of active cells: " << triangulation.n_active_cells()
 *               << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs()
 *               << std::endl;
 * 
 * @endcode
 * 
 * Setup the sparsity pattern and tangent matrix
 * 
 * @code
 *     tangent_matrix.clear();
 *     {
 *       const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
 * 
 *       BlockDynamicSparsityPattern csp(n_blocks, n_blocks);
 * 
 *       csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
 *       csp.collect_sizes();
 * 
 * @endcode
 * 
 * Naturally, for a one-field vector-valued problem, all of the
 * components of the system are coupled.
 * 
 * @code
 *       Table<2, DoFTools::Coupling> coupling(n_components, n_components);
 *       for (unsigned int ii = 0; ii < n_components; ++ii)
 *         for (unsigned int jj = 0; jj < n_components; ++jj)
 *           coupling[ii][jj] = DoFTools::always;
 *       DoFTools::make_sparsity_pattern(dof_handler_ref,
 *                                       coupling,
 *                                       csp,
 *                                       constraints,
 *                                       false);
 *       sparsity_pattern.copy_from(csp);
 *     }
 * 
 *     tangent_matrix.reinit(sparsity_pattern);
 * 
 * @endcode
 * 
 * We then set up storage vectors
 * 
 * @code
 *     system_rhs.reinit(dofs_per_block);
 *     system_rhs.collect_sizes();
 * 
 *     solution_n.reinit(dofs_per_block);
 *     solution_n.collect_sizes();
 * 
 * @endcode
 * 
 * ...and finally set up the quadrature
 * point history:
 * 
 * @code
 *     setup_qph();
 * 
 *     timer.leave_subsection();
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidsetup_qph"></a> 
 * <h4>Solid::setup_qph</h4>
 * The method used to store quadrature information is already described in
 * @ref step_18 "step-18" and @ref step_44 "step-44". Here we implement a similar setup for a SMP machine.
 * 
 
 * 
 * Firstly the actual QPH data objects are created. This must be done only
 * once the grid is refined to its finest level.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::setup_qph()
 *   {
 *     std::cout << "    Setting up quadrature point data..." << std::endl;
 * 
 *     quadrature_point_history.initialize(triangulation.begin_active(),
 *                                         triangulation.end(),
 *                                         n_q_points);
 * 
 * @endcode
 * 
 * Next we setup the initial quadrature point data. Note that when
 * the quadrature point data is retrieved, it is returned as a vector
 * of smart pointers.
 * 
 * @code
 *     for (typename Triangulation<dim>::active_cell_iterator cell =
 *            triangulation.begin_active(); cell != triangulation.end(); ++cell)
 *       {
 *         const std::vector<std::shared_ptr<PointHistory<dim,NumberType> > > lqph =
 *           quadrature_point_history.get_data(cell);
 *         Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 *         for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *           lqph[q_point]->setup_lqp(parameters);
 *       }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidsolve_nonlinear_timestep"></a> 
 * <h4>Solid::solve_nonlinear_timestep</h4>
 * 
 
 * 
 * The next function is the driver method for the Newton-Raphson scheme. At
 * its top we create a new vector to store the current Newton update step,
 * reset the error storage objects and print solver header.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void
 *   Solid<dim,NumberType>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
 *   {
 *     std::cout << std::endl << "Timestep " << time.get_timestep() << " @ "
 *               << time.current() << "s" << std::endl;
 * 
 *     BlockVector<double> newton_update(dofs_per_block);
 * 
 *     error_residual.reset();
 *     error_residual_0.reset();
 *     error_residual_norm.reset();
 *     error_update.reset();
 *     error_update_0.reset();
 *     error_update_norm.reset();
 * 
 *     print_conv_header();
 * 
 * @endcode
 * 
 * We now perform a number of Newton iterations to iteratively solve the
 * nonlinear problem.  Since the problem is fully nonlinear and we are
 * using a full Newton method, the data stored in the tangent matrix and
 * right-hand side vector is not reusable and must be cleared at each
 * Newton step.  We then initially build the right-hand side vector to
 * check for convergence (and store this value in the first iteration).
 * The unconstrained DOFs of the rhs vector hold the out-of-balance
 * forces. The building is done before assembling the system matrix as the
 * latter is an expensive operation and we can potentially avoid an extra
 * assembly process by not assembling the tangent matrix when convergence
 * is attained.
 * 
 * @code
 *     unsigned int newton_iteration = 0;
 *     for (; newton_iteration < parameters.max_iterations_NR;
 *          ++newton_iteration)
 *       {
 *         std::cout << " " << std::setw(2) << newton_iteration << " " << std::flush;
 * 
 * @endcode
 * 
 * If we have decided that we want to continue with the iteration, we
 * assemble the tangent, make and impose the Dirichlet constraints,
 * and do the solve of the linearized system:
 * 
 * @code
 *         make_constraints(newton_iteration);
 *         assemble_system(solution_delta);
 * 
 *         get_error_residual(error_residual);
 * 
 *         if (newton_iteration == 0)
 *           error_residual_0 = error_residual;
 * 
 * @endcode
 * 
 * We can now determine the normalised residual error and check for
 * solution convergence:
 * 
 * @code
 *         error_residual_norm = error_residual;
 *         error_residual_norm.normalise(error_residual_0);
 * 
 *         if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u
 *             && error_residual_norm.u <= parameters.tol_f)
 *           {
 *             std::cout << " CONVERGED! " << std::endl;
 *             print_conv_footer();
 * 
 *             break;
 *           }
 * 
 *         const std::pair<unsigned int, double>
 *         lin_solver_output = solve_linear_system(newton_update);
 * 
 *         get_error_update(newton_update, error_update);
 *         if (newton_iteration == 0)
 *           error_update_0 = error_update;
 * 
 * @endcode
 * 
 * We can now determine the normalised Newton update error, and
 * perform the actual update of the solution increment for the current
 * time step, update all quadrature point information pertaining to
 * this new displacement and stress state and continue iterating:
 * 
 * @code
 *         error_update_norm = error_update;
 *         error_update_norm.normalise(error_update_0);
 * 
 *         solution_delta += newton_update;
 * 
 *         std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
 *                   << std::scientific << lin_solver_output.first << "  "
 *                   << lin_solver_output.second << "  " << error_residual_norm.norm
 *                   << "  " << error_residual_norm.u << "  "
 *                   << "  " << error_update_norm.norm << "  " << error_update_norm.u
 *                   << "  " << std::endl;
 *       }
 * 
 * @endcode
 * 
 * At the end, if it turns out that we have in fact done more iterations
 * than the parameter file allowed, we raise an exception that can be
 * caught in the main() function. The call <code>AssertThrow(condition,
 * exc_object)</code> is in essence equivalent to <code>if (!cond) throw
 * exc_object;</code> but the former form fills certain fields in the
 * exception object that identify the location (filename and line number)
 * where the exception was raised to make it simpler to identify where the
 * problem happened.
 * 
 * @code
 *     AssertThrow (newton_iteration <= parameters.max_iterations_NR,
 *                  ExcMessage("No convergence in nonlinear solver!"));
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidprint_conv_headerSolidprint_conv_footerandSolidprint_vertical_tip_displacement"></a> 
 * <h4>Solid::print_conv_header, Solid::print_conv_footer and Solid::print_vertical_tip_displacement</h4>
 * 
 
 * 
 * This program prints out data in a nice table that is updated
 * on a per-iteration basis. The next two functions set up the table
 * header and footer:
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::print_conv_header()
 *   {
 *     static const unsigned int l_width = 87;
 * 
 *     for (unsigned int i = 0; i < l_width; ++i)
 *       std::cout << "_";
 *     std::cout << std::endl;
 * 
 *     std::cout << "    SOLVER STEP    "
 *               << " |  LIN_IT   LIN_RES    RES_NORM    "
 *               << " RES_U     NU_NORM     "
 *               << " NU_U " << std::endl;
 * 
 *     for (unsigned int i = 0; i < l_width; ++i)
 *       std::cout << "_";
 *     std::cout << std::endl;
 *   }
 * 
 * 
 * 
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::print_conv_footer()
 *   {
 *     static const unsigned int l_width = 87;
 * 
 *     for (unsigned int i = 0; i < l_width; ++i)
 *       std::cout << "_";
 *     std::cout << std::endl;
 * 
 *     std::cout << "Relative errors:" << std::endl
 *               << "Displacement:\t" << error_update.u / error_update_0.u << std::endl
 *               << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
 *               << "v / V_0:\t" << vol_current << " / " << vol_reference
 *               << std::endl;
 *   }
 * 
 * @endcode
 * 
 * At the end we also output the result that can be compared to that found in
 * the literature, namely the displacement at the upper right corner of the
 * beam.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::print_vertical_tip_displacement()
 *   {
 *     static const unsigned int l_width = 87;
 * 
 *     for (unsigned int i = 0; i < l_width; ++i)
 *       std::cout << "_";
 *     std::cout << std::endl;
 * 
 * @endcode
 * 
 * The measurement point, as stated in the reference paper, is at the midway
 * point of the surface on which the traction is applied.
 * 
 * @code
 *     const Point<dim> soln_pt = (dim == 3 ? 
 *                                 Point<dim>(48.0*parameters.scale, 52.0*parameters.scale, 0.5*parameters.scale) : 
 *                                 Point<dim>(48.0*parameters.scale, 52.0*parameters.scale));
 *     double vertical_tip_displacement = 0.0;
 *     double vertical_tip_displacement_check = 0.0;
 * 
 *     typename DoFHandler<dim>::active_cell_iterator cell =
 *       dof_handler_ref.begin_active(), endc = dof_handler_ref.end();
 *     for (; cell != endc; ++cell)
 *       {
 * @endcode
 * 
 * if (cell->point_inside(soln_pt) == true)
 * 
 * @code
 *         for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
 *           if (cell->vertex(v).distance(soln_pt) < 1e-6)
 *             {
 * @endcode
 * 
 * Extract y-component of solution at the given point
 * This point is coindicent with a vertex, so we can
 * extract it directly as we're using FE_Q finite elements
 * that have support at the vertices
 * 
 * @code
 *               vertical_tip_displacement = solution_n(cell->vertex_dof_index(v,u_dof+1));
 * 
 * @endcode
 * 
 * Sanity check using alternate method to extract the solution
 * at the given point. To do this, we must create an FEValues instance
 * to help us extract the solution value at the desired point
 * 
 * @code
 *               const MappingQ<dim> mapping (parameters.poly_degree);
 *               const Point<dim> qp_unit = mapping.transform_real_to_unit_cell(cell,soln_pt);
 *               const Quadrature<dim> soln_qrule (qp_unit);
 *               AssertThrow(soln_qrule.size() == 1, ExcInternalError());
 *               FEValues<dim> fe_values_soln (fe, soln_qrule, update_values);
 *               fe_values_soln.reinit(cell);
 * 
 * @endcode
 * 
 * Extract y-component of solution at given point
 * 
 * @code
 *               std::vector< Tensor<1,dim> > soln_values (soln_qrule.size());
 *               fe_values_soln[u_fe].get_function_values(solution_n,
 *                                                        soln_values);
 *               vertical_tip_displacement_check = soln_values[0][u_dof+1];
 * 
 *               break;
 *             }
 *       }
 *     AssertThrow(vertical_tip_displacement > 0.0, ExcMessage("Found no cell with point inside!"))
 * 
 *     std::cout << "Vertical tip displacement: " << vertical_tip_displacement
 *               << "\t Check: " << vertical_tip_displacement_check
 *               << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidget_error_residual"></a> 
 * <h4>Solid::get_error_residual</h4>
 * 
 
 * 
 * Determine the true residual error for the problem.  That is, determine the
 * error in the residual for the unconstrained degrees of freedom.  Note that to
 * do so, we need to ignore constrained DOFs by setting the residual in these
 * vector components to zero.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::get_error_residual(Errors &error_residual)
 *   {
 *     BlockVector<double> error_res(dofs_per_block);
 * 
 *     for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
 *       if (!constraints.is_constrained(i))
 *         error_res(i) = system_rhs(i);
 * 
 *     error_residual.norm = error_res.l2_norm();
 *     error_residual.u = error_res.block(u_dof).l2_norm();
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidget_error_udpate"></a> 
 * <h4>Solid::get_error_udpate</h4>
 * 
 
 * 
 * Determine the true Newton update error for the problem
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::get_error_update(const BlockVector<double> &newton_update,
 *                                                Errors &error_update)
 *   {
 *     BlockVector<double> error_ud(dofs_per_block);
 *     for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
 *       if (!constraints.is_constrained(i))
 *         error_ud(i) = newton_update(i);
 * 
 *     error_update.norm = error_ud.l2_norm();
 *     error_update.u = error_ud.block(u_dof).l2_norm();
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidget_total_solution"></a> 
 * <h4>Solid::get_total_solution</h4>
 * 
 
 * 
 * This function provides the total solution, which is valid at any Newton step.
 * This is required as, to reduce computational error, the total solution is
 * only updated at the end of the timestep.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   BlockVector<double>
 *   Solid<dim,NumberType>::get_total_solution(const BlockVector<double> &solution_delta) const
 *   {
 *     BlockVector<double> solution_total(solution_n);
 *     solution_total += solution_delta;
 *     return solution_total;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidassemble_system"></a> 
 * <h4>Solid::assemble_system</h4>
 * 
 
 * 
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   struct Assembler_Base
 *   {
 *     virtual ~Assembler_Base() {}
 * 
 * @endcode
 * 
 * Here we deal with the tangent matrix assembly structures. The
 * PerTaskData object stores local contributions.
 * 
 * @code
 *     struct PerTaskData_ASM
 *     {
 *       const Solid<dim,NumberType>          *solid;
 *       FullMatrix<double>                   cell_matrix;
 *       Vector<double>                       cell_rhs;
 *       std::vector<types::global_dof_index> local_dof_indices;
 * 
 *       PerTaskData_ASM(const Solid<dim,NumberType> *solid)
 *         :
 *         solid (solid),
 *         cell_matrix(solid->dofs_per_cell, solid->dofs_per_cell),
 *         cell_rhs(solid->dofs_per_cell),
 *         local_dof_indices(solid->dofs_per_cell)
 *       {}
 * 
 *       void reset()
 *       {
 *         cell_matrix = 0.0;
 *         cell_rhs = 0.0;
 *       }
 *     };
 * 
 * @endcode
 * 
 * On the other hand, the ScratchData object stores the larger objects such as
 * the shape-function values array (<code>Nx</code>) and a shape function
 * gradient and symmetric gradient vector which we will use during the
 * assembly.
 * 
 * @code
 *     struct ScratchData_ASM
 *     {
 *       const BlockVector<double>               &solution_total;
 *       std::vector<Tensor<2, dim,NumberType> >  solution_grads_u_total;
 * 
 *       FEValues<dim>                fe_values_ref;
 *       FEFaceValues<dim>            fe_face_values_ref;
 * 
 *       std::vector<std::vector<Tensor<2, dim,NumberType> > >         grad_Nx;
 *       std::vector<std::vector<SymmetricTensor<2,dim,NumberType> > > symm_grad_Nx;
 * 
 *       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_total)
 *         :
 *         solution_total(solution_total),
 *         solution_grads_u_total(qf_cell.size()),
 *         fe_values_ref(fe_cell, qf_cell, uf_cell),
 *         fe_face_values_ref(fe_cell, qf_face, uf_face),
 *         grad_Nx(qf_cell.size(),
 *                 std::vector<Tensor<2,dim,NumberType> >(fe_cell.dofs_per_cell)),
 *         symm_grad_Nx(qf_cell.size(),
 *                      std::vector<SymmetricTensor<2,dim,NumberType> >
 *                      (fe_cell.dofs_per_cell))
 *       {}
 * 
 *       ScratchData_ASM(const ScratchData_ASM &rhs)
 *         :
 *         solution_total (rhs.solution_total),
 *         solution_grads_u_total(rhs.solution_grads_u_total),
 *         fe_values_ref(rhs.fe_values_ref.get_fe(),
 *                       rhs.fe_values_ref.get_quadrature(),
 *                       rhs.fe_values_ref.get_update_flags()),
 *         fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
 *                            rhs.fe_face_values_ref.get_quadrature(),
 *                            rhs.fe_face_values_ref.get_update_flags()),
 *         grad_Nx(rhs.grad_Nx),
 *         symm_grad_Nx(rhs.symm_grad_Nx)
 *       {}
 * 
 *       void reset()
 *       {
 *         const unsigned int n_q_points = fe_values_ref.get_quadrature().size();
 *         const unsigned int n_dofs_per_cell = fe_values_ref.dofs_per_cell;
 *         for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *           {
 *             Assert( grad_Nx[q_point].size() == n_dofs_per_cell,
 *                     ExcInternalError());
 *             Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell,
 *                     ExcInternalError());
 * 
 *             solution_grads_u_total[q_point] = Tensor<2,dim,NumberType>();
 *             for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
 *               {
 *                 grad_Nx[q_point][k] = Tensor<2,dim,NumberType>();
 *                 symm_grad_Nx[q_point][k] = SymmetricTensor<2,dim,NumberType>();
 *               }
 *           }
 *       }
 * 
 *     };
 * 
 * @endcode
 * 
 * Of course, we still have to define how we assemble the tangent matrix
 * contribution for a single cell.
 * 
 * @code
 *     void
 *     assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                              ScratchData_ASM &scratch,
 *                              PerTaskData_ASM &data)
 *     {
 * @endcode
 * 
 * Due to the C++ specialization rules, we need one more
 * level of indirection in order to define the assembly
 * routine for all different number. The next function call
 * is specialized for each NumberType, but to prevent having
 * to specialize the whole class along with it we have inlined
 * the definition of the other functions that are common to
 * all implementations.
 * 
 * @code
 *       assemble_system_tangent_residual_one_cell(cell, scratch, data);
 *       assemble_neumann_contribution_one_cell(cell, scratch, data);
 *     }
 * 
 * @endcode
 * 
 * This function adds the local contribution to the system matrix.
 * 
 * @code
 *     void
 *     copy_local_to_global_ASM(const PerTaskData_ASM &data)
 *     {
 *       const AffineConstraints<double> &constraints = data.solid->constraints;
 *       BlockSparseMatrix<double> &tangent_matrix = const_cast<Solid<dim,NumberType> *>(data.solid)->tangent_matrix;
 *       BlockVector<double> &system_rhs =  const_cast<Solid<dim,NumberType> *>(data.solid)->system_rhs;
 * 
 *       constraints.distribute_local_to_global(
 *         data.cell_matrix, data.cell_rhs,
 *         data.local_dof_indices,
 *         tangent_matrix, system_rhs);
 *     }
 * 
 *   protected:
 * 
 * @endcode
 * 
 * This function needs to exist in the base class for
 * Workstream to work with a reference to the base class.
 * 
 * @code
 *     virtual void
 *     assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &/*cell*/,
 *                                               ScratchData_ASM &/*scratch*/,
 *                                               PerTaskData_ASM &/*data*/)
 *     {
 *       AssertThrow(false, ExcPureFunctionCalled());
 *     }
 * 
 *     void
 *     assemble_neumann_contribution_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                                            ScratchData_ASM &scratch,
 *                                            PerTaskData_ASM &data)
 *     {
 * @endcode
 * 
 * Aliases for data referenced from the Solid class
 * 
 * @code
 *       const unsigned int &n_q_points_f = data.solid->n_q_points_f;
 *       const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
 *       const Parameters::AllParameters &parameters = data.solid->parameters;
 *       const Time &time = data.solid->time;
 *       const FESystem<dim> &fe = data.solid->fe;
 *       const unsigned int &u_dof = data.solid->u_dof;
 * 
 * @endcode
 * 
 * Next we assemble the Neumann contribution. We first check to see it the
 * cell face exists on a boundary on which a traction is applied and add
 * the contribution if this is the case.
 * 
 * @code
 *       for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
 *            ++face)
 *         if (cell->face(face)->at_boundary() == true
 *             && cell->face(face)->boundary_id() == 11)
 *           {
 *             scratch.fe_face_values_ref.reinit(cell, face);
 * 
 *             for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
 *                  ++f_q_point)
 *               {
 * @endcode
 * 
 * We specify the traction in reference configuration.
 * For this problem, a defined total vertical force is applied
 * in the reference configuration.
 * The direction of the applied traction is assumed not to
 * evolve with the deformation of the domain.
 * 
 
 * 
 * Note that the contributions to the right hand side vector we
 * compute here only exist in the displacement components of the
 * vector.
 * 
 * @code
 *                 const double time_ramp = (time.current() / time.end());
 *                 const double magnitude  = (1.0/(16.0*parameters.scale*1.0*parameters.scale))*time_ramp; // (Total force) / (RHS surface area)
 *                 Tensor<1,dim> dir;
 *                 dir[1] = 1.0;
 *                 const Tensor<1, dim> traction  = magnitude*dir;
 * 
 *                 for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                   {
 *                     const unsigned int i_group =
 *                       fe.system_to_base_index(i).first.first;
 * 
 *                     if (i_group == u_dof)
 *                       {
 *                         const unsigned int component_i =
 *                           fe.system_to_component_index(i).first;
 *                         const double Ni =
 *                           scratch.fe_face_values_ref.shape_value(i,
 *                                                                  f_q_point);
 *                         const double JxW = scratch.fe_face_values_ref.JxW(
 *                                              f_q_point);
 * 
 *                         data.cell_rhs(i) += (Ni * traction[component_i])
 *                                             * JxW;
 *                       }
 *                   }
 *               }
 *           }
 *     }
 * 
 *   };
 * 
 *   template <int dim>
 *   struct Assembler<dim,double> : Assembler_Base<dim,double>
 *   {
 *     typedef double NumberType;
 *     using typename Assembler_Base<dim,NumberType>::ScratchData_ASM;
 *     using typename Assembler_Base<dim,NumberType>::PerTaskData_ASM;
 * 
 *     virtual ~Assembler() {}
 * 
 *     virtual void
 *     assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                                               ScratchData_ASM &scratch,
 *                                               PerTaskData_ASM &data)
 *     {
 * @endcode
 * 
 * Aliases for data referenced from the Solid class
 * 
 * @code
 *       const unsigned int &n_q_points = data.solid->n_q_points;
 *       const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
 *       const FESystem<dim> &fe = data.solid->fe;
 *       const unsigned int &u_dof = data.solid->u_dof;
 *       const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
 * 
 *       data.reset();
 *       scratch.reset();
 *       scratch.fe_values_ref.reinit(cell);
 *       cell->get_dof_indices(data.local_dof_indices);
 * 
 *       const std::vector<std::shared_ptr<const PointHistory<dim,NumberType> > > lqph =
 *         const_cast<const Solid<dim,NumberType> *>(data.solid)->quadrature_point_history.get_data(cell);
 *       Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 * @endcode
 * 
 * We first need to find the solution gradients at quadrature points
 * inside the current cell and then we update each local QP using the
 * displacement gradient:
 * 
 * @code
 *       scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
 *                                                          scratch.solution_grads_u_total);
 * 
 * @endcode
 * 
 * Now we build the local cell stiffness matrix. Since the global and
 * local system matrices are symmetric, we can exploit this property by
 * building only the lower half of the local matrix and copying the values
 * to the upper half.
 *       
 
 * 
 * In doing so, we first extract some configuration dependent variables
 * from our QPH history objects for the current quadrature point.
 * 
 * @code
 *       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *         {
 *           const Tensor<2,dim,NumberType> &grad_u = scratch.solution_grads_u_total[q_point];
 *           const Tensor<2,dim,NumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
 *           const NumberType               det_F = determinant(F);
 *           const Tensor<2,dim,NumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
 *           const SymmetricTensor<2,dim,NumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
 *           const Tensor<2,dim,NumberType> F_inv = invert(F);
 *           Assert(det_F > NumberType(0.0), ExcInternalError());
 * 
 *           for (unsigned int k = 0; k < dofs_per_cell; ++k)
 *             {
 *               const unsigned int k_group = fe.system_to_base_index(k).first.first;
 * 
 *               if (k_group == u_dof)
 *                 {
 *                   scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
 *                   scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
 *                 }
 *               else
 *                 Assert(k_group <= u_dof, ExcInternalError());
 *             }
 * 
 *           const SymmetricTensor<2,dim,NumberType> tau = lqph[q_point]->get_tau(det_F,b_bar);
 *           const SymmetricTensor<4,dim,NumberType> Jc  = lqph[q_point]->get_Jc(det_F,b_bar);
 *           const Tensor<2,dim,NumberType> tau_ns (tau);
 * 
 * @endcode
 * 
 * Next we define some aliases to make the assembly process easier to
 * follow
 * 
 * @code
 *           const std::vector<SymmetricTensor<2, dim> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
 *           const std::vector<Tensor<2, dim> > &grad_Nx = scratch.grad_Nx[q_point];
 *           const double JxW = scratch.fe_values_ref.JxW(q_point);
 * 
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               const unsigned int component_i = fe.system_to_component_index(i).first;
 *               const unsigned int i_group     = fe.system_to_base_index(i).first.first;
 * 
 *               if (i_group == u_dof)
 *                 data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
 *               else
 *                 Assert(i_group <= u_dof, ExcInternalError());
 * 
 *               for (unsigned int j = 0; j <= i; ++j)
 *                 {
 *                   const unsigned int component_j = fe.system_to_component_index(j).first;
 *                   const unsigned int j_group     = fe.system_to_base_index(j).first.first;
 * 
 * @endcode
 * 
 * This is the @f$\mathsf{\mathbf{k}}_{\mathbf{u} \mathbf{u}}@f$
 * contribution. It comprises a material contribution, and a
 * geometrical stress contribution which is only added along
 * the local matrix diagonals:
 * 
 * @code
 *                   if ((i_group == j_group) && (i_group == u_dof))
 *                     {
 *                       data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
 *                                                 * symm_grad_Nx[j] * JxW;
 *                       if (component_i == component_j) // geometrical stress contribution
 *                         data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau_ns
 *                                                   * grad_Nx[j][component_j] * JxW;
 *                     }
 *                   else
 *                     Assert((i_group <= u_dof) && (j_group <= u_dof),
 *                            ExcInternalError());
 *                 }
 *             }
 *         }
 * 
 * 
 * @endcode
 * 
 * Finally, we need to copy the lower half of the local matrix into the
 * upper half:
 * 
 * @code
 *       for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *         for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
 *           data.cell_matrix(i, j) = data.cell_matrix(j, i);
 *     }
 * 
 *   };
 * 
 * #ifdef ENABLE_SACADO_FORMULATION
 * 
 * 
 *   template <int dim>
 *   struct Assembler<dim,Sacado::Fad::DFad<double> > : Assembler_Base<dim,Sacado::Fad::DFad<double> >
 *   {
 *     typedef Sacado::Fad::DFad<double> ADNumberType;
 *     using typename Assembler_Base<dim,ADNumberType>::ScratchData_ASM;
 *     using typename Assembler_Base<dim,ADNumberType>::PerTaskData_ASM;
 * 
 *     virtual ~Assembler() {}
 * 
 *     virtual void
 *     assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                                               ScratchData_ASM &scratch,
 *                                               PerTaskData_ASM &data)
 *     {
 * @endcode
 * 
 * Aliases for data referenced from the Solid class
 * 
 * @code
 *       const unsigned int &n_q_points = data.solid->n_q_points;
 *       const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
 *       const FESystem<dim> &fe = data.solid->fe;
 *       const unsigned int &u_dof = data.solid->u_dof;
 *       const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
 * 
 *       data.reset();
 *       scratch.reset();
 *       scratch.fe_values_ref.reinit(cell);
 *       cell->get_dof_indices(data.local_dof_indices);
 * 
 *       const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
 *         const_cast<const Solid<dim,ADNumberType> *>(data.solid)->quadrature_point_history.get_data(cell);
 *       Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 *       const unsigned int n_independent_variables = data.local_dof_indices.size();
 *       std::vector<double> local_dof_values(n_independent_variables);
 *       cell->get_dof_values(scratch.solution_total,
 *                            local_dof_values.begin(),
 *                            local_dof_values.end());
 * 
 * @endcode
 * 
 * We now retrieve a set of degree-of-freedom values that
 * have the operations that are performed with them tracked.
 * 
 * @code
 *       std::vector<ADNumberType> local_dof_values_ad (n_independent_variables);
 *       for (unsigned int i=0; i<n_independent_variables; ++i)
 *         local_dof_values_ad[i] = ADNumberType(n_independent_variables, i, local_dof_values[i]);
 * 
 * @endcode
 * 
 * Compute all values, gradients etc. based on sensitive
 * AD degree-of-freedom values.
 * 
 * @code
 *       scratch.fe_values_ref[u_fe].get_function_gradients_from_local_dof_values(
 *         local_dof_values_ad,
 *         scratch.solution_grads_u_total);
 * 
 * @endcode
 * 
 * Accumulate the residual value for each degree of freedom.
 * Note: Its important that the vectors is initialised (zero'd) correctly.
 * 
 * @code
 *       std::vector<ADNumberType> residual_ad (dofs_per_cell, ADNumberType(0.0));
 *       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *         {
 *           const Tensor<2,dim,ADNumberType> &grad_u = scratch.solution_grads_u_total[q_point];
 *           const Tensor<2,dim,ADNumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
 *           const ADNumberType               det_F = determinant(F);
 *           const Tensor<2,dim,ADNumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
 *           const SymmetricTensor<2,dim,ADNumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
 *           const Tensor<2,dim,ADNumberType> F_inv = invert(F);
 *           Assert(det_F > ADNumberType(0.0), ExcInternalError());
 * 
 *           for (unsigned int k = 0; k < dofs_per_cell; ++k)
 *             {
 *               const unsigned int k_group = fe.system_to_base_index(k).first.first;
 * 
 *               if (k_group == u_dof)
 *                 {
 *                   scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
 *                   scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
 *                 }
 *               else
 *                 Assert(k_group <= u_dof, ExcInternalError());
 *             }
 * 
 *           const SymmetricTensor<2,dim,ADNumberType> tau = lqph[q_point]->get_tau(det_F,b_bar);
 * 
 * @endcode
 * 
 * Next we define some position-dependent aliases, again to
 * make the assembly process easier to follow.
 * 
 * @code
 *           const std::vector<SymmetricTensor<2, dim,ADNumberType> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
 *           const double JxW = scratch.fe_values_ref.JxW(q_point);
 * 
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               const unsigned int i_group     = fe.system_to_base_index(i).first.first;
 * 
 *               if (i_group == u_dof)
 *                 residual_ad[i] += (symm_grad_Nx[i] * tau) * JxW;
 *               else
 *                 Assert(i_group <= u_dof, ExcInternalError());
 *             }
 *         }
 * 
 *       for (unsigned int I=0; I<n_independent_variables; ++I)
 *         {
 *           const ADNumberType &residual_I = residual_ad[I];
 *           data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
 *           for (unsigned int J=0; J<n_independent_variables; ++J)
 *             {
 * @endcode
 * 
 * Compute the gradients of the residual entry [forward-mode]
 * 
 * @code
 *               data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
 *             }
 *         }
 *     }
 * 
 *   };
 * 
 * 
 *   template <int dim>
 *   struct Assembler<dim,Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > > : Assembler_Base<dim,Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > >
 *   {
 *     typedef Sacado::Fad::DFad<double>       ADDerivType;
 *     typedef Sacado::Rad::ADvar<ADDerivType> ADNumberType;
 *     using typename Assembler_Base<dim,ADNumberType>::ScratchData_ASM;
 *     using typename Assembler_Base<dim,ADNumberType>::PerTaskData_ASM;
 * 
 *     virtual ~Assembler() {}
 * 
 *     virtual void
 *     assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                                               ScratchData_ASM &scratch,
 *                                               PerTaskData_ASM &data)
 *     {
 * @endcode
 * 
 * Aliases for data referenced from the Solid class
 * 
 * @code
 *       const unsigned int &n_q_points = data.solid->n_q_points;
 *       const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
 * 
 *       data.reset();
 *       scratch.reset();
 *       scratch.fe_values_ref.reinit(cell);
 *       cell->get_dof_indices(data.local_dof_indices);
 * 
 *       const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
 *         data.solid->quadrature_point_history.get_data(cell);
 *       Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 *       const unsigned int n_independent_variables = data.local_dof_indices.size();
 *       std::vector<double> local_dof_values(n_independent_variables);
 *       cell->get_dof_values(scratch.solution_total,
 *                            local_dof_values.begin(),
 *                            local_dof_values.end());
 * 
 * @endcode
 * 
 * We now retrieve a set of degree-of-freedom values that
 * have the operations that are performed with them tracked.
 * 
 * @code
 *       std::vector<ADNumberType> local_dof_values_ad (n_independent_variables);
 *       for (unsigned int i=0; i<n_independent_variables; ++i)
 *         local_dof_values_ad[i] = ADDerivType(n_independent_variables, i, local_dof_values[i]);
 * 
 * @endcode
 * 
 * Compute all values, gradients etc. based on sensitive
 * AD degree-of-freedom values.
 * 
 * @code
 *       scratch.fe_values_ref[u_fe].get_function_gradients_from_local_dof_values(
 *         local_dof_values_ad,
 *         scratch.solution_grads_u_total);
 * 
 * @endcode
 * 
 * Next we compute the total potential energy of the element.
 * This is defined as follows:
 * Total energy = (internal - external) energies
 * Note: Its important that this value is initialised (zero'd) correctly.
 * 
 * @code
 *       ADNumberType cell_energy_ad = ADNumberType(0.0);
 *       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *         {
 *           const Tensor<2,dim,ADNumberType> &grad_u = scratch.solution_grads_u_total[q_point];
 *           const Tensor<2,dim,ADNumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
 *           const ADNumberType               det_F = determinant(F);
 *           const Tensor<2,dim,ADNumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
 *           const SymmetricTensor<2,dim,ADNumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
 *           Assert(det_F > ADNumberType(0.0), ExcInternalError());
 * 
 * @endcode
 * 
 * Next we define some position-dependent aliases, again to
 * make the assembly process easier to follow.
 * 
 * @code
 *           const double JxW = scratch.fe_values_ref.JxW(q_point);
 * 
 *           const ADNumberType Psi = lqph[q_point]->get_Psi(det_F,b_bar);
 * 
 * @endcode
 * 
 * We extract the configuration-dependent material energy
 * from our QPH history objects for the current quadrature point
 * and integrate its contribution to increment the total
 * cell energy.
 * 
 * @code
 *           cell_energy_ad += Psi * JxW;
 *         }
 * 
 * @endcode
 * 
 * Compute derivatives of reverse-mode AD variables
 * 
 * @code
 *       ADNumberType::Gradcomp();
 * 
 *       for (unsigned int I=0; I<n_independent_variables; ++I)
 *         {
 * @endcode
 * 
 * This computes the adjoint df/dX_{i} [reverse-mode]
 * 
 * @code
 *           const ADDerivType residual_I = local_dof_values_ad[I].adj();
 *           data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
 *           for (unsigned int J=0; J<n_independent_variables; ++J)
 *             {
 * @endcode
 * 
 * Compute the gradients of the residual entry [forward-mode]
 * 
 * @code
 *               data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
 *             }
 *         }
 *     }
 * 
 *   };
 * 
 * 
 * #endif
 * 
 * 
 * @endcode
 * 
 * Since we use TBB for assembly, we simply setup a copy of the
 * data structures required for the process and pass them, along
 * with the memory addresses of the assembly functions to the
 * WorkStream object for processing. Note that we must ensure that
 * the matrix is reset before any assembly operations can occur.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::assemble_system(const BlockVector<double> &solution_delta)
 *   {
 *     timer.enter_subsection("Assemble linear system");
 *     std::cout << " ASM " << std::flush;
 * 
 *     tangent_matrix = 0.0;
 *     system_rhs = 0.0;
 * 
 *     const UpdateFlags uf_cell(update_gradients |
 *                               update_JxW_values);
 *     const UpdateFlags uf_face(update_values |
 *                               update_JxW_values);
 * 
 *     const BlockVector<double> solution_total(get_total_solution(solution_delta));
 *     typename Assembler_Base<dim,NumberType>::PerTaskData_ASM per_task_data(this);
 *     typename Assembler_Base<dim,NumberType>::ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total);
 *     Assembler<dim,NumberType> assembler;
 * 
 *     WorkStream::run(dof_handler_ref.begin_active(),
 *                     dof_handler_ref.end(),
 *                     static_cast<Assembler_Base<dim,NumberType>&>(assembler),
 *                     &Assembler_Base<dim,NumberType>::assemble_system_one_cell,
 *                     &Assembler_Base<dim,NumberType>::copy_local_to_global_ASM,
 *                     scratch_data,
 *                     per_task_data);
 * 
 *     timer.leave_subsection();
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidmake_constraints"></a> 
 * <h4>Solid::make_constraints</h4>
 * The constraints for this problem are simple to describe.
 * However, since we are dealing with an iterative Newton method,
 * it should be noted that any displacement constraints should only
 * be specified at the zeroth iteration and subsequently no
 * additional contributions are to be made since the constraints
 * are already exactly satisfied.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::make_constraints(const int &it_nr)
 *   {
 *     std::cout << " CST " << std::flush;
 * 
 * @endcode
 * 
 * Since the constraints are different at different Newton iterations, we
 * need to clear the constraints matrix and completely rebuild
 * it. However, after the first iteration, the constraints remain the same
 * and we can simply skip the rebuilding step if we do not clear it.
 * 
 * @code
 *     if (it_nr > 1)
 *       return;
 *     const bool apply_dirichlet_bc = (it_nr == 0);
 * 
 * @endcode
 * 
 * The boundary conditions for the indentation problem are as follows: On
 * the -x face (ID = 1), we set up a zero-displacement condition, -y and +y traction 
 * free boundary condition (don't need to take care); -z and +z faces (ID = 2) are 
 * not allowed to move along z axis so that it is a plane strain problem. 
 * Finally, as described earlier, +x face (ID = 11) has an the applied 
 * distributed shear force (converted by total force per unit area) which 
 * needs to be taken care as an inhomogeneous Newmann boundary condition.
 *     
 
 * 
 * In the following, we will have to tell the function interpolation
 * boundary values which components of the solution vector should be
 * constrained (i.e., whether it's the x-, y-, z-displacements or
 * combinations thereof). This is done using ComponentMask objects (see
 * @ref GlossComponentMask) which we can get from the finite element if we
 * provide it with an extractor object for the component we wish to
 * select. To this end we first set up such extractor objects and later
 * use it when generating the relevant component masks:
 * 
 
 * 
 * 
 * @code
 *     if (apply_dirichlet_bc)
 *     {
 *       constraints.clear();
 * 
 * @endcode
 * 
 * Fixed left hand side of the beam
 * 
 * @code
 *       {
 *         const int boundary_id = 1;
 *         VectorTools::interpolate_boundary_values(dof_handler_ref,
 *                                                 boundary_id,
 *                                                 Functions::ZeroFunction<dim>(n_components),
 *                                                 constraints,
 *                                                 fe.component_mask(u_fe));
 *       }
 * 
 * @endcode
 * 
 * Zero Z-displacement through thickness direction
 * This corresponds to a plane strain condition being imposed on the beam
 * 
 * @code
 *       if (dim == 3)
 *       {
 *         const int boundary_id = 2;
 *         const FEValuesExtractors::Scalar z_displacement(2);
 *         VectorTools::interpolate_boundary_values(dof_handler_ref,
 *                                                 boundary_id,
 *                                                 Functions::ZeroFunction<dim>(n_components),
 *                                                 constraints,
 *                                                 fe.component_mask(z_displacement));
 *       }
 *     }
 *     else
 *     {
 *       if (constraints.has_inhomogeneities())
 *       {
 *         AffineConstraints<double> homogeneous_constraints(constraints);
 *         for (unsigned int dof = 0; dof != dof_handler_ref.n_dofs(); ++dof)
 *           if (homogeneous_constraints.is_inhomogeneously_constrained(dof))
 *             homogeneous_constraints.set_inhomogeneity(dof, 0.0);
 *         constraints.clear();
 *         constraints.copy_from(homogeneous_constraints);
 *       }
 *     }
 * 
 *     constraints.close();
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="Solidsolve_linear_system"></a> 
 * <h4>Solid::solve_linear_system</h4>
 * As the system is composed of a single block, defining a solution scheme
 * for the linear problem is straight-forward.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   std::pair<unsigned int, double>
 *   Solid<dim,NumberType>::solve_linear_system(BlockVector<double> &newton_update)
 *   {
 *     BlockVector<double> A(dofs_per_block);
 *     BlockVector<double> B(dofs_per_block);
 * 
 *     unsigned int lin_it = 0;
 *     double lin_res = 0.0;
 * 
 * @endcode
 * 
 * We solve for the incremental displacement @f$d\mathbf{u}@f$.
 * 
 * @code
 *     {
 *       timer.enter_subsection("Linear solver");
 *       std::cout << " SLV " << std::flush;
 *       if (parameters.type_lin == "CG")
 *         {
 *           const int solver_its = static_cast<unsigned int>(
 *                                     tangent_matrix.block(u_dof, u_dof).m()
 *                                     * parameters.max_iterations_lin);
 *           const double tol_sol = parameters.tol_lin
 *                                  * system_rhs.block(u_dof).l2_norm();
 * 
 *           SolverControl solver_control(solver_its, tol_sol);
 * 
 *           GrowingVectorMemory<Vector<double> > GVM;
 *           SolverCG<Vector<double> > solver_CG(solver_control, GVM);
 * 
 * @endcode
 * 
 * We've chosen by default a SSOR preconditioner as it appears to
 * provide the fastest solver convergence characteristics for this
 * problem on a single-thread machine.  However, for multicore
 * computing, the Jacobi preconditioner which is multithreaded may
 * converge quicker for larger linear systems.
 * 
 * @code
 *           PreconditionSelector<SparseMatrix<double>, Vector<double> >
 *           preconditioner (parameters.preconditioner_type,
 *                           parameters.preconditioner_relaxation);
 *           preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
 * 
 *           solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
 *                           newton_update.block(u_dof),
 *                           system_rhs.block(u_dof),
 *                           preconditioner);
 * 
 *           lin_it = solver_control.last_step();
 *           lin_res = solver_control.last_value();
 *         }
 *       else if (parameters.type_lin == "Direct")
 *         {
 * @endcode
 * 
 * Otherwise if the problem is small
 * enough, a direct solver can be
 * utilised.
 * 
 * @code
 *           SparseDirectUMFPACK A_direct;
 *           A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
 *           A_direct.vmult(newton_update.block(u_dof), system_rhs.block(u_dof));
 * 
 *           lin_it = 1;
 *           lin_res = 0.0;
 *         }
 *       else
 *         Assert (false, ExcMessage("Linear solver type not implemented"));
 * 
 *       timer.leave_subsection();
 *     }
 * 
 * @endcode
 * 
 * Now that we have the displacement update, distribute the constraints
 * back to the Newton update:
 * 
 * @code
 *     constraints.distribute(newton_update);
 * 
 *     return std::make_pair(lin_it, lin_res);
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="Solidoutput_results"></a> 
 * <h4>Solid::output_results</h4>
 * Here we present how the results are written to file to be viewed
 * using ParaView or Visit. The method is similar to that shown in the
 * tutorials so will not be discussed in detail.
 * 
 * @code
 *   template <int dim,typename NumberType>
 *   void Solid<dim,NumberType>::output_results() const
 *   {
 *     DataOut<dim> data_out;
 *     std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *     data_component_interpretation(dim,
 *                                   DataComponentInterpretation::component_is_part_of_vector);
 * 
 *     std::vector<std::string> solution_name(dim, "displacement");
 * 
 *     data_out.attach_dof_handler(dof_handler_ref);
 *     data_out.add_data_vector(solution_n,
 *                              solution_name,
 *                              DataOut<dim>::type_dof_data,
 *                              data_component_interpretation);
 * 
 * @endcode
 * 
 * Since we are dealing with a large deformation problem, it would be nice
 * to display the result on a displaced grid!  The MappingQEulerian class
 * linked with the DataOut class provides an interface through which this
 * can be achieved without physically moving the grid points in the
 * Triangulation object ourselves.  We first need to copy the solution to
 * a temporary vector and then create the Eulerian mapping. We also
 * specify the polynomial degree to the DataOut object in order to produce
 * a more refined output data set when higher order polynomials are used.
 * 
 * @code
 *     Vector<double> soln(solution_n.size());
 *     for (unsigned int i = 0; i < soln.size(); ++i)
 *       soln(i) = solution_n(i);
 *     MappingQEulerian<dim> q_mapping(degree, dof_handler_ref, soln);
 *     data_out.build_patches(q_mapping, degree);
 * 
 *     std::ostringstream filename;
 *     filename << "solution-" << time.get_timestep() << ".vtk";
 * 
 *     std::ofstream output(filename.str().c_str());
 *     data_out.write_vtk(output);
 *   }
 * 
 * }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Mainfunction"></a> 
 * <h3>Main function</h3>
 * Lastly we provide the main driver function which appears
 * no different to the other tutorials.
 * 
 * @code
 * int main (int argc, char *argv[])
 * {
 *   using namespace dealii;
 *   using namespace Cook_Membrane;
 * 
 *   const unsigned int dim = 3;
 *   try
 *     {
 *       deallog.depth_console(0);
 *       Parameters::AllParameters parameters("parameters.prm");
 *       if (parameters.automatic_differentiation_order == 0)
 *         {
 *           std::cout << "Assembly method: Residual and linearisation are computed manually." << std::endl;
 * 
 * @endcode
 * 
 * Allow multi-threading
 * 
 * @code
 *           Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
 *                                                               ::numbers::invalid_unsigned_int);
 * 
 *           typedef double NumberType;
 *           Solid<dim,NumberType> solid_3d(parameters);
 *           solid_3d.run();
 *         }
 * #ifdef ENABLE_SACADO_FORMULATION
 *       else if (parameters.automatic_differentiation_order == 1)
 *         {
 *           std::cout << "Assembly method: Residual computed manually; linearisation performed using AD." << std::endl;
 * 
 * @endcode
 * 
 * Allow multi-threading
 * 
 * @code
 *           Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
 *                                                               ::numbers::invalid_unsigned_int);
 * 
 *           typedef Sacado::Fad::DFad<double> NumberType;
 *           Solid<dim,NumberType> solid_3d(parameters);
 *           solid_3d.run();
 *         }
 *       else if (parameters.automatic_differentiation_order == 2)
 *         {
 *           std::cout << "Assembly method: Residual and linearisation computed using AD." << std::endl;
 * 
 * @endcode
 * 
 * Sacado Rad-Fad is not thread-safe, so disable threading.
 * Parallisation using MPI would be possible though.
 * 
 * @code
 *           Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
 *                                                               1);
 * 
 *           typedef Sacado::Rad::ADvar< Sacado::Fad::DFad<double> > NumberType;
 *           Solid<dim,NumberType> solid_3d(parameters);
 *           solid_3d.run();
 *         }
 * #endif
 *       else
 *         {
 *           AssertThrow(false,
 *                       ExcMessage("The selected assembly method is not supported. "
 *                                  "You need deal.II 9.0 and Trilinos with the Sacado package "
 *                                  "to enable assembly using automatic differentiation."));
 *         }
 *     }
 *   catch (std::exception &exc)
 *     {
 *       std::cerr << std::endl << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Exception on processing: " << std::endl << exc.what()
 *                 << std::endl << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 * 
 *       return 1;
 *     }
 *   catch (...)
 *     {
 *       std::cerr << std::endl << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Unknown exception!" << std::endl << "Aborting!"
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       return 1;
 *     }
 * 
 *   return 0;
 * }
 * @endcode
 
 
*/