]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Updated source, still not converging
authormcbride <mcbride@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 11 Feb 2012 21:01:38 +0000 (21:01 +0000)
committermcbride <mcbride@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 11 Feb 2012 21:01:38 +0000 (21:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@25035 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-44/step-44.cc

index 591243c5718ee0215ecc8730ce35c5a93ddd0049..77c4f822395b4904dc24ce6ad9d14fbbf5d91075 100644 (file)
@@ -1,4 +1,3 @@
-
 /* Authors: Jean-Paul Pelteret, University of Cape Town,            */
 /*          Andrew McBride, University of Erlangen-Nuremberg, 2010  */
 /*                                                                  */
 /*    to the file deal.II/doc/license.html for the  text  and       */
 /*    further information on this license.                          */
 
-// @sect3{Include files}
-
 // We start by including all the necessary
 // deal.II header files and some C++ related
 // ones. They have been discussed in detail
 // in previous tutorial programs, so you need
 // only refer to past tutorials for details.
-
-#include <base/function.h>
-#include <base/parameter_handler.h>
-#include <base/point.h>
-#include <base/quadrature_lib.h>
-#include <base/symmetric_tensor.h>
-#include <base/tensor.h>
-#include <base/timer.h>
-#include <base/work_stream.h>
-
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_renumbering.h>
-#include <dofs/dof_tools.h>
-
-#include <grid/grid_generator.h>
-#include <grid/grid_tools.h>
-#include <grid/grid_in.h>
-#include <grid/tria.h>
-#include <grid/tria_boundary_lib.h>
-
-#include <fe/fe_dgp_monomial.h>
-#include <fe/fe_q.h>
-#include <fe/fe_system.h>
-#include <fe/fe_tools.h>
-#include <fe/fe_values.h>
-
-#include <fe/mapping_q_eulerian.h>
-
-#include <lac/block_sparse_matrix.h>
-#include <lac/block_vector.h>
-#include <lac/compressed_sparsity_pattern.h>
-#include <lac/full_matrix.h>
-#include <lac/precondition.h>
-#include <lac/solver_cg.h>
-#include <lac/sparse_direct.h>
-
-#include <numerics/data_out.h>
-#include <numerics/vectors.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/work_stream.h>
+
+#include <deal.II/dofs/dof_constraints.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/fe/fe_dgp_monomial.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_direct.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vectors.h>
 
 #include <iostream>
 #include <fstream>
@@ -66,48 +61,49 @@ using namespace dealii;
 // @sect3{Run-time parameters}
 //
 // There are several parameters that can be set
-// so we choose to set up a parameter
-// handler object so that we can read in choices
-// at run-time.
-namespace Parameters
-{
+// in the code so we set up a ParameterHandler
+// object to read in the choices at run-time.
+namespace Parameters {
 // @sect4{Finite Element system}
-// Change the polynomial order used to approximate the solution.
-// The quadrature should be adjusted accordingly.
-struct FESystem
-{
-    int poly_degree;
-    int quad_order;
-
-    static void declare_parameters (ParameterHandler &prm);
-    void parse_parameters (ParameterHandler &prm);
+// As mentioned in the introduction, a different order
+// interpolation should be used for the displacement
+// $\mathbf{u}$ than for the pressure $p$ and
+// the dilatation $\widetilde{J}$.
+// Choosing $p$ and $\widetilde{J}$ as discontinuous (constant)
+// functions at the element level leads to the
+// mean-dilatation method. The discontinuous approximation
+// allows $p$ and $\widetilde{J}$ to be condensed out
+// and a classical displacement based method is recovered.
+// Here we specify the polynomial order used to
+// approximate the solution.
+// The quadrature order should be adjusted accordingly.
+struct FESystem {
+       int poly_degree;
+       int quad_order;
+
+       static void declare_parameters(ParameterHandler &prm);
+       void parse_parameters(ParameterHandler &prm);
 };
 
-void FESystem::declare_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Finite element system");
-    {
-        prm.declare_entry("Polynomial degree",
-                          "1",
-                          Patterns::Integer(),
-                          "Displacement system polynomial order");
-
-        prm.declare_entry("Quadrature order",
-                          "2",
-                          Patterns::Integer(),
-                          "Gauss quadrature order");
-    }
-    prm.leave_subsection();
+void FESystem::declare_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Finite element system");
+       {
+               prm.declare_entry("Polynomial degree", "1", Patterns::Integer(),
+                               "Displacement system polynomial order");
+
+               prm.declare_entry("Quadrature order", "2", Patterns::Integer(),
+                               "Gauss quadrature order");
+       }
+       prm.leave_subsection();
 }
 
-void FESystem::parse_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Finite element system");
-    {
-        poly_degree = prm.get_integer("Polynomial degree");
-        quad_order = prm.get_integer("Quadrature order");
-    }
-    prm.leave_subsection();
+void FESystem::parse_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Finite element system");
+       {
+               poly_degree = prm.get_integer("Polynomial degree");
+               quad_order = prm.get_integer("Quadrature order");
+       }
+       prm.leave_subsection();
 }
 
 // @sect4{Geometry}
@@ -115,1681 +111,1683 @@ void FESystem::parse_parameters (ParameterHandler &prm)
 // Since the problem modelled here is quite specific, the load
 // scale can be altered to specific values to attain results given
 // in the literature.
-struct Geometry
-{
-    int global_refinement;
-    double scale;
-    double p_p0;
+struct Geometry {
+       int global_refinement;
+       double scale;
+       double p_p0;
 
-    static void declare_parameters (ParameterHandler &prm);
-    void parse_parameters (ParameterHandler &prm);
+       static void declare_parameters(ParameterHandler &prm);
+       void parse_parameters(ParameterHandler &prm);
 };
 
-void Geometry::declare_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Geometry");
-    {
-        prm.declare_entry("Global refinement",
-                          "2",
-                          Patterns::Integer(),
-                          "Global refinement level");
-
-        prm.declare_entry("Grid scale",
-                          "1.0",
-                          Patterns::Double(),
-                          "Global grid scaling factor");
-
-        prm.declare_entry("Pressure ratio p/p0",
-                          "40",
-                          Patterns::Selection("20|40|60|80|100"),
-                          "Ratio of applied pressure to reference pressure");
-    }
-    prm.leave_subsection();
+void Geometry::declare_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Geometry");
+       {
+               prm.declare_entry("Global refinement", "2", Patterns::Integer(),
+                               "Global refinement level");
+
+               prm.declare_entry("Grid scale", "1.0", Patterns::Double(),
+                               "Global grid scaling factor");
+
+               prm.declare_entry("Pressure ratio p/p0", "40",
+                               Patterns::Selection("20|40|60|80|100"),
+                               "Ratio of applied pressure to reference pressure");
+       }
+       prm.leave_subsection();
 }
 
-void Geometry::parse_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Geometry");
-    {
-        global_refinement = prm.get_integer("Global refinement");
-        scale = prm.get_double("Grid scale");
-        p_p0 = prm.get_double("Pressure ratio p/p0");
-    }
-    prm.leave_subsection();
+void Geometry::parse_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Geometry");
+       {
+               global_refinement = prm.get_integer("Global refinement");
+               scale = prm.get_double("Grid scale");
+               p_p0 = prm.get_double("Pressure ratio p/p0");
+       }
+       prm.leave_subsection();
 }
 
 // @sect4{Materials}
-// Store the shear modulus and Lame constant
-// for the Neo-Hookean material
-struct Materials
-{
-    double nu;
-    double mu;
-
-    static void declare_parameters (ParameterHandler &prm);
-    void parse_parameters (ParameterHandler &prm);
+// Need the shear modulus $ \mu $
+// and Poisson ration $ \nu $
+// for the neo-Hookean material.
+struct Materials {
+       double nu;
+       double mu;
+
+       static void declare_parameters(ParameterHandler &prm);
+       void parse_parameters(ParameterHandler &prm);
 };
 
-void Materials::declare_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Material properties");
-    {
-        prm.declare_entry("Poisson's ratio",
-                          "0.49",
-                          Patterns::Double(),
-                          "Poisson's ratio");
-
-        prm.declare_entry("Shear modulus",
-                          "1.0e6",
-                          Patterns::Double(),
-                          "Shear modulus");
-    }
-    prm.leave_subsection();
+void Materials::declare_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Material properties");
+       {
+               prm.declare_entry("Poisson's ratio", "0.49", Patterns::Double(),
+                               "Poisson's ratio");
+
+               prm.declare_entry("Shear modulus", "1.0e6", Patterns::Double(),
+                               "Shear modulus");
+       }
+       prm.leave_subsection();
 }
 
-void Materials::parse_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Material properties");
-    {
-        nu = prm.get_double("Poisson's ratio");
-        mu = prm.get_double("Shear modulus");
-    }
-    prm.leave_subsection();
+void Materials::parse_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Material properties");
+       {
+               nu = prm.get_double("Poisson's ratio");
+               mu = prm.get_double("Shear modulus");
+       }
+       prm.leave_subsection();
 }
 
 // @sect4{Linear solver}
 // Choose both CG solver and SSOR preconditioner settings.
 // The default values are optimal for this particular problem.
-struct LinearSolver
-{
-    std::string type_lin;
-    double tol_lin;
-    double max_iterations_lin;
-    double ssor_relaxation;
-
-    static void declare_parameters (ParameterHandler &prm);
-    void parse_parameters (ParameterHandler &prm);
+struct LinearSolver {
+       std::string type_lin;
+       double tol_lin;
+       double max_iterations_lin;
+       double ssor_relaxation;
+
+       static void declare_parameters(ParameterHandler &prm);
+       void parse_parameters(ParameterHandler &prm);
 };
 
-void LinearSolver::declare_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Linear solver");
-    {
-        prm.declare_entry("Solver type",
-                          "CG",
-                          Patterns::Selection("CG|Direct"),
-                          "Type of solver used to solve the linear system");
-
-        prm.declare_entry("Residual",
-                          "1e-6",
-                          Patterns::Double(),
-                          "Linear solver residual (scaled by residual norm)");
-
-        prm.declare_entry("Max iteration multiplier",
-                          "2",
-                          Patterns::Double(),
-                          "Linear solver iterations (multiples of the system matrix size)");
-
-        prm.declare_entry("SSOR Relaxation",
-                          "0.6",
-                          Patterns::Double(),
-                          "SSOR preconditioner relaxation value");
-    }
-    prm.leave_subsection();
+void LinearSolver::declare_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Linear solver");
+       {
+               prm.declare_entry("Solver type", "CG", Patterns::Selection("CG|Direct"),
+                               "Type of solver used to solve the linear system");
+
+               prm.declare_entry("Residual", "1e-6", Patterns::Double(),
+                               "Linear solver residual (scaled by residual norm)");
+
+               prm.declare_entry(
+                               "Max iteration multiplier",
+                               "2",
+                               Patterns::Double(),
+                               "Linear solver iterations (multiples of the system matrix size)");
+
+               prm.declare_entry("SSOR Relaxation", "0.6", Patterns::Double(),
+                               "SSOR preconditioner relaxation value");
+       }
+       prm.leave_subsection();
 }
 
-void LinearSolver::parse_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Linear solver");
-    {
-        type_lin = prm.get("Solver type");
-        tol_lin = prm.get_double("Residual");
-        max_iterations_lin = prm.get_double("Max iteration multiplier");
-        ssor_relaxation = prm.get_double("SSOR Relaxation");
-    }
-    prm.leave_subsection();
+void LinearSolver::parse_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Linear solver");
+       {
+               type_lin = prm.get("Solver type");
+               tol_lin = prm.get_double("Residual");
+               max_iterations_lin = prm.get_double("Max iteration multiplier");
+               ssor_relaxation = prm.get_double("SSOR Relaxation");
+       }
+       prm.leave_subsection();
 }
 
 // @sect4{Nonlinear solver}
-// Define the tolerances and maximum number of iterations for the
-// Newton-Raphson nonlinear solver.
-struct NonlinearSolver
-{
-    unsigned int max_iterations_NR;
-    double tol_f;
-    double tol_u;
-
-    static void declare_parameters (ParameterHandler &prm);
-    void parse_parameters (ParameterHandler &prm);
+// A Newton-Raphson scheme is used to
+// solve the nonlinear system of governing equations.
+// Define the tolerances and the maximum number of
+// iterations for the Newton-Raphson nonlinear solver.
+struct NonlinearSolver {
+       unsigned int max_iterations_NR;
+       double tol_f;
+       double tol_u;
+
+       static void declare_parameters(ParameterHandler &prm);
+       void parse_parameters(ParameterHandler &prm);
 };
 
-void NonlinearSolver::declare_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Nonlinear solver");
-    {
-        prm.declare_entry("Max iterations Newton-Raphson",
-                          "10",
-                          Patterns::Integer(),
-                          "Number of Newton-Raphson iterations allowed");
-
-        prm.declare_entry("Tolerance force",
-                          "1.0e-9",
-                          Patterns::Double(),
-                          "Force residual tolerance");
-
-        prm.declare_entry("Tolerance displacement",
-                          "1.0e-3",
-                          Patterns::Double(),
-                          "Displacement error tolerance");
-    }
-    prm.leave_subsection();
+void NonlinearSolver::declare_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Nonlinear solver");
+       {
+               prm.declare_entry("Max iterations Newton-Raphson", "10",
+                               Patterns::Integer(),
+                               "Number of Newton-Raphson iterations allowed");
+
+               prm.declare_entry("Tolerance force", "1.0e-9", Patterns::Double(),
+                               "Force residual tolerance");
+
+               prm.declare_entry("Tolerance displacement", "1.0e-3",
+                               Patterns::Double(), "Displacement error tolerance");
+       }
+       prm.leave_subsection();
 }
 
-void NonlinearSolver::parse_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Nonlinear solver");
-    {
-        max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
-        tol_f = prm.get_double("Tolerance force");
-        tol_u = prm.get_double("Tolerance displacement");
-    }
-    prm.leave_subsection();
+void NonlinearSolver::parse_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Nonlinear solver");
+       {
+               max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
+               tol_f = prm.get_double("Tolerance force");
+               tol_u = prm.get_double("Tolerance displacement");
+       }
+       prm.leave_subsection();
 }
 
 // @sect4{Time}
-// Set the timestep size and the simulation end-time.
-struct Time
-{
-    double end_time;
-    double delta_t;
-
-    static void declare_parameters (ParameterHandler &prm);
-    void parse_parameters (ParameterHandler &prm);
+// Set the timestep size $ \varDelta t $
+// and the simulation end-time.
+struct Time {
+       double delta_t;
+       double end_time;
+
+       static void declare_parameters(ParameterHandler &prm);
+       void parse_parameters(ParameterHandler &prm);
 };
 
-void Time::declare_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Time");
-    {
-        prm.declare_entry("End time",
-                          "1",
-                          Patterns::Double(),
-                          "End time");
-
-        prm.declare_entry("Time step size",
-                          "0.1",
-                          Patterns::Double(),
-                          "Time step size");
-    }
-    prm.leave_subsection();
+void Time::declare_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Time");
+       {
+               prm.declare_entry("End time", "1", Patterns::Double(), "End time");
+
+               prm.declare_entry("Time step size", "0.1", Patterns::Double(),
+                               "Time step size");
+       }
+       prm.leave_subsection();
 }
 
-void Time::parse_parameters (ParameterHandler &prm)
-{
-    prm.enter_subsection("Time");
-    {
-        end_time  = prm.get_double("End time");
-        delta_t  = prm.get_double("Time step size");
-    }
-    prm.leave_subsection();
+void Time::parse_parameters(ParameterHandler &prm) {
+       prm.enter_subsection("Time");
+       {
+               end_time = prm.get_double("End time");
+               delta_t = prm.get_double("Time step size");
+       }
+       prm.leave_subsection();
 }
 
 // @sect4{All parameters}
 // Finally we consolidate all of the above structures into
 // a single container that holds all of our run-time selections.
-struct AllParameters
-       :
-       public FESystem,
-       public Geometry,
-       public Materials,
-       public LinearSolver,
-       public NonlinearSolver,
-       public Time
+struct AllParameters: public FESystem,
+               public Geometry,
+               public Materials,
+               public LinearSolver,
+               public NonlinearSolver,
+               public Time
 
 {
-    AllParameters (const std::string & input_file);
+       AllParameters(const std::string & input_file);
 
-    static void declare_parameters (ParameterHandler &prm);
-    void parse_parameters (ParameterHandler &prm);
+       static void declare_parameters(ParameterHandler &prm);
+       void parse_parameters(ParameterHandler &prm);
 };
 
-AllParameters::AllParameters (const std::string & input_file)
-{
-    ParameterHandler prm;
-    declare_parameters(prm);
-    prm.read_input (input_file);
-    parse_parameters(prm);
+AllParameters::AllParameters(const std::string & input_file) {
+       ParameterHandler prm;
+       declare_parameters(prm);
+       prm.read_input(input_file);
+       parse_parameters(prm);
 }
 
-void AllParameters::declare_parameters (ParameterHandler &prm)
-{
-    FESystem::declare_parameters(prm);
-    Geometry::declare_parameters(prm);
-    Materials::declare_parameters(prm);
-    LinearSolver::declare_parameters(prm);
-    NonlinearSolver::declare_parameters(prm);
-    Time::declare_parameters(prm);
+void AllParameters::declare_parameters(ParameterHandler &prm) {
+       FESystem::declare_parameters(prm);
+       Geometry::declare_parameters(prm);
+       Materials::declare_parameters(prm);
+       LinearSolver::declare_parameters(prm);
+       NonlinearSolver::declare_parameters(prm);
+       Time::declare_parameters(prm);
 }
 
-void AllParameters::parse_parameters (ParameterHandler &prm)
-{
-    FESystem::parse_parameters(prm);
-    Geometry::parse_parameters(prm);
-    Materials::parse_parameters(prm);
-    LinearSolver::parse_parameters(prm);
-    NonlinearSolver::parse_parameters(prm);
-    Time::parse_parameters(prm);
+void AllParameters::parse_parameters(ParameterHandler &prm) {
+       FESystem::parse_parameters(prm);
+       Geometry::parse_parameters(prm);
+       Materials::parse_parameters(prm);
+       LinearSolver::parse_parameters(prm);
+       NonlinearSolver::parse_parameters(prm);
+       Time::parse_parameters(prm);
+}
 }
-
-}  // End Parameters namespace
 
 // @sect3{General tools}
 // We need to perform some specific operations that are not defined
-// in the deal.II library yet. We place these common operations
-// in a seperate namespace for convenience.
-namespace AdditionalTools
-{
-// Define an operation that takes two tensors $ \mathbf{A} $ and
-// $ \mathbf{B} $ such that their outer-product
-// $ \mathbf{A} \bar{\otimes} \mathbf{B} \Rightarrow C_{ijkl} = A_{ik} B_{jl} $
-template <int dim>
-SymmetricTensor<4,dim> outer_product_T23 (const SymmetricTensor<2,dim> & A,
-                                          const SymmetricTensor<2,dim> & B)
-{
-    SymmetricTensor<4,dim> A_ik_B_jl;
-
-    for (unsigned int i=0; i<dim; ++i) {
-        for (unsigned int j=i; j<dim; ++j) {
-            for (unsigned int k=0; k<dim; ++k) {
-                for (unsigned int l=k; k<dim; ++k) {
-                    A_ik_B_jl[i][j][k][l] += A[i][k] * B[j][l];
-                }
-            }
-        }
-    }
-
-    return A_ik_B_jl;
-}
-
-// The \a extract_submatrix function takes specific entries from a \a matrix,
-// and copies them to a \a sub_matrix. The copied entries are defined by the
-// first two parameters which hold the row and column entries to be extracted.
-// The \a matrix is automatically resized to size $ r \times c $.
-template <typename MatrixType>
-void extract_submatrix(const std::vector< unsigned int > &row_index_set,
-                       const std::vector< unsigned int > &column_index_set,
-                       const MatrixType &matrix,
-                       FullMatrix< double > &sub_matrix)
-{
-
-    const unsigned int n_rows_submatrix = row_index_set.size();
-    const unsigned int n_cols_submatrix = column_index_set.size();
-
-    sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix);
-
-    for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) {
-        const unsigned int row = row_index_set[sub_row];
-        Assert (row<=matrix.m(), ExcInternalError());
+// in the deal.II library yet.
+// We place these common operations
+// in a separate namespace for convenience.
+// We also include some widely used operators
+namespace AdditionalTools {
+// Define an operation that takes two
+// symmetric second-order tensors
+// $\mathbf{A}$ and $\mathbf{B}$
+// such that their outer-product
+// $ \mathbf{A} \overline{\otimes} \mathbf{B} \Rightarrow C_{ijkl} = A_{ik} B_{jl} $
+template<int dim>
+SymmetricTensor<4, dim> outer_product_T23(const SymmetricTensor<2, dim> & A,
+const SymmetricTensor<2, dim> & B) {
+       SymmetricTensor<4, dim> A_ik_B_jl;
+
+       for (unsigned int i = 0; i < dim; ++i) {
+               for (unsigned int j = i; j < dim; ++j) {
+                       for (unsigned int k = 0; k < dim; ++k) {
+                               for (unsigned int l = k; k < dim; ++k) {
+                                       A_ik_B_jl[i][j][k][l] += A[i][k] * B[j][l];
+                               }
+                       }
+               }
+       }
 
-        for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) {
-            const unsigned int col = column_index_set[sub_col];
-            Assert (col<=matrix.n(), ExcInternalError());
+       return A_ik_B_jl;
+}
 
-            sub_matrix(sub_row,sub_col) = matrix(row, col);
-        }
-    }
+// The  extract_submatrix function
+// takes specific entries from a matrix,
+// and copies them to a  sub_matrix.
+// The copied entries are defined by the
+// first two parameters which hold the
+// row and columns to be extracted.
+// The  matrix is automatically resized
+// to size $ r \times c $.
+template<typename MatrixType>
+void extract_submatrix(const std::vector<unsigned int> &row_index_set,
+               const std::vector<unsigned int> &column_index_set,
+               const MatrixType &matrix, FullMatrix<double> &sub_matrix) {
+
+       const unsigned int n_rows_submatrix = row_index_set.size();
+       const unsigned int n_cols_submatrix = column_index_set.size();
+       // check the size of the input vectors
+       Assert(n_rows_submatrix > 0, ExcInternalError());
+       Assert(n_cols_submatrix > 0, ExcInternalError());
+
+       sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix);
+
+       for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) {
+               const unsigned int row = row_index_set[sub_row];
+               Assert(row<=matrix.m(), ExcInternalError());
+
+               for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) {
+                       const unsigned int col = column_index_set[sub_col];
+                       Assert(col<=matrix.n(), ExcInternalError());
+
+                       sub_matrix(sub_row, sub_col) = matrix(row, col);
+               }
+       }
 }
 
