step-44.h

Go to the documentation of this file.

)
 *     {}
 * 
 * @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
 *     Parameters::AllParameters parameters;
 * 
 * @endcode
 * 
 * ...the volume of the reference configuration...
 * 
 * @code
 *     double vol_reference;
 * 
 * @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;
 *     mutable 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>>
 *       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;
 *     const unsigned int               dofs_per_cell;
 *     const FEValuesExtractors::Vector u_fe;
 *     const FEValuesExtractors::Scalar p_fe;
 *     const FEValuesExtractors::Scalar J_fe;
 * 
 * @endcode
 * 
 * Description of how the block-system is arranged. There are 3 blocks,
 * the first contains a vector DOF @f$\mathbf{u}@f$ while the other two
 * describe scalar DOFs, @f$\widetilde{p}@f$ and @f$\widetilde{J}@f$.
 * 
 * @code
 *     static const unsigned int n_blocks          = 3;
 *     static const unsigned int n_components      = dim + 2;
 *     static const unsigned int first_u_component = 0;
 *     static const unsigned int p_component       = dim;
 *     static const unsigned int J_component       = dim + 1;
 * 
 *     enum
 *     {
 *       u_dof = 0,
 *       p_dof = 1,
 *       J_dof = 2
 *     };
 * 
 *     std::vector<types::global_dof_index> dofs_per_block;
 *     std::vector<types::global_dof_index> element_indices_u;
 *     std::vector<types::global_dof_index> element_indices_p;
 *     std::vector<types::global_dof_index> element_indices_J;
 * 
 * @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 an 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
 * normalization factors.
 * 
 * @code
 *     struct Errors
 *     {
 *       Errors()
 *         : norm(1.0)
 *         , u(1.0)
 *         , p(1.0)
 *         , J(1.0)
 *       {}
 * 
 *       void reset()
 *       {
 *         norm = 1.0;
 *         u    = 1.0;
 *         p    = 1.0;
 *         J    = 1.0;
 *       }
 *       void normalize(const Errors &rhs)
 *       {
 *         if (rhs.norm != 0.0)
 *           norm /= rhs.norm;
 *         if (rhs.u != 0.0)
 *           u /= rhs.u;
 *         if (rhs.p != 0.0)
 *           p /= rhs.p;
 *         if (rhs.J != 0.0)
 *           J /= rhs.J;
 *       }
 * 
 *       double norm, u, p, J;
 *     };
 * 
 *     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);
 * 
 *     std::pair<double, double> get_error_dilation() const;
 * 
 * @endcode
 * 
 * Compute the volume in the spatial configuration
 * 
 * @code
 *     double compute_vol_current() const;
 * 
 * @endcode
 * 
 * Print information to screen in a pleasing way...
 * 
 * @code
 *     static void print_conv_header();
 * 
 *     void print_conv_footer();
 *   };
 * 
 * @endcode
 * 
 * 
 * <a name="ImplementationofthecodeSolidcodeclass"></a> 
 * <h3>Implementation of the <code>Solid</code> class</h3>
 * 

 * 
 * 
 * <a name="Publicinterface"></a> 
 * <h4>Public interface</h4>
 * 

 * 
 * We initialize the Solid class using data extracted from the parameter file.
 * 
 * @code
 *   template <int dim>
 *   Solid<dim>::Solid(const std::string &input_file)
 *     : parameters(input_file)
 *     , vol_reference(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, and discontinuous pressure and dilatation DOFs. In an attempt to
 * satisfy the Babuska-Brezzi or LBB stability conditions (see Hughes
 * (2000)), we setup a @f$Q_n \times DGPM_{n-1} \times DGPM_{n-1}@f$
 * system. @f$Q_2 \times DGPM_1 \times DGPM_1@f$ elements satisfy this
 * condition, while @f$Q_1 \times DGPM_0 \times DGPM_0@f$ elements do
 * not. However, it has been shown that the latter demonstrate good
 * convergence characteristics nonetheless.
 * 
 * @code
 *     fe(FE_Q<dim>(parameters.poly_degree),
 *        dim, // displacement
 *        FE_DGPMonomial<dim>(parameters.poly_degree - 1),
 *        1, // pressure
 *        FE_DGPMonomial<dim>(parameters.poly_degree - 1),
 *        1)
 *     , // dilatation
 *     dof_handler(triangulation)
 *     , dofs_per_cell(fe.dofs_per_cell)
 *     , u_fe(first_u_component)
 *     , p_fe(p_component)
 *     , J_fe(J_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())
 *   {
 *     Assert(dim == 2 || dim == 3,
 *            ExcMessage("This problem only works in 2 or 3 space dimensions."));
 *     determine_component_extractors();
 *   }
 * 
 * 
 * @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, setting the initial dilatation
 * values, and then output the initial grid before starting the simulation
 * proper with the first time (and loading)
 * increment.
 *   

 * 
 * Care must be taken (or at least some thought given) when imposing the
 * constraint @f$\widetilde{J}=1@f$ on the initial solution field. The constraint
 * corresponds to the determinant of the deformation gradient in the
 * undeformed configuration, which is the identity tensor. We use
 * FE_DGPMonomial bases to interpolate the dilatation field, thus we can't
 * simply set the corresponding dof to unity as they correspond to the
 * monomial coefficients. Thus we use the VectorTools::project function to do
 * the work for us. The VectorTools::project function requires an argument
 * indicating the hanging node constraints. We have none in this program
 * So we have to create a constraint object. In its original state, constraint
 * objects are unsorted, and have to be sorted (using the
 * AffineConstraints::close function) before they can be used. Have a look at
 * @ref step_21 "step-21" for more information. We only need to enforce the initial condition
 * on the dilatation. In order to do this, we make use of a
 * ComponentSelectFunction which acts as a mask and sets the J_component of
 * n_components to 1. This is exactly what we want. Have a look at its usage
 * in @ref step_20 "step-20" for more information.
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::run()
 *   {
 *     make_grid();
 *     system_setup();
 *     {
 *       AffineConstraints<double> constraints;
 *       constraints.close();
 * 
 *       const ComponentSelectFunction<dim> J_mask(J_component, n_components);
 * 
 *       VectorTools::project(
 *         dof_handler, constraints, QGauss<dim>(degree + 2), J_mask, solution_n);
 *     }
 *     output_results();
 *     time.increment();
 * 
 * @endcode
 * 
 * We then declare the incremental solution update @f$\varDelta
 * \mathbf{\Xi} \dealcoloneq \{\varDelta \mathbf{u},\varDelta \widetilde{p},
 * \varDelta \widetilde{J} \}@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
 * 
 * 
 * <a name="Privateinterface"></a> 
 * <h3>Private interface</h3>
 * 

 * 
 * 
 * <a name="Threadingbuildingblocksstructures"></a> 
 * <h4>Threading-building-blocks structures</h4>
 * 

 * 
 * The first group of private member functions is related to parallelization.
 * We use the Threading Building Blocks library (TBB) to perform as many
 * computationally intensive distributed tasks as possible. In particular, we
 * assemble the tangent matrix and right hand side vector, the static
 * condensation contributions, and update data stored at the quadrature points
 * using TBB. Our main tool for this is the WorkStream class (see the @ref
 * threads module for more information).
 * 

 * 
 * Firstly we deal with the tangent matrix assembly structures.  The
 * PerTaskData object stores local contributions.
 * 
 * @code
 *   template <int dim>
 *   struct Solid<dim>::PerTaskData_K
 *   {
 *     FullMatrix<double>                   cell_matrix;
 *     std::vector<types::global_dof_index> local_dof_indices;
 * 
 *     PerTaskData_K(const unsigned int dofs_per_cell)
 *       : cell_matrix(dofs_per_cell, dofs_per_cell)
 *       , local_dof_indices(dofs_per_cell)
 *     {}
 * 
 *     void reset()
 *     {
 *       cell_matrix = 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
 *   template <int dim>
 *   struct Solid<dim>::ScratchData_K
 *   {
 *     FEValues<dim> fe_values;
 * 
 *     std::vector<std::vector<double>>                  Nx;
 *     std::vector<std::vector<Tensor<2, dim>>>          grad_Nx;
 *     std::vector<std::vector<SymmetricTensor<2, dim>>> symm_grad_Nx;
 * 
 *     ScratchData_K(const FiniteElement<dim> &fe_cell,
 *                   const QGauss<dim> &       qf_cell,
 *                   const UpdateFlags         uf_cell)
 *       : fe_values(fe_cell, qf_cell, uf_cell)
 *       , Nx(qf_cell.size(), std::vector<double>(fe_cell.dofs_per_cell))
 *       , grad_Nx(qf_cell.size(),
 *                 std::vector<Tensor<2, dim>>(fe_cell.dofs_per_cell))
 *       , symm_grad_Nx(qf_cell.size(),
 *                      std::vector<SymmetricTensor<2, dim>>(
 *                        fe_cell.dofs_per_cell))
 *     {}
 * 
 *     ScratchData_K(const ScratchData_K &rhs)
 *       : fe_values(rhs.fe_values.get_fe(),
 *                   rhs.fe_values.get_quadrature(),
 *                   rhs.fe_values.get_update_flags())
 *       , Nx(rhs.Nx)
 *       , grad_Nx(rhs.grad_Nx)
 *       , symm_grad_Nx(rhs.symm_grad_Nx)
 *     {}
 * 
 *     void reset()
 *     {
 *       const unsigned int n_q_points      = Nx.size();
 *       const unsigned int n_dofs_per_cell = Nx[0].size();
 *       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *         {
 *           Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
 *           Assert(grad_Nx[q_point].size() == n_dofs_per_cell,
 *                  ExcInternalError());
 *           Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell,
 *                  ExcInternalError());
 *           for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
 *             {
 *               Nx[q_point][k]           = 0.0;
 *               grad_Nx[q_point][k]      = 0.0;
 *               symm_grad_Nx[q_point][k] = 0.0;
 *             }
 *         }
 *     }
 *   };
 * 
 * @endcode
 * 
 * Next, the same approach is used for the right-hand side assembly.  The
 * PerTaskData object again stores local contributions and the ScratchData
 * object the shape function object and precomputed values vector:
 * 
 * @code
 *   template <int dim>
 *   struct Solid<dim>::PerTaskData_RHS
 *   {
 *     Vector<double>                       cell_rhs;
 *     std::vector<types::global_dof_index> local_dof_indices;
 * 
 *     PerTaskData_RHS(const unsigned int dofs_per_cell)
 *       : cell_rhs(dofs_per_cell)
 *       , local_dof_indices(dofs_per_cell)
 *     {}
 * 
 *     void reset()
 *     {
 *       cell_rhs = 0.0;
 *     }
 *   };
 * 
 * 
 *   template <int dim>
 *   struct Solid<dim>::ScratchData_RHS
 *   {
 *     FEValues<dim>     fe_values;
 *     FEFaceValues<dim> fe_face_values;
 * 
 *     std::vector<std::vector<double>>                  Nx;
 *     std::vector<std::vector<SymmetricTensor<2, dim>>> symm_grad_Nx;
 * 
 *     ScratchData_RHS(const FiniteElement<dim> &fe_cell,
 *                     const QGauss<dim> &       qf_cell,
 *                     const UpdateFlags         uf_cell,
 *                     const QGauss<dim - 1> &   qf_face,
 *                     const UpdateFlags         uf_face)
 *       : fe_values(fe_cell, qf_cell, uf_cell)
 *       , fe_face_values(fe_cell, qf_face, uf_face)
 *       , Nx(qf_cell.size(), std::vector<double>(fe_cell.dofs_per_cell))
 *       , symm_grad_Nx(qf_cell.size(),
 *                      std::vector<SymmetricTensor<2, dim>>(
 *                        fe_cell.dofs_per_cell))
 *     {}
 * 
 *     ScratchData_RHS(const ScratchData_RHS &rhs)
 *       : fe_values(rhs.fe_values.get_fe(),
 *                   rhs.fe_values.get_quadrature(),
 *                   rhs.fe_values.get_update_flags())
 *       , fe_face_values(rhs.fe_face_values.get_fe(),
 *                        rhs.fe_face_values.get_quadrature(),
 *                        rhs.fe_face_values.get_update_flags())
 *       , Nx(rhs.Nx)
 *       , symm_grad_Nx(rhs.symm_grad_Nx)
 *     {}
 * 
 *     void reset()
 *     {
 *       const unsigned int n_q_points      = Nx.size();
 *       const unsigned int n_dofs_per_cell = Nx[0].size();
 *       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *         {
 *           Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
 *           Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell,
 *                  ExcInternalError());
 *           for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
 *             {
 *               Nx[q_point][k]           = 0.0;
 *               symm_grad_Nx[q_point][k] = 0.0;
 *             }
 *         }
 *     }
 *   };
 * 
 * @endcode
 * 
 * Then we define structures to assemble the statically condensed tangent
 * matrix. Recall that we wish to solve for a displacement-based formulation.
 * We do the condensation at the element level as the @f$\widetilde{p}@f$ and
 * @f$\widetilde{J}@f$ fields are element-wise discontinuous.  As these operations
 * are matrix-based, we need to setup a number of matrices to store the local
 * contributions from a number of the tangent matrix sub-blocks.  We place
 * these in the PerTaskData struct.
 *   

 * 
 * We choose not to reset any data in the <code>reset()</code> function as the
 * matrix extraction and replacement tools will take care of this
 * 
 * @code
 *   template <int dim>
 *   struct Solid<dim>::PerTaskData_SC
 *   {
 *     FullMatrix<double>                   cell_matrix;
 *     std::vector<types::global_dof_index> local_dof_indices;
 * 
 *     FullMatrix<double> k_orig;
 *     FullMatrix<double> k_pu;
 *     FullMatrix<double> k_pJ;
 *     FullMatrix<double> k_JJ;
 *     FullMatrix<double> k_pJ_inv;
 *     FullMatrix<double> k_bbar;
 *     FullMatrix<double> A;
 *     FullMatrix<double> B;
 *     FullMatrix<double> C;
 * 
 *     PerTaskData_SC(const unsigned int dofs_per_cell,
 *                    const unsigned int n_u,
 *                    const unsigned int n_p,
 *                    const unsigned int n_J)
 *       : cell_matrix(dofs_per_cell, dofs_per_cell)
 *       , local_dof_indices(dofs_per_cell)
 *       , k_orig(dofs_per_cell, dofs_per_cell)
 *       , k_pu(n_p, n_u)
 *       , k_pJ(n_p, n_J)
 *       , k_JJ(n_J, n_J)
 *       , k_pJ_inv(n_p, n_J)
 *       , k_bbar(n_u, n_u)
 *       , A(n_J, n_u)
 *       , B(n_J, n_u)
 *       , C(n_p, n_u)
 *     {}
 * 
 *     void reset()
 *     {}
 *   };
 * 
 * 
 * @endcode
 * 
 * The ScratchData object for the operations we wish to perform here is empty
 * since we need no temporary data, but it still needs to be defined for the
 * current implementation of TBB in deal.II.  So we create a dummy struct for
 * this purpose.
 * 
 * @code
 *   template <int dim>
 *   struct Solid<dim>::ScratchData_SC
 *   {
 *     void reset()
 *     {}
 *   };
 * 
 * 
 * @endcode
 * 
 * And finally we define the structures to assist with updating the quadrature
 * point information. Similar to the SC assembly process, we do not need the
 * PerTaskData object (since there is nothing to store here) but must define
 * one nonetheless. Note that this is because for the operation that we have
 * here -- updating the data on quadrature points -- the operation is purely
 * local: the things we do on every cell get consumed on every cell, without
 * any global aggregation operation as is usually the case when using the
 * WorkStream class. The fact that we still have to define a per-task data
 * structure points to the fact that the WorkStream class may be ill-suited to
 * this operation (we could, in principle simply create a new task using
 * Threads::new_task for each cell) but there is not much harm done to doing
 * it this way anyway.
 * Furthermore, should there be different material models associated with a
 * quadrature point, requiring varying levels of computational expense, then
 * the method used here could be advantageous.
 * 
 * @code
 *   template <int dim>
 *   struct Solid<dim>::PerTaskData_UQPH
 *   {
 *     void reset()
 *     {}
 *   };
 * 
 * 
 * @endcode
 * 
 * The ScratchData object will be used to store an alias for the solution
 * vector so that we don't have to copy this large data structure. We then
 * define a number of vectors to extract the solution values and gradients at
 * the quadrature points.
 * 
 * @code
 *   template <int dim>
 *   struct Solid<dim>::ScratchData_UQPH
 *   {
 *     const BlockVector<double> &solution_total;
 * 
 *     std::vector<Tensor<2, dim>> solution_grads_u_total;
 *     std::vector<double>         solution_values_p_total;
 *     std::vector<double>         solution_values_J_total;
 * 
 *     FEValues<dim> fe_values;
 * 
 *     ScratchData_UQPH(const FiniteElement<dim> & fe_cell,
 *                      const QGauss<dim> &        qf_cell,
 *                      const UpdateFlags          uf_cell,
 *                      const BlockVector<double> &solution_total)
 *       : solution_total(solution_total)
 *       , solution_grads_u_total(qf_cell.size())
 *       , solution_values_p_total(qf_cell.size())
 *       , solution_values_J_total(qf_cell.size())
 *       , fe_values(fe_cell, qf_cell, uf_cell)
 *     {}
 * 
 *     ScratchData_UQPH(const ScratchData_UQPH &rhs)
 *       : solution_total(rhs.solution_total)
 *       , solution_grads_u_total(rhs.solution_grads_u_total)
 *       , solution_values_p_total(rhs.solution_values_p_total)
 *       , solution_values_J_total(rhs.solution_values_J_total)
 *       , fe_values(rhs.fe_values.get_fe(),
 *                   rhs.fe_values.get_quadrature(),
 *                   rhs.fe_values.get_update_flags())
 *     {}
 * 
 *     void reset()
 *     {
 *       const unsigned int n_q_points = solution_grads_u_total.size();
 *       for (unsigned int q = 0; q < n_q_points; ++q)
 *         {
 *           solution_grads_u_total[q]  = 0.0;
 *           solution_values_p_total[q] = 0.0;
 *           solution_values_J_total[q] = 0.0;
 *         }
 *     }
 *   };
 * 
 * 
 * @endcode
 * 
 * 
 * <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 the scaled cube with each
 * face given a boundary ID number.  The grid must be refined at least once
 * for the indentation problem.
 *   

 * 
 * We then determine the volume of the reference configuration and print it
 * for comparison:
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::make_grid()
 *   {
 *     GridGenerator::hyper_rectangle(
 *       triangulation,
 *       (dim == 3 ? Point<dim>(0.0, 0.0, 0.0) : Point<dim>(0.0, 0.0)),
 *       (dim == 3 ? Point<dim>(1.0, 1.0, 1.0) : Point<dim>(1.0, 1.0)),
 *       true);
 *     GridTools::scale(parameters.scale, triangulation);
 *     triangulation.refine_global(std::max(1U, parameters.global_refinement));
 * 
 *     vol_reference = GridTools::volume(triangulation);
 *     std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
 * 
 * @endcode
 * 
 * Since we wish to apply a Neumann BC to a patch on the top 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
 * +y surface and will get boundary ID 6 (zero through five are already
 * used when creating the six faces of the cube domain):
 * 
 * @code
 *     for (const auto &cell : triangulation.active_cell_iterators())
 *       for (const auto &face : cell->face_iterators())
 *         {
 *           if (face->at_boundary() == true &&
 *               face->center()[1] == 1.0 * parameters.scale)
 *             {
 *               if (dim == 3)
 *                 {
 *                   if (face->center()[0] < 0.5 * parameters.scale &&
 *                       face->center()[2] < 0.5 * parameters.scale)
 *                     face->set_boundary_id(6);
 *                 }
 *               else
 *                 {
 *                   if (face->center()[0] < 0.5 * parameters.scale)
 *                     face->set_boundary_id(6);
 *                 }
 *             }
 *         }
 *   }
 * 
 * 
 * @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, while the next two describe scalar
 * pressure and dilatation DOFs.
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::system_setup()
 *   {
 *     timer.enter_subsection("Setup system");
 * 
 *     std::vector<unsigned int> block_component(n_components,
 *                                               u_dof); // Displacement
 *     block_component[p_component] = p_dof;             // Pressure
 *     block_component[J_component] = J_dof;             // Dilatation
 * 
 * @endcode
 * 
 * The DOF handler is then initialized and we renumber the grid in an
 * efficient manner. We also record the number of DOFs per block.
 * 
 * @code
 *     dof_handler.distribute_dofs(fe);
 *     DoFRenumbering::Cuthill_McKee(dof_handler);
 *     DoFRenumbering::component_wise(dof_handler, block_component);
 * 
 *     dofs_per_block =
 *       DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
 * 
 *     std::cout << "Triangulation:"
 *               << "\n\t Number of active cells: "
 *               << triangulation.n_active_cells()
 *               << "\n\t Number of degrees of freedom: " << dof_handler.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];
 *       const types::global_dof_index n_dofs_p = dofs_per_block[p_dof];
 *       const types::global_dof_index n_dofs_J = dofs_per_block[J_dof];
 * 
 *       BlockDynamicSparsityPattern dsp(n_blocks, n_blocks);
 * 
 *       dsp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
 *       dsp.block(u_dof, p_dof).reinit(n_dofs_u, n_dofs_p);
 *       dsp.block(u_dof, J_dof).reinit(n_dofs_u, n_dofs_J);
 * 
 *       dsp.block(p_dof, u_dof).reinit(n_dofs_p, n_dofs_u);
 *       dsp.block(p_dof, p_dof).reinit(n_dofs_p, n_dofs_p);
 *       dsp.block(p_dof, J_dof).reinit(n_dofs_p, n_dofs_J);
 * 
 *       dsp.block(J_dof, u_dof).reinit(n_dofs_J, n_dofs_u);
 *       dsp.block(J_dof, p_dof).reinit(n_dofs_J, n_dofs_p);
 *       dsp.block(J_dof, J_dof).reinit(n_dofs_J, n_dofs_J);
 *       dsp.collect_sizes();
 * 
 * @endcode
 * 
 * The global system matrix initially has the following structure
 * @f{align*}
 * \underbrace{\begin{bmatrix}
 * \mathsf{\mathbf{K}}_{uu}  & \mathsf{\mathbf{K}}_{u\widetilde{p}} &
 * \mathbf{0}
 * \\ \mathsf{\mathbf{K}}_{\widetilde{p}u} & \mathbf{0} &
 * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}
 * \\ \mathbf{0} & \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}} &
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
 * \end{bmatrix}}_{\mathsf{\mathbf{K}}(\mathbf{\Xi}_{\textrm{i}})}
 * \underbrace{\begin{bmatrix}
 * d \mathsf{u}
 * \\  d \widetilde{\mathsf{\mathbf{p}}}
 * \\  d \widetilde{\mathsf{\mathbf{J}}}
 * \end{bmatrix}}_{d \mathbf{\Xi}}
 * =
 * \underbrace{\begin{bmatrix}
 * \mathsf{\mathbf{F}}_{u}(\mathbf{u}_{\textrm{i}})
 * \\ \mathsf{\mathbf{F}}_{\widetilde{p}}(\widetilde{p}_{\textrm{i}})
 * \\ \mathsf{\mathbf{F}}_{\widetilde{J}}(\widetilde{J}_{\textrm{i}})
 * \end{bmatrix}}_{ \mathsf{\mathbf{F}}(\mathbf{\Xi}_{\textrm{i}}) } \, .
 * @f}
 * We optimize the sparsity pattern to reflect this structure
 * and prevent unnecessary data creation for the right-diagonal
 * block components.
 * 
 * @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)
 *           if (((ii < p_component) && (jj == J_component)) ||
 *               ((ii == J_component) && (jj < p_component)) ||
 *               ((ii == p_component) && (jj == p_component)))
 *             coupling[ii][jj] = DoFTools::none;
 *           else
 *             coupling[ii][jj] = DoFTools::always;
 *       DoFTools::make_sparsity_pattern(
 *         dof_handler, coupling, dsp, constraints, false);
 *       sparsity_pattern.copy_from(dsp);
 *     }
 * 
 *     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="Soliddetermine_component_extractors"></a> 
 * <h4>Solid::determine_component_extractors</h4>
 * Next we compute some information from the FE system that describes which
 * local element DOFs are attached to which block component.  This is used
 * later to extract sub-blocks from the global matrix.
 *   

 * 
 * In essence, all we need is for the FESystem object to indicate to which
 * block component a DOF on the reference cell is attached to.  Currently, the
 * interpolation fields are setup such that 0 indicates a displacement DOF, 1
 * a pressure DOF and 2 a dilatation DOF.
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::determine_component_extractors()
 *   {
 *     element_indices_u.clear();
 *     element_indices_p.clear();
 *     element_indices_J.clear();
 * 
 *     for (unsigned int k = 0; k < fe.dofs_per_cell; ++k)
 *       {
 *         const unsigned int k_group = fe.system_to_base_index(k).first.first;
 *         if (k_group == u_dof)
 *           element_indices_u.push_back(k);
 *         else if (k_group == p_dof)
 *           element_indices_p.push_back(k);
 *         else if (k_group == J_dof)
 *           element_indices_J.push_back(k);
 *         else
 *           {
 *             Assert(k_group <= J_dof, ExcInternalError());
 *           }
 *       }
 *   }
 * 
 * @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". 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>
 *   void Solid<dim>::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 (const auto &cell : triangulation.active_cell_iterators())
 *       {
 *         const std::vector<std::shared_ptr<PointHistory<dim>>> 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="Solidupdate_qph_incremental"></a> 
 * <h4>Solid::update_qph_incremental</h4>
 * As the update of QP information occurs frequently and involves a number of
 * expensive operations, we define a multithreaded approach to distributing
 * the task across a number of CPU cores.
 *   

 * 
 * To start this, we first we need to obtain the total solution as it stands
 * at this Newton increment and then create the initial copy of the scratch
 * and copy data objects:
 * 
 * @code
 *   template <int dim>
 *   void
 *   Solid<dim>::update_qph_incremental(const BlockVector<double> &solution_delta)
 *   {
 *     timer.enter_subsection("Update QPH data");
 *     std::cout << " UQPH " << std::flush;
 * 
 *     const BlockVector<double> solution_total(
 *       get_total_solution(solution_delta));
 * 
 *     const UpdateFlags uf_UQPH(update_values | update_gradients);
 *     PerTaskData_UQPH  per_task_data_UQPH;
 *     ScratchData_UQPH  scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total);
 * 
 * @endcode
 * 
 * We then pass them and the one-cell update function to the WorkStream to
 * be processed:
 * 
 * @code
 *     WorkStream::run(dof_handler.active_cell_iterators(),
 *                     *this,
 *                     &Solid::update_qph_incremental_one_cell,
 *                     &Solid::copy_local_to_global_UQPH,
 *                     scratch_data_UQPH,
 *                     per_task_data_UQPH);
 * 
 *     timer.leave_subsection();
 *   }
 * 
 * 
 * @endcode
 * 
 * Now we describe how we extract data from the solution vector and pass it
 * along to each QP storage object for processing.
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::update_qph_incremental_one_cell(
 *     const typename DoFHandler<dim>::active_cell_iterator &cell,
 *     ScratchData_UQPH &                                    scratch,
 *     PerTaskData_UQPH & /*data*/)
 *   {
 *     const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
 *       quadrature_point_history.get_data(cell);
 *     Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 *     Assert(scratch.solution_grads_u_total.size() == n_q_points,
 *            ExcInternalError());
 *     Assert(scratch.solution_values_p_total.size() == n_q_points,
 *            ExcInternalError());
 *     Assert(scratch.solution_values_J_total.size() == n_q_points,
 *            ExcInternalError());
 * 
 *     scratch.reset();
 * 
 * @endcode
 * 
 * We first need to find the values and gradients at quadrature points
 * inside the current cell and then we update each local QP using the
 * displacement gradient and total pressure and dilatation solution
 * values:
 * 
 * @code
 *     scratch.fe_values.reinit(cell);
 *     scratch.fe_values[u_fe].get_function_gradients(
 *       scratch.solution_total, scratch.solution_grads_u_total);
 *     scratch.fe_values[p_fe].get_function_values(
 *       scratch.solution_total, scratch.solution_values_p_total);
 *     scratch.fe_values[J_fe].get_function_values(
 *       scratch.solution_total, scratch.solution_values_J_total);
 * 
 *     for (const unsigned int q_point :
 *          scratch.fe_values.quadrature_point_indices())
 *       lqph[q_point]->update_values(scratch.solution_grads_u_total[q_point],
 *                                    scratch.solution_values_p_total[q_point],
 *                                    scratch.solution_values_J_total[q_point]);
 *   }
 * 
 * 
 * @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>
 *   void Solid<dim>::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;
 * 
 *         tangent_matrix = 0.0;
 *         system_rhs     = 0.0;
 * 
 *         assemble_system_rhs();
 *         get_error_residual(error_residual);
 * 
 *         if (newton_iteration == 0)
 *           error_residual_0 = error_residual;
 * 
 * @endcode
 * 
 * We can now determine the normalized residual error and check for
 * solution convergence:
 * 
 * @code
 *         error_residual_norm = error_residual;
 *         error_residual_norm.normalize(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;
 *           }
 * 
 * @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 linearised system:
 * 
 * @code
 *         assemble_system_tangent();
 *         make_constraints(newton_iteration);
 *         constraints.condense(tangent_matrix, system_rhs);
 * 
 *         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 normalized 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.normalize(error_update_0);
 * 
 *         solution_delta += newton_update;
 *         update_qph_incremental(solution_delta);
 * 
 *         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_residual_norm.p << "  "
 *                   << error_residual_norm.J << "  " << error_update_norm.norm
 *                   << "  " << error_update_norm.u << "  " << error_update_norm.p
 *                   << "  " << error_update_norm.J << "  " << 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_headerandSolidprint_conv_footer"></a> 
 * <h4>Solid::print_conv_header and Solid::print_conv_footer</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>
 *   void Solid<dim>::print_conv_header()
 *   {
 *     static const unsigned int l_width = 155;
 * 
 *     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     RES_P      RES_J     NU_NORM     "
 *               << " NU_U       NU_P       NU_J " << std::endl;
 * 
 *     for (unsigned int i = 0; i < l_width; ++i)
 *       std::cout << "_";
 *     std::cout << std::endl;
 *   }
 * 
 * 
 * 
 *   template <int dim>
 *   void Solid<dim>::print_conv_footer()
 *   {
 *     static const unsigned int l_width = 155;
 * 
 *     for (unsigned int i = 0; i < l_width; ++i)
 *       std::cout << "_";
 *     std::cout << std::endl;
 * 
 *     const std::pair<double, double> error_dil = get_error_dilation();
 * 
 *     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
 *               << "Dilatation:\t" << error_dil.first << std::endl
 *               << "v / V_0:\t" << error_dil.second * vol_reference << " / "
 *               << vol_reference << " = " << error_dil.second << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidget_error_dilation"></a> 
 * <h4>Solid::get_error_dilation</h4>
 * 

 * 
 * Calculate the volume of the domain in the spatial configuration
 * 
 * @code
 *   template <int dim>
 *   double Solid<dim>::compute_vol_current() const
 *   {
 *     double vol_current = 0.0;
 * 
 *     FEValues<dim> fe_values(fe, qf_cell, update_JxW_values);
 * 
 *     for (const auto &cell : triangulation.active_cell_iterators())
 *       {
 *         fe_values.reinit(cell);
 * 
 * @endcode
 * 
 * In contrast to that which was previously called for,
 * in this instance the quadrature point data is specifically
 * non-modifiable since we will only be accessing data.
 * We ensure that the right get_data function is called by
 * marking this update function as constant.
 * 
 * @code
 *         const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *           quadrature_point_history.get_data(cell);
 *         Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 *         for (const unsigned int q_point : fe_values.quadrature_point_indices())
 *           {
 *             const double det_F_qp = lqph[q_point]->get_det_F();
 *             const double JxW      = fe_values.JxW(q_point);
 * 
 *             vol_current += det_F_qp * JxW;
 *           }
 *       }
 *     Assert(vol_current > 0.0, ExcInternalError());
 *     return vol_current;
 *   }
 * 
 * @endcode
 * 
 * Calculate how well the dilatation @f$\widetilde{J}@f$ agrees with @f$J
 * \dealcoloneq \textrm{det}\ \mathbf{F}@f$ from the @f$L^2@f$ error @f$ \bigl[
 * \int_{\Omega_0} {[ J - \widetilde{J}]}^{2}\textrm{d}V \bigr]^{1/2}@f$.
 * We also return the ratio of the current volume of the
 * domain to the reference volume. This is of interest for incompressible
 * media where we want to check how well the isochoric constraint has been
 * enforced.
 * 
 * @code
 *   template <int dim>
 *   std::pair<double, double> Solid<dim>::get_error_dilation() const
 *   {
 *     double dil_L2_error = 0.0;
 * 
 *     FEValues<dim> fe_values(fe, qf_cell, update_JxW_values);
 * 
 *     for (const auto &cell : triangulation.active_cell_iterators())
 *       {
 *         fe_values.reinit(cell);
 * 
 *         const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *           quadrature_point_history.get_data(cell);
 *         Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 *         for (const unsigned int q_point : fe_values.quadrature_point_indices())
 *           {
 *             const double det_F_qp   = lqph[q_point]->get_det_F();
 *             const double J_tilde_qp = lqph[q_point]->get_J_tilde();
 *             const double the_error_qp_squared =
 *               std::pow((det_F_qp - J_tilde_qp), 2);
 *             const double JxW = fe_values.JxW(q_point);
 * 
 *             dil_L2_error += the_error_qp_squared * JxW;
 *           }
 *       }
 * 
 *     return std::make_pair(std::sqrt(dil_L2_error),
 *                           compute_vol_current() / vol_reference);
 *   }
 * 
 * 
 * @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>
 *   void Solid<dim>::get_error_residual(Errors &error_residual)
 *   {
 *     BlockVector<double> error_res(dofs_per_block);
 * 
 *     for (unsigned int i = 0; i < dof_handler.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();
 *     error_residual.p    = error_res.block(p_dof).l2_norm();
 *     error_residual.J    = error_res.block(J_dof).l2_norm();
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Solidget_error_update"></a> 
 * <h4>Solid::get_error_update</h4>
 * 

 * 
 * Determine the true Newton update error for the problem
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::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.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();
 *     error_update.p    = error_ud.block(p_dof).l2_norm();
 *     error_update.J    = error_ud.block(J_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>
 *   BlockVector<double> Solid<dim>::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_tangent"></a> 
 * <h4>Solid::assemble_system_tangent</h4>
 * 

 * 
 * 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>
 *   void Solid<dim>::assemble_system_tangent()
 *   {
 *     timer.enter_subsection("Assemble tangent matrix");
 *     std::cout << " ASM_K " << std::flush;
 * 
 *     tangent_matrix = 0.0;
 * 
 *     const UpdateFlags uf_cell(update_values | update_gradients |
 *                               update_JxW_values);
 * 
 *     PerTaskData_K per_task_data(dofs_per_cell);
 *     ScratchData_K scratch_data(fe, qf_cell, uf_cell);
 * 
 * @endcode
 * 
 * The syntax used here to pass data to the WorkStream class
 * is discussed in @ref step_14 "step-14". We need to use this particular
 * call to WorkStream because assemble_system_tangent_one_cell
 * is a constant function and copy_local_to_global_K is
 * non-constant.
 * 
 * @code
 *     WorkStream::run(
 *       dof_handler.active_cell_iterators(),
 *       [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *              ScratchData_K &                                       scratch,
 *              PerTaskData_K &                                       data) {
 *         this->assemble_system_tangent_one_cell(cell, scratch, data);
 *       },
 *       [this](const PerTaskData_K &data) { this->copy_local_to_global_K(data); },
 *       scratch_data,
 *       per_task_data);
 * 
 *     timer.leave_subsection();
 *   }
 * 
 * @endcode
 * 
 * This function adds the local contribution to the system matrix.
 * Note that we choose not to use the constraint matrix to do the
 * job for us because the tangent matrix and residual processes have
 * been split up into two separate functions.
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::copy_local_to_global_K(const PerTaskData_K &data)
 *   {
 *     for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *       for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *         tangent_matrix.add(data.local_dof_indices[i],
 *                            data.local_dof_indices[j],
 *                            data.cell_matrix(i, j));
 *   }
 * 
 * @endcode
 * 
 * Of course, we still have to define how we assemble the tangent matrix
 * contribution for a single cell.  We first need to reset and initialize some
 * of the scratch data structures and retrieve some basic information
 * regarding the DOF numbering on this cell.  We can precalculate the cell
 * shape function values and gradients. Note that the shape function gradients
 * are defined with regard to the current configuration.  That is
 * @f$\textrm{grad}\ \boldsymbol{\varphi} = \textrm{Grad}\ \boldsymbol{\varphi}
 * \ \mathbf{F}^{-1}@f$.
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::assemble_system_tangent_one_cell(
 *     const typename DoFHandler<dim>::active_cell_iterator &cell,
 *     ScratchData_K &                                       scratch,
 *     PerTaskData_K &                                       data) const
 *   {
 *     data.reset();
 *     scratch.reset();
 *     scratch.fe_values.reinit(cell);
 *     cell->get_dof_indices(data.local_dof_indices);
 * 
 *     const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *       quadrature_point_history.get_data(cell);
 *     Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 *     for (const unsigned int q_point :
 *          scratch.fe_values.quadrature_point_indices())
 *       {
 *         const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
 *         for (const unsigned int k : scratch.fe_values.dof_indices())
 *           {
 *             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[u_fe].gradient(k, q_point) * F_inv;
 *                 scratch.symm_grad_Nx[q_point][k] =
 *                   symmetrize(scratch.grad_Nx[q_point][k]);
 *               }
 *             else if (k_group == p_dof)
 *               scratch.Nx[q_point][k] =
 *                 scratch.fe_values[p_fe].value(k, q_point);
 *             else if (k_group == J_dof)
 *               scratch.Nx[q_point][k] =
 *                 scratch.fe_values[J_fe].value(k, q_point);
 *             else
 *               Assert(k_group <= J_dof, ExcInternalError());
 *           }
 *       }
 * 
 * @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.  So we only assemble half of the
 * @f$\mathsf{\mathbf{k}}_{uu}@f$, @f$\mathsf{\mathbf{k}}_{\widetilde{p}
 * \widetilde{p}} = \mathbf{0}@f$, @f$\mathsf{\mathbf{k}}_{\widetilde{J}
 * \widetilde{J}}@f$ blocks, while the whole
 * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$,
 * @f$\mathsf{\mathbf{k}}_{u \widetilde{J}} = \mathbf{0}@f$,
 * @f$\mathsf{\mathbf{k}}_{u \widetilde{p}}@f$ blocks are built.
 *     

 * 
 * In doing so, we first extract some configuration dependent variables
 * from our quadrature history objects for the current quadrature point.
 * 
 * @code
 *     for (const unsigned int q_point :
 *          scratch.fe_values.quadrature_point_indices())
 *       {
 *         const Tensor<2, dim>          tau = lqph[q_point]->get_tau();
 *         const SymmetricTensor<4, dim> Jc  = lqph[q_point]->get_Jc();
 *         const double d2Psi_vol_dJ2        = lqph[q_point]->get_d2Psi_vol_dJ2();
 *         const double det_F                = lqph[q_point]->get_det_F();
 *         const SymmetricTensor<2, dim> &I =
 *           Physics::Elasticity::StandardTensors<dim>::I;
 * 
 * @endcode
 * 
 * Next we define some aliases to make the assembly process easier to
 * follow
 * 
 * @code
 *         const std::vector<double> &                 N = scratch.Nx[q_point];
 *         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.JxW(q_point);
 * 
 *         for (const unsigned int i : scratch.fe_values.dof_indices())
 *           {
 *             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;
 * 
 *             for (const unsigned int j :
 *                  scratch.fe_values.dof_indices_ending_at(i))
 *               {
 *                 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}}_{uu}@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))
 *                   {
 * @endcode
 * 
 * The material contribution:
 * 
 * @code
 *                     data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc * 
 *                                               symm_grad_Nx[j] * JxW; 
 * 
 * @endcode
 * 
 * The geometrical stress contribution:
 * 
 * @code
 *                     if (component_i == component_j)
 *                       data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau *
 *                                                 grad_Nx[j][component_j] * JxW;
 *                   }
 * @endcode
 * 
 * Next is the @f$\mathsf{\mathbf{k}}_{ \widetilde{p} u}@f$
 * contribution
 * 
 * @code
 *                 else if ((i_group == p_dof) && (j_group == u_dof))
 *                   {
 *                     data.cell_matrix(i, j) += N[i] * det_F *               
 *                                               (symm_grad_Nx[j] * I) * JxW; 
 *                   }
 * @endcode
 * 
 * and lastly the @f$\mathsf{\mathbf{k}}_{ \widetilde{J}
 * \widetilde{p}}@f$ and @f$\mathsf{\mathbf{k}}_{ \widetilde{J}
 * \widetilde{J}}@f$ contributions:
 * 
 * @code
 *                 else if ((i_group == J_dof) && (j_group == p_dof))
 *                   data.cell_matrix(i, j) -= N[i] * N[j] * JxW;
 *                 else if ((i_group == j_group) && (i_group == J_dof))
 *                   data.cell_matrix(i, j) += N[i] * d2Psi_vol_dJ2 * N[j] * JxW;
 *                 else
 *                   Assert((i_group <= J_dof) && (j_group <= J_dof),
 *                          ExcInternalError());
 *               }
 *           }
 *       }
 * 
 * @endcode
 * 
 * Finally, we need to copy the lower half of the local matrix into the
 * upper half:
 * 
 * @code
 *     for (const unsigned int i : scratch.fe_values.dof_indices())
 *       for (const unsigned int j :
 *            scratch.fe_values.dof_indices_starting_at(i + 1))
 *         data.cell_matrix(i, j) = data.cell_matrix(j, i);
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="Solidassemble_system_rhs"></a> 
 * <h4>Solid::assemble_system_rhs</h4>
 * The assembly of the right-hand side process is similar to the
 * tangent matrix, so we will not describe it in too much detail.
 * Note that since we are describing a problem with Neumann BCs,
 * we will need the face normals and so must specify this in the
 * update flags.
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::assemble_system_rhs()
 *   {
 *     timer.enter_subsection("Assemble system right-hand side");
 *     std::cout << " ASM_R " << std::flush;
 * 
 *     system_rhs = 0.0;
 * 
 *     const UpdateFlags uf_cell(update_values | update_gradients |
 *                               update_JxW_values);
 *     const UpdateFlags uf_face(update_values | update_normal_vectors |
 *                               update_JxW_values);
 * 
 *     PerTaskData_RHS per_task_data(dofs_per_cell);
 *     ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face);
 * 
 *     WorkStream::run(
 *       dof_handler.active_cell_iterators(),
 *       [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *              ScratchData_RHS &                                     scratch,
 *              PerTaskData_RHS &                                     data) {
 *         this->assemble_system_rhs_one_cell(cell, scratch, data);
 *       },
 *       [this](const PerTaskData_RHS &data) {
 *         this->copy_local_to_global_rhs(data);
 *       },
 *       scratch_data,
 *       per_task_data);
 * 
 *     timer.leave_subsection();
 *   }
 * 
 * 
 * 
 *   template <int dim>
 *   void Solid<dim>::copy_local_to_global_rhs(const PerTaskData_RHS &data)
 *   {
 *     for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *       system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);
 *   }
 * 
 * 
 * 
 *   template <int dim>
 *   void Solid<dim>::assemble_system_rhs_one_cell(
 *     const typename DoFHandler<dim>::active_cell_iterator &cell,
 *     ScratchData_RHS &                                     scratch,
 *     PerTaskData_RHS &                                     data) const
 *   {
 *     data.reset();
 *     scratch.reset();
 *     scratch.fe_values.reinit(cell);
 *     cell->get_dof_indices(data.local_dof_indices);
 * 
 *     const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
 *       quadrature_point_history.get_data(cell);
 *     Assert(lqph.size() == n_q_points, ExcInternalError());
 * 
 *     for (const unsigned int q_point :
 *          scratch.fe_values.quadrature_point_indices())
 *       {
 *         const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
 * 
 *         for (const unsigned int k : scratch.fe_values.dof_indices())
 *           {
 *             const unsigned int k_group = fe.system_to_base_index(k).first.first;
 * 
 *             if (k_group == u_dof)
 *               scratch.symm_grad_Nx[q_point][k] = symmetrize(
 *                 scratch.fe_values[u_fe].gradient(k, q_point) * F_inv);
 *             else if (k_group == p_dof)
 *               scratch.Nx[q_point][k] =
 *                 scratch.fe_values[p_fe].value(k, q_point);
 *             else if (k_group == J_dof)
 *               scratch.Nx[q_point][k] =
 *                 scratch.fe_values[J_fe].value(k, q_point);
 *             else
 *               Assert(k_group <= J_dof, ExcInternalError());
 *           }
 *       }
 * 
 *     for (const unsigned int q_point :
 *          scratch.fe_values.quadrature_point_indices())
 *       {
 *         const SymmetricTensor<2, dim> tau     = lqph[q_point]->get_tau();
 *         const double                  det_F   = lqph[q_point]->get_det_F();
 *         const double                  J_tilde = lqph[q_point]->get_J_tilde();
 *         const double                  p_tilde = lqph[q_point]->get_p_tilde();
 *         const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ();
 * 
 *         const std::vector<double> &                 N = scratch.Nx[q_point];
 *         const std::vector<SymmetricTensor<2, dim>> &symm_grad_Nx =
 *           scratch.symm_grad_Nx[q_point];
 *         const double JxW = scratch.fe_values.JxW(q_point);
 * 
 * @endcode
 * 
 * We first compute the contributions
 * from the internal forces.  Note, by
 * definition of the rhs as the negative
 * of the residual, these contributions
 * are subtracted.
 * 
 * @code
 *         for (const unsigned int i : scratch.fe_values.dof_indices())
 *           {
 *             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 if (i_group == p_dof)
 *               data.cell_rhs(i) -= N[i] * (det_F - J_tilde) * JxW;
 *             else if (i_group == J_dof)
 *               data.cell_rhs(i) -= N[i] * (dPsi_vol_dJ - p_tilde) * JxW;
 *             else
 *               Assert(i_group <= J_dof, ExcInternalError());
 *           }
 *       }
 * 
 * @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 (const auto &face : cell->face_iterators())
 *       if (face->at_boundary() && face->boundary_id() == 6)
 *         {
 *           scratch.fe_face_values.reinit(cell, face);
 * 
 *           for (const unsigned int f_q_point :
 *                scratch.fe_face_values.quadrature_point_indices())
 *             {
 *               const Tensor<1, dim> &N =
 *                 scratch.fe_face_values.normal_vector(f_q_point);
 * 
 * @endcode
 * 
 * Using the face normal at this quadrature point we specify the
 * traction in reference configuration. For this problem, a
 * defined pressure is applied in the reference configuration.
 * The direction of the applied traction is assumed not to
 * evolve with the deformation of the domain. The traction is
 * defined using the first Piola-Kirchhoff stress is simply
 * @f$\mathbf{t} = \mathbf{P}\mathbf{N} = [p_0 \mathbf{I}]
 * \mathbf{N} = p_0 \mathbf{N}@f$ We use the time variable to
 * linearly ramp up the pressure load.
 *               

 * 
 * Note that the contributions to the right hand side vector we
 * compute here only exist in the displacement components of the
 * vector.
 * 
 * @code
 *               static const double p0 =
 *                 -4.0 / (parameters.scale * parameters.scale);
 *               const double         time_ramp = (time.current() / time.end());
 *               const double         pressure  = p0 * parameters.p_p0 * time_ramp;
 *               const Tensor<1, dim> traction  = pressure * N;
 * 
 *               for (const unsigned int i : scratch.fe_values.dof_indices())
 *                 {
 *                   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.shape_value(i, f_q_point);
 *                       const double JxW = scratch.fe_face_values.JxW(f_q_point);
 * 
 *                       data.cell_rhs(i) += (Ni * traction[component_i]) * JxW;
 *                     }
 *                 }
 *             }
 *         }
 *   }
 * 
 * @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>
 *   void Solid<dim>::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;
 *     constraints.clear();
 *     const bool apply_dirichlet_bc = (it_nr == 0);
 * 
 * @endcode
 * 
 * In this particular example, the boundary values will be calculated for
 * the two first iterations of Newton's algorithm. In general, one would
 * build non-homogeneous constraints in the zeroth iteration (that is, when
 * `apply_dirichlet_bc == true`) and build only the corresponding
 * homogeneous constraints in the following step. While the current
 * example has only homogeneous constraints, previous experiences have
 * shown that a common error is forgetting to add the extra condition when
 * refactoring the code to specific uses. This could lead errors that are
 * hard to debug. In this spirit, we choose to make the code more verbose
 * in terms of what operations are performed at each Newton step.
 *     

 * 
 * The boundary conditions for the indentation problem are as follows: On
 * the -x, -y and -z faces (IDs 0,2,4) we set up a symmetry condition to
 * allow only planar movement while the +x and +z faces (IDs 1,5) are
 * traction free. In this contrived problem, part of the +y face (ID 3) is
 * set to have no motion in the x- and z-component. Finally, as described
 * earlier, the other part of the +y face has an the applied pressure but
 * is also constrained in the x- and z-directions.
 *     

 * 
 * 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
 *     const FEValuesExtractors::Scalar x_displacement(0);
 *     const FEValuesExtractors::Scalar y_displacement(1);
 * 
 *     {
 *       const int boundary_id = 0;
 * 
 *       if (apply_dirichlet_bc == true)
 *         VectorTools::interpolate_boundary_values(
 *           dof_handler,
 *           boundary_id,
 *           Functions::ZeroFunction<dim>(n_components),
 *           constraints,
 *           fe.component_mask(x_displacement));
 *       else
 *         VectorTools::interpolate_boundary_values(
 *           dof_handler,
 *           boundary_id,
 *           Functions::ZeroFunction<dim>(n_components),
 *           constraints,
 *           fe.component_mask(x_displacement));
 *     }
 *     {
 *       const int boundary_id = 2;
 * 
 *       if (apply_dirichlet_bc == true)
 *         VectorTools::interpolate_boundary_values(
 *           dof_handler,
 *           boundary_id,
 *           Functions::ZeroFunction<dim>(n_components),
 *           constraints,
 *           fe.component_mask(y_displacement));
 *       else
 *         VectorTools::interpolate_boundary_values(
 *           dof_handler,
 *           boundary_id,
 *           Functions::ZeroFunction<dim>(n_components),
 *           constraints,
 *           fe.component_mask(y_displacement));
 *     }
 * 
 *     if (dim == 3)
 *       {
 *         const FEValuesExtractors::Scalar z_displacement(2);
 * 
 *         {
 *           const int boundary_id = 3;
 * 
 *           if (apply_dirichlet_bc == true)
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               (fe.component_mask(x_displacement) |
 *                fe.component_mask(z_displacement)));
 *           else
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               (fe.component_mask(x_displacement) |
 *                fe.component_mask(z_displacement)));
 *         }
 *         {
 *           const int boundary_id = 4;
 * 
 *           if (apply_dirichlet_bc == true)
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               fe.component_mask(z_displacement));
 *           else
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               fe.component_mask(z_displacement));
 *         }
 * 
 *         {
 *           const int boundary_id = 6;
 * 
 *           if (apply_dirichlet_bc == true)
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               (fe.component_mask(x_displacement) |
 *                fe.component_mask(z_displacement)));
 *           else
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               (fe.component_mask(x_displacement) |
 *                fe.component_mask(z_displacement)));
 *         }
 *       }
 *     else
 *       {
 *         {
 *           const int boundary_id = 3;
 * 
 *           if (apply_dirichlet_bc == true)
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               (fe.component_mask(x_displacement)));
 *           else
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               (fe.component_mask(x_displacement)));
 *         }
 *         {
 *           const int boundary_id = 6;
 * 
 *           if (apply_dirichlet_bc == true)
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               (fe.component_mask(x_displacement)));
 *           else
 *             VectorTools::interpolate_boundary_values(
 *               dof_handler,
 *               boundary_id,
 *               Functions::ZeroFunction<dim>(n_components),
 *               constraints,
 *               (fe.component_mask(x_displacement)));
 *         }
 *       }
 * 
 *     constraints.close();
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="Solidassemble_sc"></a> 
 * <h4>Solid::assemble_sc</h4>
 * Solving the entire block system is a bit problematic as there are no
 * contributions to the @f$\mathsf{\mathbf{K}}_{ \widetilde{J} \widetilde{J}}@f$
 * block, rendering it noninvertible (when using an iterative solver).
 * Since the pressure and dilatation variables DOFs are discontinuous, we can
 * condense them out to form a smaller displacement-only system which
 * we will then solve and subsequently post-process to retrieve the
 * pressure and dilatation solutions.
 * 

 * 
 * The static condensation process could be performed at a global level but we
 * need the inverse of one of the blocks. However, since the pressure and
 * dilatation variables are discontinuous, the static condensation (SC)
 * operation can also be done on a per-cell basis and we can produce the
 * inverse of the block-diagonal
 * @f$\mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}@f$ block by inverting the
 * local blocks. We can again use TBB to do this since each operation will be
 * independent of one another.
 *   

 * 
 * Using the TBB via the WorkStream class, we assemble the contributions to
 * form
 * @f$
 * \mathsf{\mathbf{K}}_{\textrm{con}}
 * = \bigl[ \mathsf{\mathbf{K}}_{uu} +
 * \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]
 * @f$
 * from each element's contributions. These contributions are then added to
 * the global stiffness matrix. Given this description, the following two
 * functions should be clear:
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::assemble_sc()
 *   {
 *     timer.enter_subsection("Perform static condensation");
 *     std::cout << " ASM_SC " << std::flush;
 * 
 *     PerTaskData_SC per_task_data(dofs_per_cell,
 *                                  element_indices_u.size(),
 *                                  element_indices_p.size(),
 *                                  element_indices_J.size());
 *     ScratchData_SC scratch_data;
 * 
 *     WorkStream::run(dof_handler.active_cell_iterators(),
 *                     *this,
 *                     &Solid::assemble_sc_one_cell,
 *                     &Solid::copy_local_to_global_sc,
 *                     scratch_data,
 *                     per_task_data);
 * 
 *     timer.leave_subsection();
 *   }
 * 
 * 
 *   template <int dim>
 *   void Solid<dim>::copy_local_to_global_sc(const PerTaskData_SC &data)
 *   {
 *     for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *       for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *         tangent_matrix.add(data.local_dof_indices[i],
 *                            data.local_dof_indices[j],
 *                            data.cell_matrix(i, j));
 *   }
 * 
 * 
 * @endcode
 * 
 * Now we describe the static condensation process. As per usual, we must
 * first find out which global numbers the degrees of freedom on this cell
 * have and reset some data structures:
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::assemble_sc_one_cell(
 *     const typename DoFHandler<dim>::active_cell_iterator &cell,
 *     ScratchData_SC &                                      scratch,
 *     PerTaskData_SC &                                      data)
 *   {
 *     data.reset();
 *     scratch.reset();
 *     cell->get_dof_indices(data.local_dof_indices);
 * 
 * @endcode
 * 
 * We now extract the contribution of the dofs associated with the current
 * cell to the global stiffness matrix.  The discontinuous nature of the
 * @f$\widetilde{p}@f$ and @f$\widetilde{J}@f$ interpolations mean that their is
 * no coupling of the local contributions at the global level. This is not
 * the case with the @f$\mathbf{u}@f$ dof.  In other words,
 * @f$\mathsf{\mathbf{k}}_{\widetilde{J} \widetilde{p}}@f$,
 * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{p}}@f$ and
 * @f$\mathsf{\mathbf{k}}_{\widetilde{J} \widetilde{p}}@f$, when extracted
 * from the global stiffness matrix are the element contributions.  This
 * is not the case for @f$\mathsf{\mathbf{k}}_{uu}@f$.
 *     

 * 
 * Note: A lower-case symbol is used to denote element stiffness matrices.
 * 

 * 
 * Currently the matrix corresponding to
 * the dof associated with the current element
 * (denoted somewhat loosely as @f$\mathsf{\mathbf{k}}@f$)
 * is of the form:
 * @f{align*}
 * \begin{bmatrix}
 * \mathsf{\mathbf{k}}_{uu}  &  \mathsf{\mathbf{k}}_{u\widetilde{p}}
 * & \mathbf{0}
 * \\ \mathsf{\mathbf{k}}_{\widetilde{p}u} & \mathbf{0}  &
 * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}
 * \\ \mathbf{0}  &  \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}}  &
 * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}} \end{bmatrix}
 * @f}
 *     

 * 
 * We now need to modify it such that it appear as
 * @f{align*}
 * \begin{bmatrix}
 * \mathsf{\mathbf{k}}_{\textrm{con}}   &
 * \mathsf{\mathbf{k}}_{u\widetilde{p}}    & \mathbf{0}
 * \\ \mathsf{\mathbf{k}}_{\widetilde{p}u} & \mathbf{0} &
 * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}^{-1}
 * \\ \mathbf{0} & \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}} &
 * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}} \end{bmatrix}
 * @f}
 * with @f$\mathsf{\mathbf{k}}_{\textrm{con}} = \bigl[
 * \mathsf{\mathbf{k}}_{uu} +\overline{\overline{\mathsf{\mathbf{k}}}}~
 * \bigr]@f$ where @f$               \overline{\overline{\mathsf{\mathbf{k}}}}
 * \dealcoloneq \mathsf{\mathbf{k}}_{u\widetilde{p}}
 * \overline{\mathsf{\mathbf{k}}} \mathsf{\mathbf{k}}_{\widetilde{p}u}
 * @f$
 * and
 * @f$
 * \overline{\mathsf{\mathbf{k}}} =
 * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}}^{-1}
 * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}}
 * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}^{-1}
 * @f$.
 *     

 * 
 * At this point, we need to take note of
 * the fact that global data already exists
 * in the @f$\mathsf{\mathbf{K}}_{uu}@f$,
 * @f$\mathsf{\mathbf{K}}_{\widetilde{p} \widetilde{J}}@f$
 * and
 * @f$\mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{p}}@f$
 * sub-blocks.  So if we are to modify them, we must account for the data
 * that is already there (i.e. simply add to it or remove it if
 * necessary).  Since the copy_local_to_global operation is a "+="
 * operation, we need to take this into account
 *     

 * 
 * For the @f$\mathsf{\mathbf{K}}_{uu}@f$ block in particular, this means that
 * contributions have been added from the surrounding cells, so we need to
 * be careful when we manipulate this block.  We can't just erase the
 * sub-blocks.
 *     

 * 
 * This is the strategy we will employ to get the sub-blocks we want:
 *     

 * 
 * - @f$ {\mathsf{\mathbf{k}}}_{\textrm{store}}@f$:
 * Since we don't have access to @f$\mathsf{\mathbf{k}}_{uu}@f$,
 * but we know its contribution is added to
 * the global @f$\mathsf{\mathbf{K}}_{uu}@f$ matrix, we just want
 * to add the element wise
 * static-condensation @f$\overline{\overline{\mathsf{\mathbf{k}}}}@f$.
 *     

 * 
 * - @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}@f$:
 * Similarly, @f$\mathsf{\mathbf{k}}_{\widetilde{p}
 * \widetilde{J}}@f$ exists in
 * the subblock. Since the copy
 * operation is a += operation, we
 * need to subtract the existing
 * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$
 * submatrix in addition to
 * "adding" that which we wish to
 * replace it with.
 *     

 * 
 * - @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}@f$:
 * Since the global matrix
 * is symmetric, this block is the
 * same as the one above and we
 * can simply use
 * @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}@f$
 * as a substitute for this one.
 *     

 * 
 * We first extract element data from the
 * system matrix. So first we get the
 * entire subblock for the cell, then
 * extract @f$\mathsf{\mathbf{k}}@f$
 * for the dofs associated with
 * the current element
 * 
 * @code
 *     data.k_orig.extract_submatrix_from(tangent_matrix,
 *                                        data.local_dof_indices,
 *                                        data.local_dof_indices);
 * @endcode
 * 
 * and next the local matrices for
 * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} u}@f$
 * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} \widetilde{J}}@f$
 * and
 * @f$\mathsf{\mathbf{k}}_{ \widetilde{J} \widetilde{J}}@f$:
 * 
 * @code
 *     data.k_pu.extract_submatrix_from(data.k_orig,
 *                                      element_indices_p,
 *                                      element_indices_u);
 *     data.k_pJ.extract_submatrix_from(data.k_orig,
 *                                      element_indices_p,
 *                                      element_indices_J);
 *     data.k_JJ.extract_submatrix_from(data.k_orig,
 *                                      element_indices_J,
 *                                      element_indices_J);
 * 
 * @endcode
 * 
 * To get the inverse of @f$\mathsf{\mathbf{k}}_{\widetilde{p}
 * \widetilde{J}}@f$, we invert it directly.  This operation is relatively
 * inexpensive since @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$
 * since block-diagonal.
 * 
 * @code
 *     data.k_pJ_inv.invert(data.k_pJ);
 * 
 * @endcode
 * 
 * Now we can make condensation terms to
 * add to the @f$\mathsf{\mathbf{k}}_{uu}@f$
 * block and put them in
 * the cell local matrix
 * @f$
 * \mathsf{\mathbf{A}}
 * =
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{k}}_{\widetilde{p} u}
 * @f$:
 * 
 * @code
 *     data.k_pJ_inv.mmult(data.A, data.k_pu);
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{B}}
 * =
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{k}}_{\widetilde{p} u}
 * @f$
 * 
 * @code
 *     data.k_JJ.mmult(data.B, data.A);
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{C}}
 * =
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{k}}_{\widetilde{p} u}
 * @f$
 * 
 * @code
 *     data.k_pJ_inv.Tmmult(data.C, data.B);
 * @endcode
 * 
 * @f$
 * \overline{\overline{\mathsf{\mathbf{k}}}}
 * =
 * \mathsf{\mathbf{k}}_{u \widetilde{p}}
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{k}}_{\widetilde{p} u}
 * @f$
 * 
 * @code
 *     data.k_pu.Tmmult(data.k_bbar, data.C);
 *     data.k_bbar.scatter_matrix_to(element_indices_u,
 *                                   element_indices_u,
 *                                   data.cell_matrix);
 * 
 * @endcode
 * 
 * Next we place
 * @f$\mathsf{\mathbf{k}}^{-1}_{ \widetilde{p} \widetilde{J}}@f$
 * in the
 * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} \widetilde{J}}@f$
 * block for post-processing.  Note again
 * that we need to remove the
 * contribution that already exists there.
 * 
 * @code
 *     data.k_pJ_inv.add(-1.0, data.k_pJ);
 *     data.k_pJ_inv.scatter_matrix_to(element_indices_p,
 *                                     element_indices_J,
 *                                     data.cell_matrix);
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="Solidsolve_linear_system"></a> 
 * <h4>Solid::solve_linear_system</h4>
 * We now have all of the necessary components to use one of two possible
 * methods to solve the linearised system. The first is to perform static
 * condensation on an element level, which requires some alterations
 * to the tangent matrix and RHS vector. Alternatively, the full block
 * system can be solved by performing condensation on a global level.
 * Below we implement both approaches.
 * 
 * @code
 *   template <int dim>
 *   std::pair<unsigned int, double>
 *   Solid<dim>::solve_linear_system(BlockVector<double> &newton_update)
 *   {
 *     unsigned int lin_it  = 0;
 *     double       lin_res = 0.0;
 * 
 *     if (parameters.use_static_condensation == true)
 *       {
 * @endcode
 * 
 * Firstly, here is the approach using the (permanent) augmentation of
 * the tangent matrix. For the following, recall that
 * @f{align*}
 * \mathsf{\mathbf{K}}_{\textrm{store}}
 * \dealcoloneq
 * \begin{bmatrix}
 * \mathsf{\mathbf{K}}_{\textrm{con}}      &
 * \mathsf{\mathbf{K}}_{u\widetilde{p}}    & \mathbf{0}
 * \\  \mathsf{\mathbf{K}}_{\widetilde{p}u}    &       \mathbf{0} &
 * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}
 * \\  \mathbf{0}      &
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}                &
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}} \end{bmatrix} \, .
 * @f}
 * and
 * @f{align*}
 * d \widetilde{\mathsf{\mathbf{p}}}
 * & =
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
 * \bigl[
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * -
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
 * d \widetilde{\mathsf{\mathbf{J}}} \bigr]
 * \\ d \widetilde{\mathsf{\mathbf{J}}}
 * & =
 * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}
 * \bigl[
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * - \mathsf{\mathbf{K}}_{\widetilde{p}u} d
 * \mathsf{\mathbf{u}} \bigr]
 * \\ \Rightarrow d \widetilde{\mathsf{\mathbf{p}}}
 * &= \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * -
 * \underbrace{\bigl[\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
 * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}\bigr]}_{\overline{\mathsf{\mathbf{K}}}}\bigl[
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * - \mathsf{\mathbf{K}}_{\widetilde{p}u} d
 * \mathsf{\mathbf{u}} \bigr]
 * @f}
 * and thus
 * @f[
 * \underbrace{\bigl[ \mathsf{\mathbf{K}}_{uu} +
 * \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]
 * }_{\mathsf{\mathbf{K}}_{\textrm{con}}} d
 * \mathsf{\mathbf{u}}
 * =
 * \underbrace{
 * \Bigl[
 * \mathsf{\mathbf{F}}_{u}
 * - \mathsf{\mathbf{K}}_{u\widetilde{p}} \bigl[
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * -
 * \overline{\mathsf{\mathbf{K}}}\mathsf{\mathbf{F}}_{\widetilde{p}}
 * \bigr]
 * \Bigr]}_{\mathsf{\mathbf{F}}_{\textrm{con}}}
 * @f]
 * where
 * @f[
 * \overline{\overline{\mathsf{\mathbf{K}}}} \dealcoloneq
 * \mathsf{\mathbf{K}}_{u\widetilde{p}}
 * \overline{\mathsf{\mathbf{K}}}
 * \mathsf{\mathbf{K}}_{\widetilde{p}u} \, .
 * @f]
 * 

 * 
 * At the top, we allocate two temporary vectors to help with the
 * static condensation, and variables to store the number of
 * linear solver iterations and the (hopefully converged) residual.
 * 
 * @code
 *         BlockVector<double> A(dofs_per_block);
 *         BlockVector<double> B(dofs_per_block);
 * 
 * 
 * @endcode
 * 
 * In the first step of this function, we solve for the incremental
 * displacement @f$d\mathbf{u}@f$.  To this end, we perform static
 * condensation to make
 * @f$\mathsf{\mathbf{K}}_{\textrm{con}}
 * = \bigl[ \mathsf{\mathbf{K}}_{uu} +
 * \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]@f$
 * and put
 * @f$\mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}@f$
 * in the original @f$\mathsf{\mathbf{K}}_{\widetilde{p} \widetilde{J}}@f$
 * block. That is, we make @f$\mathsf{\mathbf{K}}_{\textrm{store}}@f$.
 * 
 * @code
 *         {
 *           assemble_sc();
 * 
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{J}}
 * =
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * @f$
 * 
 * @code
 *           tangent_matrix.block(p_dof, J_dof)
 *             .vmult(A.block(J_dof), system_rhs.block(p_dof));
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{B}}_{\widetilde{J}}
 * =
 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * @f$
 * 
 * @code
 *           tangent_matrix.block(J_dof, J_dof)
 *             .vmult(B.block(J_dof), A.block(J_dof));
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{J}}
 * =
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * -
 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * @f$
 * 
 * @code
 *           A.block(J_dof) = system_rhs.block(J_dof);
 *           A.block(J_dof) -= B.block(J_dof);
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{J}}
 * =
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
 * [
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * -
 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * ]
 * @f$
 * 
 * @code
 *           tangent_matrix.block(p_dof, J_dof)
 *             .Tvmult(A.block(p_dof), A.block(J_dof));
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{u}
 * =
 * \mathsf{\mathbf{K}}_{u \widetilde{p}}
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
 * [
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * -
 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * ]
 * @f$
 * 
 * @code
 *           tangent_matrix.block(u_dof, p_dof)
 *             .vmult(A.block(u_dof), A.block(p_dof));
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{F}}_{\text{con}}
 * =
 * \mathsf{\mathbf{F}}_{u}
 * -
 * \mathsf{\mathbf{K}}_{u \widetilde{p}}
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
 * [
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * -
 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * ]
 * @f$
 * 
 * @code
 *           system_rhs.block(u_dof) -= A.block(u_dof);
 * 
 *           timer.enter_subsection("Linear solver");
 *           std::cout << " SLV " << std::flush;
 *           if (parameters.type_lin == "CG")
 *             {
 *               const auto 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, this might not be
 * true for different problem sizes.
 * 
 * @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);
 * 
 *         timer.enter_subsection("Linear solver postprocessing");
 *         std::cout << " PP " << std::flush;
 * 
 * @endcode
 * 
 * The next step after solving the displacement
 * problem is to post-process to get the
 * dilatation solution from the
 * substitution:
 * @f$
 * d \widetilde{\mathsf{\mathbf{J}}}
 * = \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1} \bigl[
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * - \mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
 * \bigr]
 * @f$
 * 
 * @code
 *         {
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{p}}
 * =
 * \mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
 * @f$
 * 
 * @code
 *           tangent_matrix.block(p_dof, u_dof)
 *             .vmult(A.block(p_dof), newton_update.block(u_dof));
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{p}}
 * =
 * -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
 * @f$
 * 
 * @code
 *           A.block(p_dof) *= -1.0;
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{p}}
 * =
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
 * @f$
 * 
 * @code
 *           A.block(p_dof) += system_rhs.block(p_dof);
 * @endcode
 * 
 * @f$
 * d\mathsf{\mathbf{\widetilde{J}}}
 * =
 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p}\widetilde{J}}
 * [
 * \mathsf{\mathbf{F}}_{\widetilde{p}}
 * -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
 * ]
 * @f$
 * 
 * @code
 *           tangent_matrix.block(p_dof, J_dof)
 *             .vmult(newton_update.block(J_dof), A.block(p_dof));
 *         }
 * 
 * @endcode
 * 
 * we ensure here that any Dirichlet constraints
 * are distributed on the updated solution:
 * 
 * @code
 *         constraints.distribute(newton_update);
 * 
 * @endcode
 * 
 * Finally we solve for the pressure
 * update with the substitution:
 * @f$
 * d \widetilde{\mathsf{\mathbf{p}}}
 * =
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
 * \bigl[
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * - \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
 * d \widetilde{\mathsf{\mathbf{J}}}
 * \bigr]
 * @f$
 * 
 * @code
 *         {
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{J}}
 * =
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
 * d \widetilde{\mathsf{\mathbf{J}}}
 * @f$
 * 
 * @code
 *           tangent_matrix.block(J_dof, J_dof)
 *             .vmult(A.block(J_dof), newton_update.block(J_dof));
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{J}}
 * =
 * -\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
 * d \widetilde{\mathsf{\mathbf{J}}}
 * @f$
 * 
 * @code
 *           A.block(J_dof) *= -1.0;
 * @endcode
 * 
 * @f$
 * \mathsf{\mathbf{A}}_{\widetilde{J}}
 * =
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * -
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
 * d \widetilde{\mathsf{\mathbf{J}}}
 * @f$
 * 
 * @code
 *           A.block(J_dof) += system_rhs.block(J_dof);
 * @endcode
 * 
 * and finally....
 * @f$
 * d \widetilde{\mathsf{\mathbf{p}}}
 * =
 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
 * \bigl[
 * \mathsf{\mathbf{F}}_{\widetilde{J}}
 * - \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
 * d \widetilde{\mathsf{\mathbf{J}}}
 * \bigr]
 * @f$
 * 
 * @code
 *           tangent_matrix.block(p_dof, J_dof)
 *             .Tvmult(newton_update.block(p_dof), A.block(J_dof));
 *         }
 * 
 * @endcode
 * 
 * We are now at the end, so we distribute all
 * constrained dofs back to the Newton
 * update:
 * 
 * @code
 *         constraints.distribute(newton_update);
 * 
 *         timer.leave_subsection();
 *       }
 *     else
 *       {
 *         std::cout << " ------ " << std::flush;
 * 
 *         timer.enter_subsection("Linear solver");
 *         std::cout << " SLV " << std::flush;
 * 
 *         if (parameters.type_lin == "CG")
 *           {
 * @endcode
 * 
 * Manual condensation of the dilatation and pressure fields on
 * a local level, and subsequent post-processing, took quite a
 * bit of effort to achieve. To recap, we had to produce the
 * inverse matrix
 * @f$\mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}@f$, which
 * was permanently written into the global tangent matrix. We then
 * permanently modified @f$\mathsf{\mathbf{K}}_{uu}@f$ to produce
 * @f$\mathsf{\mathbf{K}}_{\textrm{con}}@f$. This involved the
 * extraction and manipulation of local sub-blocks of the tangent
 * matrix. After solving for the displacement, the individual
 * matrix-vector operations required to solve for dilatation and
 * pressure were carefully implemented. Contrast these many sequence
 * of steps to the much simpler and transparent implementation using
 * functionality provided by the LinearOperator class.
 * 

 * 
 * For ease of later use, we define some aliases for
 * blocks in the RHS vector
 * 
 * @code
 *             const Vector<double> &f_u = system_rhs.block(u_dof);
 *             const Vector<double> &f_p = system_rhs.block(p_dof);
 *             const Vector<double> &f_J = system_rhs.block(J_dof);
 * 
 * @endcode
 * 
 * ... and for blocks in the Newton update vector.
 * 
 * @code
 *             Vector<double> &d_u = newton_update.block(u_dof);
 *             Vector<double> &d_p = newton_update.block(p_dof);
 *             Vector<double> &d_J = newton_update.block(J_dof);
 * 
 * @endcode
 * 
 * We next define some linear operators for the tangent matrix
 * sub-blocks We will exploit the symmetry of the system, so not all
 * blocks are required.
 * 
 * @code
 *             const auto K_uu =
 *               linear_operator(tangent_matrix.block(u_dof, u_dof));
 *             const auto K_up =
 *               linear_operator(tangent_matrix.block(u_dof, p_dof));
 *             const auto K_pu =
 *               linear_operator(tangent_matrix.block(p_dof, u_dof));
 *             const auto K_Jp =
 *               linear_operator(tangent_matrix.block(J_dof, p_dof));
 *             const auto K_JJ =
 *               linear_operator(tangent_matrix.block(J_dof, J_dof));
 * 
 * @endcode
 * 
 * We then construct a LinearOperator that represents the inverse of
 * (square block)
 * @f$\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}@f$. Since it is
 * diagonal (or, when a higher order ansatz it used, nearly
 * diagonal), a Jacobi preconditioner is suitable.
 * 
 * @code
 *             PreconditionSelector<SparseMatrix<double>, Vector<double>>
 *               preconditioner_K_Jp_inv("jacobi");
 *             preconditioner_K_Jp_inv.use_matrix(
 *               tangent_matrix.block(J_dof, p_dof));
 *             ReductionControl solver_control_K_Jp_inv(
 *               static_cast<unsigned int>(tangent_matrix.block(J_dof, p_dof).m() *
 *                                         parameters.max_iterations_lin),
 *               1.0e-30,
 *               parameters.tol_lin);
 *             SolverSelector<Vector<double>> solver_K_Jp_inv;
 *             solver_K_Jp_inv.select("cg");
 *             solver_K_Jp_inv.set_control(solver_control_K_Jp_inv);
 *             const auto K_Jp_inv =
 *               inverse_operator(K_Jp, solver_K_Jp_inv, preconditioner_K_Jp_inv);
 * 
 * @endcode
 * 
 * Now we can construct that transpose of
 * @f$\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}@f$ and a
 * linear operator that represents the condensed operations
 * @f$\overline{\mathsf{\mathbf{K}}}@f$ and
 * @f$\overline{\overline{\mathsf{\mathbf{K}}}}@f$ and the final
 * augmented matrix
 * @f$\mathsf{\mathbf{K}}_{\textrm{con}}@f$.
 * Note that the schur_complement() operator could also be of use
 * here, but for clarity and the purpose of demonstrating the
 * similarities between the formulation and implementation of the
 * linear solution scheme, we will perform these operations
 * manually.
 * 
 * @code
 *             const auto K_pJ_inv     = transpose_operator(K_Jp_inv);
 *             const auto K_pp_bar     = K_Jp_inv * K_JJ * K_pJ_inv;
 *             const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu;
 *             const auto K_uu_con     = K_uu + K_uu_bar_bar;
 * 
 * @endcode
 * 
 * Lastly, we define an operator for inverse of augmented stiffness
 * matrix, namely @f$\mathsf{\mathbf{K}}_{\textrm{con}}^{-1}@f$. Note
 * that the preconditioner for the augmented stiffness matrix is
 * different to the case when we use static condensation. In this
 * instance, the preconditioner is based on a non-modified
 * @f$\mathsf{\mathbf{K}}_{uu}@f$, while with the first approach we
 * actually modified the entries of this sub-block. However, since
 * @f$\mathsf{\mathbf{K}}_{\textrm{con}}@f$ and
 * @f$\mathsf{\mathbf{K}}_{uu}@f$ operate on the same space, it remains
 * adequate for this problem.
 * 
 * @code
 *             PreconditionSelector<SparseMatrix<double>, Vector<double>>
 *               preconditioner_K_con_inv(parameters.preconditioner_type,
 *                                        parameters.preconditioner_relaxation);
 *             preconditioner_K_con_inv.use_matrix(
 *               tangent_matrix.block(u_dof, u_dof));
 *             ReductionControl solver_control_K_con_inv(
 *               static_cast<unsigned int>(tangent_matrix.block(u_dof, u_dof).m() *
 *                                         parameters.max_iterations_lin),
 *               1.0e-30,
 *               parameters.tol_lin);
 *             SolverSelector<Vector<double>> solver_K_con_inv;
 *             solver_K_con_inv.select("cg");
 *             solver_K_con_inv.set_control(solver_control_K_con_inv);
 *             const auto K_uu_con_inv =
 *               inverse_operator(K_uu_con,
 *                                solver_K_con_inv,
 *                                preconditioner_K_con_inv);
 * 
 * @endcode
 * 
 * Now we are in a position to solve for the displacement field.
 * We can nest the linear operations, and the result is immediately
 * written to the Newton update vector.
 * It is clear that the implementation closely mimics the derivation
 * stated in the introduction.
 * 
 * @code
 *             d_u =
 *               K_uu_con_inv * (f_u - K_up * (K_Jp_inv * f_J - K_pp_bar * f_p));
 * 
 *             timer.leave_subsection();
 * 
 * @endcode
 * 
 * The operations need to post-process for the dilatation and
 * pressure fields are just as easy to express.
 * 
 * @code
 *             timer.enter_subsection("Linear solver postprocessing");
 *             std::cout << " PP " << std::flush;
 * 
 *             d_J = K_pJ_inv * (f_p - K_pu * d_u);
 *             d_p = K_Jp_inv * (f_J - K_JJ * d_J);
 * 
 *             lin_it  = solver_control_K_con_inv.last_step();
 *             lin_res = solver_control_K_con_inv.last_value();
 *           }
 *         else if (parameters.type_lin == "Direct")
 *           {
 * @endcode
 * 
 * Solve the full block system with
 * a direct solver. As it is relatively
 * robust, it may be immune to problem
 * arising from the presence of the zero
 * @f$\mathsf{\mathbf{K}}_{ \widetilde{J} \widetilde{J}}@f$
 * block.
 * 
 * @code
 *             SparseDirectUMFPACK A_direct;
 *             A_direct.initialize(tangent_matrix);
 *             A_direct.vmult(newton_update, system_rhs);
 * 
 *             lin_it  = 1;
 *             lin_res = 0.0;
 * 
 *             std::cout << " -- " << std::flush;
 *           }
 *         else
 *           Assert(false, ExcMessage("Linear solver type not implemented"));
 * 
 *         timer.leave_subsection();
 * 
 * @endcode
 * 
 * Finally, we again ensure here that any Dirichlet
 * constraints are distributed on the updated solution:
 * 
 * @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 previous
 * tutorials so will not be discussed in detail.
 * 
 * @code
 *   template <int dim>
 *   void Solid<dim>::output_results() const
 *   {
 *     DataOut<dim> data_out;
 *     std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *       data_component_interpretation(
 *         dim, DataComponentInterpretation::component_is_part_of_vector);
 *     data_component_interpretation.push_back(
 *       DataComponentInterpretation::component_is_scalar);
 *     data_component_interpretation.push_back(
 *       DataComponentInterpretation::component_is_scalar);
 * 
 *     std::vector<std::string> solution_name(dim, "displacement");
 *     solution_name.emplace_back("pressure");
 *     solution_name.emplace_back("dilatation");
 * 
 *     DataOutBase::VtkFlags output_flags;
 *     output_flags.write_higher_order_cells = true;
 *     data_out.set_flags(output_flags);
 * 
 *     data_out.attach_dof_handler(dof_handler);
 *     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, soln);
 *     data_out.build_patches(q_mapping, degree);
 * 
 *     std::ofstream output("solution-" + std::to_string(dim) + "d-" +
 *                          std::to_string(time.get_timestep()) + ".vtu");
 *     data_out.write_vtu(output);
 *   }
 * 
 * } // namespace Step44
 * 
 * 
 * @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()
 * {
 *   using namespace Step44;
 * 
 *   try
 *     {
 *       const unsigned int dim = 3;
 *       Solid<dim>         solid("parameters.prm");
 *       solid.run();
 *     }
 *   catch (std::exception &exc)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Exception on processing: " << std::endl
 *                 << exc.what() << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 * 
 *       return 1;
 *     }
 *   catch (...)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Unknown exception!" << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       return 1;
 *     }
 * 
 *   return 0;
 * }
 * @endcode
