From: pelteret Date: Sat, 24 Dec 2011 05:32:33 +0000 (+0000) Subject: Step-44: Made some corrections, added a convergence table and startedsome more in... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=852faab1bade520b6046f93609b36a8e10864cb7;p=dealii-svn.git Step-44: Made some corrections, added a convergence table and startedsome more in-depth documentation. git-svn-id: https://svn.dealii.org/trunk@24855 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-44/parameters.prm b/deal.II/examples/step-44/parameters.prm index 17f5a04b4d..2932d03d22 100644 --- a/deal.II/examples/step-44/parameters.prm +++ b/deal.II/examples/step-44/parameters.prm @@ -2,10 +2,10 @@ # --------------------- subsection Finite element system # Displacement system polynomial order - set Polynomial degree = 1 + set Polynomial degree = 2 # Gauss quadrature order - set Quadrature order = 2 + set Quadrature order = 3 end diff --git a/deal.II/examples/step-44/step-44.cc b/deal.II/examples/step-44/step-44.cc index 426a571184..15b7c82cdc 100644 --- a/deal.II/examples/step-44/step-44.cc +++ b/deal.II/examples/step-44/step-44.cc @@ -2,7 +2,7 @@ /* Authors: Jean-Paul Pelteret, University of Cape Town, */ /* Andrew McBride, University of Erlangen-Nuremberg, 2010 */ /* */ -/* Copyright (C) 2010, 2011 by the deal.II authors */ +/* Copyright (C) 2010 by the deal.II authors */ /* & Jean-Paul Pelteret and Andrew McBride */ /* */ /* This file is subject to QPL and may not be distributed */ @@ -10,310 +10,334 @@ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include +// @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 +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + #include #include -#include - -namespace Step44 -{ - using namespace dealii; +// Next we import all the deal.II +// function and class names to the global namespace +using namespace dealii; // @sect3{Run-time parameters} - namespace Parameters - { -// Finite Element system - struct FESystem - { - int poly_degree; - int quad_order; +// +// 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 +{ +// @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); - }; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); +}; - void FESystem::declare_parameters (ParameterHandler &prm) +void FESystem::declare_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Finite element system"); { - 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(); + 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) +void FESystem::parse_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Finite element system"); { - prm.enter_subsection("Finite element system"); - { - poly_degree = prm.get_integer("Polynomial degree"); - quad_order = prm.get_integer("Quadrature order"); - } - prm.leave_subsection(); + poly_degree = prm.get_integer("Polynomial degree"); + quad_order = prm.get_integer("Quadrature order"); } + prm.leave_subsection(); +} -// Geometry - struct Geometry - { - int global_refinement; - double scale; - double p_p0; +// @sect4{Geometry} +// Make adjustments to the problem geometry and the applied load. +// 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; - 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) +void Geometry::declare_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Geometry"); { - 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(); + 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) +void Geometry::parse_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Geometry"); { - 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(); + 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(); +} -// Materials - struct Materials - { - double nu; - double mu; +// @sect{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); - }; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); +}; - void Materials::declare_parameters (ParameterHandler &prm) +void Materials::declare_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Material properties"); { - 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(); + 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) +void Materials::parse_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Material properties"); { - prm.enter_subsection("Material properties"); - { - nu = prm.get_double("Poisson's ratio"); - mu = prm.get_double("Shear modulus"); - } - prm.leave_subsection(); + nu = prm.get_double("Poisson's ratio"); + mu = prm.get_double("Shear modulus"); } + prm.leave_subsection(); +} -// Linear solver - struct LinearSolver - { - std::string type_lin; - double tol_lin; - double max_iterations_lin; - double ssor_relaxation; +// @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); - }; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); +}; - void LinearSolver::declare_parameters (ParameterHandler &prm) +void LinearSolver::declare_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Linear solver"); { - 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(); + 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) +void LinearSolver::parse_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Linear solver"); { - 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(); + 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(); +} // Nonlinear solver - struct NonlinearSolver - { - unsigned int max_iterations_NR; - double tol_f; - double tol_u; +// Define the tolerances and maximum number of iterations for the +// Newton-Raphson nono-linear 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); - }; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); +}; - void NonlinearSolver::declare_parameters (ParameterHandler &prm) +void NonlinearSolver::declare_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Nonlinear solver"); { - 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(); + 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) +void NonlinearSolver::parse_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Nonlinear solver"); { - 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(); + 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(); +} -// Time - struct Time - { - double end_time; - double delta_t; +// @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); - }; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); +}; - void Time::declare_parameters (ParameterHandler &prm) +void Time::declare_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Time"); { - 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(); + 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) +void Time::parse_parameters (ParameterHandler &prm) +{ + prm.enter_subsection("Time"); { - prm.enter_subsection("Time"); - { - end_time = prm.get_double("End time"); - delta_t = prm.get_double("Time step size"); - } - prm.leave_subsection(); + end_time = prm.get_double("End time"); + delta_t = prm.get_double("Time step size"); } + prm.leave_subsection(); +} -// All parameters - struct AllParameters - : +// 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, @@ -321,912 +345,1153 @@ namespace Step44 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} - namespace AdditionalTools - { - template - 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 ) - { +// 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 \f$ \mathbf{A} \f$ and +// \f$ \mathbf{B} \f$ such that their outer-product +// \f$ \mathbf{A} \bar{\otimes} \mathbf{B} \Rightarrow C_{ijkl} = A_{ik} B_{jl} \f$ +template +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 +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) +{ - 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())); + const unsigned int n_rows_submatrix = row_index_set.size(); + const unsigned int n_cols_submatrix = column_index_set.size(); - 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.reinit(n_rows_submatrix, n_cols_submatrix); - 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()); + + 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 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(), 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())); +template <> +void extract_submatrix < dealii::BlockSparseMatrix >(const std::vector< unsigned int > &row_index_set, + const std::vector< unsigned int > &column_index_set, + const dealii::BlockSparseMatrix &matrix, + FullMatrix< double > &sub_matrix) +{ - 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())); + const unsigned int n_rows_submatrix = row_index_set.size(); + const unsigned int n_cols_submatrix = column_index_set.size(); - 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.reinit(n_rows_submatrix, n_cols_submatrix); - matrix(row, col) = sub_matrix(sub_row, sub_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()); - } - } - } + 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; - } + sub_matrix(sub_row,sub_col) = matrix(row, col); + } + } +} -// @sect3{Time class} - 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;} - - private: - unsigned int timestep; - double time_current; - const double time_end; - const double delta_t; - }; +// 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 +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()); -// @sect3{Neo-Hookean material} - template - class Material_NH - { - public: - /// \brief Class constructor - Material_NH (const double & lambda, - const double & mu) - : - lambda_0 (lambda), - mu_0 (mu), - 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) - { - const double dW_dJ = get_dU_dtheta (J); - return mu_0*B + dW_dJ*J*I; - } + 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()); - virtual 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 ); - } + 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()); - // 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));} + matrix(row, col) = sub_matrix(sub_row, sub_col); - 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; - }; +} - template SymmetricTensor<2, dim> const Material_NH::I = SymmetricTensor<2, dim> (unit_symmetric_tensor ()); - template SymmetricTensor<4, dim> const Material_NH::IxI = SymmetricTensor<4, dim> (outer_product (I, I)); - template SymmetricTensor<4, dim> const Material_NH::II = SymmetricTensor<4, dim> (identity_tensor ()); +// @sect3{Time class} +// A simple class to store time data is created. Its +// functioning is transparent so no discussion is +// necessary. +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;} + +private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; +}; -// @sect3{Quadrature point history} - template - class PointHistory - { - public: - PointHistory (void) - : - material (NULL), - dilatation_n (1.0), - 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 (lambda, - parameters.mu); +// @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 +// can be described by a strain-energy function (SEF) +// \f$ \phi = \phi_{B} + \phi_{V} \f$ +// where the bulk deformation is given by +// \f$ \phi_{B} = C_{1} \left( I_{1} - 3 \right) \f$ +// where \f$ C_{1} - \frac{\mu}{2} \f$ and $I_{1}$ is the first +// invariant of the left- or right- Cauchy deformation tensors. +// In this example the SEF that governs the volumetric +// response is defined as +// \f$ \phi_{V} = \kappa \left( \frac{1}{2} \left( \theta^{2} - 1 \right) - ln \left( \theta \right) \right) \f$ +// where $\kappa$ is the bulk modulus. +template +class Material_NH +{ +public: + /// \brief Class constructor + 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 + // \f$ \mathbf{T} = \mathbf{F} \frac{\partial \hat{\phi}}{\partial \mathbf{C}} \mathbf{F}^{T} = \frac{\partial \phi}{\partial \mathbf{B}} \mathbf{B} \f$ + 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; + } - // Initialise all tensors correctly - update_values (Tensor <2,dim> (), 0.0, 1.0); - } + // The tangent matrix for this material is also calculated from the SEF by + // \f$ JC_{ijkl} = F_{iA} F_{jB} H_{ABCD} F_{kC} F_{lD}\f$ + // with + // \f$ \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 ); + } - // Total Variables - void update_values (const Tensor<2, dim> & grad_u_n, - const double & pressure, - const double & dilatation) - { - // Calculated variables from displacement, displacement gradients - const Tensor <2,dim> F = static_cast > (unit_symmetric_tensor ()) + grad_u_n; - 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) - T_bar = material->get_T (get_J(), get_B_bar()); - T_iso = dev_P*get_T_bar(); // Note: T_iso depends on T_bar - T_vol = get_pressure()*get_J()*I; - } + // 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));} - // 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;} +protected: + // Material properties + const double lambda_0; // Lame constant + const double mu_0; // Shear modulus + const double kappa_0; // Bulk modulus - // Volumetric terms - double get_dU_dtheta (void) { - return material->get_dU_dtheta(get_dilatation()); - } + // 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; +}; - double get_d2U_dtheta2 (void) { - return material->get_d2U_dtheta2(get_dilatation()); - } +template SymmetricTensor<2, dim> const Material_NH::I = SymmetricTensor<2, dim> (unit_symmetric_tensor ()); +template SymmetricTensor<4, dim> const Material_NH::IxI = SymmetricTensor<4, dim> (outer_product (I, I)); +template SymmetricTensor<4, dim> const Material_NH::II = SymmetricTensor<4, dim> (identity_tensor ()); - // 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;}; +// @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: +// \f$ \mathbf{F} = \hat{\mathbf{F}} \bar{\mathbf{F}} \f$ +// where the volumetric part is +// \f$ \hat{\mathbf{F}} = J^{\frac{1}{3}} \mathbf{I} \f$ +// and the isochoric part is given by +// \f$ \bar{\mathbf{F}} = J^{-\frac{1}{3}} \mathbf{F} \f$ +// . From this, the deviatoric left Cauchy-Green deformation +// tensor can be defined as +// \f$ \bar{\mathbf{B}} = \bar{\mathbf{F}} \bar{\mathbf{F}}^{T} = J^{-\frac{2}{3}} \mathbf{F} \mathbf{F}^{T} \f$ +// +// Here we also introduce an additive volumetric-deviatoric split +// in the material reponse. We can express the governing SEF as +// \f$ \phi = \phi_{V} + \phi_{I} \f$ +// with the result that the Kirchhoff stress is additively +// decomposed into +// \f$ \mathbf{\tau} = \mathbf{\tau}_{V} + \mathbf{\tau}_{I} \f$ +// as is the tangent matrix +// \f$ J\mathbf{C} = J\mathbf{C}_{V} + J\mathbf{C}_{I} \f$. +// +// These quantities are calculated as +// \f$ \mathbf{\tau}_{I} = pJ\mathbf{I} \f$ +// \f$ \mathbf{\tau}_{V} = \mathcal{P}:\bar{\mathbf{\tau}} \f$ +// with \f$ \bar{\mathbf{\tau}} = \mathbf{\tau} \vert_{\mathbf{B} = \bar{\mathbf{B}}} \f$ +// and the deviatoric tensor \f$ \mathcal{P} = \mathcal{I} - \mathbf{I} \otimes \mathbf{I} \f$ +// \f$ J\mathbf{C}_{I} = pJ(\mathbf{I} \otimes \mathbf{I} - 2 \mathcal{I}) \f$ +// \f$ 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} \f$ +// with \f$ \bar{\mathcal{C}} = \mathcal{C} \vert_{\mathbf{B} = \bar{\mathbf{B}}} \f$ +template +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 (lambda, + parameters.mu); + + // Initialise all tensors correctly + update_values (Tensor <2,dim> (), + 0.0, + 1.0); + } - // 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(); + // 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 > (unit_symmetric_tensor ()); + 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 + // Stress update can only update the stresses once the + // dilatation has been set as p = p(d). + // 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; + } - 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); + // 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(); - 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; - } + 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); - SymmetricTensor <4,dim> get_C_vol(void) - { - const double & p = get_pressure(); - const double & J = get_J(); - return p*J*(IxI - 2.0*II); - } + 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); + } - private: - // === MATERIAL === - Material_NH * 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; - static SymmetricTensor<4, dim> const II; - static SymmetricTensor<4, dim> const dev_P; - }; - - template SymmetricTensor<2,dim> const PointHistory::I - = SymmetricTensor<2,dim> (unit_symmetric_tensor ()); - template SymmetricTensor<4,dim> const PointHistory::IxI - = SymmetricTensor<4,dim> (outer_product (I, I)); - template SymmetricTensor<4,dim> const PointHistory::II - = SymmetricTensor<4,dim> (identity_tensor ()); - template SymmetricTensor<4,dim> const PointHistory::dev_P - = SymmetricTensor<4,dim> (II - 1.0/3.0*IxI); +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 * 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; +}; + +template SymmetricTensor<2,dim> const PointHistory::I += SymmetricTensor<2,dim> (unit_symmetric_tensor ()); +template SymmetricTensor<4,dim> const PointHistory::IxI += SymmetricTensor<4,dim> (outer_product (I, I)); +template SymmetricTensor<4,dim> const PointHistory::II += SymmetricTensor<4,dim> (identity_tensor ()); +template SymmetricTensor<4,dim> const PointHistory::dev_P += SymmetricTensor<4,dim> (II - 1.0/3.0*IxI); // @sect3{Quasi-static quasi-incompressible finite-strain solid} - template - class Solid - { - public: - Solid (const std::string & input_file); - virtual ~Solid (void); - void run (void); - - private: - - // === DATA STRUCTS === - - struct PerTaskData_K - { - FullMatrix cell_matrix; - std::vector 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 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 & fe_cell, - const QGauss & 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; - } - } - } - - }; - - struct PerTaskData_F - { - Vector cell_rhs; - std::vector 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 fe_values_ref; - FEFaceValues fe_face_values_ref; - - std::vector < std::vector< double > > Nx; - std::vector < std::vector< SymmetricTensor<2, dim> > > symm_grad_Nx; - std::vector< Vector > rhs_values; - - // Solution data - std::vector< std::vector > > solution_grads; - - ScratchData_F ( const FiniteElement & fe_cell, - const QGauss & qf_cell, - const UpdateFlags uf_cell, - const QGauss & 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)), - rhs_values (qf_cell.size(), - Vector(dim)) - { } - - 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), - 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) { - Nx[q_point][k] = 0.0; - symm_grad_Nx[q_point][k] = 0.0; - rhs_values[q_point] = 0.0; - } - } - } - - }; - - struct PerTaskData_SC - { - FullMatrix cell_matrix; - std::vector local_dof_indices; - - // Calculation matrices (auto resized) - FullMatrix K_orig; - FullMatrix K_pu; - FullMatrix K_pt; - FullMatrix K_tt; - // Calculation matrices (manual resized) - FullMatrix K_pt_inv; - FullMatrix K_tt_inv; - FullMatrix K_con; - FullMatrix A; - FullMatrix B; - FullMatrix 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 - // The matrix extraction and replacement tools will take care of this - void reset(void) { } - }; - - // Dummy struct for TBB - struct ScratchData_SC - { - ScratchData_SC (void) { } - 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 fe_values_ref; - std::vector< Tensor< 2, dim> > solution_grads_u_total; - std::vector solution_values_p_total; - std::vector solution_values_t_total; - const BlockVector & solution_total; - - ScratchData_UQPH (const FiniteElement & fe_cell, - const QGauss & qf_cell, - const UpdateFlags uf_cell, - const BlockVector & solution_total) - : - fe_values_ref (fe_cell, - qf_cell, - uf_cell), - solution_grads_u_total (qf_cell.size()), - solution_values_p_total (qf_cell.size()), - 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(), - rhs.fe_values_ref.get_quadrature(), - rhs.fe_values_ref.get_update_flags()), - 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), - solution_total (rhs.solution_total) - { } - - 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; - } - } - }; - - // === 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::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::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::active_cell_iterator & cell, - ScratchData_SC & scratch, - PerTaskData_SC & data); - void copy_local_to_global_SC (const PerTaskData_SC & data); - /// \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 - void update_qph_incremental ( const BlockVector & solution_delta ); - void update_qph_incremental_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_UQPH & scratch, - PerTaskData_UQPH & data); - void copy_local_to_global_UQPH (const PerTaskData_UQPH & data) {} - /// \brief Solve for the displacement using a Newton-Rhapson method - void solve_nonlinear_timestep (BlockVector & solution_delta); - void solve_linear_system (BlockVector & newton_update); - - /// \brief Error measurement - void get_error_res (const BlockVector & residual, BlockVector & error_res); - void get_error_update (const BlockVector & newton_update, BlockVector & error_update); - double get_error_dil (void); - - // Solution - BlockVector get_solution_total (const BlockVector & solution_delta); - - // Postprocessing - void output_results(void); - - // === ATTRIBUTES === - // Parameters - Parameters::AllParameters parameters; - - // Geometry - Triangulation triangulation; // Describes the triangulation - - // Time - Time time; - TimerOutput timer; - - // === Quadrature points === - std::vector< PointHistory > quadrature_point_history; // Quadrature point history - - // === Finite element system === - DoFHandler dof_handler_ref; // Describes the degrees of freedom - const unsigned int degree; - const FESystem 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 dofs_per_block; - std::vector element_indices_u; - std::vector element_indices_p; - std::vector element_indices_t; - - // === Quadrature === - QGauss qf_cell; // Cell quadrature formula - QGauss 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 - BlockSparseMatrix tangent_matrix; // Global stiffness matrix - BlockVector residual; // Holds the residual vector - BlockVector solution_n; // Holds the solution vector: Total displacement over all time-steps - }; +template +class Solid +{ +public: + 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::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::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::active_cell_iterator & cell, + ScratchData_SC & scratch, + PerTaskData_SC & data); + void copy_local_to_global_SC (const PerTaskData_SC & data); + /// \brief 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 & solution_delta ); + void update_qph_incremental_one_cell (const typename DoFHandler::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 & solution_delta); + std::pair solve_linear_system (BlockVector & newton_update); + + // Solution retrieval + BlockVector get_solution_total (const BlockVector & 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 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 > 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 fe; + DoFHandler 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 dofs_per_block; + std::vector element_indices_u; + std::vector element_indices_p; + std::vector 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 qf_cell; + QGauss 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 tangent_matrix; + BlockVector residual; + BlockVector 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 & 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 Solid class} // @sect4{Public interface} - template - Solid::Solid (const std::string & input_file) - : - parameters (input_file), - triangulation (Triangulation::maximum_smoothing), - time (parameters.end_time, parameters.delta_t), - timer (std::cout, - TimerOutput::summary, - TimerOutput::wall_times), - dof_handler_ref (triangulation), - degree (parameters.poly_degree), - fe (FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree-1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree-1), 1), // dilatation - 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) - { +// We initialise the the solid class using data extracted +// from the parameter file. +template +Solid::Solid (const std::string & input_file) + : + parameters (input_file), + triangulation (Triangulation::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(parameters.poly_degree), dim, + FE_DGPMonomial(parameters.poly_degree-1), 1, + FE_DGPMonomial(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 - Solid::~Solid (void) - { +// The class destructor simply needs to clear the data held by the DOFHandler +template +Solid::~Solid (void) +{ dof_handler_ref.clear (); - } +} - template - void Solid::run (void) - { - // Pre-processing +// In solving the quasti-static problem, the time +// becomes a loading parameter. We choose to increment +// time linearly using a constant timestep size. +template +void Solid::run (void) +{ + // After preprocessing, we output the initial grid + // before starting the simulation proper. make_grid (); system_setup (); - output_results (); // Output initial grid position + output_results (); time.increment(); BlockVector solution_delta (dofs_per_block); solution_delta.collect_sizes (); - while (time.current() <= time.end()) { - solution_delta = 0.0; + 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(); + } +} + +// @sect3{Private interface} + +// @sect4{Threaded-building-blocks structures} +// We choose to 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. + +// Firstly we deal with the tangent matrix assembly structures. +// The PerTaskData object stores local contributions. +template +struct Solid::PerTaskData_K +{ + FullMatrix cell_matrix; + std::vector 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 +struct Solid::ScratchData_K +{ + FEValues 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 & fe_cell, + const QGauss & 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. +// The PerTaskData object again stores local contributions +template +struct Solid::PerTaskData_F +{ + Vector cell_rhs; + std::vector 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; } +}; +// and the ScratchData object the shape function object +// and precomputed values vector +template +struct Solid::ScratchData_F +{ + FEValues fe_values_ref; + FEFaceValues 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 > > solution_grads; + + ScratchData_F ( const FiniteElement & fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell, + const QGauss & 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; + } + } + } - // Solve step and update total solution vector - solve_nonlinear_timestep (solution_delta); - solution_n += solution_delta; +}; - output_results (); - time.increment(); +// 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. +// We place these in the PerTaskData struct. +template +struct Solid::PerTaskData_SC +{ + FullMatrix cell_matrix; + std::vector local_dof_indices; + + // Calculation matrices (auto resized) + FullMatrix K_orig; + FullMatrix K_pu; + FullMatrix K_pt; + FullMatrix K_tt; + // Calculation matrices (manual resized) + FullMatrix K_pt_inv; + FullMatrix K_tt_inv; + FullMatrix K_con; + FullMatrix A; + FullMatrix B; + FullMatrix 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) { } +}; +// 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 +struct Solid::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 +struct Solid::PerTaskData_UQPH +{ + PerTaskData_UQPH (void) { } + void reset(void) { } +}; +// The ScratchData object will be used to store a alias fort 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 +struct Solid::ScratchData_UQPH +{ + const BlockVector & solution_total; + + std::vector< Tensor< 2, dim> > solution_grads_u_total; + std::vector solution_values_p_total; + std::vector solution_values_t_total; + + FEValues fe_values_ref; + + ScratchData_UQPH (const FiniteElement & fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell, + const BlockVector & 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; + } } - } +}; // @sect4{Solid::make_grid} - template - void Solid::make_grid (void) - { +// Here we create the grid on which the minimisation problem is to be solved. +template +void Solid::make_grid (void) +{ + // Create a unit cube with each face given a boundary ID number GridGenerator::hyper_rectangle ( triangulation, Point (0.0, 0.0, 0.0), Point (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 + 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::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); + cell = triangulation.begin_active(), + endc = triangulation.end(); for (; cell!=endc; ++cell) - { + { if (cell->at_boundary() == true) { - for (unsigned int face=0; face < GeometryInfo::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 - } - } - } + for (unsigned int face=0; face < GeometryInfo::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} - template - void Solid::system_setup (void) - { +// Next we describe how the FE system is setup. +template +void Solid::system_setup (void) +{ timer.enter_subsection ("Setup system"); - // Number of components per block + // 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 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 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); - // Count number of dofs per block - DoFTools::count_dofs_per_block (dof_handler_ref, dofs_per_block, block_component); + 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; - - // 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 | - // | 0 | K_tp | K_tt | | dU_t | | dR_t | - // reflect this structure in the sparsity pattern - Table<2,DoFTools::Coupling> coupling (n_components, n_components); - for (unsigned int ii = 0; ii < n_components; ++ii) { - for (unsigned int jj = ii; jj < n_components; ++jj) { - if ((ii < p_component) && (jj == t_component)) { - coupling[jj][ii] = DoFTools::none; - coupling[ii][jj] = DoFTools::none; - } - else { - coupling[ii][jj] = DoFTools::always; - coupling[jj][ii] = DoFTools::always; - } - } - } + << "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 system matrix + // 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(); - - 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); + 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 + // 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; // Dilatation is 1 in the initial configuration + solution_n.block(t_dof) = 1.0; - // Set up the quadrature point history + // and finally set up the quadrature point history setup_qph (); timer.leave_subsection(); - } +} -// A way to extract subblocks from the matrix - template - void Solid::determine_component_extractors(void) - { +// 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 +void Solid::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; - 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()); - } + // 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()); + } } - } +} // @sect4{Solid::setup_qph} - template - void Solid::setup_qph (void) - { +// The method used to store quadrature information is already described in +// tutorial 18. Here we implement a similar setup for a SMP machine. +template +void Solid::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. { - typename Triangulation::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 > 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()); + quadrature_point_history = std::vector< PointHistory > (triangulation.n_active_cells() * n_q_points); + + unsigned int history_index = 0; + typename Triangulation::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()); } - // Setup initial data + // Next we setup the initial QP data typename DoFHandler::active_cell_iterator - cell = dof_handler_ref.begin_active(), - endc = dof_handler_ref.end(); + cell = dof_handler_ref.begin_active(), + endc = dof_handler_ref.end(); for (; cell != endc; ++cell) { - PointHistory* lqph = reinterpret_cast*> (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 ); - } + PointHistory* lqph = reinterpret_cast*> (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} - template - void Solid::update_qph_incremental (const BlockVector & solution_delta) - { +// 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 +void Solid::update_qph_incremental (const BlockVector & solution_delta) +{ timer.enter_subsection("Update QPH data"); - std::cout << "Update QPH data..."<< std::endl; + std::cout << " UQPH "<< std::flush; - // Get total solution as it stands at this update increment + // Firstly we need to attain the total solution as it stands + // at this Newton increment const BlockVector 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, @@ -1234,6 +1499,7 @@ namespace Step44 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, @@ -1243,13 +1509,15 @@ namespace Step44 per_task_data_UQPH); timer.leave_subsection(); - } +} - template - void Solid::update_qph_incremental_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_UQPH & scratch, - PerTaskData_UQPH & data) - { +// Now we describe how we extract data from the solution vector and pass it +// along to each QP storage object for processing. +template +void Solid::update_qph_incremental_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_UQPH & scratch, + PerTaskData_UQPH & data) +{ PointHistory* lqph = reinterpret_cast*> (cell->user_pointer()); Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); Assert(lqph < &quadrature_point_history.back(), ExcInternalError()); @@ -1258,823 +1526,1043 @@ namespace Step44 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 + // 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); - // === UPDATE DATA AT EACH GAUSS POINT === - // Update displacement and deformation gradient at all quadrature points + // 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]); + 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]); } - } +} // @sect4{Solid::solve_nonlinear_timestep} - template - void Solid::solve_nonlinear_timestep (BlockVector & solution_delta) - { - // timer.enter_subsection("Nonlinear solver"); +template +void Solid::solve_nonlinear_timestep (BlockVector & solution_delta) +{ + // timer.enter_subsection("Nonlinear solver"); std::cout - << "Timestep " << time.get_timestep() - << std::endl; + << std::endl + << "Timestep " << time.get_timestep() + << " @ " << time.current() << "s" + << std::endl; - // Newton update vector + // We create a new vector to store the current Newton update step BlockVector newton_update (dofs_per_block); newton_update.collect_sizes (); - // Solution error vectors - BlockVector soln_error_res (dofs_per_block); // Holds the true residual vector - BlockVector soln_error_update (dofs_per_block); // Holds the update error vector - soln_error_res.collect_sizes (); - soln_error_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(); - double res_u = 0.0, res_f = 0.0; - double res_u_0 = 1.0, res_f_0 = 1.0; + // 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::endl - << "Newton iteration: " << it_nr - << std::endl; - + << " " + << 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; - // Check residual - make_constraints (it_nr, constraints); // Make boundary conditions + // 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_res(residual, soln_error_res); - // Residual scaling factors - res_f = soln_error_res.block(u_dof).l2_norm(); - if (it_nr == 0) res_f_0 = res_f; + get_error_residual(error_residual); - // Check for solution convergence + // 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 - && res_u/res_u_0 <= parameters.tol_u - && res_f/res_f_0 <= parameters.tol_f) - { + && error_update_norm.u <= parameters.tol_u + && error_residual_norm.u <= parameters.tol_f) + { std::cout - << std::endl - << "Solution for timestep " << time.get_timestep() - << " converged on Newton iteration " << it_nr-1 << "." - << std::endl - << "Relative displacement error: " << res_u/res_u_0 - << "\t Relative force error: " << res_f/res_f_0 - << "\t Dilatation error: " << get_error_dil() - << std::endl << std::endl; - - // timer.leave_subsection(); - return; - } + << " CONVERGED! " + << std::endl; - // No convergence -> continue with calculations - // Assemble stiffness matrix - assemble_system_K (); + print_conv_footer(); - // 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(); + // timer.leave_subsection(); + return; + } - // 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(); + assemble_system_K (); // Assemble stiffness matrix + make_constraints (it_nr, constraints); // Make boundary conditions + constraints.condense (tangent_matrix, + residual); // Apply BC's - constraints.condense (tangent_matrix, residual); // Apply BC's - solve_linear_system (newton_update); + const std::pair lin_solver_output = 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(); + 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); - // Residual scaling factors - if (it_nr == 0) res_u_0 = res_u; std::cout - << "Nonlinear system error: " - << std::endl << std::scientific - << " Solution update \t ||dU||: " << soln_error_update.l2_norm() - << "\t ||dU_u||: " << soln_error_update.block(u_dof).l2_norm() - << "\t ||dU_p||: " << soln_error_update.block(p_dof).l2_norm() - << "\t ||dU_t||: " << soln_error_update.block(t_dof).l2_norm() - << std::endl; - std::cout << std::scientific - << " Residual \t ||dF||: " << soln_error_res.l2_norm() - << "\t ||dR_u||: " << soln_error_res.block(u_dof).l2_norm() - << "\t ||dR_p||: " << soln_error_res.block(p_dof).l2_norm() - << "\t ||dR_t||: " << soln_error_res.block(t_dof).l2_norm() - << std::endl; - std::cout << std::scientific - << " Relative displacement error: " << res_u/res_u_0 - << "\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 - } + << " | " + << 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 - void Solid::get_error_res (const BlockVector & residual, BlockVector & error_res) - { - for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i) - if (!constraints.is_constrained(i)) - error_res(i) = residual(i); - } +// We print out data in a nice table that is updated +// on a per-iteration basis. Here we set up the table +// header +template +void Solid::print_conv_header (void) +{ + static const unsigned int l_width = 155; - template - void Solid::get_error_update (const BlockVector & newton_update, BlockVector & error_update) - { - for (unsigned int i=0; i < dof_handler_ref.n_dofs(); ++i) - if (!constraints.is_constrained(i)) - error_update(i) = newton_update(i); - } + 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 +void Solid::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; - template - double Solid::get_error_dil (void) - { + + 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 +double Solid::get_error_dil (void) +{ double v_e = 0.0; // Volume in current configuration double V_e = 0.0; // Volume in reference configuration FEValues fe_values_ref (fe, qf_cell, update_JxW_values); typename DoFHandler::active_cell_iterator - cell = dof_handler_ref.begin_active(), - endc = dof_handler_ref.end(); + cell = dof_handler_ref.begin_active(), + endc = dof_handler_ref.end(); for (; cell != endc; ++cell) { - fe_values_ref.reinit (cell); - PointHistory* lqph = reinterpret_cast*> (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); - } + fe_values_ref.reinit (cell); + PointHistory* lqph = reinterpret_cast*> (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 - } - -// Solution (valid at any Newton step) - template - BlockVector Solid::get_solution_total (const BlockVector & solution_delta) - { - BlockVector solution_total (solution_n); - solution_total += solution_delta; - - return solution_total; - } - -// @sect4{Solid::solve_linear_system} - template - void Solid::solve_linear_system (BlockVector & newton_update) - { - std::cout << "Solve linear system..." << std::endl; - - BlockVector A (dofs_per_block); - BlockVector 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 - { - - // 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)); - 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"); - 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(); +// Determine the true residual error for the problem +template +void Solid::get_error_residual (Errors & error_residual) +{ + BlockVector error_res (dofs_per_block); + error_res.collect_sizes (); - SolverControl solver_control (solver_its , tol_sol); + // 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); - GrowingVectorMemory < Vector > GVM; - SolverCG < Vector > solver_CG (solver_control, GVM); + 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(); +} - // SSOR -> much better than Jacobi for symmetric systems - PreconditionSSOR > preconditioner; - preconditioner.initialize (tangent_matrix.block(u_dof, u_dof), parameters.ssor_relaxation); +// Determine the true Newton update error for the problem +template +void Solid::get_error_update (const BlockVector & newton_update, + Errors & error_update) +{ + BlockVector error_ud (dofs_per_block); + error_ud.collect_sizes (); - solver_CG.solve (tangent_matrix.block(u_dof, u_dof), - newton_update.block(u_dof), - residual.block(u_dof), - preconditioner); + // 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); - std::cout - << "\t Iterations: " << solver_control.last_step() - << "\n\t Residual: " << solver_control.last_value() - << std::endl; - } - else if (parameters.type_lin == "Direct") - { - 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)); - } - else throw (ExcMessage("Linear solver type not implemented")); - timer.leave_subsection(); - } + 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(); +} - timer.enter_subsection("Linear solver postprocessing"); - // Postprocess for dp - { - // dp = Ktp^{-1} ( Rt − Ktt Kpt^{-1} (Rp − Kpu du) ) - tangent_matrix.block(p_dof, u_dof).vmult (A.block(p_dof), newton_update.block(u_dof)); - B.block(p_dof).equ(1.0, residual.block(p_dof), -1.0, A.block(p_dof)); - tangent_matrix.block(p_dof, t_dof).vmult(A.block(t_dof), B.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 (newton_update.block(p_dof), A.block(t_dof)); - } +// 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 +BlockVector Solid::get_solution_total (const BlockVector & solution_delta) +{ + BlockVector solution_total (solution_n); + solution_total += solution_delta; - // Postprocess for dt - { - // dt = Ktt^{-1} (Rt − Ktp dp) - tangent_matrix.block(t_dof, p_dof).vmult (A.block(t_dof), newton_update.block(p_dof)); - residual.block(t_dof) -= A.block(t_dof); - tangent_matrix.block(p_dof, p_dof).vmult (newton_update.block(t_dof), residual.block(t_dof)); - } - timer.leave_subsection(); - } + return solution_total; +} // @sect4{Solid::assemble_system_K} - template - void Solid::assemble_system_K (void) - { - timer.enter_subsection("Assemble system matrix"); - std::cout << "Assemble system matrix..."<< std::endl; +// Since we use TBB for assembly, we simply setup a copy of the +// data structures required for the process and pass them, along +// with the memory addresses of the assembly functions to the +// WorkStream object for processing. Note that we must ensure that +// the matrix is reset before any assembly operations can occur. +template +void Solid::assemble_system_K (void) +{ + timer.enter_subsection("Assemble tangent matrix"); + std::cout << " ASM_K " << std::flush; - tangent_matrix = 0.0; // Clear the matrix + 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); // Initialise members of per_task_data to the correct sizes. - 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); + 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(); - } +} - template - void Solid::copy_local_to_global_K (const PerTaskData_K & data) - { - // Add the local contribution to the system matrix +// 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 +void Solid::copy_local_to_global_K (const PerTaskData_K & data) +{ for (unsigned int i=0; i - void Solid::assemble_system_K_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_K & scratch, - PerTaskData_K & data) - { - data.reset(); // Reset data in the PerTaskData_K storage unit - scratch.reset(); // Reset data in the Scratch storage unit + for (unsigned int j=0; j +void Solid::assemble_system_K_one_cell (const typename DoFHandler::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); // Find out which global numbers the degrees of freedom on this cell have + cell->get_dof_indices (data.local_dof_indices); PointHistory *lqph = reinterpret_cast*>(cell->user_pointer()); - // Set up cell shape function gradients + // 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 (); 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(); + 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; + 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()); + 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()); + } } - } } - // Build cell stiffness matrix - // Global and local system matrices are symmetric - // => Take advantage of this: Build only the lower half of the local matrix - // Only assemble 1/2 of the K_uu, K_pp = 0, K_tt blocks and the whole K_pt, K_ut, K_up blocks + // 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) { - 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(); - - const std::vector & 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 - + ( component_i == component_j ? - B[i][component_i] * T * B[j][component_j] : - 0.0 ) // Geometric stiffness. Only add this along local diagonals - ) * JxW; // K_uu - } - else if ( (i_group == p_dof) && (j_group == u_dof) ) { - data.cell_matrix(i,j) += N[i]*J*(symm_B[j]*I)*JxW; // K_pu - } - else if ( (i_group == t_dof) && (j_group == p_dof) ) { - data.cell_matrix(i,j) -= N[i]*N[j]*JxW; // K_tp - } - else if ( (i_group == j_group) && (i_group == t_dof) ) { - data.cell_matrix(i,j) += N[i]*C_v*N[j]*JxW; // K_tt - } - 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 + // 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 & 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; + } + // 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; + } + // 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; + } + 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 - void Solid::assemble_system_F (void) - { - timer.enter_subsection("Assemble system RHS"); - std::cout << "Assemble system RHS..."<< std::endl; +// The setup of the residual assembly 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 +void Solid::assemble_system_F (void) +{ + timer.enter_subsection("Assemble residual"); + std::cout << " ASM_R "<< std::flush; - residual = 0.0; // Clear the vector + 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); + 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. + PerTaskData_F per_task_data (dofs_per_cell); ScratchData_F scratch_data (fe, - qf_cell, - uf_cell, - qf_face, - uf_face); + 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 ); + 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 - void Solid::copy_local_to_global_F (const PerTaskData_F & data) - { - // Add the local contribution to the system RHS vector +template +void Solid::copy_local_to_global_F (const PerTaskData_F & data) +{ for (unsigned int i=0; i - void Solid::assemble_system_F_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_F & scratch, - PerTaskData_F & data) - { - data.reset(); // Reset data in the PerTaskData_K storage unit - scratch.reset(); // Reset data in the ScratchData_F storage unit +} + +template +void Solid::assemble_system_F_one_cell (const typename DoFHandler::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); // Find out which global numbers the degrees of freedom on this cell have + cell->get_dof_indices (data.local_dof_indices); PointHistory *lqph = reinterpret_cast*>(cell->user_pointer()); - // Precompute some data + // 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(); + const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv(); - for (unsigned int k=0; k 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; iat_boundary() == true) { + for (unsigned int face=0; face < GeometryInfo::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; + + 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; + } + } + } + } } - 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(); - 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(); +// @sect4{Solid::make_constraints} +// The constraints for this problem are simple to describe. +// However, since we are dealing with an iterative Newton method, +// it should be noted that any displacement constraints should only +// be specified at the zeroth iteration and subsequently no +// additional contributions are to be made since the constraints +// are already exactly satisfied. So we describe this process for +// 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 +void Solid::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; - 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); + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + components ); + } + } + { + const int boundary_id = 2; - for (unsigned int i=0; i components (n_components, false); + components[1] = true; - if (i_group == u_dof) { - data.cell_rhs(i) -= ( symm_B[i]*T )*JxW; // R_u + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + components ); } - else if (i_group == p_dof ) { - data.cell_rhs(i) -= N[i]*(J - D)*JxW; // R_p + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + components ); } - else if ( i_group == t_dof) { - data.cell_rhs(i) -= N[i]*(p_star-p)*JxW; // R_t + } + { + 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(n_components), + constraints, + components ); } - 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 () ); - - for (unsigned int face=0; face < GeometryInfo::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); + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, + boundary_id, + ZeroFunction(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(n_components), + constraints, + components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, + boundary_id, + ZeroFunction(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(n_components), + constraints, + components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + components ); + } + } - // 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; + constraints.close(); +} - 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; - } - } // END i LOOP - } // END face q_point LOOP - } // END at boundary check LOOP - - } // END face LOOP - } - } +// @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. +// 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 +std::pair Solid::solve_linear_system (BlockVector & newton_update) +{ + // Need to create two temporary vectors so that the static condensation operation can be performed + BlockVector A (dofs_per_block); + BlockVector 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") + { + 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 > GVM; + SolverCG < Vector > 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 > 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); + + 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), + residual.block(u_dof)); + + lin_it = 1; + lin_res = 0.0; + } + else throw (ExcMessage("Linear solver type not implemented")); + 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); +} // @sect4{Solid::assemble_system_SC} - template - void Solid::assemble_SC (void) - { +// 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 +// 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 +void Solid::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. + 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 ); + 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 - void Solid::copy_local_to_global_SC (const PerTaskData_SC & data) - { - // Add the local contribution to the system matrix +// We need to describe how to add the local contribution to the tangent matrix. +template +void Solid::copy_local_to_global_SC (const PerTaskData_SC & data) +{ for (unsigned int i=0; i - void Solid::assemble_SC_one_cell (const typename DoFHandler::active_cell_iterator & cell, - ScratchData_SC & scratch, - PerTaskData_SC & data) - { + for (unsigned int j=0; j +void Solid::assemble_SC_one_cell (const typename DoFHandler::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); // 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 | - // | 0 | K_tp | K_tt | - // - // We are going to exploit the zeros for post-processing as: - // | K'_uu | K_up | 0 | - // | 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 - // - // 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 and - // Additionally the copy_local_to_global operation is a "+=" operation -> need to take this - // into account - // - // 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 | K_tt^-1 | K_pt^-1 - K_pt | - // | 0 | 0 | 0 | - // - // Strategy 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: 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 - + 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); + 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); + element_indices_u, + data.K_orig, + data.K_pu); AdditionalTools::extract_submatrix(element_indices_p, - element_indices_t, - data.K_orig, - data.K_pt); + element_indices_t, + data.K_orig, + data.K_pt); AdditionalTools::extract_submatrix(element_indices_t, - element_indices_t, - data.K_orig, - data.K_tt); + element_indices_t, + data.K_orig, + data.K_tt); - // Place K_pt^-1 in the K_pt block + // 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); - 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); - // 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 + // 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); // Symmetric matrix - data.K_pu.Tmmult(data.K_con, data.C); // Symmetric matrix + 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); - } - -// @sect4{Solid::make_constraints} - template - void Solid::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; + element_indices_u, + data.K_con, + data.cell_matrix); - 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(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(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(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(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(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(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(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(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(n_components), constraints, components ); - } - else { - VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); - } - } - - constraints.close(); - } + // 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); +} // @sect4{Solid::output_results} - template - void Solid::output_results(void) - { +// Here we present how the results are written to file to be viewed +// using Paraview. The method is similar to that shown in previous +// tutorials so will not be discussed in detail. +template +void Solid::output_results(void) +{ DataOut data_out; - - std::vector data_component_interpretation (dim, DataComponentInterpretation::component_is_part_of_vector); + std::vector 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 solution_name (dim, "displacement"); + std::vector 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::type_dof_data, data_component_interpretation); - // MappingQEulerian q_mapping (degree, solution_n.block(u_dof), dof_handler_ref); - // MappingQEulerian q_mapping (degree, solution_n, dof_handler_ref); - Vector soln; - soln.reinit(solution_n.size()); - for (unsigned int i=0; i < soln.size(); ++i) soln(i) = solution_n(i); - MappingQEulerian q_mapping (degree, soln, dof_handler_ref); - data_out.build_patches (q_mapping,degree); + DataOut::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 soln (solution_n.size()); + for (unsigned int i=0; i < soln.size(); ++i) + soln(i) = solution_n(i); + MappingQEulerian q_mapping (degree, + soln, + dof_handler_ref); + data_out.build_patches (q_mapping, + degree); std::ostringstream filename; filename << "solution-" - << time.get_timestep() - << ".vtk"; + << time.get_timestep() + << ".vtk"; std::ofstream output (filename.str().c_str()); data_out.write_vtk (output); - } } - // @sect3{Main function} -int main () +// Lastly we provide the main driver function which appears +// no different to the other tutorials. +int main (void) { - try + try { - using namespace dealii; - using namespace Step44; - - deallog.depth_console (0); + deallog.depth_console (0); - Solid<3> solid_3d ("parameters.prm"); - solid_3d.run(); + Solid<3> solid_3d ("parameters.prm"); + solid_3d.run(); } - catch (std::exception &exc) + 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; + 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 (...) + catch (...) { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; } - return 0; + return 0; } -