]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix accidental checkin of a merge commit.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 22 Apr 2011 21:41:27 +0000 (21:41 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 22 Apr 2011 21:41:27 +0000 (21:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@23625 0785d39b-7218-0410-832d-ea1e28bc413d

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

index e8c203822ee154cf27968d69e53f83e3729c8220..55cf8cdbce25002bd6f27f43a51d747e56404d41 100644 (file)
@@ -2,7 +2,7 @@
 /* Authors: Jean-Paul Pelteret, University of Cape Town,            */
 /*          Andrew McBride, University of Erlangen-Nuremberg, 2010  */
 /*                                                                  */
-/*    Copyright (C) 2010 by the deal.II authors                     */
+/*    Copyright (C) 2010, 2011 by the deal.II authors                     */
 /*                        & Jean-Paul Pelteret and Andrew McBride   */
 /*                                                                  */
 /*    This file is subject to QPL and may not be  distributed       */
@@ -63,7 +63,7 @@ struct FESystem
 {
     int poly_degree;
     int quad_order;
-    
+
     static void declare_parameters (ParameterHandler &prm);
     void parse_parameters (ParameterHandler &prm);
 };
@@ -76,7 +76,7 @@ void FESystem::declare_parameters (ParameterHandler &prm)
                          "1",
                          Patterns::Integer(),
                          "Displacement system polynomial order");
-       
+
        prm.declare_entry("Quadrature order",
                          "2",
                          Patterns::Integer(),
@@ -101,7 +101,7 @@ struct Geometry
     int global_refinement;
     double scale;
     double p_p0;
-    
+
     static void declare_parameters (ParameterHandler &prm);
     void parse_parameters (ParameterHandler &prm);
 };
@@ -114,12 +114,12 @@ void Geometry::declare_parameters (ParameterHandler &prm)
                          "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"),
@@ -144,7 +144,7 @@ struct Materials
 {
     double nu;
     double mu;
-    
+
     static void declare_parameters (ParameterHandler &prm);
     void parse_parameters (ParameterHandler &prm);
 };
@@ -157,7 +157,7 @@ void Materials::declare_parameters (ParameterHandler &prm)
                          "0.49",
                          Patterns::Double(),
                          "Poisson's ratio");
-       
+
        prm.declare_entry("Shear modulus",
                          "1.0e6",
                          Patterns::Double(),
@@ -183,7 +183,7 @@ struct LinearSolver
     double tol_lin;
     double max_iterations_lin;
     double ssor_relaxation;
-    
+
     static void declare_parameters (ParameterHandler &prm);
     void parse_parameters (ParameterHandler &prm);
 };
@@ -196,17 +196,17 @@ void LinearSolver::declare_parameters (ParameterHandler &prm)
                          "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(),
@@ -233,7 +233,7 @@ struct NonlinearSolver
     unsigned int max_iterations_NR;
     double tol_f;
     double tol_u;
-    
+
     static void declare_parameters (ParameterHandler &prm);
     void parse_parameters (ParameterHandler &prm);
 };
@@ -246,12 +246,12 @@ void NonlinearSolver::declare_parameters (ParameterHandler &prm)
                          "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(),
@@ -276,7 +276,7 @@ struct Time
 {
     double end_time;
     double delta_t;
-    
+
     static void declare_parameters (ParameterHandler &prm);
     void parse_parameters (ParameterHandler &prm);
 };
@@ -289,7 +289,7 @@ void Time::declare_parameters (ParameterHandler &prm)
                          "1",
                          Patterns::Double(),
                          "End time");
-       
+
        prm.declare_entry("Time step size",
                          "0.1",
                          Patterns::Double(),
@@ -317,10 +317,10 @@ struct AllParameters
        public LinearSolver,
        public NonlinearSolver,
        public Time
-       
+
 {
     AllParameters (const std::string & input_file);
-    
+
     static void declare_parameters (ParameterHandler &prm);
     void parse_parameters (ParameterHandler &prm);
 };
@@ -363,20 +363,20 @@ void extract_submatrix(const std::vector< unsigned int > &row_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(), ExcIndexRange(row, 0, matrix.m()));
-       
+
        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(), ExcIndexRange(col, 0, matrix.n()));
-           
+
            sub_matrix(sub_row,sub_col) = matrix(row, col);
        }
     }
@@ -392,17 +392,17 @@ void replace_submatrix(const std::vector< unsigned int > &row_index_set,
     Assert (n_rows_submatrix<=sub_matrix.m(), ExcIndexRange(n_rows_submatrix, 0, sub_matrix.m()));
     const unsigned int n_cols_submatrix = column_index_set.size();
     Assert (n_cols_submatrix<=sub_matrix.n(), ExcIndexRange(n_cols_submatrix, 0, sub_matrix.n()));
-    
+
     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(), ExcIndexRange(row, 0, matrix.m()));
