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