<a name="Results"></a><h1>Results</h1>


Firstly, we present a comparison of a series of 3-d results with those
in the literature (see Reese et al (2000)) to demonstrate that the program works as expected.

We begin with a comparison of the convergence with mesh refinement for the @f$Q_1-DGPM_0-DGPM_0@f$ and
@f$Q_2-DGPM_1-DGPM_1@f$ formulations, as summarised in the figure below.
The vertical displacement of the midpoint of the upper surface of the block is used to assess convergence.
Both schemes demonstrate good convergence properties for varying values of the load parameter @f$p/p_0@f$.
The results agree with those in the literature.
The lower-order formulation typically overestimates the displacement for low levels of refinement,
while the higher-order interpolation scheme underestimates it, but be a lesser degree.
This benchmark, and a series of others not shown here, give us confidence that the code is working
as it should.

<table align="center" class="tutorial" cellspacing="3" cellpadding="3">
  <tr>
     <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Q1-P0_convergence.png" alt="">
        <p align="center">
        Convergence of the @f$Q_1-DGPM_0-DGPM_0@f$ formulation in 3-d.
        </p>
    </td>
    <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Q2-P1_convergence.png" alt="">
        <p align="center">
        Convergence of the @f$Q_2-DGPM_1-DGPM_1@f$ formulation in 3-d.
        </p>
    </td>
  </tr>