-       
+
        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(), ExcIndexRange(col, 0, matrix.n()));
-           
+
            matrix(row, col) = sub_matrix(sub_row, sub_col);
-           
+
        }
     }
 }
@@ -421,13 +421,13 @@ public:
          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;
@@ -449,7 +449,7 @@ public:
          kappa_0 (lambda + 2.0/3.0*mu)
     { }
     virtual ~Material_NH (void) {};
-    
+
     // Stress and constitutive tensors
     virtual SymmetricTensor<2, dim> get_T (const double & J,
                                           const SymmetricTensor <2, dim> & B)
@@ -457,7 +457,7 @@ public:
        const double dW_dJ  = get_dU_dtheta (J);
        return mu_0*B + dW_dJ*J*I;
     }
-    
+
     virtual SymmetricTensor<4, dim> get_JC (const double & J,
                                            const SymmetricTensor <2, dim> & B)
     {
@@ -465,17 +465,17 @@ public:
        const double d2W_dJ2 = get_d2U_dtheta2 (J);
        return  J*(  (dW_dJ + J*d2W_dJ2)*IxI - (2.0*dW_dJ)*II  );
     }
-    
+
     // Volumetric quantities methods
     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));}
-    
+
 protected:
     // Material properties
     const double lambda_0; // Lame modulus
     const double mu_0;     // Shear modulus
     const double kappa_0;  // Bulk modulus
-    
+
     static SymmetricTensor<2, dim> const I;
     static SymmetricTensor<4, dim> const IxI;
     static SymmetricTensor<4, dim> const II;
@@ -497,17 +497,17 @@ public:
          pressure_n (0.0)
     { }
     virtual ~PointHistory (void) {delete material;}
-    
+
     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);
     }
-    
+
     // Total Variables
     void update_values (const Tensor<2, dim> & grad_u_n,
                        const double & pressure,
@@ -518,11 +518,11 @@ public:
        J     = determinant(F);
        F_inv = invert(F);
        B_bar = std::pow(get_J(), -2.0/3.0) * symmetrize ( F* transpose (F) );
-       
+
         // Precalculated pressure, dilatation
        pressure_n = pressure;
        dilatation_n = dilatation;
-       
+
         // Now that all the necessary variables are set, we can update the stress tensors
         // Stress update can only update the stresses once the
         // dilatation has been set as p = p(d)
@@ -530,68 +530,68 @@ public:
         T_iso = dev_P*get_T_bar(); // Note: T_iso depends on T_bar
         T_vol = get_pressure()*get_J()*I;
     }
-    
+
     // Displacement and strain
     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;}
     const SymmetricTensor <2,dim> & get_B_bar (void) const {return B_bar;}
-    
+
     // Volumetric terms
     double get_dU_dtheta (void) {
        return material->get_dU_dtheta(get_dilatation());
     }
-    
+
     double get_d2U_dtheta2 (void) {
        return material->get_d2U_dtheta2(get_dilatation());
     }
-    
+
     // Stress
     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;};
-    
+
     // Tangent matrices
     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> CC_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*CC_bar*dev_P;
     }
-    
+
     SymmetricTensor <4,dim> get_C_vol(void)
     {
        const double & p = get_pressure();
        const double & J = get_J();
        return p*J*(IxI - 2.0*II);
     }
-    
+
 private:
     // === MATERIAL ===
     Material_NH <dim>* material;
-    
+
     // ==== VOLUME, DISPLACEMENT AND STRAIN VARIABLES ====
     double                  dilatation_n;   // Current dilatation
     double                  J;
     Tensor <2,dim>         F_inv;
     SymmetricTensor <2,dim> B_bar;
     SymmetricTensor <2,dim> E;
-    
+
     // ==== STRESS VARIABLES ====
     double                  pressure_n; // Current pressure
     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;}
-    
+
     // Basis tensors
     static SymmetricTensor<2, dim> const I;
     static SymmetricTensor<4, dim> const IxI;
@@ -617,36 +617,36 @@ public:
     Solid (const std::string & input_file);
     virtual ~Solid (void);
     void run (void);
-    
+
 private:
-    
+
     // === DATA STRUCTS ===
-    
+
     struct 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;
         }
     };
-    
+
     struct 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)
@@ -661,7 +661,7 @@ private:
              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(),
@@ -670,7 +670,7 @@ private:
            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) {
@@ -680,35 +680,35 @@ private:
                }
            }
        }
-       
+
     };
-    
+
     struct PerTaskData_F
     {
        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)
        { }
-       
+
        void reset (void) { cell_rhs = 0.0; }
     };
-    
+
     struct 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;
        std::vector< Vector<double> > rhs_values;
-       
+
        // 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,