-template <>
-void extract_submatrix < dealii::BlockSparseMatrix <double> >(const std::vector< unsigned int > &row_index_set,
-                                                             const std::vector< unsigned int > &column_index_set,
-                                                             const dealii::BlockSparseMatrix <double> &matrix,
-                                                             FullMatrix< double > &sub_matrix)
-{
+// As above, but to extract entries from
+// a <code> BlockSparseMatrix </code>.
+template<>
+void extract_submatrix<dealii::BlockSparseMatrix<double> >(
+               const std::vector<unsigned int> &row_index_set,
+               const std::vector<unsigned int> &column_index_set,
+               const dealii::BlockSparseMatrix<double> &matrix,
+               FullMatrix<double> &sub_matrix) {
 
-    const unsigned int n_rows_submatrix = row_index_set.size();
-    const unsigned int n_cols_submatrix = column_index_set.size();
+       const unsigned int n_rows_submatrix = row_index_set.size();
+       const unsigned int n_cols_submatrix = column_index_set.size();
 
-    sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix);
+       // check the size of the input vectors
+       Assert(n_rows_submatrix > 0, ExcInternalError());
+       Assert(n_cols_submatrix > 0, ExcInternalError());
 
-    for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) {
-        const unsigned int row = row_index_set[sub_row];
-        Assert (row<=matrix.m(), ExcInternalError());
+       sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix);
 
-        for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) {
-            const unsigned int col = column_index_set[sub_col];
-            Assert (col<=matrix.n(), ExcInternalError());
-            if (matrix.get_sparsity_pattern().exists(row, col) == false) continue;
+       for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) {
+               const unsigned int row = row_index_set[sub_row];
+               Assert(row<=matrix.m(), ExcInternalError());
 
-            sub_matrix(sub_row,sub_col) = matrix(row, col);
-        }
-    }
-}
+               for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) {
+                       const unsigned int col = column_index_set[sub_col];
+                       Assert(col<=matrix.n(), ExcInternalError());
+                       if (matrix.get_sparsity_pattern().exists(row, col) == false)
+                               continue;
 
-// The \a replace_submatrix function takes specific entries from a \a matrix,
-// and copies them to a \a sub_matrix. The copied entries are defined by the
-// first two parameters which hold the row and column entries to be replaced.
-// The \a matrix expected to be of the correct size.
-template <typename MatrixType>
-void replace_submatrix(const std::vector< unsigned int > &row_index_set,
-                       const std::vector< unsigned int > &column_index_set,
-                       const MatrixType &sub_matrix,
-                       FullMatrix< double >  &matrix)
-{
-    const unsigned int n_rows_submatrix = row_index_set.size();
-    Assert (n_rows_submatrix<=sub_matrix.m(), ExcInternalError());
-    const unsigned int n_cols_submatrix = column_index_set.size();
-    Assert (n_cols_submatrix<=sub_matrix.n(), ExcInternalError());
+                       sub_matrix(sub_row, sub_col) = matrix(row, col);
+               }
+       }
+}
 
-    for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) {
-        const unsigned int row = row_index_set[sub_row];
-        Assert (row<=matrix.m(), ExcInternalError());
+// The replace_submatrix function takes
+// specific entries from a  sub_matrix,
+// and copies them into a  matrix.
+// The copied entries are defined by the
+// first two parameters which hold the
+// row and column entries to be replaced.
+// The matrix expected to be of the correct size.
+template<typename MatrixType>
+void replace_submatrix(const std::vector<unsigned int> &row_index_set,
+               const std::vector<unsigned int> &column_index_set,
+               const MatrixType &sub_matrix, FullMatrix<double> &matrix) {
+       const unsigned int n_rows_submatrix = row_index_set.size();
+       Assert(n_rows_submatrix<=sub_matrix.m(), ExcInternalError());
+       const unsigned int n_cols_submatrix = column_index_set.size();
+       Assert(n_cols_submatrix<=sub_matrix.n(), ExcInternalError());
+
+       for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) {
+               const unsigned int row = row_index_set[sub_row];
+               Assert(row<=matrix.m(), ExcInternalError());
+
+               for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) {
+                       const unsigned int col = column_index_set[sub_col];
+                       Assert(col<=matrix.n(), ExcInternalError());
+
+                       matrix(row, col) = sub_matrix(sub_row, sub_col);
 
-        for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) {
-            const unsigned int col = column_index_set[sub_col];
-            Assert (col<=matrix.n(), ExcInternalError());
+               }
+       }
+}
 
-            matrix(row, col) = sub_matrix(sub_row, sub_col);
+// Define some frequently used
+// second and fourth-order tensors:
+template<int dim>
+class StandardTensors {
+public:
 
-        }
-    }
-}
+       // $\mathbf{I}$
+       static SymmetricTensor<2, dim> const I;
+       // $\mathbf{I} \otimes \mathbf{I}$
+       static SymmetricTensor<4, dim> const IxI;
+       // $\mathcal{S}$, note that as we only use
+       // this fourth-order unit tensor to operate
+       // on symmetric second-order tensors.
+       // To maintain notation consistent with Holzapfel (2001)
+       // we name the tensor $\mathcal{I}$
+       static SymmetricTensor<4, dim> const II;
+       // Fourth-order deviatoric such that
+       // $\textrm{dev}(\bullet) = (\bullet) - (1/3)[(\bullet):\mathbf{I}]\mathbf{I}$
+       static SymmetricTensor<4, dim> const dev_P;
+};
 
+template<int dim>
+SymmetricTensor<2, dim> const StandardTensors<dim>::I = SymmetricTensor<2, dim>(
+               unit_symmetric_tensor<dim>());
+template<int dim>
+SymmetricTensor<4, dim> const StandardTensors<dim>::IxI =
+               SymmetricTensor<4, dim>(outer_product(I, I));
+template<int dim>
+SymmetricTensor<4, dim> const StandardTensors<dim>::II =
+               SymmetricTensor<4, dim>(identity_tensor<dim>());
+template<int dim>
+SymmetricTensor<4, dim> const StandardTensors<dim>::dev_P = (II
+               - 1.0 / 3.0 * IxI);
 }
 
 // @sect3{Time class}
-// A simple class to store time data is created. Its
+// A simple class to store time data. Its
 // functioning is transparent so no discussion is
-// necessary.
+// necessary. For simplicity we assume a constant
+// time step size.
 class Time {
 public:
-    Time (const double & time_end,
-          const double & delta_t)
-       :
-          timestep (0),
-          time_current (0.0),
-          time_end (time_end),
-          delta_t (delta_t)
-    {}
-    virtual ~Time (void) {}
-
-    const double & current (void) const {return time_current;}
-    const double & end (void) const {return time_end;}
-    const double & get_delta_t (void) const {return delta_t;}
-    const unsigned int & get_timestep (void) const {return timestep;}
-    void increment (void) {time_current += delta_t; ++timestep;}
+       Time(const double & time_end, const double & delta_t) :
+                       timestep(0), time_current(0.0), time_end(time_end), delta_t(delta_t) {
+       }
+       virtual ~Time(void) {
+       }
+
+       const double current(void) const {
+               return time_current;
+       }
+       const double end(void) const {
+               return time_end;
+       }
+       const double get_delta_t(void) const {
+               return delta_t;
+       }
+       const unsigned int get_timestep(void) const {
+               return timestep;
+       }
+       void increment(void) {
+               time_current += delta_t;
+               ++timestep;
+       }
 
 private:
-    unsigned int timestep;
-    double time_current;
-    const double time_end;
-    const double delta_t;
+       unsigned int timestep;
+       double time_current;
+       const double time_end;
+       const double delta_t;
 };
 
-// @sect3{Neo-Hookean material}
-// The entire domain is to made of a Neo-Hookean material
-// with constant properties throughout. This class defines
-// the behaviour of this material. Neo-Hookean materials
+// @sect3{Compressible neo-Hookean material}
+
+// As discussed in the Introduction,
+// Neo-Hookean materials are a
+// type of hyperelastic materials.
+// The entire domain is assumed 
+// to be composed of a compressible neo-Hookean material. 
+// This class defines
+// the behaviour of this material. 
+// Compressible neo-Hookean materials
 // can be described by a strain-energy function (SEF)
-// $ \phi = \phi_{B} + \phi_{V}  $
-// where the bulk deformation is given by
-// $ \phi_{B} = C_{1} \left( I_{1} - 3 \right)  $
-// where $ C_{1} - \frac{\mu}{2} $ and $I_{1}$ is the first
-// invariant of the left- or right- Cauchy deformation tensors.
+// $ \Psi = \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(J) $.
+//
+// The isochoric response is given by
+// $ \Psi_{\text{iso}}(\mathbf{b}) = c_{1} [\overline{I}_{1} - 3]  $
+// where $ c_{1} = \frac{\mu}{2} $ and $\overline{I}_{1}$ is the first
+// invariant of the left- or right- isochoric Cauchy-Green deformation tensors.
+// That is $\overline{I}_1 :=\textrm{tr}(\overline{\mathbf{b}})$.
 // In this example the SEF that governs the volumetric
 // response is defined as
-// $ \phi_{V} = \kappa \left( \frac{1}{2} \left( \theta^{2} - 1 \right) - ln \left( \theta \right) \right)  $
-// where $\kappa$ is the bulk modulus.
-template <int dim>
-class Material_NH
-{
+// $ \Psi_{\text{vol}}(\widetilde{J})  = \kappa \bigl[ \frac{1}{2} [ \widetilde{J}^{2} - 1 ] - \textrm{ln}( \widetilde{J}) ] \bigr]  $
+// where $\kappa:= \lambda + 2/3 \mu$ is the bulk modulus and
+// $\lambda$ is a Lame moduli.
+template<int dim>
+class Material_Compressilbe_Neo_Hook_Uncoupled {
 public:
-    Material_NH (const double & lambda,
-                const double & mu)
-        :
-          lambda_0 (lambda),
-          mu_0 (mu),
-          kappa_0 (lambda + 2.0/3.0*mu)
-    { }
-    ~Material_NH (void) {}
-
-    // The Kirchhoff stress tensor is required in the formulation
-    // used in this work. This is obtained from the SEF by
-    // $ \mathbf{T} = \mathbf{F} \frac{\partial \hat{\phi}}{\partial \mathbf{C}} \mathbf{F}^{T} = \frac{\partial \phi}{\partial \mathbf{B}} \mathbf{B} $
-    SymmetricTensor<2, dim> get_T (const double & J,
-                                   const SymmetricTensor <2, dim> & B)
-    {
-       const double dW_dJ  = get_dU_dtheta (J);
-       return mu_0*B + dW_dJ*J*I;
-    }
-
-    // The tangent matrix for this material is also calculated from the SEF by
-    // $ JC_{ijkl} = F_{iA} F_{jB} H_{ABCD} F_{kC} F_{lD}$
-    // with
-    // $ \mathbf{H} = \frac{\partial^{2} \hat{\phi}}{\partial \mathbf{C} \partial \mathbf{C}} $
-    SymmetricTensor<4, dim> get_JC (const double & J,
-                                    const SymmetricTensor <2, dim> & B)
-    {
-       const double dW_dJ   = get_dU_dtheta (J);
-       const double d2W_dJ2 = get_d2U_dtheta2 (J);
-       return  J*(  (dW_dJ + J*d2W_dJ2)*IxI - (2.0*dW_dJ)*II  );
-    }
-
-    // From the volumetric strain-energy function we calculate the
-    // first and second derivatives with respect to the dilatation
-    double get_dU_dtheta    (const double & d) {return kappa_0*(d - 1.0/d);}
-    double get_d2U_dtheta2  (const double & d) {return kappa_0*(1.0 + 1.0/(d*d));}
+       Material_Compressilbe_Neo_Hook_Uncoupled(const double mu, const double nu) :
+                       kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))), c_1(
+                                       mu / 2.0), det_F(1.0), J_tilde(1.0), b_bar(
+                                       AdditionalTools::StandardTensors<dim>::I) {
+               Assert(kappa > 0, ExcInternalError());
+       }
+       ~Material_Compressilbe_Neo_Hook_Uncoupled(void) {
+       }
+
+       // The Kirchhoff stress tensor $\boldsymbol{\tau}$ is
+       // the chosen stress measure.
+       // Recall that
+       // $\boldsymbol{\tau} = \chi_{*}(\mathbf{S})$, i.e.
+       // $\boldsymbol{\tau} = \mathbf{F} \mathbf{S} \mathbf{F}^{T}$.
+       // Furthermore,
+       // $\boldsymbol{\tau} = 2 \mathbf{F} \frac{\partial \Psi(\mathbf{C})}{\partial \mathbf{C}}  \mathbf{F}^{T} = 2 \mathbf{b} \frac{\partial \Psi(\mathbf{b})}{\partial \mathbf{b}}$.
+       // Therefore,
+       // $\boldsymbol{\tau} = 2 \mathbf{b} \bigl[ \frac{\partial \Psi_{\text{iso}}(\mathbf{b})}{\partial \mathbf{b}} + \frac{\partial \Psi_{\text{vol}}(J)}{\partial J}\frac{\partial J}{\partial \mathbf{b}} \bigr] =  2 \mathbf{b} \frac{\partial \Psi_{\text{iso}}(\mathbf{b})}{\partial \mathbf{b}} + J\frac{\partial \Psi_{\text{vol}}(J)}{\partial J}\mathbf{I} $
+
+       // We update the material model with various deformation
+       // dependent data based on F
+       void update_material_data(const Tensor<2, dim> & F,
+       const double J_tilde_in) {
+               det_F = determinant(F);
+               b_bar = std::pow(det_F, -2.0 / 3.0) * symmetrize(F * transpose(F));
+               J_tilde = J_tilde_in;
+
+               // include a coupled of checks on the input data
+               Assert(det_F > 0, ExcInternalError());
+               Assert(J_tilde > 0, ExcInternalError());
+
+       }
+
+       // Determine the Kirchhoff stress
+       // $\boldsymbol{\tau} = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}$
+       SymmetricTensor<2, dim> get_tau(void) {
+               return get_tau_iso() + get_tau_vol();
+       }
+
+       // The fourth-order elasticity tensor in the spatial setting
+       // $\mathfrak{c}$ is calculated from the SEF $\Psi$ as
+       // $ J \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}$
+       // where
+       // $ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial \mathbf{C} \partial \mathbf{C}}$
+       SymmetricTensor<4, dim> get_Jc(void) const {
+               return get_Jc_vol() + get_Jc_iso();
+       }
+
+       // Derivative of the volumetric free energy wrt $\widetilde{J}$
+       // return $\frac{\partial \Psi_{\text{vol}}(\widetilde{J})}{\partial \widetilde{J}}$
+       double get_dPsi_vol_dJ(void) const {
+               return kappa * (J_tilde - 1.0 / J_tilde);
+       }
+
+       // Second derivative of the volumetric free energy wrt $\widetilde{J}$
+       // We need the following computation explicitly in the tangent
+       // so we make it public.
+       // calculate
+       // $\frac{\partial^2 \Psi_{\textrm{vol}}(\widetilde{J})}{\partial \widetilde{J} \partial \widetilde{J}}$
+       double get_d2Psi_vol_dJ2(void) const {
+               return kappa * (1.0 + 1.0 / (J_tilde * J_tilde));
+       }
 
 protected:
-    // Material properties
-    const double lambda_0; // Lame constant
-    const double mu_0;     // Shear modulus
-    const double kappa_0;  // Bulk modulus
-
-    // We also choose to precalculate and store some frequently used
-    // second and fourth-order tensors.
-    static SymmetricTensor<2, dim> const I;
-    static SymmetricTensor<4, dim> const IxI;
-    static SymmetricTensor<4, dim> const II;
-};
+       // Model properties $\kappa$ and $c_1$
+       const double kappa; // Bulk modulus
+       const double c_1; // neo-Hookean model parameter
+
+       // Model specific data that is convenient to store with the material
+       double det_F;
+       double J_tilde;
+       SymmetricTensor<2, dim> b_bar;
+
+       // Determine the volumetric Kirchhoff stress
+       // $\boldsymbol{\tau}_{\textrm{vol}}$
+       SymmetricTensor<2, dim> get_tau_vol(void) const {
+               // calculate
+               // $\frac{\partial \Psi_{\text{vol}}(\widetilde{J})}{\partial \widetilde{J}}$
+               const double dPsi_vol_dJ = get_dPsi_vol_dJ();
+               // $\boldsymbol{\tau} = J \frac{\partial \Psi_{\textrm{vol}}}{\partial J} \mathbf{I}$
+               return det_F * dPsi_vol_dJ * AdditionalTools::StandardTensors<dim>::I;
+       }
+
+       // Determine the isochoric Kirchhoff stress
+       // $\boldsymbol{\tau}_{\textrm{iso}} = \mathcal{P}:\overline{\boldsymbol{\tau}}$
+       SymmetricTensor<2, dim> get_tau_iso(void) const {
+               return AdditionalTools::StandardTensors<dim>::dev_P * get_tau_bar();
+       }
+
+       // Determine the fictitious Kirchhoff stress
+       SymmetricTensor<2, dim> get_tau_bar(void) const {
+               return 2.0 * c_1 * b_bar;
+       }
+
+       // Calculate the volumetric part of the tangent $J \mathfrak{c}_\textrm{vol}$
+       SymmetricTensor<4, dim> get_Jc_vol(void) const {
+               // now get 
+               // $ \frac{\partial p}{\partial J} = \frac{\partial^2 \Psi_{\textrm{vol}}(J)}{\partial J \partial J}$
+               const double d2Psi_vol_dJ2 = get_d2Psi_vol_dJ2();
+               const double dPsi_vol_dJ = get_dPsi_vol_dJ();
+               const double p_tilde = dPsi_vol_dJ + det_F * d2Psi_vol_dJ2;
+
+               return det_F
+                               * (p_tilde * AdditionalTools::StandardTensors<dim>::IxI
+                                               - (2.0 * dPsi_vol_dJ)
+                                                               * AdditionalTools::StandardTensors<dim>::II);
+       }
+
+       // Calculate the isochoric part of the tangent $J \mathfrak{c}_\textrm{iso}$
+       SymmetricTensor<4, dim> get_Jc_iso(void) const {
+               const SymmetricTensor<2, dim> tau_bar = get_tau_bar();
+               const SymmetricTensor<2, dim> tau_iso = get_tau_iso();
+               const SymmetricTensor<4, dim> tau_iso_x_I = outer_product(tau_iso,
+                               AdditionalTools::StandardTensors<dim>::I);
+               const SymmetricTensor<4, dim> I_x_tau_iso = outer_product(
+                               AdditionalTools::StandardTensors<dim>::I, tau_iso);
+               const SymmetricTensor<4, dim> c_bar = get_c_bar();
+
+               return (2.0 / 3.0) * trace(tau_bar)
+                               * AdditionalTools::StandardTensors<dim>::dev_P
+                               - (2.0 / 3.0) * (tau_iso_x_I + I_x_tau_iso)
+                               + AdditionalTools::StandardTensors<dim>::dev_P * c_bar
+                                               * AdditionalTools::StandardTensors<dim>::dev_P;
+       }
 