</table>


A typical screen output generated by running the problem is shown below.
The particular case demonstrated is that of the @f$Q_2-DGPM_1-DGPM_1@f$ formulation.
It is clear that, using the Newton-Raphson method, quadratic convergence of the solution is obtained.
Solution convergence is achieved within 5 Newton increments for all time-steps.
The converged displacement's @f$L_2@f$-norm is several orders of magnitude less than the geometry scale.

@code
Grid:
         Reference volume: 1e-09
Triangulation:
         Number of active cells: 64
         Number of degrees of freedom: 2699
    Setting up quadrature point data...

Timestep 1 @ 0.1s
___________________________________________________________________________________________________________________________________________________________
                 SOLVER STEP                   |  LIN_IT   LIN_RES    RES_NORM     RES_U     RES_P      RES_J     NU_NORM      NU_U       NU_P       NU_J
___________________________________________________________________________________________________________________________________________________________
  0  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     786  2.118e-06  1.000e+00  1.000e+00  0.000e+00  0.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00
  1  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     552  1.031e-03  8.563e-02  8.563e-02  9.200e-13  3.929e-08  1.060e-01  3.816e-02  1.060e-01  1.060e-01
  2  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     667  5.602e-06  2.482e-03  2.482e-03  3.373e-15  2.982e-10  2.936e-03  2.053e-04  2.936e-03  2.936e-03
  3  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     856  6.469e-10  2.129e-06  2.129e-06  2.245e-19  1.244e-13  1.887e-06  7.289e-07  1.887e-06  1.887e-06
  4  ASM_R  CONVERGED!