@@ -728,7 +728,7 @@ private:
              rhs_values   (qf_cell.size(),
                  Vector<double>(dim))
        {  }
-       
+
        ScratchData_F ( const ScratchData_F & rhs )
            :
              fe_values_ref ( rhs.fe_values_ref.get_fe(),
@@ -741,7 +741,7 @@ private:
              symm_grad_Nx (rhs.symm_grad_Nx),
              rhs_values (rhs.rhs_values)
        {  }
-       
+
        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) {
@@ -751,14 +751,14 @@ private:
                }
            }
        }
-       
+
     };
-    
+
     struct 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;
@@ -771,7 +771,7 @@ private:
        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,
@@ -787,12 +787,12 @@ private:
              B (n_t, n_u),
              C (n_p, n_u)
        {  }
-       
+
        // Choose not to reset any data
        // The matrix extraction and replacement tools will take care of this
        void reset(void) { }
     };
-    
+
     // Dummy struct for TBB
     struct ScratchData_SC
     {
@@ -800,14 +800,14 @@ private:
        ScratchData_SC (const ScratchData_SC & rhs) { }
        void reset (void) { }
     };
-    
+
     // Dummy struct for TBB
     struct PerTaskData_UQPH
     {
        PerTaskData_UQPH (void) { }
        void reset(void) { }
     };
-    
+
     struct ScratchData_UQPH
     {
        FEValues<dim> fe_values_ref;
@@ -815,7 +815,7 @@ private:
        std::vector <double> solution_values_p_total;
        std::vector <double> solution_values_t_total;
        const BlockVector <double> & solution_total;
-       
+
        ScratchData_UQPH (const FiniteElement <dim> & fe_cell,
                          const QGauss <dim> & qf_cell,
                          const UpdateFlags uf_cell,
@@ -829,7 +829,7 @@ private:
              solution_values_t_total (qf_cell.size()),
              solution_total (solution_total)
        { }
-       
+
        ScratchData_UQPH (const ScratchData_UQPH & rhs)
            :
              fe_values_ref (rhs.fe_values_ref.get_fe(),
@@ -840,7 +840,7 @@ private:
              solution_values_t_total (rhs.solution_values_t_total),
              solution_total (rhs.solution_total)
        { }
-       
+
        void reset (void)
        {
            // Is this necessary? Won't the call to fe_values.get_gradient overwrite this data?
@@ -852,15 +852,15 @@ private:
            }
        }
     };
-    
+
     // === METHODS ===
-    
+
     /// \brief Print out a greeting for the user
     void make_grid (void);
     /// \brief Setup the Finite Element system to be solved
     void system_setup (void);
     void determine_component_extractors(void);
-    
+
     /// \brief 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,
@@ -880,7 +880,7 @@ private:
     /// \brief Apply Dirichlet boundary values
     void make_constraints (const int & it_nr,
                           ConstraintMatrix & constraints);
-    
+
     //    /// \brief Setup the quadrature point history for each cell
     void setup_qph(void);
     //    /// \brief Update the quadrature points stress and strain values, and fibre directions
@@ -892,61 +892,61 @@ private:
     /// \brief Solve for the displacement using a Newton-Rhapson method
     void solve_nonlinear_timestep (BlockVector <double> & solution_delta);
     void solve_linear_system (BlockVector <double> & newton_update);
-    
+
     /// \brief Error measurement
     void get_error_res (const BlockVector <double> & residual, BlockVector <double> & error_res);
     void get_error_update (const BlockVector <double> & newton_update, BlockVector <double> & error_update);
     double get_error_dil (void);
-    
+
     // Solution
     BlockVector <double> get_solution_total (const BlockVector <double> & solution_delta);
-    
+
     // Postprocessing
     void output_results(void);
-    
+
     // === ATTRIBUTES ===
     // Parameters
     Parameters::AllParameters parameters;
-    
+
     // Geometry
     Triangulation<dim> triangulation; // Describes the triangulation
-    
+
     // Time
     Time time;
     TimerOutput timer;
-    
+
     // === Quadrature points ===
     std::vector< PointHistory <dim> > quadrature_point_history; // Quadrature point history
-    
+
     // === Finite element system ===
     DoFHandler<dim>     dof_handler_ref; // Describes the degrees of freedom
     const unsigned int  degree;
     const FESystem<dim> fe; // Describes the global FE system
-    
+
     unsigned int dofs_per_cell; // Number of degrees of freedom on each cell
     const FEValuesExtractors::Vector u_fe;
     const FEValuesExtractors::Scalar p_fe;
     const FEValuesExtractors::Scalar t_fe;
-    
+
     // Block description
     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;
-    
+
     // === Quadrature ===
     QGauss<dim> qf_cell; // Cell quadrature formula
     QGauss<dim-1> qf_face; // Face quadrature formula
     unsigned int n_q_points; // Number of quadrature points in a cell
     unsigned int n_q_points_f; // Number of quadrature points in a face
-    
+
     // === Stiffness matrix setup ====
     ConstraintMatrix constraints; // Matrix to keep track of all constraints
     BlockSparsityPattern sparsity_pattern; // Sparsity pattern for the stiffness matrix
@@ -999,22 +999,17 @@ void Solid<dim>::run (void)
     system_setup ();
     output_results (); // Output initial grid position
     time.increment();
-    
+
     BlockVector <double> solution_delta (dofs_per_block);
     solution_delta.collect_sizes ();
-<<<<<<< .mine
 
     while (time.current() <= time.end()) {
-=======
-    
-    while (time.current() <= time.end()) {
->>>>>>> .r23597
        solution_delta = 0.0;
-       
+
        // Solve step and update total solution vector
        solve_nonlinear_timestep (solution_delta);
        solution_n += solution_delta;
-       
+
        output_results ();
        time.increment();
     }
@@ -1029,11 +1024,11 @@ void Solid<dim>::make_grid (void)
                                    Point<dim> (1.0, 1.0, 1.0),
                                    true );
     GridTools::scale (parameters.scale, triangulation);
-    
+
     // Need to refine at least once for the indentation problem
     if (parameters.global_refinement == 0) triangulation.refine_global (1);
     else triangulation.refine_global (parameters.global_refinement);
-    
+
     // Apply different BC's to a patch on the top surface
     typename Triangulation<dim>::active_cell_iterator
            cell = triangulation.begin_active(),
@@ -1062,25 +1057,25 @@ template <int dim>
 void Solid<dim>::system_setup (void)
 {
     timer.enter_subsection ("Setup system");
-    
+
     // Number of components per block
     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
-    
+
     // Setup DOF handler
     dof_handler_ref.distribute_dofs (fe);
     DoFRenumbering::Cuthill_McKee (dof_handler_ref);
     DoFRenumbering::component_wise (dof_handler_ref, block_component);
     // Count number of dofs per block
     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;
-    
+
     // the global system matrix will have the following structure
     //      | K'_uu |   K_up    |     0     |         | dU_u |         | dR_u |
     // K =  | K_pu  |   K_tt^-1 |   K_pt^-1 | , dU =  | dU_p | , dR =  | dR_p |
@@ -1096,52 +1091,52 @@ void Solid<dim>::system_setup (void)
             else {
                coupling[ii][jj] = DoFTools::always;
                coupling[jj][ii] = DoFTools::always;
-            }  
+            }
         }
     }