-template <int dim> SymmetricTensor<2, dim> const Material_NH<dim>::I   = SymmetricTensor<2, dim> (unit_symmetric_tensor <dim> ());
-template <int dim> SymmetricTensor<4, dim> const Material_NH<dim>::IxI = SymmetricTensor<4, dim> (outer_product (I, I));
-template <int dim> SymmetricTensor<4, dim> const Material_NH<dim>::II  = SymmetricTensor<4, dim> (identity_tensor <dim> ());
+       // Calculate the fictitious elasticity tensor $\overline{\mathfrak{c}}$
+       SymmetricTensor<4, dim> get_c_bar() const {
+               SymmetricTensor<4, dim> c_bar;
+               c_bar = 0.0;
+               return c_bar;
+       }
+};
 
 // @sect3{Quadrature point history}
-// As seen in step-18, the point history class offers
-// a method of storing data defined at the quadrature points.
-// As this method requires the nonlinear stress and
-// material tangents to be evaluated at these points,
-// we used this class to perform these operations.
-//
-// We introduce the multiplicative decomposition of the
-// deformation gradient into a volume-preserving and volume
-// changing component:
-// $ \mathbf{F} = \hat{\mathbf{F}} \bar{\mathbf{F}} $
-// where the volumetric part is
-// $ \hat{\mathbf{F}} = J^{\frac{1}{3}} \mathbf{I} $
-// and the isochoric part is given by
-// $ \bar{\mathbf{F}} = J^{-\frac{1}{3}} \mathbf{F} $
-// . From this, the deviatoric left Cauchy-Green deformation
-// tensor can be defined as
-// $ \bar{\mathbf{B}} = \bar{\mathbf{F}} \bar{\mathbf{F}}^{T} = J^{-\frac{2}{3}} \mathbf{F} \mathbf{F}^{T} $
-//
-// Here we also introduce an additive volumetric-deviatoric split
-// in the material reponse. We can express the governing SEF as
-// $ \phi = \phi_{V} + \phi_{I} $
-// with the result that the Kirchhoff stress is additively
-// decomposed into
-// $ \mathbf{\tau} = \mathbf{\tau}_{V} + \mathbf{\tau}_{I} $
-// as is the tangent matrix
-// $ J\mathbf{C} = J\mathbf{C}_{V} + J\mathbf{C}_{I} $.
-//
-// These quantities are calculated as
-// $  \mathbf{\tau}_{I} = pJ\mathbf{I} $
-// $  \mathbf{\tau}_{V} = \mathcal{P}:\bar{\mathbf{\tau}} $
-// with $ \bar{\mathbf{\tau}} = \mathbf{\tau} \vert_{\mathbf{B} = \bar{\mathbf{B}}} $
-// and the deviatoric tensor $ \mathcal{P} = \mathcal{I} - \frac{1}{3} \mathbf{I} \otimes \mathbf{I} $
-// $  J\mathbf{C}_{I} = pJ(\mathbf{I} \otimes \mathbf{I} - 2 \mathcal{I}) $
-// $  J\mathbf{C}_{V} = \frac{2}{3} tr\left(\bar{\mathbf{\tau}}\right) \mathcal{P} - \frac{2}{3} \left(\mathbf{\tau}_{I}\otimes\mathbf{I} + \mathbf{I}\otimes\mathbf{\tau}_{I} \right) +  \mathcal{P}:\bar{\mathcal{C}}:\mathcal{P} $
-// with $ \bar{\mathcal{C}} = \mathcal{C} \vert_{\mathbf{B} = \bar{\mathbf{B}}} $
-template <int dim>
-class PointHistory
-{
+// As seen in step-18, the <code> PointHistory </code> class offers
+// a method for storing data at the quadrature points.
+// We need to evaluate the Kirchhoff stress $\boldsymbol{\tau}$ and
+// the tangent $J\mathfrak{c}$ at the quadrature points.
+
+template<int dim>
+class PointHistory {
 public:
-    PointHistory (void)
-       :
-          material (NULL),
-          dilatation_n (1.0),
-          pressure_n (0.0)
-    { }
-    virtual ~PointHistory (void) {delete material;}
-
-    // We first create a material object based on the data sent in.
-    // This object could potentially be shared amoung QPH objects
-    // but this could cause data-race issues when assembling the system matrix.
-    void setup_lqp ( Parameters::AllParameters & parameters )
-    {
-       const double lambda = 2.0*parameters.mu*parameters.nu / (1.0-2.0*parameters.nu);
-       material = new Material_NH<dim> (lambda,
-                                        parameters.mu);
-
-        // Initialise all tensors correctly
-        update_values (Tensor <2,dim> (),
-                       0.0,
-                       1.0);
-    }
-
-    // We can update the stored values and stresses based on the current
-    // deformation configuration and pressure and dilation field values
-    void update_values (const Tensor<2, dim> & grad_u_n,
-                       const double & pressure,
-                       const double & dilatation)
-    {
-        // Deformation variables calculated from displacement, displacement gradients
-        static const Tensor < 2, dim> I =  static_cast <Tensor < 2, dim> > (unit_symmetric_tensor <dim> ());
-        const Tensor <2,dim>  F = I + grad_u_n;
-       J     = determinant(F);
-       F_inv = invert(F);
-       B_bar = std::pow(get_J(), -2.0/3.0) * symmetrize ( F* transpose (F) );
-
-       // Store the precalculated pressure and dilatation
-       pressure_n = pressure;
-       dilatation_n = dilatation;
-
-        // Now that all the necessary variables are set, we can update the stress tensors.
-        // Note that T_iso depends on T_bar so it must be calculated afterwards.
-        T_bar = material->get_T (get_J(), get_B_bar());
-        T_iso = dev_P*get_T_bar();
-        T_vol =-get_pressure()*get_J()*I;
-    }
-
-    // We offer and interface to retrieve certain data.
-    // Here are the displacement and strain variables
-    const double & get_dilatation(void) const {return dilatation_n;}
-    const double & get_J (void) const {return J;}
-    const Tensor <2,dim> & get_F_inv (void) const {return F_inv;}
-
-    //, the volumetric SEF quantities
-    double get_dU_dtheta (void) { return material->get_dU_dtheta(get_dilatation()); }
-    double get_d2U_dtheta2 (void) { return material->get_d2U_dtheta2(get_dilatation()); }
-
-    // and stress-based variables. These are used in the material and global
-    // tangent matrix and residual assembly operations so we compute these and
-    // store them.
-    double get_pressure(void) {return pressure_n;}
-    const SymmetricTensor<2, dim> & get_T_iso (void) const {return T_iso;}
-    const SymmetricTensor<2, dim> & get_T_vol (void) const {return T_vol;}
-
-    // Here we provide the local material tangent matrix contribution.
-    // Since they are only used in the tangent matrix assembly process
-    // we compute them as required.
-    // This is the isochoric contribution
-    SymmetricTensor <4,dim> get_C_iso(void)
-    {
-        const double & J = get_J();
-        const SymmetricTensor<2, dim> & B_bar = get_B_bar();
-        const SymmetricTensor<2, dim> & T_iso = get_T_iso();
-
-        const SymmetricTensor <4,dim> T_iso_x_I = outer_product(T_iso, I);
-        const SymmetricTensor <4,dim> I_x_T_iso = outer_product(I, T_iso);
-        const SymmetricTensor <4,dim> C_bar = material->get_JC (J, B_bar);
-
-       return     2.0/3.0*trace(get_T_bar())*dev_P
-               -  2.0/3.0*(T_iso_x_I + I_x_T_iso)
-               +  dev_P*C_bar*dev_P;
-    }
-    // and the volumetric contribution
-    SymmetricTensor <4,dim> get_C_vol(void)
-    {
-        const double & p = get_pressure();
-       const double & J = get_J();
-       return -p*J*(IxI - 2.0*II);
-    }
+       PointHistory(void) :
+                       material(NULL), J_tilde_n(1.0), det_F(1.0), F_inv(
+                                       AdditionalTools::StandardTensors<dim>::I), p_n(0.0), d2Psi_vol_dJ2(
+                                       0.0), dPsi_vol_dJ(0.0) {
+       }
+       virtual ~PointHistory(void) {
+               delete material;
+               material = NULL;
+       }
 
-private:
-    // We specify that each QP has a copy of a material
-    // type in case different materials are used
-    // in different regions of the domain. This also
-    // deals with the issue of preventing data-races during
-    // multi-threading operations when using shared objects.
-    Material_NH <dim>* material;
-
-    // These are all the volume, displacement and strain variables
-    double                  dilatation_n;
-    double                  J;
-    Tensor <2,dim>         F_inv;
-    SymmetricTensor <2,dim> B_bar;
-    SymmetricTensor <2,dim> E;
-    const SymmetricTensor <2,dim> & get_B_bar (void) const {return B_bar;}
-
-    // and the stress-type variables
-    double                  pressure_n;
-    SymmetricTensor<2, dim> T_bar;
-    SymmetricTensor<2, dim> T_iso;
-    SymmetricTensor<2, dim> T_vol;
-    const SymmetricTensor<2, dim> & get_T_bar (void) const {return T_bar;}
-
-    // Some higher-order tensors are frequently used but
-    // remain unchanged. We calculate these once-off
-    // and store them such that they are shared between
-    // all QPH objects.
-    static SymmetricTensor<2, dim> const I;
-    static SymmetricTensor<4, dim> const IxI;
-    static SymmetricTensor<4, dim> const II;
-    static SymmetricTensor<4, dim> const dev_P;
-};
+       // We first create a material object.
+       // This object could, potentially, be shared among QPH objects
+       // but this could cause data-race issues when assembling the system matrix.
+       // ToDo: This issue of the data race needs to be clarified
+       void setup_lqp(Parameters::AllParameters & parameters) {
+
+               // Create an instance of a neo-Hookean material
+               material = new Material_Compressilbe_Neo_Hook_Uncoupled<dim>(
+                               parameters.mu, parameters.nu);
+
+               // Initialise all tensors correctly
+               update_values(Tensor<2, dim>(), 0.0, 1.0);
+       }
+
+       // Update the stored values and stresses based on the current
+       // deformation configuration, pressure $p$ and
+       // dilation $\widetilde{J}$ field values.
+       // The input is the material gradient of the displacement
+       // $\textrm{Grad}\mathbf{u}_{\textrm{n}}$
+       void update_values(const Tensor<2, dim> & Grad_u_n, const double p
+                       ,const double J_tilde) {
+               // Store the calculated pressure $p$
+               // and dilatation $\widetilde{J}$
+               p_n = p;
+               J_tilde_n = J_tilde;
+
+               // Various deformation gradient $\mathbf{F}$ from the
+               // displacement gradient $\textrm{Grad}\mathbf{u}$, i.e.
+               // $\mathbf{F}(\mathbf{u}) = \mathbf{I} + \textrm{Grad} \mathbf{u}$
+               static const Tensor<2, dim> I =
+                               static_cast<Tensor<2, dim> >(AdditionalTools::StandardTensors<
+                                               dim>::I);
+               const Tensor<2, dim> F = I + Grad_u_n;
+
+
+
+               // We use the inverse of $\mathbf{F}$ frequently so we store it
+               F_inv = invert(F);
+               // as well as the determinant $\textrm{det}\mathbf{F}$
+               det_F = determinant(F);
+
+               std::cout << det_F << "\t" << J_tilde << std::endl;
+
+               // Now we update the material model with the new deformation measures
+               material->update_material_data(F, J_tilde);
 
-template <int dim> SymmetricTensor<2,dim> const PointHistory<dim>::I
-= SymmetricTensor<2,dim> (unit_symmetric_tensor <dim> ());
-template <int dim> SymmetricTensor<4,dim> const PointHistory<dim>::IxI
-= SymmetricTensor<4,dim> (outer_product (I, I));
-template <int dim> SymmetricTensor<4,dim> const PointHistory<dim>::II
-= SymmetricTensor<4,dim> (identity_tensor <dim> ());
-template <int dim> SymmetricTensor<4,dim> const PointHistory<dim>::dev_P
-= SymmetricTensor<4,dim> (II - 1.0/3.0*IxI);
+               // The material has been updated so we now calculate the
+               // Kirchhoff stress $\mathbf{\tau}$ and the tangent $J\mathfrak{c}$
+               tau = material->get_tau();
 
+               Jc = material->get_Jc();
+               dPsi_vol_dJ = material->get_dPsi_vol_dJ();
+               d2Psi_vol_dJ2 = material->get_d2Psi_vol_dJ2();
+
+       }
+
+       // We offer an interface to retrieve certain data.
+       // Here are the kinematic variables
+       double get_J_tilde(void) const {
+               return J_tilde_n;
+       }
+       double get_det_F(void) const {
+               return det_F;
+       }
+       Tensor<2, dim> get_F_inv(void) const {
+               return F_inv;
+       }
+
+       // and the kinetic variables.
+       // These are used in the material and global
+       // tangent matrix and residual assembly operations
+       // so we compute these and store them.
+       double get_p(void) const {
+               return p_n;
+       }
+       SymmetricTensor<2, dim> get_tau(void) const {
+               return tau;
+       }
+
+       double get_dPsi_vol_dJ(void) const {
+               return dPsi_vol_dJ;
+       }
+
+       double get_d2Psi_vol_dJ2(void) const {
+               return d2Psi_vol_dJ2;
+       }
+
+       // and finally the tangent
+       SymmetricTensor<4, dim> get_Jc(void) const {
+               return Jc;
+       }
+
+private:
+       // We specify that each QP has a copy of a material
+       // type in case different materials are used
+       // in different regions of the domain.
+       // This also
+       // deals with the issue of preventing data-races during
+       // multi-threading operations when using shared objects.
+       Material_Compressilbe_Neo_Hook_Uncoupled<dim>* material;
+
+       // These are all the volume, displacement and strain variables
+       double J_tilde_n;
+       double det_F;
+       Tensor<2, dim> F_inv;
+
+       // and the stress-type variables
+       double p_n;
+       SymmetricTensor<2, dim> tau;
+       double d2Psi_vol_dJ2;
+       double dPsi_vol_dJ;
+
+       // and the tangent
+       SymmetricTensor<4, dim> Jc;
+};
 
 // @sect3{Quasi-static quasi-incompressible finite-strain solid}