___________________________________________________________________________________________________________________________________________________________
Relative errors:
Displacement:   7.289e-07
Force:          2.451e-10
Dilatation:     1.353e-07
v / V_0:        1.000e-09 / 1.000e-09 = 1.000e+00


[...]

Timestep 10 @ 1.000e+00s
___________________________________________________________________________________________________________________________________________________________
                 SOLVER STEP                   |  LIN_IT   LIN_RES    RES_NORM     RES_U     RES_P      RES_J     NU_NORM      NU_U       NU_P       NU_J
___________________________________________________________________________________________________________________________________________________________
  0  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     874  2.358e-06  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00
  1  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     658  2.942e-04  1.544e-01  1.544e-01  1.208e+13  1.855e+06  6.014e-02  7.398e-02  6.014e-02  6.014e-02
  2  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     790  2.206e-06  2.908e-03  2.908e-03  7.302e+10  2.067e+03  2.716e-03  1.433e-03  2.716e-03  2.717e-03
  3  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     893  2.374e-09  1.919e-06  1.919e-06  4.527e+07  4.100e+00  1.672e-06  6.842e-07  1.672e-06  1.672e-06
  4  ASM_R  CONVERGED!
___________________________________________________________________________________________________________________________________________________________
Relative errors:
Displacement:   6.842e-07
Force:          8.995e-10
Dilatation:     1.528e-06
v / V_0:        1.000e-09 / 1.000e-09 = 1.000e+00
@endcode