-    
+
     // Setup system 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();
-       
+
        DoFTools::make_sparsity_pattern (dof_handler_ref, csp);
        //        DoFTools::make_sparsity_pattern (dof_handler_ref, csp, constraints, false);
        //        DoFTools::make_sparsity_pattern (dof_handler_ref, coupling, csp, constraints, false);
         sparsity_pattern.copy_from (csp);
     }
-    
-    
+
+
     tangent_matrix.reinit (sparsity_pattern);
-    
+
     // Setup storage vectors
     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; // Dilatation is 1 in the initial configuration
-    
+
     // Set up the quadrature point history
     setup_qph ();
-    
+
     timer.leave_subsection();
 }
 
@@ -1152,7 +1147,7 @@ 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) {
        // 0 = u, 1 = p, 2 = dilatation interpolation fields
        const unsigned int k_group = fe.system_to_base_index(k).first.first;
@@ -1176,34 +1171,34 @@ template <int dim>
 void Solid<dim>::setup_qph (void)
 {
     std::cout << "    Setting up quadrature point data..." << std::endl;
-    
+
     {
        typename Triangulation<dim>::active_cell_iterator
                cell = triangulation.begin_active(),
                endc = triangulation.end();
-       
+
        unsigned int our_cells = 0;
        for (; cell != endc; ++cell) {
            cell->clear_user_pointer();
            ++our_cells;
        }
-       
+
        {
            std::vector<PointHistory <dim> > tmp;
            tmp.swap(quadrature_point_history);
        }
-       
+
        quadrature_point_history.resize(our_cells * n_q_points);
-       
+
        unsigned int history_index = 0;
        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());
     }
-    
+
     // Setup initial data
     typename DoFHandler<dim>::active_cell_iterator
            cell = dof_handler_ref.begin_active(),
@@ -1212,7 +1207,7 @@ void Solid<dim>::setup_qph (void)
        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 );