-template <int dim>
-class Solid
-{
+template<int dim>
+class Solid {
 public:
-    Solid (const std::string & input_file);
-    virtual ~Solid (void);
-    void run (void);
+       Solid(const std::string & input_file);
+       virtual ~Solid(void);
+       void run(void);
 
 private:
 
-    // Threaded building-blocks data structures
-    struct PerTaskData_K;
-    struct ScratchData_K;
-    struct PerTaskData_F;
-    struct ScratchData_F;
-    struct PerTaskData_SC;
-    struct ScratchData_SC;
-    struct PerTaskData_UQPH;
-    struct ScratchData_UQPH;
-
-    // Build the grid
-    void make_grid (void);
-
-    // Setup the Finite Element system to be solved
-    void system_setup (void);
-    void determine_component_extractors(void);
-
-    // Assemble the system and right hand side matrices using multi-threading
-    void assemble_system_K          (void);
-    void assemble_system_K_one_cell (const typename DoFHandler<dim>::active_cell_iterator & cell,
-                                    ScratchData_K & scratch,
-                                    PerTaskData_K & data);
-    void copy_local_to_global_K     (const PerTaskData_K & data);
-    void assemble_system_F          (void);
-    void assemble_system_F_one_cell (const typename DoFHandler<dim>::active_cell_iterator & cell,
-                                    ScratchData_F & scratch,
-                                    PerTaskData_F & data);
-    void copy_local_to_global_F     (const PerTaskData_F & data);
-    void assemble_SC                (void);
-    void assemble_SC_one_cell       (const typename DoFHandler<dim>::active_cell_iterator & cell,
-                                     ScratchData_SC & scratch,
-                                     PerTaskData_SC & data);
-    void copy_local_to_global_SC    (const PerTaskData_SC & data);
-    // Apply Dirichlet boundary values
-    void make_constraints (const int & it_nr,
-                          ConstraintMatrix & constraints);
-
-    // Create and update the quadrature points stress and strain values
-    void setup_qph(void);
-    void update_qph_incremental ( const BlockVector <double> & solution_delta );
-    void update_qph_incremental_one_cell (const typename DoFHandler<dim>::active_cell_iterator & cell,
-                                         ScratchData_UQPH & scratch,
-                                         PerTaskData_UQPH & data);
-    void copy_local_to_global_UQPH    (const PerTaskData_UQPH & data) {}
-
-    // Solve for the displacement using a Newton-Rhapson method
-    void solve_nonlinear_timestep (BlockVector <double> & solution_delta);
-    std::pair <unsigned int, double> solve_linear_system (BlockVector <double> & newton_update);
-
-    // Solution retrieval
-    BlockVector <double> get_solution_total (const BlockVector <double> & solution_delta);
-
-    // Postprocessing and writing data to file
-    void output_results(void);
-
-    // A collection of the parameters used to describe the problem setup
-    Parameters::AllParameters parameters;
-
-    // Description of the geometry on which the problem is solved
-    Triangulation<dim> triangulation;
-
-    // Keep track of the current time and the time spent evaluating certain functions
-    Time time;
-    TimerOutput timer;
-
-    // A storage object for quadrature point information
-    std::vector< PointHistory <dim> > quadrature_point_history;
-
-    // A desciption of the finite-element system including the displacement polynomial degree,
-    // the degree-of-freedom handler, number of dof's per cell and the extractor objects used
-    // to retrieve information from the solution vectors
-    const unsigned int  degree;
-    const FESystem<dim> fe;
-    DoFHandler<dim>     dof_handler_ref;
-    unsigned int dofs_per_cell;
-    const FEValuesExtractors::Vector u_fe;
-    const FEValuesExtractors::Scalar p_fe;
-    const FEValuesExtractors::Scalar t_fe;
-
-    // Description of how the block-system is arranged
-    // There are 3 blocks, the first contains a vector DOF
-    // while the other two describe scalar DOFs.
-    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 t_component = dim + 1;
-
-    enum {u_dof=0 , p_dof, t_dof};
-    std::vector<unsigned int> dofs_per_block;
-    std::vector<unsigned int> element_indices_u;
-    std::vector<unsigned int> element_indices_p;
-    std::vector<unsigned int> element_indices_t;
-
-    // Rules for gauss-quadrature on both the cell and faces. The
-    // number of quadrature points on both cells and faces is
-    // recorded.
-    QGauss<dim> qf_cell;
-    QGauss<dim-1> qf_face;
-    unsigned int n_q_points;
-    unsigned int n_q_points_f;
-
-    // Objects that store the converged solution and residual vectors,
-    // as well as the tangent matrix. There is a ConstraintMatrix object
-    // used to keep track of constraints for the nonlinear problem.
-    ConstraintMatrix constraints;
-    BlockSparsityPattern sparsity_pattern;
-    BlockSparseMatrix <double> tangent_matrix;
-    BlockVector <double> residual;
-    BlockVector <double> solution_n;
-
-    // Then define a number of variables to store residual and update
-    // norms and normalisation factors.
-    struct Errors
-    {
-        Errors (void) : norm(1.0), u (1.0), p(1.0), t(1.0) {}
-        double norm,u,p,t;
-        void reset (void) {norm = 1.0; u = 1.0; p = 1.0; t = 1.0;}
-        void normalise (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.t != 0.0) t /= rhs.t;
-        }
-    }
-    error_residual, error_residual_0, error_residual_norm,
-    error_update, error_update_0, error_update_norm;
-
-    // Methods to calculate error measures
-    void get_error_residual (Errors & error_residual);
-    void get_error_update (const BlockVector <double> & newton_update,
-                           Errors & error_update);
-    double get_error_dil (void);
-
-    // Print information to screen
-    void print_conv_header (void);
-    void print_conv_footer (void);
+       // Threaded building-blocks data structures:
+       // for the tangent matrix
+       struct PerTaskData_K;
+       struct ScratchData_K;
+       // for the right-hand side
+       struct PerTaskData_RHS;
+       struct ScratchData_RHS;
+       // for the static-condensation
+       struct PerTaskData_SC;
+       struct ScratchData_SC;
+       // for the updating of the quadrature points
+       struct PerTaskData_UQPH;
+       struct ScratchData_UQPH;
+
+       // Build the grid
+       void make_grid(void);
+
+       // Setup the Finite Element system to be solved
+       void system_setup(void);
+       void determine_component_extractors(void);
+
+       // Assemble the system and right hand side matrices using multi-threading
+       void assemble_system_K(void);
+       void assemble_system_K_one_cell(
+                       const typename DoFHandler<dim>::active_cell_iterator & cell,
+                       ScratchData_K & scratch, PerTaskData_K & data);
+       void copy_local_to_global_K(const PerTaskData_K & data);
+       void assemble_system_rhs(void);
+       void assemble_system_rhs_one_cell(
+                       const typename DoFHandler<dim>::active_cell_iterator & cell,
+                       ScratchData_RHS & scratch, PerTaskData_RHS & data);
+       void copy_local_to_global_rhs(const PerTaskData_RHS & data);
+       void assemble_sc(void);
+       void assemble_sc_one_cell(
+                       const typename DoFHandler<dim>::active_cell_iterator & cell,
+                       ScratchData_SC & scratch, PerTaskData_SC & data);
+       void copy_local_to_global_sc(const PerTaskData_SC & data);
+       // Apply Dirichlet boundary values
+       void make_constraints(const int & it_nr, ConstraintMatrix & constraints);
+
+       // Create and update the quadrature points stress and strain values
+       void setup_qph(void);
+       void update_qph_incremental(const BlockVector<double> & solution_delta);
+       void update_qph_incremental_one_cell(
+                       const typename DoFHandler<dim>::active_cell_iterator & cell,
+                       ScratchData_UQPH & scratch, PerTaskData_UQPH & data);
+       void copy_local_to_global_UQPH(const PerTaskData_UQPH & data) {
+       }
+
+       // Solve for the displacement using a Newton-Rhapson method
+       void solve_nonlinear_timestep(BlockVector<double> & solution_delta);
+       std::pair<unsigned int, double> solve_linear_system(
+                       BlockVector<double> & newton_update);
+
+       // Solution retrieval
+       BlockVector<double> get_solution_total(
+                       const BlockVector<double> & solution_delta);
+
+       // Post-processing and writing data to file
+       void output_results(void);
+
+       // A collection of the parameters used to describe the problem setup
+       Parameters::AllParameters parameters;
+
+       // Description of the geometry on which the problem is solved
+       Triangulation<dim> triangulation;
+
+       // Keep track of the current time and the time spent evaluating certain functions
+       Time time;
+       TimerOutput timer;
+
+       // A storage object for quadrature point information
+       std::vector<PointHistory<dim> > quadrature_point_history;
+
+       // A description of the finite-element system including the displacement polynomial degree,
+       // the degree-of-freedom handler, number of dof's per cell and the extractor objects used
+       // to retrieve information from the solution vectors
+       const unsigned int degree;
+       const FESystem<dim> fe;
+       DoFHandler<dim> dof_handler_ref;
+       unsigned int dofs_per_cell;
+       const FEValuesExtractors::Vector u_fe;
+       const FEValuesExtractors::Scalar p_fe;
+       const FEValuesExtractors::Scalar J_fe;
+
+       // Description of how the block-system is arranged
+       // There are 3 blocks, the first contains a vector DOF $\mathbf{u}$
+       // while the other two describe scalar DOFs, $p$ and $\widetilde{J}$.
+       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, J_dof
+       };
+       std::vector<unsigned int> dofs_per_block;
+       std::vector<unsigned int> element_indices_u;
+       std::vector<unsigned int> element_indices_p;
+       std::vector<unsigned int> element_indices_J;
+
+       // Rules for Gauss-quadrature on both the cell and faces. The
+       // number of quadrature points on both cells and faces is
+       // recorded.
+       QGauss<dim> qf_cell;
+       QGauss<dim - 1> qf_face;
+       unsigned int n_q_points;
+       unsigned int n_q_points_f;
+
+       // Objects that store the converged solution and right-hand side vectors,
+       // as well as the tangent matrix. There is a ConstraintMatrix object
+       // used to keep track of constraints.
+       ConstraintMatrix constraints;
+       BlockSparsityPattern sparsity_pattern;
+       BlockSparseMatrix<double> tangent_matrix;
+       BlockVector<double> system_rhs;
+       BlockVector<double> solution_n;
+
+       // Then define a number of variables to store norms and update
+       // norms and normalisation factors.
+       struct Errors {
+               Errors(void) :
+                               norm(1.0), u(1.0), p(1.0), J(1.0) {
+               }
+               double norm, u, p, J;
+               void reset(void) {
+                       norm = 1.0;
+                       u = 1.0;
+                       p = 1.0;
+                       J = 1.0;
+               }
+               void normalise(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;
+               }
+       } error_residual, error_residual_0, error_residual_norm, error_update,
+                       error_update_0, error_update_norm;
+
+       // Methods to calculate error measures
+       void get_error_residual(Errors & error_residual);
+       void get_error_update(const BlockVector<double> & newton_update,
+                       Errors & error_update);
+       double get_error_dil(void);
+
+       // Print information to screen
+       void print_conv_header(void);
+       void print_conv_footer(void);
 };
 
 // @sect3{Implementation of the <code>Solid</code> class}
 
 // @sect4{Public interface}
-// We initialise the the solid class using data extracted
+// We initialise the Solid class using data extracted
 // from the parameter file.
-template <int dim>
-Solid<dim>::Solid (const std::string & input_file)
-    :
-      parameters (input_file),
-      triangulation (Triangulation<dim>::maximum_smoothing),
-      time (parameters.end_time,
-            parameters.delta_t),
-      timer (std::cout,
-          TimerOutput::summary,
-          TimerOutput::wall_times),
-      degree (parameters.poly_degree),
-      // The Finite Element System is composed of dim continuous
-      // displacment DOFs and linear discontinuous pressure and
-      // dilatation DOFs. In an attempt to satisfy the LBB conditions,
-      // we setup a Q(n)-P(n-1)-P(n-1) system. Q2-P1 element satisfy
-      // this condition, while Q1-P0 elements do not. However, it
-      // has been shown that they demonstrate good convergence
-      // characteristics nonetheless.
-      fe (FE_Q<dim>(parameters.poly_degree), dim,
-          FE_DGPMonomial<dim>(parameters.poly_degree-1), 1,
-          FE_DGPMonomial<dim>(parameters.poly_degree-1), 1),
-      dof_handler_ref (triangulation),
-      u_fe (first_u_component),
-      p_fe (p_component),
-      t_fe (t_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();
-    dofs_per_cell = fe.dofs_per_cell;
-    determine_component_extractors();
+template<int dim>
+Solid<dim>::Solid(const std::string & input_file) :
+               parameters(input_file), triangulation(
+                               Triangulation<dim>::maximum_smoothing), time(
+                               parameters.end_time, parameters.delta_t), timer(std::cout,
+                               TimerOutput::summary, TimerOutput::wall_times), degree(
+                               parameters.poly_degree),
+               // The Finite Element System is composed of dim continuous
+               // displacement DOFs, and discontinuous pressure and
+               // dilatation DOFs. In an attempt to satisfy the LBB conditions,
+               // we setup a Q(n)-P(n-1)-P(n-1) system. Q2-P1-P1 elements satisfy
+               // this condition, while Q1-P0-P0 elements do not. However, it
+               // has been shown that the latter demonstrate good convergence
+               // characteristics nonetheless.
+               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_ref(triangulation), 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();
+       dofs_per_cell = fe.dofs_per_cell;
+       determine_component_extractors();
 }
 