Using the Timer class, we can discern which parts of the code require the highest computational expense.
For a case with a large number of degrees-of-freedom (i.e. a high level of refinement), a typical output of the Timer is given below.
Much of the code in the tutorial has been developed based on the optimizations described,
discussed and demonstrated in @ref step_18 "step-18" and others.
With over 93% of the time being spent in the linear solver, it is obvious that it may be necessary
to invest in a better solver for large three-dimensional problems.
The SSOR preconditioner is not multithreaded but is effective for this class of solid problems.
It may be beneficial to investigate the use of another solver such as those available through the Trilinos library.


@code
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    | 9.874e+02s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble system right-hand side |        53 | 1.727e+00s |  1.75e-01% |
| Assemble tangent matrix         |        43 | 2.707e+01s |  2.74e+00% |
| Linear solver                   |        43 | 9.248e+02s |  9.37e+01% |
| Linear solver postprocessing    |        43 | 2.743e-02s |  2.78e-03% |
| Perform static condensation     |        43 | 1.437e+01s |  1.46e+00% |
| Setup system                    |         1 | 3.897e-01s |  3.95e-02% |
| Update QPH data                 |        43 | 5.770e-01s |  5.84e-02% |
+---------------------------------+-----------+------------+------------+
@endcode


We then used ParaView to visualize the results for two cases.
The first was for the coarsest grid and the lowest-order interpolation method: @f$Q_1-DGPM_0-DGPM_0@f$.
The second was on a refined grid using a @f$Q_2-DGPM_1-DGPM_1@f$ formulation.
The vertical component of the displacement, the pressure @f$\widetilde{p}@f$ and the dilatation @f$\widetilde{J}@f$ fields
are shown below.