@@ -1226,7 +1221,7 @@ void Solid<dim>::update_qph_incremental (const BlockVector <double> & solution_d
 {
     timer.enter_subsection("Update QPH data");
     std::cout << "Update QPH data..."<< std::endl;
-    
+
     // Get total solution as it stands at this update increment
     const BlockVector <double> solution_total = get_solution_total(solution_delta);
     const UpdateFlags uf_UQPH ( update_values | update_gradients );
@@ -1235,7 +1230,7 @@ void Solid<dim>::update_qph_incremental (const BlockVector <double> & solution_d
                                        qf_cell,
                                        uf_UQPH,
                                        solution_total);
-    
+
     WorkStream::run (  dof_handler_ref.begin_active(),
                     dof_handler_ref.end(),
                     *this,
@@ -1243,7 +1238,7 @@ void Solid<dim>::update_qph_incremental (const BlockVector <double> & solution_d
                     &Solid::copy_local_to_global_UQPH,
                     scratch_data_UQPH,
                     per_task_data_UQPH);
-    
+
     timer.leave_subsection();
 }
 
@@ -1255,17 +1250,17 @@ void Solid<dim>::update_qph_incremental_one_cell (const typename DoFHandler<dim>
     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());
-    
+
     // 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);
-    
+
     // === UPDATE DATA AT EACH GAUSS POINT ===
     // Update displacement and deformation gradient at all quadrature points
     for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
@@ -1283,17 +1278,17 @@ void Solid<dim>::solve_nonlinear_timestep (BlockVector <double> & solution_delta
     std::cout
            << "Timestep " << time.get_timestep()
            << std::endl;
-    
+
     // Newton update vector
     BlockVector <double> newton_update (dofs_per_block);
     newton_update.collect_sizes ();
-    
+
     // Solution error vectors
     BlockVector <double> soln_error_res (dofs_per_block); // Holds the true residual vector
     BlockVector <double> soln_error_update (dofs_per_block); // Holds the update error vector
     soln_error_res.collect_sizes ();
     soln_error_update .collect_sizes ();
-    
+
     double res_u = 0.0, res_f = 0.0;
     double res_u_0 = 1.0, res_f_0 = 1.0;
     for (unsigned int it_nr=0; it_nr < parameters.max_iterations_NR; ++ it_nr)
@@ -1302,10 +1297,10 @@ void Solid<dim>::solve_nonlinear_timestep (BlockVector <double> & solution_delta
                << std::endl
                << "Newton iteration: " << it_nr
                << std::endl;
-       
+
        tangent_matrix = 0.0;
        residual = 0.0;
-       
+
        // Check residual
        make_constraints (it_nr, constraints); // Make boundary conditions
        assemble_system_F (); // Assemble RHS
@@ -1313,7 +1308,7 @@ void Solid<dim>::solve_nonlinear_timestep (BlockVector <double> & solution_delta
        // Residual scaling factors
        res_f = soln_error_res.block(u_dof).l2_norm();
        if (it_nr == 0) res_f_0 = res_f;
-       
+
        // Check for solution convergence
        if (   it_nr > 0
                && res_u/res_u_0 <= parameters.tol_u
@@ -1328,19 +1323,19 @@ void Solid<dim>::solve_nonlinear_timestep (BlockVector <double> & solution_delta
                    << "\t Relative force error: " << res_f/res_f_0
                    << "\t Dilatation error: " << get_error_dil()
                    << std::endl << std::endl;
-           
+
            //      timer.leave_subsection();
            return;
        }
-       
+
        // No convergence -> continue with calculations
        // Assemble stiffness matrix
-       assemble_system_K (); 
-       
+       assemble_system_K ();
+
        // Do the static condensation to make K'_uu, and put K_pt^{-1}
        // in the K_pt block and K_tt^{-1} in the K_pp block
         assemble_SC();
-       
+
 
         // Do the static condensation to make K'_uu, and put K_pt^{-1}
         // in the K_pt block and K_tt^{-1} in the K_pp block
@@ -1349,11 +1344,11 @@ void Solid<dim>::solve_nonlinear_timestep (BlockVector <double> & solution_delta
        constraints.condense (tangent_matrix, residual); // Apply BC's
        solve_linear_system (newton_update);
        constraints.distribute(newton_update); // Populate the constrained DOF's with their values
-       
+
        // Newton update error
        get_error_update(newton_update, soln_error_update);
        res_u = soln_error_update.block(u_dof).l2_norm();
-       
+
        // Residual scaling factors
        if (it_nr == 0) res_u_0 = res_u;
        std::cout
@@ -1375,28 +1370,28 @@ void Solid<dim>::solve_nonlinear_timestep (BlockVector <double> & solution_delta
                  << "\t Relative force error: " << res_f/res_f_0
                  << "\t Dilatation error: " << get_error_dil()
                  << std::endl;
-       
+
        // Update and continue iterating
        solution_delta += newton_update; // Update current solution
        update_qph_incremental (solution_delta); // Update quadrature point information
     }
-    
+
     throw(ExcMessage("No convergence in nonlinear solver!"));
 }
 
 template <int dim>
 void Solid<dim>::get_error_res (const BlockVector <double> & residual, BlockVector <double> & error_res)
-{   
+{
     for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i)
-       if (!constraints.is_constrained(i)) 
+       if (!constraints.is_constrained(i))
            error_res(i) = residual(i);
 }
 
 template <int dim>
 void Solid<dim>::get_error_update (const BlockVector <double> & newton_update, BlockVector <double> & error_update)
-{   
+{
     for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i)
-       if (!constraints.is_constrained(i)) 
+       if (!constraints.is_constrained(i))
            error_update(i) = newton_update(i);
 }
 
@@ -1405,9 +1400,9 @@ 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();
@@ -1416,13 +1411,13 @@ double Solid<dim>::get_error_dil (void)
         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
 }
 
@@ -1432,7 +1427,7 @@ BlockVector <double> Solid<dim>::get_solution_total (const BlockVector <double>
 {
     BlockVector <double> solution_total (solution_n);
     solution_total += solution_delta;
-    
+
     return solution_total;
 }
 
@@ -1441,22 +1436,19 @@ template <int dim>
 void Solid<dim>::solve_linear_system (BlockVector <double> & newton_update)
 {
     std::cout << "Solve linear system..." << std::endl;
-    
+
     BlockVector <double> A (dofs_per_block);
     BlockVector <double> B (dofs_per_block);
     A.collect_sizes ();
     B.collect_sizes ();
-    
+
     //      | K'_uu |   K_up    |     0     |         | dU_u |         | dR_u |
     // K =  | K_pu  |   K_tt^-1 |   K_pt^-1 | , dU =  | dU_p | , dR =  | dR_p |
     //      |   0   |   K_tp    |   K_tt    |         | dU_t |         | dR_t |
-    
+
     // Solve for du
     {
-<<<<<<< .mine
 
-=======
->>>>>>> .r23597
        // K'uu du = Ru âˆ’ Kup Ktp^-1 (Rt âˆ’ Ktt Kpt^{-1} Rp)
        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));
@@ -1464,27 +1456,27 @@ void Solid<dim>::solve_linear_system (BlockVector <double> & newton_update)
        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");
        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 * 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);
-           
+
            // SSOR -> much better than Jacobi for symmetric systems
            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),
                             residual.block(u_dof),
                             preconditioner);
-           
+
            std::cout
                    << "\t Iterations: " << solver_control.last_step()
                    << "\n\t Residual: " << solver_control.last_value()
@@ -1500,7 +1492,7 @@ void Solid<dim>::solve_linear_system (BlockVector <double> & newton_update)
        else throw (ExcMessage("Linear solver type not implemented"));
        timer.leave_subsection();
     }
-    
+
     timer.enter_subsection("Linear solver postprocessing");
     // Postprocess for dp
     {
@@ -1512,7 +1504,7 @@ void Solid<dim>::solve_linear_system (BlockVector <double> & newton_update)
        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 (newton_update.block(p_dof), A.block(t_dof));
     }
-    
+
     // Postprocess for dt
     {
        // dt = Ktt^{-1} (Rt âˆ’ Ktp dp)
@@ -1529,14 +1521,14 @@ void Solid<dim>::assemble_system_K (void)
 {
     timer.enter_subsection("Assemble system matrix");
     std::cout << "Assemble system matrix..."<< std::endl;
-    
+
     tangent_matrix = 0.0; // Clear the matrix
-    
+
     const UpdateFlags uf_cell ( update_values | update_gradients | update_JxW_values  );
-    
+
     PerTaskData_K per_task_data (dofs_per_cell); // Initialise members of per_task_data to the correct sizes.
     ScratchData_K scratch_data (fe, qf_cell, uf_cell);
-    
+
     WorkStream::run (  dof_handler_ref.begin_active(),
                     dof_handler_ref.end(),
                     *this,
@@ -1544,7 +1536,7 @@ void Solid<dim>::assemble_system_K (void)
                     &Solid::copy_local_to_global_K,
                     scratch_data,
                     per_task_data);
-    
+
     timer.leave_subsection();
 }
 
@@ -1569,15 +1561,15 @@ void Solid<dim>::assemble_system_K_one_cell (const typename DoFHandler<dim>::act
     scratch.fe_values_ref.reinit (cell);
     cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have
     PointHistory<dim> *lqph = reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
-    
+
     // Set up cell shape function gradients
     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]);
@@ -1593,7 +1585,7 @@ void Solid<dim>::assemble_system_K_one_cell (const typename DoFHandler<dim>::act
            }
        }
     }