-// The class destructor simply needs to clear the data held by the DOFHandler
-template <int dim>
-Solid<dim>::~Solid (void)
-{
-    dof_handler_ref.clear ();
+// The class destructor simply clears the data held by the DOFHandler
+template<int dim>
+Solid<dim>::~Solid(void) {
+       dof_handler_ref.clear();
 }
 
-// In solving the quasti-static problem, the time
+// In solving the quasi-static problem, the time
 // becomes a loading parameter. We choose to increment
-// time linearly using a constant timestep size.
-template <int dim>
-void Solid<dim>::run (void)
-{
-    // After preprocessing, we output the initial grid
-    // before starting the simulation proper.
-    make_grid ();
-    system_setup ();
-    output_results ();
-    time.increment();
-
-    BlockVector <double> solution_delta (dofs_per_block);
-    solution_delta.collect_sizes ();
-
-    while (time.current() < time.end()) {
-        // We need to reset the solution update
-        // for this timestep
-       solution_delta = 0.0;
-
-       // Solve the current timestep and update total
-       // solution vector
-       solve_nonlinear_timestep (solution_delta);
-       solution_n += solution_delta;
-       output_results ();
-
-        time.increment();
-    }
+// time linearly using a constant time step size.
+template<int dim>
+void Solid<dim>::run(void) {
+       // After preprocessing, we output the initial grid
+       // before starting the simulation proper.
+       make_grid();
+       system_setup();
+       output_results();
+       time.increment();
+
+       // Here we define 
+       // $\varDelta \mathbf{\Xi}:= \{\varDelta \mathbf{u},\varDelta p, \varDelta \widetilde{J} \}$.
+       BlockVector<double> solution_delta(dofs_per_block);
+       solution_delta.collect_sizes();
+
+       // Now we loop over the time domain
+       while (time.current() < time.end()) {
+               // We need to reset the solution update
+               // for this time step
+               solution_delta = 0.0;
+
+               // Solve the current time step and update total
+               // solution vector
+               solve_nonlinear_timestep(solution_delta);
+               // $\varDelta \mathbf{\Xi}_{\textrm{n}} = \varDelta \mathbf{\Xi}_{\textrm{n-1}} + \varDelta \mathbf{\Xi}$
+               solution_n += solution_delta;
+               output_results();
+
+               time.increment();
+       }
 }
 
 // @sect3{Private interface}
 
 // @sect4{Threaded-building-blocks structures}
-// We choose to use TBB to perform as many computationally intensive
+// We use TBB to perform as many computationally intensive
 // distributed tasks as possible. In particular, we assemble the
-// tangent matrix and residual vector, assemble the static
-// condensation contributions and update data stored
-// at the quadrature points.
+// tangent matrix and residual vector, the static
+// condensation contributions, and update data stored
+// at the quadrature points using TBB.
 
 // Firstly we deal with the tangent matrix assembly structures.
-// The PerTaskData object stores local contributions.
-template <int dim>
-struct Solid<dim>::PerTaskData_K
-{
-    FullMatrix<double>          cell_matrix;
-    std::vector<unsigned int>   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 (void) {
-        cell_matrix = 0.0;
-    }
+// The PerTaskData object stores local contributions. 
+template<int dim>
+struct Solid<dim>::PerTaskData_K {
+       FullMatrix<double> cell_matrix;
+       std::vector<unsigned int> 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(void) {
+               cell_matrix = 0.0;
+       }
 };
 // while the ScratchData object stores the larger objects
 // such as the shape-function values object and a shape function
-// values and gradient vector which we will precompute later.
-template <int dim>
-struct Solid<dim>::ScratchData_K
-{
-    FEValues <dim> fe_values_ref;
-
-    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_ref   (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_ref ( rhs.fe_values_ref.get_fe(),
-            rhs.fe_values_ref.get_quadrature(),
-            rhs.fe_values_ref.get_update_flags() ),
-        Nx (rhs.Nx),
-        grad_Nx (rhs.grad_Nx),
-        symm_grad_Nx (rhs.symm_grad_Nx)
-    {  }
-
-    void reset (void) {
-        for (unsigned int q_point=0; q_point < grad_Nx.size(); ++q_point) {
-            for (unsigned int k=0; k < Nx.size(); ++k) {
-                Nx[q_point][k] = 0.0;
-                grad_Nx[q_point][k] = 0.0;
-                symm_grad_Nx[q_point][k] = 0.0;
-            }
-        }
-    }
+// gradient and symmetric gradient vector which we will precompute later.
+template<int dim>
+struct Solid<dim>::ScratchData_K {
+       FEValues<dim> fe_values_ref;
+
+       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_ref(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_ref(rhs.fe_values_ref.get_fe(),
+                                       rhs.fe_values_ref.get_quadrature(),
+                                       rhs.fe_values_ref.get_update_flags()), Nx(rhs.Nx), grad_Nx(
+                                       rhs.grad_Nx), symm_grad_Nx(rhs.symm_grad_Nx) {
+       }
+
+       void reset(void) {
+               for (unsigned int q_point = 0; q_point < grad_Nx.size(); ++q_point) {
+                       for (unsigned int k = 0; k < Nx.size(); ++k) {
+                               Nx[q_point][k] = 0.0;
+                               grad_Nx[q_point][k] = 0.0;
+                               symm_grad_Nx[q_point][k] = 0.0;
+                       }
+               }
+       }
 
 };
 
-// Next are the same data structures used for the residual assembly.
+// Next are the same data structures used for the
+// right-hand side assembly.
 // The PerTaskData object again stores local contributions
-template <int dim>
-struct Solid<dim>::PerTaskData_F
-{
-    Vector<double>              cell_rhs;
-    std::vector<unsigned int>   local_dof_indices;
+template<int dim>
+struct Solid<dim>::PerTaskData_RHS {
+       Vector<double> cell_rhs;
+       std::vector<unsigned int> local_dof_indices;
 
-    PerTaskData_F (const unsigned int dofs_per_cell)
-        :
-          cell_rhs           (dofs_per_cell),
-          local_dof_indices  (dofs_per_cell)
-    { }
+       PerTaskData_RHS(const unsigned int dofs_per_cell) :
+                       cell_rhs(dofs_per_cell), local_dof_indices(dofs_per_cell) {
+       }
 
-    void reset (void) { cell_rhs = 0.0; }
+       void reset(void) {
+               cell_rhs = 0.0;
+       }
 };
 // and the ScratchData object the shape function object
 // and precomputed values vector
-template <int dim>
-struct Solid<dim>::ScratchData_F
-{
-    FEValues <dim>     fe_values_ref;
-    FEFaceValues <dim> fe_face_values_ref;
-
-    std::vector < std::vector< double > > Nx;
-    std::vector < std::vector< SymmetricTensor<2, dim> > > symm_grad_Nx;
-
-    // Solution data
-    std::vector< std::vector<Tensor <1,dim> > > solution_grads;
-
-    ScratchData_F ( 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_ref   (fe_cell,
-              qf_cell,
-              uf_cell),
-          fe_face_values_ref   (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_F ( const ScratchData_F & rhs )
-        :
-          fe_values_ref ( rhs.fe_values_ref.get_fe(),
-              rhs.fe_values_ref.get_quadrature(),
-              rhs.fe_values_ref.get_update_flags() ),
-          fe_face_values_ref ( rhs.fe_face_values_ref.get_fe(),
-              rhs.fe_face_values_ref.get_quadrature(),
-              rhs.fe_face_values_ref.get_update_flags() ),
-          Nx (rhs.Nx),
-          symm_grad_Nx (rhs.symm_grad_Nx)
-    {  }
-
-    void reset (void) {
-        for (unsigned int q_point=0; q_point < symm_grad_Nx.size(); ++q_point) {
-            for (unsigned int k=0; k < symm_grad_Nx[q_point].size(); ++k) {
-                Nx[q_point][k] = 0.0;
-                symm_grad_Nx[q_point][k] = 0.0;
-            }
-        }
-    }
+template<int dim>
+struct Solid<dim>::ScratchData_RHS {
+       FEValues<dim> fe_values_ref;
+       FEFaceValues<dim> fe_face_values_ref;
+
+       std::vector<std::vector<double> > Nx;
+       std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
+
+       // Solution data
+       std::vector<std::vector<Tensor<1, dim> > > solution_grads;
+
+       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_ref(fe_cell, qf_cell, uf_cell), fe_face_values_ref(
+                                       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_ref(rhs.fe_values_ref.get_fe(),
+                                       rhs.fe_values_ref.get_quadrature(),
+                                       rhs.fe_values_ref.get_update_flags()), fe_face_values_ref(
+                                       rhs.fe_face_values_ref.get_fe(),
+                                       rhs.fe_face_values_ref.get_quadrature(),
+                                       rhs.fe_face_values_ref.get_update_flags()), Nx(rhs.Nx), symm_grad_Nx(
+                                       rhs.symm_grad_Nx) {
+       }
+
+       void reset(void) {
+               for (unsigned int q_point = 0; q_point < symm_grad_Nx.size();
+                               ++q_point) {
+                       for (unsigned int k = 0; k < symm_grad_Nx[q_point].size(); ++k) {
+                               Nx[q_point][k] = 0.0;
+                               symm_grad_Nx[q_point][k] = 0.0;
+                       }
+               }
+       }
 
 };
 
-// Here we define structures to assemble the static condensation contributions.
-// As the operations are matrix-based, we need to setup a number of matrices
-// to store the local contributions from a number of the tangent matrix subblocks.
+// Here 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 $p$ and $\widetilde{J}$
+// 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.
-template <int dim>
-struct Solid<dim>::PerTaskData_SC
-{
-    FullMatrix<double>          cell_matrix;
-    std::vector<unsigned int>   local_dof_indices;
-
-    // Calculation matrices (auto resized)
-    FullMatrix<double> K_orig;
-    FullMatrix<double> K_pu;
-    FullMatrix<double> K_pt;
-    FullMatrix<double> K_tt;
-    // Calculation matrices (manual resized)
-    FullMatrix<double> K_pt_inv;
-    FullMatrix<double> K_tt_inv;
-    FullMatrix<double> K_con;
-    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_t)
-        :
-          cell_matrix (dofs_per_cell,
-                       dofs_per_cell),
-          local_dof_indices  (dofs_per_cell),
-          K_pt_inv (n_t, n_p),
-          K_tt_inv (n_t, n_t),
-          K_con (n_u, n_u),
-          A (n_t, n_u),
-          B (n_t, n_u),
-          C (n_p, n_u)
-    {  }
-
-    // Choose not to reset any data as the matrix extraction and
-    // replacement tools will take care of this
-    void reset(void) { }
+template<int dim>
+struct Solid<dim>::PerTaskData_SC {
+       FullMatrix<double> cell_matrix;
+       std::vector<unsigned int> local_dof_indices;
+
+       // Calculation matrices (auto resized)
+       FullMatrix<double> k_orig;
+       FullMatrix<double> k_pu;
+       FullMatrix<double> k_pJ;
+       FullMatrix<double> k_JJ;
+       // Calculation matrices (manual resized)
+       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_pJ_inv(n_J, n_p), k_bbar(n_u, n_u), A(n_J,
+                                       n_u), B(n_J, n_u), C(n_p, n_u) {
+       }
+
+       // Choose not to reset any data as the matrix extraction and
+       // replacement tools will take care of this
+       void reset(void) {
+       }
 };
 // The ScratchData object is not strictly necessary for the
 // operations we wish to perform, but it still needs to be defined for the
-// current implementation of TBB in deal.II.So we creatre a dummy struct for this purpose.
-template <int dim>
-struct Solid<dim>::ScratchData_SC
-{
-    ScratchData_SC (void) { }
-    ScratchData_SC (const ScratchData_SC & rhs) { }
-    void reset (void) { }
+// current implementation of TBB in deal.II. 
+// So we create a dummy struct for this purpose.
+template<int dim>
+struct Solid<dim>::ScratchData_SC {
+       ScratchData_SC(void) {
+       }
+       ScratchData_SC(const ScratchData_SC & rhs) {
+       }
+       void reset(void) {
+       }
 };
 
 // And finally we define the structures to assist with updating the quadrature
 // point information. Similar to the SC assembly process, we choose not to use
 // the PerTaskData object to store any information but must define one nonetheless.
-template <int dim>
-struct Solid<dim>::PerTaskData_UQPH
-{
-    PerTaskData_UQPH (void) { }
-    void reset(void) { }
+template<int dim>
+struct Solid<dim>::PerTaskData_UQPH {
+       PerTaskData_UQPH(void) {
+       }
+       void reset(void) {
+       }
 };
-// The ScratchData object will be used to store a alias fort the solution vector
+// 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.
-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_t_total;
-
-    FEValues<dim> fe_values_ref;
-
-    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_t_total (qf_cell.size()),
-          fe_values_ref (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_t_total (rhs.solution_values_t_total),
-          fe_values_ref (rhs.fe_values_ref.get_fe(),
-                            rhs.fe_values_ref.get_quadrature(),
-                            rhs.fe_values_ref.get_update_flags())
-    { }
-
-    void reset (void)
-    {
-        // Is this necessary? Won't the call to fe_values.get_gradient overwrite this data?
-        for (unsigned int q=0; q < qf_cell.size(); ++q)
-        {
-            solution_grads_u_total[q] = 0.0;
-            solution_values_p_total[q] = 0.0;
-            solution_values_t_total[q] = 0.0;
-        }
-    }
+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_ref;
+
+       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_ref(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_ref(
+                                       rhs.fe_values_ref.get_fe(),
+                                       rhs.fe_values_ref.get_quadrature(),
+                                       rhs.fe_values_ref.get_update_flags()) {
+       }
+
+       void reset(void) {
+               // ToDo: Is this necessary? Won't the call to fe_values.get_gradient overwrite this data?
+               for (unsigned int q = 0; q < qf_cell.size(); ++q) {
+                       solution_grads_u_total[q] = 0.0;
+                       solution_values_p_total[q] = 0.0;
+                       solution_values_J_total[q] = 0.0;
+               }
+       }
 };
 
 // @sect4{Solid::make_grid}
-// Here we create the grid on which the minimisation problem is to be solved.
-template <int dim>
-void Solid<dim>::make_grid (void)
-{
-    // Create a unit cube with each face given a boundary ID number
-    GridGenerator::hyper_rectangle ( triangulation,
-                                    Point<dim> (0.0, 0.0, 0.0),
-                                    Point<dim> (1.0, 1.0, 1.0),
-                                    true );
-    GridTools::scale (parameters.scale,
-                      triangulation);
-
-    // The grid must be refined at least once for the indentation problem
-    if (parameters.global_refinement == 0)
-        triangulation.refine_global (1);
-    else
-        triangulation.refine_global (parameters.global_refinement);
-
-    // 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
-    typename Triangulation<dim>::active_cell_iterator
-           cell = triangulation.begin_active(),
-           endc = triangulation.end();
-    for (; cell!=endc; ++cell)
-    {
-        if (cell->at_boundary() == true) {
-           for (unsigned int face=0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
-               // Find faces on the +y surface
-               if (   cell->face(face)->at_boundary() == true
-                      && cell->face(face)->center()[2] == 1.0*parameters.scale)
-               {
-                   if (   cell->face(face)->center()[0] < 0.5*parameters.scale
-                          && cell->face(face)->center()[1] < 0.5*parameters.scale)
-                   {
-                       cell->face(face)->set_boundary_indicator (6); // Set a new boundary id on a patch
-                   }
+// Here we create the triangulation of the domain
+template<int dim>
+void Solid<dim>::make_grid(void) {
+       // Create a unit cube with each face given a boundary ID number
+       GridGenerator::hyper_rectangle(triangulation, Point<dim>(0.0, 0.0, 0.0),
+                       Point<dim>(1.0, 1.0, 1.0), true);
+       GridTools::scale(parameters.scale, triangulation);
+
+       // The grid must be refined at least once for the indentation problem
+       if (parameters.global_refinement == 0)
+               triangulation.refine_global(1);
+       else
+               triangulation.refine_global(parameters.global_refinement);
+
+       // 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
+       typename Triangulation<dim>::active_cell_iterator cell =
+                       triangulation.begin_active(), endc = triangulation.end();
+       for (; cell != endc; ++cell) {
+               if (cell->at_boundary() == true) {
+                       for (unsigned int face = 0;
+                                       face < GeometryInfo<dim>::faces_per_cell; ++face) {
+                               // Find faces on the +y surface
+                               if (cell->face(face)->at_boundary() == true
+                                               && cell->face(face)->center()[2]
+                                                               == 1.0 * parameters.scale) {
+                                       if (cell->face(face)->center()[0] < 0.5 * parameters.scale
+                                                       && cell->face(face)->center()[1]
+                                                                       < 0.5 * parameters.scale) {
+                                               cell->face(face)->set_boundary_indicator(6); // Set a new boundary id on a patch
+                                       }
+                               }
+                       }
                }
-           }
        }
-    }
 }
 
 // @sect4{Solid::system_setup}
 // Next we describe how the FE system is setup.
-template <int dim>
-void Solid<dim>::system_setup (void)
-{
-    timer.enter_subsection ("Setup system");
-
-    // We first describe 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.
-    std::vector<unsigned int> block_component (n_components, u_dof); // Displacement
-    block_component[p_component] = p_dof; // Pressure
-    block_component[t_component] = t_dof; // Dilatation
-
-    // DOF handler is then initialised and we renumber the grid in an
-    // efficient manner. We also record the number of DOF's per block.
-    dof_handler_ref.distribute_dofs (fe);
-    DoFRenumbering::Cuthill_McKee (dof_handler_ref);
-    DoFRenumbering::component_wise (dof_handler_ref,
-                                    block_component);
-    DoFTools::count_dofs_per_block (dof_handler_ref,
-                                    dofs_per_block,
-                                    block_component);
-
-    std::cout
-           << "Triangulation:"
-           << "\n\t Number of active cells: " << triangulation.n_active_cells()
-           << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs()
-           << std::endl;
-
-    // Setup the sparsity pattern and tangent matrix
-    tangent_matrix.clear ();
-    {
-       const unsigned int n_dofs_u = dofs_per_block[u_dof];
-       const unsigned int n_dofs_p = dofs_per_block[p_dof];
-       const unsigned int n_dofs_t = dofs_per_block[t_dof];
-
-        BlockCompressedSimpleSparsityPattern csp (n_blocks,
-                                                  n_blocks);
-
-        csp.block(u_dof,u_dof).reinit (n_dofs_u, n_dofs_u);
-        csp.block(u_dof,p_dof).reinit (n_dofs_u, n_dofs_p);
-        csp.block(u_dof,t_dof).reinit (n_dofs_u, n_dofs_t);
-
-        csp.block(p_dof,u_dof).reinit (n_dofs_p, n_dofs_u);
-        csp.block(p_dof,p_dof).reinit (n_dofs_p, n_dofs_p);
-        csp.block(p_dof,t_dof).reinit (n_dofs_p, n_dofs_t);
-
-        csp.block(t_dof,u_dof).reinit (n_dofs_t, n_dofs_u);
-        csp.block(t_dof,p_dof).reinit (n_dofs_t, n_dofs_p);
-        csp.block(t_dof,t_dof).reinit (n_dofs_t, n_dofs_t);
-        csp.collect_sizes();
-
-        // The global system matrix will have the following structure
-        //      | K'_uu |  K_up  |     0     |         | dU_u |         | dR_u |
-        // K =  | K_pu  |    0   |   K_pt^-1 | , dU =  | dU_p | , dR =  | dR_p |
-        //      |   0   |  K_tp  |   K_tt    |         | dU_t |         | dR_t |
-        // We optimise the sparsity pattern to reflect this structure
-        // and prevent unnecessary data creation for the right-diagonal
-        // block components.
-        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 == t_component) )
-                     || ( (ii == t_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_ref, coupling, csp, constraints, false);
-        sparsity_pattern.copy_from (csp);
-    }
-    
-    tangent_matrix.reinit (sparsity_pattern);
-
-    // Setup storage vectors noting that the dilatation is unity
-    // in the reference configuration
-    residual.reinit (dofs_per_block);
-    residual.collect_sizes ();
-
-    solution_n.reinit (dofs_per_block);
-    solution_n.collect_sizes ();
-    solution_n.block(t_dof) = 1.0;
-
-    // and finally set up the quadrature point history
-    setup_qph ();
-
-    timer.leave_subsection();
+template<int dim>
+void Solid<dim>::system_setup(void) {
+       timer.enter_subsection("Setup system");
+
+       // We first describe 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.
+       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
+
+       // DOF handler is then initialised and we renumber the grid in an
+       // efficient manner. We also record the number of DOF's per block.
+       dof_handler_ref.distribute_dofs(fe);
+       DoFRenumbering::Cuthill_McKee(dof_handler_ref);
+       DoFRenumbering::component_wise(dof_handler_ref, block_component);
+       DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block,
+                       block_component);
+
+       std::cout << "Triangulation:" << "\n\t Number of active cells: "
+                       << triangulation.n_active_cells()
+                       << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs()
+                       << std::endl;
+
+       // Setup the sparsity pattern and tangent matrix
+       tangent_matrix.clear();
+       {
+               const unsigned int n_dofs_u = dofs_per_block[u_dof];
+               const unsigned int n_dofs_p = dofs_per_block[p_dof];
+               const unsigned int n_dofs_J = dofs_per_block[J_dof];
+
+               BlockCompressedSimpleSparsityPattern csp(n_blocks, n_blocks);
+
+               csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
+               csp.block(u_dof, p_dof).reinit(n_dofs_u, n_dofs_p);
+               csp.block(u_dof, J_dof).reinit(n_dofs_u, n_dofs_J);
+
+               csp.block(p_dof, u_dof).reinit(n_dofs_p, n_dofs_u);
+               csp.block(p_dof, p_dof).reinit(n_dofs_p, n_dofs_p);
+               csp.block(p_dof, J_dof).reinit(n_dofs_p, n_dofs_J);
+
+               csp.block(J_dof, u_dof).reinit(n_dofs_J, n_dofs_u);
+               csp.block(J_dof, p_dof).reinit(n_dofs_J, n_dofs_p);
+               csp.block(J_dof, J_dof).reinit(n_dofs_J, n_dofs_J);
+               csp.collect_sizes();
+
+               // In order to perform the static condensation efficiently,
+               // we choose to exploit the symmetry of the the system matrix.
+               // The global system matrix has the following structure
+               //      | K_con |  K_up  |     0     |         | dU_u |         | R_u |
+               // K =  | K_pu  |    0   |   K_pJ^-1 | , dU =  | dU_p | , R =   | R_p |
+               //      |   0   |  K_Jp  |   K_JJ    |         | dU_J |         | R_J |
+               // We optimise the sparsity pattern to reflect this structure
+               // and prevent unnecessary data creation for the right-diagonal
+               // block components.
+               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_ref, coupling, csp,
+                               constraints, false);
+               sparsity_pattern.copy_from(csp);
+       }
+
+       tangent_matrix.reinit(sparsity_pattern);
+
+       // Setup storage vectors noting that the dilatation is unity
+       // (i.e. $\widetilde{J} = 1$)
+       // in the undeformed configuration
+       system_rhs.reinit(dofs_per_block);
+       system_rhs.collect_sizes();
+
+       solution_n.reinit(dofs_per_block);
+       solution_n.collect_sizes();
+       solution_n.block(J_dof) = 1.0;
+
+       // and finally set up the quadrature point history
+       setup_qph();
+
+       timer.leave_subsection();
 }
 
 // We next get information from the FE system
 // that describes which local element DOFs are
 // attached to which block component.
-// This is used later to extract subblocks from the global matrix.
-template <int dim>
-void Solid<dim>::determine_component_extractors(void)
-{
-    element_indices_u.clear();
-    element_indices_p.clear();
-    element_indices_t.clear();
-
-    for (unsigned int k=0; k < fe.dofs_per_cell; ++k) {
-        // The next call has the FE System indicate to which block component
-        // the current DOF is attached to.
-        // Currently, the interpotation fields are setup such that
-        // 0 indicates a displacement DOF, 1 a pressure DOF and 2 a dilatation DOF.
-       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 == t_dof) {
-           element_indices_t.push_back(k);
-       }
-       else {
-           Assert (k_group <= t_dof, ExcInternalError());
+// This is used later to extract sub-blocks from the global matrix.
+template<int dim>
+void Solid<dim>::determine_component_extractors(void) {
+       element_indices_u.clear();
+       element_indices_p.clear();
+       element_indices_J.clear();
+
+       for (unsigned int k = 0; k < fe.dofs_per_cell; ++k) {
+               // The next call has the FE System indicate to which block component
+               // the current DOF 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.
+               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());
+               }
        }
-    }
 }
 
 // @sect4{Solid::setup_qph}
 // The method used to store quadrature information is already described in
-// tutorial 18. Here we implement a similar setup for a SMP machine.
-template <int dim>
-void Solid<dim>::setup_qph (void)
-{
-    std::cout << "    Setting up quadrature point data..." << std::endl;
-
-    // Firstly the actual QPH data objects are created. This must be done
-    // only once the grid is refined to its finest level.
-    {
-        quadrature_point_history = std::vector< PointHistory <dim> > (triangulation.n_active_cells() * n_q_points);
-
-       unsigned int history_index = 0;
-        typename Triangulation<dim>::active_cell_iterator
-                        cell = triangulation.begin_active(),
-                        endc = triangulation.end();
-       for (cell = triangulation.begin_active(); cell != endc; ++cell) {
-           cell->set_user_pointer(&quadrature_point_history[history_index]);
-           history_index += n_q_points;
-       }
-
-       Assert(history_index == quadrature_point_history.size(), ExcInternalError());
-    }
-
-    // Next we setup the initial QP data
-    typename DoFHandler<dim>::active_cell_iterator
-           cell = dof_handler_ref.begin_active(),
-           endc = dof_handler_ref.end();
-    for (; cell != endc; ++cell) {
-       PointHistory<dim>* lqph = reinterpret_cast<PointHistory<dim>*> (cell->user_pointer());
-       Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
-       Assert(lqph < &quadrature_point_history.back(), ExcInternalError());
-
-       // Setup any initial information at displacement gauss points
-       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
-           lqph[q_point].setup_lqp( parameters );
-       }
-    }
+// step-18. Here we implement a similar setup for a SMP machine.
+template<int dim>
+void Solid<dim>::setup_qph(void) {
+       std::cout << "    Setting up quadrature point data..." << std::endl;
+
+       // Firstly the actual QPH data objects are created. This must be done
+       // only once the grid is refined to its finest level.
+       {
+               quadrature_point_history = std::vector<PointHistory<dim> >(
+                               triangulation.n_active_cells() * n_q_points);
+
+               unsigned int history_index = 0;
+               typename Triangulation<dim>::active_cell_iterator cell =
+                               triangulation.begin_active(), endc = triangulation.end();
+               for (cell = triangulation.begin_active(); cell != endc; ++cell) {
+                       cell->set_user_pointer(&quadrature_point_history[history_index]);
+                       history_index += n_q_points;
+               }
+
+               Assert(history_index == quadrature_point_history.size(),
+                               ExcInternalError());
+       }
+
+       // Next we setup the initial QP data
+       typename DoFHandler<dim>::active_cell_iterator cell =
+                       dof_handler_ref.begin_active(), endc = dof_handler_ref.end();
+       for (; cell != endc; ++cell) {
+               PointHistory<dim>* lqph =
+                               reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+               Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
+               Assert(lqph < &quadrature_point_history.back(), ExcInternalError());
+
+               // Setup any initial information at displacement Gauss points
+               for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
+                       lqph[q_point].setup_lqp(parameters);
+               }
+       }
 }
 
 // @sect4{Solid::update_qph_incremental}
 // As the update of QP information occurs frequently and involves a number of
 // expensive operations, we define a multi-threaded approach to distributing
 // the task across a number of CPU cores.
-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;
-
-    // Firstly we need to attain the total solution as it stands
-    // at this Newton increment
-    const BlockVector <double> solution_total = get_solution_total(solution_delta);
-
-    // Next we create the initial copy of TBB objects
-    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);
-
-    // and pass them and the one-cell update function to the workstream to be processed
-    WorkStream::run (  dof_handler_ref.begin_active(),
-                      dof_handler_ref.end(),
-                      *this,
-                      &Solid::update_qph_incremental_one_cell,
-                      &Solid::copy_local_to_global_UQPH,
-                      scratch_data_UQPH,
-                      per_task_data_UQPH);
-
-    timer.leave_subsection();
+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;
+
+       // Firstly we need to obtain the total solution as it stands
+       // at this Newton increment
+       const BlockVector<double> solution_total = get_solution_total(
+                       solution_delta);
+
+       // Next we create the initial copy of TBB objects
+       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);
+
+       // and pass them and the one-cell update function to the workstream to be processed
+       WorkStream::run(dof_handler_ref.begin_active(), dof_handler_ref.end(),
+                       *this, &Solid::update_qph_incremental_one_cell,
+                       &Solid::copy_local_to_global_UQPH, scratch_data_UQPH,
+                       per_task_data_UQPH);
+
+       timer.leave_subsection();
 }
 
 // Now we describe how we extract data from the solution vector and pass it
 // along to each QP storage object for processing.
-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)
-{
-    PointHistory<dim>* lqph = reinterpret_cast<PointHistory<dim>*> (cell->user_pointer());
-    Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
-    Assert(lqph < &quadrature_point_history.back(), 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_t_total.size() == n_q_points, ExcInternalError());
-
-    // Firstly we need to find the values and gradients at quadrature points
-    // inside the current cell
-    scratch.fe_values_ref.reinit(cell);
-    scratch.fe_values_ref[u_fe].get_function_gradients (scratch.solution_total, scratch.solution_grads_u_total);
-    scratch.fe_values_ref[p_fe].get_function_values (scratch.solution_total, scratch.solution_values_p_total);
-    scratch.fe_values_ref[t_fe].get_function_values (scratch.solution_total,scratch. solution_values_t_total);
-
-    // and then we update the each local QP using the displacment deformation gradient
-    // and total pressure and dilatation solution values.
-    for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
-       lqph[q_point].update_values (scratch.solution_grads_u_total [q_point],
-                                    scratch.solution_values_p_total[q_point],
-                                    scratch.solution_values_t_total[q_point]);
-    }
+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) {
+       PointHistory<dim>* lqph =
+                       reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+       Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
+       Assert(lqph < &quadrature_point_history.back(), 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());
+
+       // Firstly we need to find the values and gradients at quadrature points
+       // inside the current cell
+       scratch.fe_values_ref.reinit(cell);
+       scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
+                       scratch.solution_grads_u_total);
+       scratch.fe_values_ref[p_fe].get_function_values(scratch.solution_total,
+                       scratch.solution_values_p_total);
+       scratch.fe_values_ref[J_fe].get_function_values(scratch.solution_total,
+                       scratch.solution_values_J_total);
+
+       // and then we update each local QP
+       // using the displacement gradient
+       // and total pressure and dilatation solution values.
+       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
+               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]);
+       }
 }
 
 // @sect4{Solid::solve_nonlinear_timestep}
-template <int dim>
-void Solid<dim>::solve_nonlinear_timestep (BlockVector <double> & solution_delta)
-{
-    //    timer.enter_subsection("Nonlinear solver");
-    std::cout
-            << std::endl
-           << "Timestep " << time.get_timestep()
-           << " @ " << time.current() << "s"
-           << std::endl;
-
-    // We create a new vector to store the current Newton update step
-    BlockVector <double> newton_update (dofs_per_block);
-    newton_update.collect_sizes ();
-
-    // Reset the error storage objects
-    error_residual.reset();
-    error_residual_0.reset();
-    error_residual_norm.reset();
-    error_update.reset();
-    error_update_0.reset();
-    error_update_norm.reset();
-
-    // Print solver header
-    print_conv_header();
-
-    // We now perform a number of Newton iterations to iteratively solve
-    // the nonlinear problem.
-    for (unsigned int it_nr=0; it_nr < parameters.max_iterations_NR; ++ it_nr)
-    {
-        // Print Newton iteration
-       std::cout
-               << " "
-               << std::setw(2)
-               << it_nr
-               << " "
-               << std::flush;
-
-       // Since the problem is fully nonlinear and we are using a
-       // full Newton method, the data stored in the tangent matrix
-       // and residual vector is not reusable and must be cleared
-       // at each Newton step.
-       tangent_matrix = 0.0;
-       residual = 0.0;
-
-       // We initially build the residual vector to check for convergence.
-       // The unconstrained DOF's of the residual vector hold the out-of-balance
-       // forces. This 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.
-       assemble_system_F (); // Assemble RHS
-       get_error_residual(error_residual);
-
-       // We store the residual errors after the first iteration
-       // in order to normalise by their value
-       if (it_nr == 0) error_residual_0 = error_residual;
-
-       // We can now determine the normalised residual error
-       error_residual_norm = error_residual;
-       error_residual_norm.normalise(error_residual_0);
-
-       // Check for solution convergence
-       if (   it_nr > 0
-              && error_update_norm.u <= parameters.tol_u
-              && error_residual_norm.u <= parameters.tol_f)
-       {
-           std::cout
-                   << " CONVERGED! "
-                   << std::endl;
-
-           print_conv_footer();
-
-           //      timer.leave_subsection();
-           return;
-       }
-
-
-       assemble_system_K (); // Assemble stiffness matrix
-       make_constraints (it_nr, constraints); // Make boundary conditions
-       constraints.condense (tangent_matrix,
-                             residual); // Apply BC's
-
-       const std::pair <unsigned int, double> lin_solver_output = solve_linear_system (newton_update);
-       constraints.distribute(newton_update); // Populate the constrained DOF's with their values
-
-       get_error_update(newton_update,
-                        error_update);
-       if (it_nr == 0) error_update_0 = error_update;
-       // We can now determine the normalised newton update error
-       error_update_norm = error_update;
-       error_update_norm.normalise(error_update_0);
-
-       // The current solution state unacceptable, so we need to update
-       // the solution increment for this timestep, update all quadrature
-       // point inforation pertaining to this new displacment and stress state
-       // and continue iterating.
-       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.t << "  "
-               << error_update_norm.norm << "  "
-               << error_update_norm.u << "  "
-               << error_update_norm.p << "  "
-               << error_update_norm.t << "  "
-               << std::endl;
-    }
-
-    throw(ExcMessage("No convergence in nonlinear solver!"));
+template<int dim>
+void Solid<dim>::solve_nonlinear_timestep(
+               BlockVector<double> & solution_delta) {
+       //    timer.enter_subsection("Nonlinear solver");
+       std::cout << std::endl << "Timestep " << time.get_timestep() << " @ "
+                       << time.current() << "s" << std::endl;
+
+       // We create a new vector to store the current Newton update step
+       BlockVector<double> newton_update(dofs_per_block);
+       newton_update.collect_sizes();
+
+       // Reset the error storage objects
+       error_residual.reset();
+       error_residual_0.reset();
+       error_residual_norm.reset();
+       error_update.reset();
+       error_update_0.reset();
+       error_update_norm.reset();
+
+       // Print solver header
+       print_conv_header();
+
+       // We now perform a number of Newton iterations to iteratively solve
+       // the nonlinear problem.
+       for (unsigned int it_nr = 0; it_nr < parameters.max_iterations_NR;
+                       ++it_nr) {
+               // Print Newton iteration
+               std::cout << " " << std::setw(2) << it_nr << " " << std::flush;
+
+               // 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.
+               tangent_matrix = 0.0;
+               system_rhs = 0.0;
+
+               // We initially build the right-hand side vector to check for convergence.
+               // The unconstrained DOF's 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.
+               assemble_system_rhs(); // Assemble RHS
+               get_error_residual(error_residual);
+
+               // We store the residual errors after the first iteration
+               // in order to normalise by their value
+               if (it_nr == 0)
+                       error_residual_0 = error_residual;
+
+               // We can now determine the normalised residual error
+               error_residual_norm = error_residual;
+               error_residual_norm.normalise(error_residual_0);
+
+               // Check for solution convergence
+               if (it_nr > 0 && error_update_norm.u <= parameters.tol_u
+                               && error_residual_norm.u <= parameters.tol_f) {
+                       std::cout << " CONVERGED! " << std::endl;
+                       print_conv_footer();
+                       return;
+               }
+
+               assemble_system_K(); // Assemble stiffness matrix
+               make_constraints(it_nr, constraints); // Make boundary conditions
+               constraints.condense(tangent_matrix, system_rhs); // Apply BC's
+
+               const std::pair<unsigned int, double> lin_solver_output =
+                               solve_linear_system(newton_update);
+
+               get_error_update(newton_update, error_update);
+               if (it_nr == 0)
+                       error_update_0 = error_update;
+
+               // We can now determine the normalised Newton update error
+               error_update_norm = error_update;
+               error_update_norm.normalise(error_update_0);
+
+               // The current solution state is unacceptable, so we need to update
+               // the solution increment for this time step, update all quadrature
+               // point information pertaining to this new displacement and stress state
+               // and continue iterating.
+               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;
+       }
+
+       throw(ExcMessage("No convergence in nonlinear solver!"));
 }
 
 // We print out data in a nice table that is updated
 // on a per-iteration basis. Here we set up the table
 // header
-template <int dim>
-void Solid<dim>::print_conv_header (void)
-{
-    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_T    "
-            << " NU_NORM     "
-            << " NU_U      "
-            << " NU_P      "
-            << " NU_T "
-            << 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_header(void) {
+       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_T    " << " NU_NORM     "
+                       << " NU_U      " << " NU_P      " << " NU_T " << std::endl;
+
+       for (unsigned int i = 0; i < l_width; ++i)
+               std::cout << "_";
+       std::cout << std::endl;
 }
 // and here the footer
-template <int dim>
-void Solid<dim>::print_conv_footer (void)
-{
-    static const unsigned int l_width = 155;
-
-    for (unsigned int i=0; i < l_width; ++i)
-        std::cout << "_";
-    std::cout << std::endl;
-
-
-    std::cout
-            << "Relative errors:" << std::endl
-            << "Displacement:\t" << error_update.u/error_update_0.u << std::endl
-            << "Force: \t\t" << error_residual.u/error_residual_0.u << std::endl
-            << "Dilatation:\t" << get_error_dil()
-            << std::endl;
+template<int dim>
+void Solid<dim>::print_conv_footer(void) {
+       static const unsigned int l_width = 155;
+
+       for (unsigned int i = 0; i < l_width; ++i)
+               std::cout << "_";
+       std::cout << std::endl;
+
+       std::cout << "Relative errors:" << std::endl << "Displacement:\t"
+                       << error_update.u / error_update_0.u << std::endl << "Force: \t\t"
+                       << error_residual.u / error_residual_0.u << std::endl
+                       << "Dilatation:\t" << get_error_dil() << std::endl;
 }
 
-// Calculate the ratio of the volume of the domain in the
-// current configuration and the reference configuration
-template <int dim>
-double Solid<dim>::get_error_dil (void)
-{
-    double v_e = 0.0; // Volume in current configuration
-    double V_e = 0.0; // Volume in reference configuration
-
-    FEValues<dim> fe_values_ref (fe, qf_cell, update_JxW_values);
-
-    typename DoFHandler<dim>::active_cell_iterator
-            cell = dof_handler_ref.begin_active(),
-            endc = dof_handler_ref.end();
-    for (; cell != endc; ++cell) {
-        fe_values_ref.reinit (cell);
-        PointHistory<dim>* lqph = reinterpret_cast<PointHistory<dim>*> (cell->user_pointer());
-       Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
-       Assert(lqph < &quadrature_point_history.back(), ExcInternalError());
-
-        for (unsigned int q_point=0; q_point < n_q_points; ++q_point) {
-            v_e += lqph[q_point].get_dilatation() * fe_values_ref.JxW(q_point);
-            V_e += fe_values_ref.JxW(q_point);
-        }
-    }
-
-    return std::abs((v_e - V_e)/V_e); // Difference between initial and current volume
+// Calculate how well the dilatation $\widetilde{J}$ 
+// agrees with $J := \textrm{det}\mathbf{F}$
+// from the $L^2$ error
+// $ \bigl[ \int_{\Omega_0} {[ J - \widetilde{J}]}^{2}\textrm{d}V \bigr]^{1/2}$ 
+// which is then normalised by the current volume
+// $\int_{\Omega_0}  J ~\textrm{d}V = \int_\Omega  ~\textrm{d}v$.
+template<int dim>
+double Solid<dim>::get_error_dil(void) {
+       double vol = 0.0; // Volume of current configuration
+       double dil_L2_error = 0.0;
+
+       FEValues<dim> fe_values_ref(fe, qf_cell, update_JxW_values);
+
+       typename DoFHandler<dim>::active_cell_iterator cell =
+                       dof_handler_ref.begin_active(), endc = dof_handler_ref.end();
+       for (; cell != endc; ++cell) {
+               fe_values_ref.reinit(cell);
+               PointHistory<dim>* lqph =
+                               reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+               Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
+               Assert(lqph < &quadrature_point_history.back(), ExcInternalError());
+
+               for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
+
+                       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_ref.JxW(q_point);
+
+                       dil_L2_error += the_error_qp_squared * JxW;
+                       vol += det_F_qp * JxW;
+               }
+       }
+       Assert(vol >= 0, ExcInternalError());
+       return std::sqrt(dil_L2_error) / vol;
 }
 
-// Determine the true residual error for the problem
-template <int dim>
-void Solid<dim>::get_error_residual (Errors & error_residual)
-{
-    BlockVector <double> error_res (dofs_per_block);
-    error_res.collect_sizes ();
-
-    // Need to ignore constrained DOFs
-    for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i)
-        if (!constraints.is_constrained(i))
-            error_res(i) = residual(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.t = error_res.block(t_dof).l2_norm();
+// Determine the true residual error for the problem. 
+// That is, determine the error in the residual for
+// unconstrained dof.
+template<int dim>
+void Solid<dim>::get_error_residual(Errors & error_residual) {
+       BlockVector<double> error_res(dofs_per_block);
+       error_res.collect_sizes();
+
+       // Need to ignore constrained DOFs
+       for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
+               if (!constraints.is_constrained(i))
+                       error_res(i) = system_rhs(i);
+
+       error_residual.norm = error_res.l2_norm();
+       error_residual.u = error_res.block(u_dof).l2_norm();
+       error_residual.p = error_res.block(p_dof).l2_norm();
+       error_residual.J = error_res.block(J_dof).l2_norm();
 }
 
 // Determine the true Newton update error for the problem
-template <int dim>
-void Solid<dim>::get_error_update (const BlockVector <double> & newton_update,
-                                   Errors & error_update)
-{
-    BlockVector <double> error_ud (dofs_per_block);
-    error_ud.collect_sizes ();
-
-    // Need to ignore constrained DOFs as they have a prescribed
-    // value
-    for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i)
-        if (!constraints.is_constrained(i))
-            error_ud(i) = newton_update(i);
-
-    error_update.norm = error_ud.l2_norm();
-    error_update.u = error_ud.block(u_dof).l2_norm();
-    error_update.p = error_ud.block(p_dof).l2_norm();
-    error_update.t = error_ud.block(t_dof).l2_norm();
+template<int dim>
+void Solid<dim>::get_error_update(const BlockVector<double> & newton_update,
+               Errors & error_update) {
+       BlockVector<double> error_ud(dofs_per_block);
+       error_ud.collect_sizes();
+
+       // Need to ignore constrained DOFs as they have a prescribed
+       // value
+       for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
+               if (!constraints.is_constrained(i))
+                       error_ud(i) = newton_update(i);
+
+       error_update.norm = error_ud.l2_norm();
+       error_update.u = error_ud.block(u_dof).l2_norm();
+       error_update.p = error_ud.block(p_dof).l2_norm();
+       error_update.J = error_ud.block(J_dof).l2_norm();
 }
 
 // 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.
-template <int dim>
-BlockVector <double> Solid<dim>::get_solution_total (const BlockVector <double> & solution_delta)
-{
-    BlockVector <double> solution_total (solution_n);
-    solution_total += solution_delta;
+template<int dim>
+BlockVector<double> Solid<dim>::get_solution_total(
+               const BlockVector<double> & solution_delta) {
+       BlockVector<double> solution_total(solution_n);
+       solution_total += solution_delta;
+       return solution_total;
 
-    return solution_total;
 }
 
 // @sect4{Solid::assemble_system_K}
@@ -1798,298 +1796,318 @@ BlockVector <double> Solid<dim>::get_solution_total (const BlockVector <double>
 // 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.
-template <int dim>
-void Solid<dim>::assemble_system_K (void)
-{
-    timer.enter_subsection("Assemble tangent matrix");
-    std::cout << " ASM_K " << std::flush;
+template<int dim>
+void Solid<dim>::assemble_system_K(void) {
+       timer.enter_subsection("Assemble tangent matrix");
+       std::cout << " ASM_K " << std::flush;
 
-    tangent_matrix = 0.0;
+       tangent_matrix = 0.0;
 
-    const UpdateFlags uf_cell (update_values | update_gradients | update_JxW_values);
+       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);
+       PerTaskData_K per_task_data(dofs_per_cell);
+       ScratchData_K scratch_data(fe, qf_cell, uf_cell);
 
-    WorkStream::run (  dof_handler_ref.begin_active(),
-                       dof_handler_ref.end(),
-                       *this,
-                       &Solid::assemble_system_K_one_cell,
-                       &Solid::copy_local_to_global_K,
-                       scratch_data,
-                       per_task_data);
+       WorkStream::run(dof_handler_ref.begin_active(), dof_handler_ref.end(),
+                       *this, &Solid::assemble_system_K_one_cell,
+                       &Solid::copy_local_to_global_K, scratch_data, per_task_data);
 
-    timer.leave_subsection();
+       timer.leave_subsection();
 }
 
 // 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 seperate functions.
-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));
+// been split up into two separate functions.
+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));
 }
 
 // Here we define how we assemble the tangent matrix contribution for a
 // single cell.
-template <int dim>
-void Solid<dim>::assemble_system_K_one_cell (const typename DoFHandler<dim>::active_cell_iterator & cell,
-                                             ScratchData_K & scratch,
-                                             PerTaskData_K & data)
-{
-    // We first need to reset and initialise some of the data structures and retrieve some
-    // basic information regarding the DOF numbering on this cell
-    data.reset();
-    scratch.reset();
-    scratch.fe_values_ref.reinit (cell);
-    cell->get_dof_indices (data.local_dof_indices);
-    PointHistory<dim> *lqph = reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
-
-    // We can precalculate the cell shape function values and gradients. Note that the
-    // shape function gradients are defined in the current configuration for this problem.
-    static const SymmetricTensor<2, dim> I = unit_symmetric_tensor <dim> ();
-    for (unsigned int q_point=0; q_point < n_q_points; ++q_point) {
-        const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
-
-       for (unsigned int k=0; k< dofs_per_cell; ++k) {
-           const unsigned int k_group = fe.system_to_base_index(k).first.first;
-
-           if (k_group == u_dof) {
-               scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
-               scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
-           }
-           else if (k_group == p_dof) {
-               scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k, q_point);
-           }
-           else if (k_group == t_dof) {
-               scratch.Nx[q_point][k] = scratch.fe_values_ref[t_fe].value(k, q_point);
-           }
-           else {
-               Assert (k_group <= t_dof, ExcInternalError());
-           }
-       }
-    }
-
-    // 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 those values to the upper half.
-    // So we only assemble half of the K_uu, K_pp (= 0), K_tt blocks, while the whole
-    // K_pt, K_ut, K_up blocks are built.
-    for (unsigned int q_point=0; q_point < n_q_points; ++q_point) {
-        // We first extract some configuration dependent variables from our
-        // QPH history objects that remain constant at each QP.
-        const Tensor <2,dim>          T   = static_cast < Tensor<2, dim> > (lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol());
-        const SymmetricTensor <4,dim> C   = lqph[q_point].get_C_iso() + lqph[q_point].get_C_vol();
-        const double                  C_v = lqph[q_point].get_d2U_dtheta2();
-        const double                  J   = lqph[q_point].get_J();
-
-       // Next we define some aliases to make the assembly process easier to follow
-       const std::vector<double> & N = scratch.Nx[q_point];
-       const std::vector< SymmetricTensor <2,dim> > & symm_B = scratch.symm_grad_Nx[q_point];
-       const std::vector< Tensor <2,dim> > & B = scratch.grad_Nx[q_point];
-       const double & JxW = scratch.fe_values_ref.JxW(q_point);
-
-       for (unsigned int i=0; i < dofs_per_cell; ++i) {
-           const unsigned int component_i = fe.system_to_component_index(i).first;
-           // Determine the dimensional component that matches the dof component (i.e. i % dim)
-           const unsigned int i_group = fe.system_to_base_index(i).first.first;
-
-           for (unsigned int j=0; j <= i; ++j) {
-               const unsigned int component_j = fe.system_to_component_index(j).first;
-               const unsigned int j_group = fe.system_to_base_index(j).first.first;
-
-               // This is the K_{uu} contribution. It comprises of a material stiffness
-               // contribution and a geometric stiffness contribution which is only
-               // added along the local matrix diagonals
-               if (   (i_group == j_group) && (i_group == u_dof ) ) {
-                   data.cell_matrix(i,j) += symm_B[i] * C * symm_B[j] * JxW;
-                   if (component_i == component_j)
-                       data.cell_matrix(i,j) += B[i][component_i] * T * B[j][component_j] * JxW;
-               }
-               // Next is the K_{pu} contibution
-               else if ( (i_group == p_dof) && (j_group == u_dof) ) {
-                   data.cell_matrix(i,j) -= N[i]*J*(symm_B[j]*I)*JxW;
+template<int dim>
+void Solid<dim>::assemble_system_K_one_cell(
+               const typename DoFHandler<dim>::active_cell_iterator & cell,
+               ScratchData_K & scratch, PerTaskData_K & data) {
+       // We first need to reset and initialise some
+       // of the data structures and retrieve some
+       // basic information regarding the DOF numbering on this cell
+       data.reset();
+       scratch.reset();
+       scratch.fe_values_ref.reinit(cell);
+       cell->get_dof_indices(data.local_dof_indices);
+       PointHistory<dim> *lqph =
+                       reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+
+       // We can precalculate the cell shape function values and gradients. Note that the
+       // shape function gradients are defined wrt the current configuration.
+       // That is
+       // $\textrm{grad}\boldsymbol{\varphi} = \textrm{Grad}\boldsymbol{\varphi} \mathbf{F}^{-1}$
+       static const SymmetricTensor<2, dim> I = AdditionalTools::StandardTensors<
+                       dim>::I;
+       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
+               const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
+               for (unsigned int k = 0; k < dofs_per_cell; ++k) {
+                       const unsigned int k_group = fe.system_to_base_index(k).first.first;
+
+                       if (k_group == u_dof) {
+                               scratch.grad_Nx[q_point][k] =
+                                               scratch.fe_values_ref[u_fe].gradient(k, q_point)
+                                                               * F_inv;
+                               scratch.symm_grad_Nx[q_point][k] = symmetrize(
+                                               scratch.grad_Nx[q_point][k]);
+                       } else if (k_group == p_dof) {
+                               scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
+                                               q_point);
+                       } else if (k_group == J_dof) {
+                               scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
+                                               q_point);
+                       } else {
+                               Assert(k_group <= J_dof, ExcInternalError());
+                       }
                }
-               // and the K_{tp} contibution
-               else if ( (i_group == t_dof) && (j_group == p_dof) ) {
-                   data.cell_matrix(i,j) += N[i]*N[j]*JxW;
+       }
+
+       // 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 K_uu, K_pp (= 0), K_JJ blocks, while the whole
+       // K_pJ, K_uJ (=0), K_up blocks are built.
+       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
+               // We first extract some configuration dependent variables from our
+               // QPH history objects that for the current q_point.
+               // Get the current stress state $\boldsymbol{\tau}$
+               const Tensor<2, dim> tau =
+                               static_cast<Tensor<2, dim> >(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();
+
+               // Next we define some aliases to make the assembly process easier to follow
+               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_ref.JxW(q_point);
+
+               for (unsigned int i = 0; i < dofs_per_cell; ++i) {
+                       const unsigned int component_i =
+                                       fe.system_to_component_index(i).first;
+                       // Determine the dimensional component that matches the dof component (i.e. i % dim)
+                       const unsigned int i_group = fe.system_to_base_index(i).first.first;
+
+                       for (unsigned int j = 0; j <= i; ++j) {
+                               const unsigned int component_j =
+                                               fe.system_to_component_index(j).first;
+                               const unsigned int j_group =
+                                               fe.system_to_base_index(j).first.first;
+
+                               // This is the K_{uu} contribution. It comprises of a material
+                               // contribution and a geometrical stress contribution which is only
+                               // added along the local matrix diagonals
+                               if ((i_group == j_group) && (i_group == u_dof)) {
+                                       // The material contribution:
+                                       data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc
+                                                       * symm_grad_Nx[j] * JxW;
+                                       if (component_i == component_j) // geometrical stress contribution
+                                               data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau
+                                                               * grad_Nx[j][component_j] * JxW;
+                               }
+                               // Next is the K_{pu} contribution
+                               else if ((i_group == p_dof) && (j_group == u_dof)) {
+                                       data.cell_matrix(i, j) += N[i] * det_F
+                                                       * (symm_grad_Nx[j]
+                                                                       * AdditionalTools::StandardTensors<dim>::I)
+                                                       * JxW;
+                               }
+                               // and the K_{Jp} contribution
+                               else if ((i_group == J_dof) && (j_group == p_dof)) {
+                                       data.cell_matrix(i, j) -= N[i] * N[j] * JxW;
+                               }
+                               // and lastly the K_{JJ} contribution
+                               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());
+                       }
                }
-               // and lastly the K_{tt} contibution
-               else if ( (i_group == j_group) && (i_group == t_dof)  ) {
-                   data.cell_matrix(i,j) -=  N[i]*C_v*N[j]*JxW;
+       }
+
+       // Here we copy the lower half of the local matrix in the upper
+       // half of the local matrix
+       for (unsigned int i = 0; i < dofs_per_cell; ++i) {
+               for (unsigned int j = i + 1; j < dofs_per_cell; ++j) {
+                       data.cell_matrix(i, j) = data.cell_matrix(j, i);
                }
-               else Assert ((i_group <= t_dof) && (j_group <= t_dof), ExcInternalError());
-           }
-       }
-    }
-
-    // Here we copy the lower half of the local matrix in the upper
-    // half of the local matrix
-    for (unsigned int i=0; i<dofs_per_cell; ++i) {
-        for (unsigned int j=i+1; j<dofs_per_cell; ++j) {
-            data.cell_matrix(i,j) = data.cell_matrix(j,i);
-        }
-    }
+       }
 }
 
-// @sect4{Solid::assemble_system_F}
-// The setup of the residual assembly process is similar to the
+// @sect4{Solid::assemble_system_rhs}
+// 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.
-template <int dim>
-void Solid<dim>::assemble_system_F (void)
-{
-    timer.enter_subsection("Assemble residual");
-    std::cout << " ASM_R "<< std::flush;
-
-    residual  = 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_F per_task_data (dofs_per_cell);
-    ScratchData_F scratch_data (fe,
-                                qf_cell,
-                                uf_cell,
-                                qf_face,
-                                uf_face);
-
-    WorkStream::run ( dof_handler_ref.begin_active(),
-                      dof_handler_ref.end(),
-                      *this,
-                      &Solid::assemble_system_F_one_cell,
-                      &Solid::copy_local_to_global_F,
-                      scratch_data,
-                      per_task_data );
-
-    timer.leave_subsection();
+template<int dim>
+void Solid<dim>::assemble_system_rhs(void) {
+       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_ref.begin_active(), dof_handler_ref.end(),
+                       *this, &Solid::assemble_system_rhs_one_cell,
+                       &Solid::copy_local_to_global_rhs, scratch_data, per_task_data);
+
+       timer.leave_subsection();
 }
 
-template <int dim>
-void Solid<dim>::copy_local_to_global_F (const PerTaskData_F & data)
-{
-    for (unsigned int i=0; i<dofs_per_cell; ++i) {
-        residual(data.local_dof_indices[i]) += data.cell_rhs(i);
-    }
+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_F_one_cell (const typename DoFHandler<dim>::active_cell_iterator & cell,
-                                             ScratchData_F & scratch,
-                                             PerTaskData_F & data)
-{
-    // Again we reset the data structures
-    data.reset();
-    scratch.reset();
-    scratch.fe_values_ref.reinit (cell);
-    cell->get_dof_indices (data.local_dof_indices);
-    PointHistory<dim> *lqph = reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
-
-    // and then precompute some shape function data
-    for (unsigned int q_point=0; q_point < n_q_points; ++q_point) {
-        const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
-
-       for (unsigned int k=0; k<dofs_per_cell; ++k) {
-           const unsigned int k_group = fe.system_to_base_index(k).first.first;
-
-           if (k_group == u_dof) {
-               scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv);
-           }
-           else if (k_group == p_dof) {
-               scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k, q_point);
-           }
-           else if (k_group == t_dof) {
-               scratch.Nx[q_point][k] = scratch.fe_values_ref[t_fe].value(k, q_point);
-           }
-           else Assert (k_group <= t_dof, ExcInternalError());
-       }
-    }
-
-    // and can now assemble the residual contribution
-    for (unsigned int q_point=0; q_point < n_q_points; ++q_point) {
-        // We fist retrieve data that remains constant a QP
-        const SymmetricTensor <2,dim>  T = lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol();
-        const double  J = lqph[q_point].get_J();
-        const double  D = lqph[q_point].get_dilatation();
-        const double  p = lqph[q_point].get_pressure();
-        const double  p_star = lqph[q_point].get_dU_dtheta();
-
-       // define some shortcuts
-       const std::vector< double > & N = scratch.Nx[q_point];
-       const std::vector< SymmetricTensor <2,dim> > & symm_B = scratch.symm_grad_Nx[q_point];
-       const double  JxW = scratch.fe_values_ref.JxW(q_point);
-
-       for (unsigned int i=0; i<dofs_per_cell; ++i) {
-           const unsigned int i_group = fe.system_to_base_index(i).first.first;
-           // Add the contribution to the R_{u} block
-           if (i_group == u_dof) {
-               data.cell_rhs(i) -= ( symm_B[i]*T )*JxW;
-           }
-           // the R_{p} block
-           else if (i_group == p_dof ) {
-               data.cell_rhs(i) += N[i]*(J - D)*JxW;
-           }
-           // and finally the R_{t} block
-           else if ( i_group == t_dof) {
-               data.cell_rhs(i) += N[i]*(p_star-p)*JxW;
-           }
-           else Assert (i_group <= t_dof, ExcInternalError());
-       }
-    }
-
-    // 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.
-    if (cell->at_boundary() == true) {
-        for (unsigned int face=0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
-            if (    cell->face(face)->at_boundary() == true
-                    &&  cell->face(face)->boundary_indicator() == 6 ) {
-                scratch.fe_face_values_ref.reinit (cell, face);
-
-               for (unsigned int f_q_point=0; f_q_point < n_q_points_f; ++f_q_point) {
-                   // We retrieve the face normal at this QP
-                   const Tensor <1, dim> & N = scratch.fe_face_values_ref.normal_vector(f_q_point);
-
-                   // and specify the traction in reference configuration. For this problem,
-                   // a defined pressure is applied in the reference configuration. so the
-                   // traction defined using the first Piola-Kirchhoff stress is simply
-                   // t_0 = P*N = (pI)*N = p*N
-                   // We choose to use the time variable to linearly ramp up the pressure
-                   // load.
-                   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 (unsigned int i=0; i < dofs_per_cell; ++i) {
-                       const unsigned int i_group = fe.system_to_base_index(i).first.first;
+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) {
+       // Again we reset the data structures
+       data.reset();
+       scratch.reset();
+       scratch.fe_values_ref.reinit(cell);
+       cell->get_dof_indices(data.local_dof_indices);
+       PointHistory<dim> *lqph =
+                       reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+
+       // and then precompute some shape function data
+       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
+               const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
+
+               for (unsigned int k = 0; k < dofs_per_cell; ++k) {
+                       const unsigned int k_group = fe.system_to_base_index(k).first.first;
+
+                       if (k_group == u_dof) {
+                               scratch.symm_grad_Nx[q_point][k] = symmetrize(
+                                               scratch.fe_values_ref[u_fe].gradient(k, q_point)
+                                                               * F_inv);
+                       } else if (k_group == p_dof) {
+                               scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
+                                               q_point);
+                       } else if (k_group == J_dof) {
+                               scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
+                                               q_point);
+                       } else
+                               Assert(k_group <= J_dof, ExcInternalError());
+               }
+       }
 
+       // and can now assemble the right-hand side contribution
+       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
+               // We fist retrieve data stored at the qp
+               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 = lqph[q_point].get_p();
+               const double dPsi_vol_dJ = lqph[q_point].get_dPsi_vol_dJ();
+
+               // define some shortcuts
+               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_ref.JxW(q_point);
+
+               // 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.
+               for (unsigned int i = 0; i < dofs_per_cell; ++i) {
+                       const unsigned int i_group = fe.system_to_base_index(i).first.first;
+                       // Add the contribution to the F_u block
                        if (i_group == u_dof) {
-                           // More shortcuts being assigned
-                           const unsigned int component_i = fe.system_to_component_index(i).first;
-                           const double & Ni = scratch.fe_face_values_ref.shape_value(i,f_q_point);
-                           const double & JxW = scratch.fe_face_values_ref.JxW(f_q_point);
-
-                           // And finally we can add the traction vector contribution to
-                           // the local RHS vector. Note that this contribution is present
-                           // on displacement DOFs only.
-                           data.cell_rhs(i) += (Ni * traction[component_i]) * JxW;
+                               data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
+                       }
+                       // the F_p block
+                       else if (i_group == p_dof) {
+                               data.cell_rhs(i) -= N[i] * (det_F - J_tilde) * JxW;
+                       }
+                       // and finally the F_J block
+                       else if (i_group == J_dof) {
+                               data.cell_rhs(i) -= N[i] * (dPsi_vol_dJ - p) * JxW;
+                       } else
+                               Assert(i_group <= J_dof, ExcInternalError());
+               }
+       }
+
+       // 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.
+       if (cell->at_boundary() == true) {
+               for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+                               ++face) {
+                       if (cell->face(face)->at_boundary() == true
+                                       && cell->face(face)->boundary_indicator() == 6) {
+                               scratch.fe_face_values_ref.reinit(cell, face);
+
+                               for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
+                                               ++f_q_point) {
+                                       // We retrieve the face normal at this QP
+                                       const Tensor<1, dim> & N =
+                                                       scratch.fe_face_values_ref.normal_vector(f_q_point);
+
+                                       // and 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
+                                       // t_0 = P*N = (pI)*N = p*N
+                                       // We choose to use the time variable to linearly ramp up the pressure
+                                       // load.
+                                       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 (unsigned int i = 0; i < dofs_per_cell; ++i) {
+                                               const unsigned int i_group =
+                                                               fe.system_to_base_index(i).first.first;
+
+                                               if (i_group == u_dof) {
+                                                       // More shortcuts being assigned
+                                                       const unsigned int component_i =
+                                                                       fe.system_to_component_index(i).first;
+                                                       const double Ni =
+                                                                       scratch.fe_face_values_ref.shape_value(i,
+                                                                                       f_q_point);
+                                                       const double JxW = scratch.fe_face_values_ref.JxW(
+                                                                       f_q_point);
+
+                                                       // And finally we can add the traction vector contribution to
+                                                       // the local RHS vector. Note that this contribution is present
+                                                       // on displacement DOFs only.
+                                                       data.cell_rhs(i) += (Ni * traction[component_i])
+                                                                       * JxW;
+                                               }
+                                       }
+                               }
                        }
-                   }
                }
-           }
        }
-    }
 }
 
 // @sect4{Solid::make_constraints}