For the first case it is clear that the coarse spatial discretization coupled with large displacements leads to a low quality solution
(the loading ratio is  @f$p/p_0=80@f$).
Additionally, the pressure difference between elements is very large.
The constant pressure field on the element means that the large pressure gradient is not captured.
However, it should be noted that locking, which would be present in a standard @f$Q_1@f$ displacement formulation does not arise
even in this poorly discretised case.
The final vertical displacement of the tracked node on the top surface of the block is still within 12.5% of the converged solution.
The pressure solution is very coarse and has large jumps between adjacent cells.
It is clear that the volume nearest to the applied traction undergoes compression while the outer extents
of the domain are in a state of expansion.
The dilatation solution field and pressure field are clearly linked,
with positive dilatation indicating regions of positive pressure and negative showing regions placed in compression.
As discussed in the Introduction, a compressive pressure has a negative sign
while an expansive pressure takes a positive sign.
This stems from the definition of the volumetric strain energy function
and is opposite to the physically realistic interpretation of pressure.


<table align="center" class="tutorial" cellspacing="3" cellpadding="3">
  <tr>
    <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-displacement.png" alt="">
        <p align="center">
        Z-displacement solution for the 3-d problem.
        </p>
    </td>
    <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-pressure.png" alt="">
        <p align="center">
        Discontinuous piece-wise constant pressure field.
        </p>
    </td>
     <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-dilatation.png" alt="">
        <p align="center">
        Discontinuous piece-wise constant dilatation field.
        </p>
    </td>
  </tr>