-    
+
     // Build cell stiffness matrix
     // Global and local system matrices are symmetric
     //  => Take advantage of this:  Build only the lower half of the local matrix
@@ -1603,23 +1595,23 @@ void Solid<dim>::assemble_system_K_one_cell (const typename DoFHandler<dim>::act
        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();
-       
+
        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;
            const unsigned int i_group = fe.system_to_base_index(i).first.first;
-           
+
            // Only assemble the lower diagonal part of the local matrix
            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;
-               
+
                if (   (i_group == j_group) && (i_group == u_dof ) ) {
                    data.cell_matrix(i,j)
                            += ( symm_B[i] * C * symm_B[j]   // Material stiffness
@@ -1640,9 +1632,9 @@ void Solid<dim>::assemble_system_K_one_cell (const typename DoFHandler<dim>::act
                else Assert ((i_group <= t_dof) && (j_group <= t_dof), ExcInternalError());
            } // END j LOOP
        } // END i LOOP
-       
+
     } // END q_point LOOP
-    
+
     // Global and local system matrices are symmetric
     // => Copy the upper half of the local matrix in the bottom  half of the local matrix
     for (unsigned int i=0; i<dofs_per_cell; ++i) {
@@ -1658,19 +1650,19 @@ void Solid<dim>::assemble_system_F (void)
 {
     timer.enter_subsection("Assemble system RHS");
     std::cout << "Assemble system RHS..."<< std::endl;
-    
+
     residual  = 0.0; // Clear the vector
-    
+
     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); // Initialise members of per_task_data to the correct sizes.
     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,
@@ -1678,7 +1670,7 @@ void Solid<dim>::assemble_system_F (void)
                     &Solid::copy_local_to_global_F,
                     scratch_data,
                     per_task_data );
-    
+
     timer.leave_subsection();
 }
 