@@ -2102,465 +2120,462 @@ void Solid<dim>::assemble_system_F_one_cell (const typename DoFHandler<dim>::act
 // completeness although for this problem the constraints are
 // trivial and it would not have made a difference if this had
 // not been accounted for in this problem.
-template <int dim>
-void Solid<dim>::make_constraints (const int & it_nr,
-                                   ConstraintMatrix & constraints)
-{
-    std::cout << " CST "<< std::flush;
-
-    // Since the constraints are different at 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.
-    if (it_nr > 1) return;
-    constraints.clear();
-    const bool apply_dirichlet_bc = (it_nr == 0);
-
-    // The boundary conditions for the indentation problem are as follows:
-    // On the -x, -y and -z faces (ID's 0,2,4) we set up a symmetry condition
-    // to allow only planar movement while the +x and +y faces (ID's 1,3) are
-    // traction free. In this contrived problem, part of the +z face (ID 5) is
-    // set to have no motion in the x- and y-component. Finally, as described
-    // earlier, the other part of the +z face has an the applied pressure but
-    // is also constrained in the x- and y-directions.
-    {
-        const int boundary_id = 0;
-
-       std::vector< bool > components (n_components, false);
-       components[0] = true;
-
-       if (apply_dirichlet_bc == true) {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-       else {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-    }
-    {
-        const int boundary_id = 2;
-
-       std::vector< bool > components (n_components, false);
-       components[1] = true;
-
-       if (apply_dirichlet_bc == true) {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-       else {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-    }
-    {
-        const int boundary_id = 4;
-        std::vector< bool > components (n_components, false);
-        components[2] = true;
-
-       if (apply_dirichlet_bc == true) {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-       else {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-    }
-    {
-        const int boundary_id = 5;
-        std::vector< bool > components (n_components, true);
-        components[2] = false;
-
-       if (apply_dirichlet_bc == true) {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-       else {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-    }
-    {
-        const int boundary_id = 6;
-        std::vector< bool > components (n_components, true);
-        components[2] = false;
-
-       if (apply_dirichlet_bc == true) {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-       else {
-           VectorTools::interpolate_boundary_values ( dof_handler_ref,
-                                                      boundary_id,
-                                                      ZeroFunction<dim>(n_components),
-                                                      constraints,
-                                                      components );
-       }
-    }
-
-    constraints.close();
+template<int dim>
+void Solid<dim>::make_constraints(const int & it_nr,
+               ConstraintMatrix & constraints) {
+       std::cout << " CST " << std::flush;
+
+       // Since the constraints are different at 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.
+       if (it_nr > 1)
+               return;
+       constraints.clear();
+       const bool apply_dirichlet_bc = (it_nr == 0);
+
+       // The boundary conditions for the indentation problem are as follows:
+       // On the -x, -y and -z faces (ID's 0,2,4) we set up a symmetry condition
+       // to allow only planar movement while the +x and +y faces (ID's 1,3) are
+       // traction free. In this contrived problem, part of the +z face (ID 5) is
+       // set to have no motion in the x- and y-component. Finally, as described
+       // earlier, the other part of the +z face has an the applied pressure but
+       // is also constrained in the x- and y-directions.
+       {
+               const int boundary_id = 0;
+
+               std::vector<bool> components(n_components, false);
+               components[0] = true;
+
+               if (apply_dirichlet_bc == true) {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               } else {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               }
+       }
+       {
+               const int boundary_id = 2;
+
+               std::vector<bool> components(n_components, false);
+               components[1] = true;
+
+               if (apply_dirichlet_bc == true) {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               } else {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               }
+       }
+       {
+               const int boundary_id = 4;
+               std::vector<bool> components(n_components, false);
+               components[2] = true;
+
+               if (apply_dirichlet_bc == true) {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               } else {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               }
+       }
+       {
+               const int boundary_id = 5;
+               std::vector<bool> components(n_components, true);
+               components[2] = false;
+
+               if (apply_dirichlet_bc == true) {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               } else {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               }
+       }
+       {
+               const int boundary_id = 6;
+               std::vector<bool> components(n_components, true);
+               components[2] = false;
+
+               if (apply_dirichlet_bc == true) {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               } else {
+                       VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                       boundary_id, ZeroFunction<dim>(n_components), constraints,
+                                       components);
+               }
+       }
+
+       constraints.close();
 }
 
 // @sect4{Solid::solve_linear_system}
 // Solving the entire block system is a bit problematic as there are no
-// contributions to the K_{pp} block, rendering it non-invertable.
+// contributions to the K_{JJ} block, rendering it non-invertible.
 // 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.
-template <int dim>
-std::pair <unsigned int, double> Solid<dim>::solve_linear_system (BlockVector <double> & newton_update)
-{
-    // Need to create two temporary vectors so that the static condensation operation can be performed
-    BlockVector <double> A (dofs_per_block);
-    BlockVector <double> B (dofs_per_block);
-    A.collect_sizes ();
-    B.collect_sizes ();
-
-    // Store the number of linear solver iterations and residual
-    unsigned int lin_it = 0;
-    double lin_res = 0.0;
-
-    //      | K'_uu |   K_up  |     0     |         | dU_u |         | dR_u |
-    // K =  | K_pu  |     0   |   K_pt^-1 | , dU =  | dU_p | , dR =  | dR_p |
-    //      |   0   |   K_tp  |   K_tt    |         | dU_t |         | dR_t |
-
-    // Solve for du
-    {
-        // Do the static condensation to make K'_uu,
-        // and put K_pt^{-1} in the K_pt block
-        assemble_SC();
-
-       // K'uu du = Ru'
-       // with Ru' = Ru âˆ’ Kup Ktp^-1 (Rt âˆ’ Ktt Kpt^{-1} Rp)
-       // Assemble the RHS vector to solve for du
-       tangent_matrix.block(p_dof, t_dof).vmult(A.block(t_dof), residual.block(p_dof));
-       tangent_matrix.block(t_dof, t_dof).vmult (B.block(t_dof), A.block(t_dof));
-       A.block(t_dof).equ(1.0, residual.block(t_dof), -1.0, B.block(t_dof));
-       tangent_matrix.block(p_dof, t_dof).Tvmult(A.block(p_dof), A.block(t_dof));
-       tangent_matrix.block(u_dof, p_dof).vmult(A.block(u_dof), A.block(p_dof));
-       residual.block(u_dof) -= A.block(u_dof);
-
-       timer.enter_subsection("Linear solver");
-       std::cout << " SLV " << std::flush;
-       if (parameters.type_lin == "CG")
+template<int dim>
+std::pair<unsigned int, double> Solid<dim>::solve_linear_system(
+               BlockVector<double> & newton_update) {
+       // Need to create two temporary vectors to help
+       // with the static condensation.
+       BlockVector<double> A(dofs_per_block);
+       BlockVector<double> B(dofs_per_block);
+       A.collect_sizes();
+       B.collect_sizes();
+
+       // Store the number of linear solver iterations and residuals
+       unsigned int lin_it = 0;
+       double lin_res = 0.0;
+
+       //      | K_con |   K_up  |     0     |         | du |         | F_u |
+       // K =  | K_pu  |     0   |   K_pJ^-1 | , dXi = | dp | , R =   | F_p |
+       //      |   0   |   K_Jp  |   K_JJ    |         | dJ |         | F_J |
+
+       // Solve for du
        {
-           const int solver_its = tangent_matrix.block(u_dof, u_dof).m() * parameters.max_iterations_lin;
-           const double tol_sol = parameters.tol_lin * residual.block(u_dof).l2_norm();
-
-           SolverControl solver_control (solver_its , tol_sol);
-
-           GrowingVectorMemory < Vector<double> > GVM;
-           SolverCG < Vector<double> >  solver_CG (solver_control, GVM);
+               // Perform static condensation to make K_con,
+               // and put K_pJ^{-1} in the original K_pJ block.
+               // That is, we make K_store.
+               assemble_sc();
+
+               // K_con du = F_con
+               // with F_con = F_u + K_up [- K_Jp^-1 F_j + K_bar F_p]
+               // Assemble the RHS vector to solve for du
+               // A_J = K_pJ^-1 F_p
+               tangent_matrix.block(p_dof, J_dof).vmult(A.block(J_dof),
+                               system_rhs.block(p_dof));
+               // B_J = K_JJ  K_pJ^-1  F_p
+               tangent_matrix.block(J_dof, J_dof).vmult(B.block(J_dof),
+                               A.block(J_dof));
+               // A_J = F_J - K_JJ  K_pJ^-1  F_p
+               A.block(J_dof).equ(1.0, system_rhs.block(J_dof), -1.0, B.block(J_dof));
+               // A_p = K_Jp^-1 [  F_J - K_JJ  K_pJ^-1  F_p ]
+               tangent_matrix.block(p_dof, J_dof).Tvmult(A.block(p_dof),
+                               A.block(J_dof));
+               // A_u = K_up  K_Jp^-1 [  F_J - K_JJ  K_pJ^-1  F_p ]
+               tangent_matrix.block(u_dof, p_dof).vmult(A.block(u_dof),
+                               A.block(p_dof));
+               // F_con = F_u -  K_up  K_Jp^-1 [  F_J - K_JJ  K_pJ^-1  F_p ]
+               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 int solver_its = 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);
+
+                       // We've chosen a SSOR preconditioner as it appears to provide
+                       // the fastest solver convergence characteristics for this problem.
+                       PreconditionSSOR<SparseMatrix<double> > preconditioner;
+                       preconditioner.initialize(tangent_matrix.block(u_dof, u_dof),
+                                       parameters.ssor_relaxation);
+
+                       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") {
+                       // Otherwise if the problem is small enough, a direct solver
+                       // can be utilised.
+                       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
+                       throw(ExcMessage("Linear solver type not implemented"));
+               timer.leave_subsection();
+       }
 
-           // We've chosen a SSOR preconditioner as it appears to provide
-           // the fastest solver convergence characteristics for this problem.
-           PreconditionSSOR <SparseMatrix<double> > preconditioner;
-           preconditioner.initialize (tangent_matrix.block(u_dof, u_dof), parameters.ssor_relaxation);
+       // distribute the constrained dof back to the Newton update
+       constraints.distribute(newton_update);
 
-           solver_CG.solve (tangent_matrix.block(u_dof, u_dof),
-                            newton_update.block(u_dof),
-                            residual.block(u_dof),
-                            preconditioner);
+       timer.enter_subsection("Linear solver postprocessing");
+       std::cout << " PP " << std::flush;
 
-           lin_it = solver_control.last_step();
-           lin_res = solver_control.last_value();
+       // Now that we've solved the displacement problem, we can post-process
+       // to get the dilatation solution from the substitution
+       // dJ = KpJ^{-1} (F_p - K_pu du )
+       {
+               // A_p  = K_pu du
+               tangent_matrix.block(p_dof, u_dof).vmult(A.block(p_dof),
+                               newton_update.block(u_dof));
+               // A_p  = -K_pu du
+               A.block(p_dof) *= -1.0;
+               // A_p  = F_p - K_pu du
+               A.block(p_dof) += system_rhs.block(p_dof);
+               // d_J = K_pJ^{-1} [ F_p - K_pu du ]
+               tangent_matrix.block(p_dof, J_dof).vmult(newton_update.block(J_dof),
+                               A.block(p_dof));
        }
-       else if (parameters.type_lin == "Direct")
+       // and finally we solve for the pressure update with the substitution
+       // dp = KJp^{-1} ( R_J - K_JJ dJ )
        {
-           // Otherwise if the problem is small enough, a direct solver
-           // can be utilised.
-           SparseDirectUMFPACK  A_direct;
-           A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
-           A_direct.vmult (newton_update.block(u_dof),
-                           residual.block(u_dof));
-
-           lin_it = 1;
-           lin_res = 0.0;
+               // A_J = K_JJ dJ
+               tangent_matrix.block(J_dof, J_dof).vmult(A.block(J_dof),
+                               newton_update.block(J_dof));
+               // A_J = -K_JJ dJ
+               A.block(J_dof) *= -1.0;
+               // A_J = F_J - K_JJ dJ
+               A.block(J_dof) += system_rhs.block(J_dof);
+               // dp = K_Jp^{-1}   [F_J - K_JJ dJ]
+               tangent_matrix.block(p_dof, J_dof).Tvmult(newton_update.block(p_dof),
+                               A.block(J_dof));
        }
-       else throw (ExcMessage("Linear solver type not implemented"));
+
+       // distribute the constrained dof back to the Newton update
+       constraints.distribute(newton_update);
+
        timer.leave_subsection();
-    }
-
-    timer.enter_subsection("Linear solver postprocessing");
-    std::cout << " PP " << std::flush;
-    // Now that we've solved the displacement problem, we can post-process
-    // to get the dilatation solution from the substitution
-    // dt = Kpt^{-1} ( Rp - Kpu du )
-    {
-        tangent_matrix.block(p_dof, u_dof).vmult (A.block(p_dof), newton_update.block(u_dof));
-        A.block(p_dof) *= -1.0;
-        A.block(p_dof) += residual.block(p_dof);
-        tangent_matrix.block(p_dof, t_dof).Tvmult (newton_update.block(t_dof), A.block(p_dof));
-    }
-    // and finally we solve for the pressure update with the substitution
-    // dp = Ktp^{-1} ( Rt - Ktt dt )
-    {
-        tangent_matrix.block(t_dof, t_dof).vmult (A.block(t_dof), newton_update.block(t_dof));
-       A.block(t_dof) *= -1.0;
-        A.block(t_dof) += residual.block(t_dof);
-        tangent_matrix.block(p_dof, t_dof).vmult (newton_update.block(p_dof), A.block(t_dof));
-    }
-    timer.leave_subsection();
-
-    return std::make_pair(lin_it, lin_res);
+
+       return std::make_pair(lin_it, lin_res);
 }
 
 // @sect4{Solid::assemble_system_SC}
 // 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 discontinous, the SC operation
+// pressure and dilatation variables are discontinuous, the SC operation
 // can be done on a per-cell basis and we can produce the inverse of the
 // block-diagonal K_{pt} block by inverting the local blocks. We can
 // again use TBB to do this since each operation will be independent of
 // one another.
-template <int dim>
-void Solid<dim>::assemble_SC  (void)
-{
-    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_t.size()); // Initialise members of per_task_data to the correct sizes.
-    ScratchData_SC scratch_data;
-
-    WorkStream::run (  dof_handler_ref.begin_active(),
-                       dof_handler_ref.end(),
-                       *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>::assemble_sc(void) {
+       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()); // Initialise members of per_task_data to the correct sizes.
+       ScratchData_SC scratch_data;
+
+       // Using TBB, we assemble the contributions to add to
+       // K_uu to form K_con from each elements contributions.
+       // These contributions are then added to the glabal stiffness
+       // matrix.
+       WorkStream::run(dof_handler_ref.begin_active(), dof_handler_ref.end(),
+                       *this, &Solid::assemble_sc_one_cell,
+                       &Solid::copy_local_to_global_sc, scratch_data, per_task_data);
+
+       timer.leave_subsection();
 }
 
-// We need to describe how to add the local contribution to the tangent matrix.
-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));
+// We need to describe how to add the local contributions
+// to K to form K_store
+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));
 }
 
 // Now we describe the static condensation process.
-template <int dim>
-void Solid<dim>::assemble_SC_one_cell (const typename DoFHandler<dim>::active_cell_iterator & cell,
-                                       ScratchData_SC & scratch,
-                                       PerTaskData_SC & data)
-{
-    // As per usual, we must first find out which global numbers the
-    // degrees of freedom on this cell have and reset some data structures
-    data.reset();
-    scratch.reset();
-    cell->get_dof_indices (data.local_dof_indices);
-
-    // Currently the the local stifness matrix K_e is of the form
-    //  | K_uu  |   K_up   |   0   |
-    //  | K_pu  |     0    |  K_pt |
-    //  |   0   |   K_tp   |  K_tt |
-    //
-    // We now need to modify it such that it appear as
-    //  | K'_uu |   K_up   |     0     |
-    //  | K_pu  |     0    |   K_pt^-1 |
-    //  |   0   |   K_tp   |   K_tt    |
-    // with K'_uu = K_uu + Kup Ktp^{-1} Ktt Kpt^{-1} Kpu
-    //
-    // At this point, we need to take note of the fact that
-    // global data already exists in the K_uu, K_pt, K_tp subblocks.
-    // 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 K_uu 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 subblocks.
-    //
-    // So the intermediate matrix that we need to get from what we have in K_uu and what we
-    // are actually wanting is:
-    //  | K'_uu - K_uu |   0   |        0        |
-    //  |       0      |   0   |  K_pt^-1 - K_pt |
-    //  |       0      |   0   |        0        |
-    //
-    // This is the strategy we will employ to get the subblocks we want:
-    // K'_{uu}: Since we don't have access to K_{uu}^h, but we know its contribution is added to the global
-    //        K_{uu} matrix, we just want to add the element wise static-condensation
-    //        K'_{uu}^h = K_{uu}^h + K_{up}^h K_{tp}^{-1, h} K_{tt}^h K_{pt}^{-1, h} K_{pu}^h
-    //        Since we already have K_uu^h in the system matrix, we just need to do the following
-    //        K'_{uu}^h == (K_{uu}^h += K_{up}^h K_{tp}^{-1}^h K_{tt}^h K_{pt}^{-1, h} K_{pu}^h)
-    // K_{pt}^-1: Similarly, K_pt exists in the subblock. Since the copy operation is a += operation, we need
-    //          to subtract the existing K_pt submatrix in addition to "adding" that which we wish to
-    //          replace it with.
-    // K_{tp}^-1: Since the global matrix is symmetric, this block is the same as the one above
-    //          and we can simply use K_pt^-1 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
-    AdditionalTools::extract_submatrix(data.local_dof_indices,
-                                       data.local_dof_indices,
-                                       tangent_matrix,
-                                       data.K_orig);
-    // and next the local matrices for K_{pu}, K_{pt} and K_{tt}
-    AdditionalTools::extract_submatrix(element_indices_p,
-                                       element_indices_u,
-                                       data.K_orig,
-                                       data.K_pu);
-    AdditionalTools::extract_submatrix(element_indices_p,
-                                       element_indices_t,
-                                       data.K_orig,
-                                       data.K_pt);
-    AdditionalTools::extract_submatrix(element_indices_t,
-                                       element_indices_t,
-                                       data.K_orig,
-                                       data.K_tt);
-
-    // To get the inverse of K_{pt}, we invert it directly.
-    // This operation is relatively inexpensive since
-    // K_{pt} is block-diagonal.
-    data.K_pt_inv.invert(data.K_pt);
-
-    // Now we can make condensation terms to add to the
-    // K_{uu} block and put them in the cell local matrix
-    data.K_pt_inv.mmult(data.A, data.K_pu);
-    data.K_tt.mmult(data.B, data.A);
-    data.K_pt_inv.Tmmult(data.C, data.B);
-    data.K_pu.Tmmult(data.K_con, data.C);
-    AdditionalTools::replace_submatrix(element_indices_u,
-                                       element_indices_u,
-                                       data.K_con,
-                                       data.cell_matrix);
-
-    // Next we place K_{pt}^-1 in the K_{pt} block for post-processing
-    // Note again that we need to remove the K_pt contribution that
-    // already exists there.
-    data.K_pt_inv.add (-1.0, data.K_pt);
-    AdditionalTools::replace_submatrix(element_indices_p,
-                                       element_indices_t,
-                                       data.K_pt_inv,
-                                       data.cell_matrix);
+template<int dim>
+void Solid<dim>::assemble_sc_one_cell(
+               const typename DoFHandler<dim>::active_cell_iterator & cell,
+               ScratchData_SC & scratch, PerTaskData_SC & data) {
+       // As per usual, we must first find out which global numbers the
+       // degrees of freedom on this cell have and reset some data structures
+       data.reset();
+       scratch.reset();
+       cell->get_dof_indices(data.local_dof_indices);
+
+       // We now extract the contribution of
+       // the  dof associated with the current cell
+       // to the global stiffness matrix.
+       // The discontinuous nature of the p and J
+       // interpolations mean that their is no
+       // coupling of the local contributions at the
+       // global level. This is not the case with the u dof.
+       // In other words, k_Jp, k_pJ and k_JJ, when extracted
+       // from the global stiffness matrix are the element
+       // contributions. This is not the case for k_uu.
+
+       // Currently the matrix corresponding to
+       // the dof associated with the current element
+       // (denoted somewhat loosely as k) is of the form
+       //  | k_uu  |   k_up   |   0   |
+       //  | k_pu  |     0    |  k_pJ |
+       //  |   0   |   k_Jp   |  k_JJ |
+       //
+       // We now need to modify it such that it appear as
+       //  | k_con |   k_up   |     0     |
+       //  | k_pu  |     0    |   k_pJ^-1 |
+       //  |   0   |   k_Jp   |   k_JJ    |
+       // with k_con = k_uu + k_bbar
+       // where
+       // k_bbar = k_up k_bar k_pu
+       // and
+       // k_bar = k_Jp^{-1} k_JJ kpJ^{-1}
+       //
+       // At this point, we need to take note of the fact that
+       // global data already exists in the K_uu, K_pt, K_tp subblocks.
+       // 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 K_uu 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 subblocks.
+       //
+       // So the intermediate matrix that we need to get from what we have in K_uu and what we
+       // are actually wanting is:
+       //  | K'_uu - K_uu |   0   |        0        |
+       //  |       0      |   0   |  K_pt^-1 - K_pt |
+       //  |       0      |   0   |        0        |
+       //
+       // This is the strategy we will employ to get the subblocks we want:
+       // K'_{uu}: Since we don't have access to K_{uu}^h, but we know its contribution is added to the global
+       //        K_{uu} matrix, we just want to add the element wise static-condensation
+       //        K'_{uu}^h = K_{uu}^h + K_{up}^h K_{tp}^{-1, h} K_{tt}^h K_{pt}^{-1, h} K_{pu}^h
+       //        Since we already have K_uu^h in the system matrix, we just need to do the following
+       //        K'_{uu}^h == (K_{uu}^h += K_{up}^h K_{tp}^{-1}^h K_{tt}^h K_{pt}^{-1, h} K_{pu}^h)
+       // K_{pt}^-1: Similarly, K_pt exists in the subblock. Since the copy operation is a += operation, we need
+       //          to subtract the existing K_pt submatrix in addition to "adding" that which we wish to
+       //          replace it with.
+       // K_{tp}^-1: Since the global matrix is symmetric, this block is the same as the one above
+       //          and we can simply use K_pt^-1 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
+
+       // extract k for the dof associated with the current element
+       AdditionalTools::extract_submatrix(data.local_dof_indices,
+                       data.local_dof_indices, tangent_matrix, data.k_orig);
+       // and next the local matrices for k_pu, k_pJ and k_JJ
+       AdditionalTools::extract_submatrix(element_indices_p, element_indices_u,
+                       data.k_orig, data.k_pu);
+       AdditionalTools::extract_submatrix(element_indices_p, element_indices_J,
+                       data.k_orig, data.k_pJ);
+       AdditionalTools::extract_submatrix(element_indices_J, element_indices_J,
+                       data.k_orig, data.k_JJ);
+
+       // To get the inverse of k_pJ, we invert it directly.
+       // This operation is relatively inexpensive since
+       // k_pJ is block-diagonal.
+       data.k_pJ_inv.invert(data.k_pJ);
+
+       // Now we can make condensation terms to add to the
+       // k_uu block and put them in the cell local matrix
+       // A = k_pJ^-1 k_pu
+       data.k_pJ_inv.mmult(data.A, data.k_pu);
+       // B = k_JJ k_pJ^-1 k_pu
+       data.k_JJ.mmult(data.B, data.A);
+       // C = k_Jp^-1 k_JJ k_pJ^-1 k_pu
+       data.k_pJ_inv.Tmmult(data.C, data.B);
+       // k_bbar = k_up k_Jp^-1 k_JJ k_pJ^-1 k_pu
+       data.k_pu.Tmmult(data.k_bbar, data.C);
+       AdditionalTools::replace_submatrix(element_indices_u, element_indices_u,
+                       data.k_bbar, data.cell_matrix);
+
+       // Next we place k_{pJ}^-1 in the k_{pJ} block for post-processing.
+       // Note again that we need to remove the k_pJ contribution that
+       // already exists there.
+       data.k_pJ_inv.add(-1.0, data.k_pJ);
+       AdditionalTools::replace_submatrix(element_indices_p, element_indices_J,
+                       data.k_pJ_inv, data.cell_matrix);
 }
 
 // @sect4{Solid::output_results}
 // Here we present how the results are written to file to be viewed
-// using Paraview. The method is similar to that shown in previous
+// using ParaView. The method is similar to that shown in previous
 // tutorials so will not be discussed in detail.
-template <int dim>
-void Solid<dim>::output_results(void)
-{
-    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.push_back ("pressure");
-    solution_name.push_back ("dilatation");
-
-    data_out.attach_dof_handler (dof_handler_ref);
-    data_out.add_data_vector (solution_n,
-                             solution_name,
-                             DataOut<dim>::type_dof_data,
-                             data_component_interpretation);
-
-    // 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 ourselves.
-    // We first need to copy the solution to a temporary vector and then
-    // create the Eularian mapping. We also specify the polynomial degree
-    // to the DataOut object in order to produce a more refined output dataset
-    // when higher order polynomials are used.
-    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,
-                                     soln,
-                                     dof_handler_ref);
-    data_out.build_patches (q_mapping,
-                            degree);
-
-    std::ostringstream filename;
-    filename << "solution-"
-             << time.get_timestep()
-             << ".vtk";
-
-    std::ofstream output (filename.str().c_str());
-    data_out.write_vtk (output);
+template<int dim>
+void Solid<dim>::output_results(void) {
+       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.push_back("pressure");
+       solution_name.push_back("dilatation");
+
+       data_out.attach_dof_handler(dof_handler_ref);
+       data_out.add_data_vector(solution_n, solution_name,
+                       DataOut<dim>::type_dof_data, data_component_interpretation);
+
+       // 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 ourselves.
+       // We first need to copy the solution to a temporary vector and then
+       // create the Eularian mapping. We also specify the polynomial degree
+       // to the DataOut object in order to produce a more refined output dataset
+       // when higher order polynomials are used.
+       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, soln, dof_handler_ref);
+       data_out.build_patches(q_mapping, degree);
+
+       std::ostringstream filename;
+       filename << "solution-" << time.get_timestep() << ".vtk";
+
+       std::ofstream output(filename.str().c_str());
+       data_out.write_vtk(output);
 }
 
 // @sect3{Main function}
 // Lastly we provide the main driver function which appears
 // no different to the other tutorials.
-int main (void)
-{
-    try
-    {
-       deallog.depth_console (0);
-
-       Solid<3> solid_3d ("parameters.prm");
-       solid_3d.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;
+int main(void) {
+       try {
+               deallog.depth_console(0);
+
+               Solid<3> solid_3d("parameters.prm");
+               solid_3d.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;
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.