</table>

Combining spatial refinement and a higher-order interpolation scheme results in a high-quality solution.
Three grid refinements coupled with a @f$Q_2-DGPM_1-DGPM_1@f$ formulation produces
a result that clearly captures the mechanics of the problem.
The deformation of the traction surface is well resolved.
We can now observe the actual extent of the applied traction, with the maximum force being applied
at the central point of the surface causing the largest compression.
Even though very high strains are experienced in the domain,
especially at the boundary of the region of applied traction,
the solution remains accurate.
The pressure field is captured in far greater detail than before.
There is a clear distinction and transition between regions of compression and expansion,
and the linear approximation of the pressure field allows a refined visualization
of the pressure at the sub-element scale.
It should however be noted that the pressure field remains discontinuous
and could be smoothed on a continuous grid for the post-processing purposes.



<table align="center" class="tutorial" cellspacing="3" cellpadding="3">
  <tr>
    <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-displacement.png" alt="">
        <p align="center">
        Z-displacement solution for the 3-d problem.
        </p>
    </td>
    <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-pressure.png" alt="">
        <p align="center">
        Discontinuous linear pressure field.
        </p>
    </td>
    <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-dilatation.png" alt="">
        <p align="center">
        Discontinuous linear dilatation field.
        </p>
    </td>
  </tr>