@@ -1701,14 +1693,14 @@ void Solid<dim>::assemble_system_F_one_cell (const typename DoFHandler<dim>::act
     scratch.fe_values_ref.reinit (cell);
     cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have
     PointHistory<dim> *lqph = reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
-    
+
     // Precompute some 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);
            }
@@ -1721,7 +1713,7 @@ void Solid<dim>::assemble_system_F_one_cell (const typename DoFHandler<dim>::act
            else Assert (k_group <= t_dof, ExcInternalError());
        }
     }
-    
+
     // Assembly for residual contribution
     for (unsigned int q_point=0; q_point < n_q_points; ++q_point) {
        const SymmetricTensor <2,dim>  T = lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol();
@@ -1729,14 +1721,14 @@ void Solid<dim>::assemble_system_F_one_cell (const typename DoFHandler<dim>::act
        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();
-       
+
        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;
-           
+
            if (i_group == u_dof) {
                data.cell_rhs(i) -= ( symm_B[i]*T )*JxW; // R_u
            }
@@ -1749,39 +1741,39 @@ void Solid<dim>::assemble_system_F_one_cell (const typename DoFHandler<dim>::act
            else Assert (i_group <= t_dof, ExcInternalError());
        } // END i LOOP
     } // END q_point LOOP
-    
+
     // Assembly for Neumann RHS contribution
     if (cell->at_boundary() == true)
     {
        static const Tensor <2, dim> I = static_cast < Tensor <2, dim> > ( unit_symmetric_tensor <dim> () );
-       
+
        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)
                {
                    const Tensor <1, dim> & N = scratch.fe_face_values_ref.normal_vector(f_q_point);
-                   
+
                    // Traction in reference configuration
                    // t_0 = p*N
                    static const double p0 = -4.0/(parameters.scale*parameters.scale); // Reference pressure of 4 Pa
                    const double time_ramp = (time.current() / time.end()); // Linearly ramp up the pressure with time
                    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) {
                        // 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;
-                       
+
                        if (i_group == u_dof) {
                            const unsigned int component_i = fe.system_to_component_index(i).first;
                            const double & Ni = scratch.fe_face_values_ref.shape_value(i,f_q_point);
                            const double & JxW = scratch.fe_face_values_ref.JxW(f_q_point);
-                           
+
                            // Add traction vector contribution to the local RHS vector (displacement dofs only)
                            data.cell_rhs(i) += (Ni * traction[component_i])  // Contribution from external forces
                                    * JxW;
@@ -1789,7 +1781,7 @@ void Solid<dim>::assemble_system_F_one_cell (const typename DoFHandler<dim>::act
                    } // END i LOOP
                } // END face q_point LOOP
            } // END at boundary check LOOP
-           
+
        } // END face LOOP
     }
 }
@@ -1799,13 +1791,13 @@ template <int dim>
 void Solid<dim>::assemble_SC  (void)
 {
     timer.enter_subsection("Perform static condensation");
-    
+
     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,
@@ -1813,7 +1805,7 @@ void Solid<dim>::assemble_SC  (void)
                     &Solid::copy_local_to_global_SC,
                     scratch_data,
                     per_task_data  );