</table>

This brief analysis of the results demonstrates that the three-field formulation is effective
in circumventing volumetric locking for highly-incompressible media.
The mixed formulation is able to accurately simulate the displacement of a
near-incompressible block under compression.
The command-line output indicates that the volumetric change under extreme compression resulted in
less than 0.01% volume change for a Poisson's ratio of 0.4999.

In terms of run-time, the @f$Q_2-DGPM_1-DGPM_1@f$ formulation tends to be more computationally expensive
than the @f$Q_1-DGPM_0-DGPM_0@f$ for a similar number of degrees-of-freedom
(produced by adding an extra grid refinement level for the lower-order interpolation).
This is shown in the graph below for a batch of tests run consecutively on a single 4-core (8-thread) machine.
The increase in computational time for the higher-order method is likely due to
the increased band-width required for the higher-order elements.
As previously mentioned, the use of a better solver and preconditioner may mitigate the
expense of using a higher-order formulation.
It was observed that for the given problem using the multithreaded Jacobi preconditioner can reduce the
computational runtime by up to 72% (for the worst case being a higher-order formulation with a large number
of degrees-of-freedom) in comparison to the single-thread SSOR preconditioner.
However, it is the author's experience that the Jacobi method of preconditioning may not be suitable for
some finite-strain problems involving alternative constitutive models.


<table align="center" class="tutorial" cellspacing="3" cellpadding="3">
  <tr>
     <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.Normalised_runtime.png" alt="">
        <p align="center">
        Runtime on a 4-core machine, normalised against the lowest grid resolution @f$Q_1-DGPM_0-DGPM_0@f$ solution that utilised a SSOR preconditioner.
        </p>
    </td>
  </tr>
</table>


Lastly, results for the displacement solution for the 2-d problem are showcased below for
two different levels of grid refinement.
It is clear that due to the extra constraints imposed by simulating in 2-d that the resulting
displacement field, although qualitatively similar, is different to that of the 3-d case.


<table align="center" class="tutorial" cellspacing="3" cellpadding="3">
  <tr>
    <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.2d-gr_2.png" alt="">
        <p align="center">
        Y-displacement solution in 2-d for 2 global grid refinement levels.
        </p>
    </td>
    <td align="center">
        <img src="https://www.dealii.org/images/steps/developer/step-44.2d-gr_5.png" alt="">
        <p align="center">
        Y-displacement solution in 2-d for 5 global grid refinement levels.
        </p>
    </td>
  </tr>
</table>

<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>


There are a number of obvious extensions for this work:

- Firstly, an additional constraint could be added to the free-energy
  function in order to enforce a high degree of incompressibility in
  materials. An additional Lagrange multiplier would be introduced,
  but this could most easily be dealt with using the principle of
  augmented Lagrange multipliers. This is demonstrated in <em>Simo and
  Taylor (1991) </em>.
- The constitutive relationship used in this
  model is relatively basic. It may be beneficial to split the material
  class into two separate classes, one dealing with the volumetric
  response and the other the isochoric response, and produce a generic
  materials class (i.e. having abstract virtual functions that derived
  classes have to implement) that would allow for the addition of more complex
  material models. Such models could include other hyperelastic
  materials, plasticity and viscoelastic materials and others.
- The program has been developed for solving problems on single-node
  multicore machines. With a little effort, the program could be
  extended to a large-scale computing environment through the use of
  Petsc or Trilinos, using a similar technique to that demonstrated in
  @ref step_40 "step-40". This would mostly involve changes to the setup, assembly,
  <code>PointHistory</code> and linear solver routines.
- As this program assumes quasi-static equilibrium, extensions to
  include dynamic effects would be necessary to study problems where
  inertial effects are important, e.g. problems involving impact.
- Load and solution limiting procedures may be necessary for highly
  nonlinear problems. It is possible to add a linesearch algorithm to
  limit the step size within a Newton increment to ensure optimum
  convergence. It may also be necessary to use a load limiting method,
  such as the Riks method, to solve unstable problems involving
  geometric instability such as buckling and snap-through.
- Many physical problems involve contact. It is possible to include
  the effect of frictional or frictionless contact between objects
  into this program. This would involve the addition of an extra term
  in the free-energy functional and therefore an addition to the
  assembly routine. One would also need to manage the contact problem
  (detection and stress calculations) itself. An alternative to
  additional penalty terms in the free-energy functional would be to
  use active set methods such as the one used in @ref step_41 "step-41".
- The complete condensation procedure using LinearOperators has been
  coded into the linear solver routine. This could also have been
  achieved through the application of the schur_complement()
  operator to condense out one or more of the fields in a more
  automated manner.
- Finally, adaptive mesh refinement, as demonstrated in @ref step_6 "step-6" and
  @ref step_18 "step-18", could provide additional solution accuracy.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-44.cc"
*/