-    
+
     timer.leave_subsection();
 }
 
@@ -1836,7 +1828,7 @@ void Solid<dim>::assemble_SC_one_cell (const typename DoFHandler<dim>::active_ce
     data.reset();
     scratch.reset();
     cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have
-    
+
     // The local stifness matrix K_e is:
     //  | K_uu  |   K_up   |   0   |
     //  | K_pu  |     0    |  K_pt |
@@ -1847,7 +1839,7 @@ void Solid<dim>::assemble_SC_one_cell (const typename DoFHandler<dim>::active_ce
     //  | K_pu  |   K_tt^-1 |   K_pt^-1 |
     //  |   0   |   K_tp    |   K_tt    |
     // with K'_uu = K_uu + Kup Ktp^{-1} Ktt Kpt^{-1} Kpu
-    
+
     // NOTE:
     // GLOBAL Data already exists in the K_uu, K_pt, K_tp subblocks
     //
@@ -1874,9 +1866,9 @@ void Solid<dim>::assemble_SC_one_cell (const typename DoFHandler<dim>::active_ce
     //          replace it with.
     // K_tp^-1: Same as above
     // K_tt^-1: Nothing exists in the original K_pp subblock, so we can just add this contribution as is.
-    
+
     // Extract element data from the system matrix
-    
+
     AdditionalTools::extract_submatrix(data.local_dof_indices,
                                       data.local_dof_indices,
                                       tangent_matrix,
@@ -1893,7 +1885,7 @@ void Solid<dim>::assemble_SC_one_cell (const typename DoFHandler<dim>::active_ce
                                       element_indices_t,
                                       data.K_orig,
                                       data.K_tt);
-    
+
     // Place K_pt^-1 in the K_pt block
     data.K_pt_inv.invert(data.K_pt);
     data.K_pt_inv.add (-1.0, data.K_pt);
@@ -1901,14 +1893,14 @@ void Solid<dim>::assemble_SC_one_cell (const typename DoFHandler<dim>::active_ce
                                       element_indices_t,
                                       data.K_pt_inv,
                                       data.cell_matrix);
-    
+
     // Place K_tt^-1 in the K_pp block
     data.K_tt_inv.invert(data.K_tt);
     AdditionalTools::replace_submatrix(element_indices_p,
                                       element_indices_p,
                                       data.K_tt_inv,
                                       data.cell_matrix);
-    
+
     // Make condensation terms to add to the K_uu block
     data.K_pt_inv.mmult(data.A, data.K_pu);
     data.K_tt.mmult(data.B, data.A);
@@ -1926,25 +1918,25 @@ void Solid<dim>::make_constraints (const int & it_nr,
                                   ConstraintMatrix & constraints)
 {
     std::cout << "Make constraints..."<< std::endl;
-    
+
     constraints.clear();
     const bool apply_dirichlet_bc = (it_nr == 0);
-    
+
     // Boundary conditions:
     // b_id 0: -x face: Zero x-component of displacement : Symmetry plane
     // b_id 2: -y face: Zero y-component of displacement : Symmetry plane
     // b_id 4: -z face: Zero z-component of displacement : Symmetry plane
-    
+
     // b_id 5: +z face: Zero x-component and Zero y-component
     // b_id 6: Applied pressure face: Zero x-component and Zero y-component
     // b_id 1: +x face: Traction free
     // b_id 3: +y face: Traction free
     {
        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 );
        }
@@ -1954,10 +1946,10 @@ void Solid<dim>::make_constraints (const int & it_nr,
     }
     {
        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 );
        }
@@ -1969,7 +1961,7 @@ void Solid<dim>::make_constraints (const int & it_nr,
        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 );
        }
@@ -1981,7 +1973,7 @@ void Solid<dim>::make_constraints (const int & it_nr,
        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 );
        }
@@ -1993,7 +1985,7 @@ void Solid<dim>::make_constraints (const int & it_nr,
        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 );
        }
@@ -2001,7 +1993,7 @@ void Solid<dim>::make_constraints (const int & it_nr,
            VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction<dim>(n_components), constraints, components );
        }
     }
-    
+
     constraints.close();
 }
 
@@ -2010,15 +2002,15 @@ 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,
@@ -2030,12 +2022,12 @@ void Solid<dim>::output_results(void)
     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);
 }
@@ -2046,7 +2038,7 @@ int main ()
     try
     {
        deallog.depth_console (0);
-       
+
        Solid<3> solid_3d ("parameters.prm");
        solid_3d.run();
     }
@@ -2060,7 +2052,7 @@ int main ()
                  << "Aborting!" << std::endl
                  << "----------------------------------------------------"
                  << std::endl;
-       
+
        return 1;
     }
     catch (...)
@@ -2074,7 +2066,7 @@ int main ()
                  << 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.