From: bangerth Date: Mon, 20 Feb 2012 09:13:13 +0000 (+0000) Subject: Much reorganization but no functional changes. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=75b816734662e9d19ea48ba661d9e6abd8fca458;p=dealii-svn.git Much reorganization but no functional changes. git-svn-id: https://svn.dealii.org/trunk@25117 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-44/step-44.cc b/deal.II/examples/step-44/step-44.cc index 506d58121c..4a880e65ca 100644 --- a/deal.II/examples/step-44/step-44.cc +++ b/deal.II/examples/step-44/step-44.cc @@ -1,7 +1,7 @@ /* Authors: Jean-Paul Pelteret, University of Cape Town, */ /* Andrew McBride, University of Erlangen-Nuremberg, 2010 */ /* */ -/* Copyright (C) 2010 by the deal.II authors */ +/* Copyright (C) 2010, 2011, 2012 by the deal.II authors */ /* & Jean-Paul Pelteret and Andrew McBride */ /* */ /* This file is subject to QPL and may not be distributed */ @@ -9,11 +9,11 @@ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ -// 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. + // 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 @@ -54,16 +54,22 @@ #include #include -// Next we import all the deal.II -// function and class names to the global namespace -using namespace dealii; + + // We then stick everything that relates to + // this tutorial program into a namespace of + // its own, and import all the deal.II + // function and class names into it: +namespace Step44 +{ + using namespace dealii; // @sect3{Run-time parameters} // // There are several parameters that can be set // in the code so we set up a ParameterHandler // object to read in the choices at run-time. -namespace Parameters { + namespace Parameters + { // @sect4{Finite Element system} // As mentioned in the introduction, a different order // interpolation should be used for the displacement @@ -77,289 +83,331 @@ namespace Parameters { // 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; + struct FESystem + { + unsigned int poly_degree; + unsigned 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", "2", Patterns::Integer(0), - "Displacement system polynomial order"); - - prm.declare_entry("Quadrature order", "3", Patterns::Integer(0), - "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::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + prm.declare_entry("Polynomial degree", "2", + Patterns::Integer(0), + "Displacement system polynomial order"); + + prm.declare_entry("Quadrature order", "3", + Patterns::Integer(0), + "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(); + } // @sect4{Geometry} // Make adjustments to the problem geometry and the applied load. // Since the problem modelled here is quite specific, the load // scale can be altered to specific values to attain results given // in the literature. -struct Geometry { - int global_refinement; + struct Geometry + { + int global_refinement; double scale; double p_p0; 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(0), - "Global refinement level"); - - prm.declare_entry("Grid scale", "1e-3", Patterns::Double(0.0), - "Global grid scaling factor"); - - prm.declare_entry("Pressure ratio p/p0", "100", - 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::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + prm.declare_entry("Global refinement", "2", + Patterns::Integer(0), + "Global refinement level"); + + prm.declare_entry("Grid scale", "1e-3", + Patterns::Double(0.0), + "Global grid scaling factor"); + + prm.declare_entry("Pressure ratio p/p0", "100", + 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(); + } // @sect4{Materials} -// Need the shear modulus $ \mu $ +// We also need the shear modulus $ \mu $ // and Poisson ration $ \nu $ // for the neo-Hookean material. -struct Materials { + 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.4999", Patterns::Double(-1.0,0.5), - "Poisson's ratio"); - - prm.declare_entry("Shear modulus", "80.194e6", 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::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Material properties"); + { + prm.declare_entry("Poisson's ratio", "0.4999", + Patterns::Double(-1.0,0.5), + "Poisson's ratio"); + + prm.declare_entry("Shear modulus", "80.194e6", + 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(); + } // @sect4{Linear solver} -// Choose both solver and preconditioner settings. +// Next, choose both solver and preconditioner settings. // The use of an effective preconditioner is critical to ensure // convergence when a large nonlinear motion occurs // in a Newton increment. // ToDo: explain // The default values are optimal for single-thread conditions this particular problem. -struct LinearSolver { + struct LinearSolver + { std::string type_lin; - double tol_lin; - double max_iterations_lin; + double tol_lin; + double max_iterations_lin; std::string preconditioner_type; - double preconditioner_relaxation; + double preconditioner_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(0.0), - "Linear solver residual (scaled by residual norm)"); - - prm.declare_entry( - "Max iteration multiplier", - "1", - Patterns::Double(0.0), - "Linear solver iterations (multiples of the system matrix size)"); - - prm.declare_entry("Preconditioner type", "ssor", Patterns::Selection("jacobi|ssor"), - "Type of preconditioner"); - - prm.declare_entry("Preconditioner relaxation", "0.65", Patterns::Double(0.0), - "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"); - preconditioner_type = prm.get("Preconditioner type"); - preconditioner_relaxation = prm.get_double("Preconditioner relaxation"); - } - 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(0.0), + "Linear solver residual (scaled by residual norm)"); + + prm.declare_entry("Max iteration multiplier", "1", + Patterns::Double(0.0), + "Linear solver iterations (multiples of the system matrix size)"); + + prm.declare_entry("Preconditioner type", "ssor", + Patterns::Selection("jacobi|ssor"), + "Type of preconditioner"); + + prm.declare_entry("Preconditioner relaxation", "0.65", + Patterns::Double(0.0), + "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"); + preconditioner_type = prm.get("Preconditioner type"); + preconditioner_relaxation = prm.get_double("Preconditioner relaxation"); + } + prm.leave_subsection(); + } // @sect4{Nonlinear solver} // 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 { + struct NonlinearSolver + { unsigned int max_iterations_NR; - double tol_f; - double tol_u; + 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(0), - "Number of Newton-Raphson iterations allowed"); - - prm.declare_entry("Tolerance force", "1.0e-9", Patterns::Double(0.0), - "Force residual tolerance"); - - prm.declare_entry("Tolerance displacement", "1.0e-6", - Patterns::Double(0.0), "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::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Nonlinear solver"); + { + prm.declare_entry("Max iterations Newton-Raphson", "10", + Patterns::Integer(0), + "Number of Newton-Raphson iterations allowed"); + + prm.declare_entry("Tolerance force", "1.0e-9", + Patterns::Double(0.0), + "Force residual tolerance"); + + prm.declare_entry("Tolerance displacement", "1.0e-6", + Patterns::Double(0.0), + "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(); + } // @sect4{Time} // Set the timestep size $ \varDelta t $ // and the simulation end-time. -struct 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::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::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(); + } // @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); 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); -} - -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); -} -} + }; + + 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::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); + } + } // @sect3{General tools} // We need to perform some specific operations that are not defined @@ -367,29 +415,8 @@ void AllParameters::parse_parameters(ParameterHandler &prm) { // 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]; - } - } - } - } - - return A_ik_B_jl; -} + namespace AdditionalTools + { // The extract_submatrix function // takes specific entries from a matrix, @@ -398,65 +425,72 @@ SymmetricTensor<4, dim> outer_product_T23(const SymmetricTensor<2, dim> & A, // 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()); +// to size $ r \times c $. At the beginning we +// check the size of the input vectors + 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(); + 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()); + 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); - } + sub_matrix(sub_row, sub_col) = matrix(row, col); + } } -} + } // 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) { + 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(); - // check the size of the input vectors - Assert(n_rows_submatrix > 0, ExcInternalError()); - Assert(n_cols_submatrix > 0, ExcInternalError()); + Assert(n_rows_submatrix > 0, ExcInternalError()); + Assert(n_cols_submatrix > 0, ExcInternalError()); - sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix); + 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_row = 0; sub_row < n_rows_submatrix; ++sub_row) + { + const unsigned int row = row_index_set[sub_row]; + Assert(row<=matrix.m(), ExcInternalError()); - for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) { - const unsigned int col = column_index_set[sub_col]; - Assert(col<=matrix.n(), ExcInternalError()); - if (matrix.get_sparsity_pattern().exists(row, col) == false) - continue; + for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) + { + const unsigned int col = column_index_set[sub_col]; + Assert(col<=matrix.n(), ExcInternalError()); + if (matrix.get_sparsity_pattern().exists(row, col) == false) + continue; - sub_matrix(sub_row, sub_col) = matrix(row, col); - } + sub_matrix(sub_row, sub_col) = matrix(row, col); + } } -} + } // The replace_submatrix function takes // specific entries from a sub_matrix, @@ -465,670 +499,912 @@ void extract_submatrix >( // 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()); + 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()); + 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); + 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 / dim) * IxI); -} + template + class StandardTensors + { + public: + + // $\mathbf{I}$ + static const SymmetricTensor<2, dim> I; + // $\mathbf{I} \otimes \mathbf{I}$ + static const SymmetricTensor<4, dim> 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 const SymmetricTensor<4, dim> II; + // Fourth-order deviatoric such that + // $\textrm{dev}(\bullet) = (\bullet) - (1/3)[(\bullet):\mathbf{I}]\mathbf{I}$ + static const SymmetricTensor<4, dim> dev_P; + }; + + template + const SymmetricTensor<2, dim> + StandardTensors::I = unit_symmetric_tensor(); + + template + const SymmetricTensor<4, dim> + StandardTensors::IxI = outer_product(I, I); + + template + const SymmetricTensor<4, dim> + StandardTensors::II = identity_tensor(); + + template + const SymmetricTensor<4, dim> + StandardTensors::dev_P = (II - (1.0 / dim) * IxI); + } // @sect3{Time class} // A simple class to store time data. Its // functioning is transparent so no discussion is // 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) { - } - - double current(void) const { - return time_current; + 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() + {} + + double current() const + { + return time_current; } - double end(void) const { - return time_end; + double end() const + { + return time_end; } - double get_delta_t(void) const { - return delta_t; + double get_delta_t() const + { + return delta_t; } - unsigned int get_timestep(void) const { - return timestep; + unsigned int get_timestep() const + { + return timestep; } - void increment(void) { - time_current += delta_t; - ++timestep; + void increment() + { + time_current += delta_t; + ++timestep; } -private: - unsigned int timestep; - double time_current; - const double time_end; - const double delta_t; -}; + private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; + }; // @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 within a three-field formulation. -// Compressible neo-Hookean materials -// can be described by a strain-energy function (SEF) -// $ \Psi = \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(\widetilde{J}) $. +// 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 within a three-field formulation. Compressible neo-Hookean +// materials can be described by a strain-energy function (SEF) $ \Psi = +// \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(\widetilde{J}) +// $. // -// The isochoric response is given by -// $ \Psi_{\text{iso}}(\overline{\mathbf{b}}) = c_{1} [\overline{I}_{1} - 3] $ +// The isochoric response is given by $ +// \Psi_{\text{iso}}(\overline{\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 -// $ \Psi_{\text{vol}}(\widetilde{J}) = \kappa \frac{1}{4} [ \widetilde{J}^2 - 1 - 2\textrm{ln}\widetilde{J} ]$. -// where $\kappa:= \lambda + 2/3 \mu$ is the bulk modulus and -// $\lambda$ is Lame's first parameter. -template -class Material_Compressible_Neo_Hook_Three_Field { -public: - Material_Compressible_Neo_Hook_Three_Field(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), - p_tilde(0.0), - J_tilde(1.0), - b_bar(AdditionalTools::StandardTensors::I) { - Assert(kappa > 0, ExcInternalError()); - } - - ~Material_Compressible_Neo_Hook_Three_Field(void) { - } +// 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 $ +// \Psi_{\text{vol}}(\widetilde{J}) = \kappa \frac{1}{4} [ \widetilde{J}^2 - 1 +// - 2\textrm{ln}\; \widetilde{J} ]$. where $\kappa:= \lambda + 2/3 \mu$ is +// the bulk modulus +// and $\lambda$ is Lame's first +// parameter. +// +// The following class will be used to characterize the material we work with, +// and provides a central point that one would need to modify if one were to +// implement a different material model. For it to work, we will store one +// object of this type per quadrature point, and in each of these objects +// store the current state (characterized by the values of the three fields) +// so that we can compute the elastic coefficients linearized around the +// current state. + template + class Material_Compressible_Neo_Hook_Three_Field + { + public: + Material_Compressible_Neo_Hook_Three_Field(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), + p_tilde(0.0), + J_tilde(1.0), + b_bar(AdditionalTools::StandardTensors::I) + { + Assert(kappa > 0, ExcInternalError()); + } + + ~Material_Compressible_Neo_Hook_Three_Field() + {} + + // 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$ and at the end of the + // function include a safety check for + // internal consistency: + void update_material_data(const Tensor<2, dim> & F, + const double p_tilde_in, + const double J_tilde_in) + { + det_F = determinant(F); + b_bar = std::pow(det_F, -2.0 / 3.0) * symmetrize(F * transpose(F)); + p_tilde = p_tilde_in; + J_tilde = J_tilde_in; - // 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 p_tilde_in, - const double J_tilde_in - ) { - det_F = determinant(F); - b_bar = std::pow(det_F, -2.0 / 3.0) * symmetrize(F * transpose(F)); - p_tilde = p_tilde_in; - J_tilde = J_tilde_in; - - // include a coupled of checks on the input data - Assert(det_F > 0, ExcInternalError()); + Assert(det_F > 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 second function determines the + // Kirchhoff stress $\boldsymbol{\tau} + // = \boldsymbol{\tau}_{\textrm{iso}} + + // \boldsymbol{\tau}_{\textrm{vol}}$ + SymmetricTensor<2, dim> get_tau() + { + 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() const + { + return get_Jc_vol() + get_Jc_iso(); } - // 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() const + { + return (kappa / 2.0) * (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. We calculate + // $\frac{\partial^2 + // \Psi_{\textrm{vol}}(\widetilde{J})}{\partial + // \widetilde{J} \partial + // \widetilde{J}}$ + double get_d2Psi_vol_dJ2() const + { + return ( (kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde))); } - // 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 / 2.0) * (J_tilde - 1.0 / J_tilde); + // The next few functions return + // various data that we choose to store + // with the material: + double get_det_F() const + { + return det_F; } - // Second derivative of the volumetric free energy wrt $\widetilde{J}$ - // We need the following computation explicitly in the tangent - // so we make it public. - // We calculate - // $\frac{\partial^2 \Psi_{\textrm{vol}}(\widetilde{J})}{\partial \widetilde{J} \partial \widetilde{J}}$ - double get_d2Psi_vol_dJ2(void) const { - return ( (kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde))); + double get_p_tilde() const + { + return p_tilde; } - // Return various data that we choose to store with the material - double get_det_F(void) const { - return det_F; + double get_J_tilde() const + { + return J_tilde; + } + + protected: + // Define constitutive model paramaters + // $\kappa$ and the neo-Hookean model + // parameter $c_1$: + const double kappa; + const double c_1; + + // Model specific data that is + // convenient to store with the + // material: + double det_F; + double p_tilde; + double J_tilde; + SymmetricTensor<2, dim> b_bar; + + // The following functions are used + // internally in determining the result + // of some of the public functions + // above. The first one determines the + // volumetric Kirchhoff stress + // $\boldsymbol{\tau}_{\textrm{vol}}$: + SymmetricTensor<2, dim> get_tau_vol() const + { + return p_tilde * det_F * AdditionalTools::StandardTensors::I; } - double get_p_tilde(void) const { - return p_tilde; + // Next, determine the isochoric + // Kirchhoff stress + // $\boldsymbol{\tau}_{\textrm{iso}} = + // \mathcal{P}:\overline{\boldsymbol{\tau}}$: + SymmetricTensor<2, dim> get_tau_iso() const + { + return AdditionalTools::StandardTensors::dev_P * get_tau_bar(); } - double get_J_tilde(void) const { - return J_tilde; + // Then, tetermine the fictitious + // Kirchhoff stress + // $\overline{\boldsymbol{\tau}}$: + SymmetricTensor<2, dim> get_tau_bar() const + { + return 2.0 * c_1 * b_bar; } -protected: - // Define constitutive model paramaters $\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 p_tilde; - 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 { - return p_tilde * det_F * AdditionalTools::StandardTensors::I; - } + // Calculate the volumetric part of the + // tangent $J + // \mathfrak{c}_\textrm{vol}$: + SymmetricTensor<4, dim> get_Jc_vol() const + { - // 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(); + return p_tilde * det_F + * ( AdditionalTools::StandardTensors::IxI + - (2.0 * AdditionalTools::StandardTensors::II) ); } - // Determine the fictitious Kirchhoff stress $\overline{\boldsymbol{\tau}}$ - SymmetricTensor<2, dim> get_tau_bar(void) const { - return 2.0 * c_1 * b_bar; + // Calculate the isochoric part of the + // tangent $J + // \mathfrak{c}_\textrm{iso}$: + SymmetricTensor<4, dim> get_Jc_iso() 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; + } + + // Calculate the fictitious elasticity + // tensor $\overline{\mathfrak{c}}$. + // For the material model chosen this + // is simply zero: + SymmetricTensor<4, dim> get_c_bar() const + { + return SymmetricTensor<4, dim>(); } + }; - // Calculate the volumetric part of the tangent $J \mathfrak{c}_\textrm{vol}$ - SymmetricTensor<4, dim> get_Jc_vol(void) const { - - return p_tilde * det_F - * ( AdditionalTools::StandardTensors::IxI - - (2.0 * AdditionalTools::StandardTensors::II) ); - } +// @sect3{Quadrature point history} - // 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; +// As seen in step-18, the +// PointHistory class offers a method for storing data at the +// quadrature points. Here each quadrature point holds a pointer to a +// material description. Thus, different material models can be used in +// different regions of the domain. Among other data, we choose to store the +// Kirchhoff stress $\boldsymbol{\tau}$ and the tangent $J\mathfrak{c}$ for +// the quadrature points. + template + class PointHistory + { + public: + PointHistory() + : + material(NULL), + F_inv(AdditionalTools::StandardTensors::I), + tau(SymmetricTensor<2, dim>()), + d2Psi_vol_dJ2(0.0), + dPsi_vol_dJ(0.0), + Jc(SymmetricTensor<4, dim>()) + {} + + virtual ~PointHistory() + { + delete material; + material = NULL; } - // Calculate the fictitious elasticity tensor $\overline{\mathfrak{c}}$. - // For the material model chosen this is simply zero. - SymmetricTensor<4, dim> get_c_bar() const { - SymmetricTensor<4, dim> c_bar; - c_bar = 0.0; - return c_bar; + // The first function is used to create + // a material object and to initialize + // all tensors correctly: + void setup_lqp (const Parameters::AllParameters & parameters) + { + material = new Material_Compressible_Neo_Hook_Three_Field(parameters.mu, + parameters.nu); + update_values(Tensor<2, dim>(), 0.0, 1.0); + } + + // The second one updates the stored + // values and stresses based on the + // current deformation measure + // $\textrm{Grad}\mathbf{u}_{\textrm{n}}$, + // pressure $\widetilde{p}$ and + // dilation $\widetilde{J}$ field + // values. + // + // To this end, we calculate the + // deformation gradient $\mathbf{F}$ + // from the displacement gradient + // $\textrm{Grad}\ \mathbf{u}$, i.e. + // $\mathbf{F}(\mathbf{u}) = \mathbf{I} + // + \textrm{Grad}\ \mathbf{u}$ and + // then let the material model + // associated with this quadrature + // point update itself. + void update_values (const Tensor<2, dim> & Grad_u_n, + const double p_tilde, + const double J_tilde) + { + const Tensor<2, dim> F = AdditionalTools::StandardTensors::I + + Grad_u_n; + material->update_material_data(F, p_tilde, J_tilde); + + // 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(); + + // Finally, we store the inverse of + // the deformation gradient since + // we frequently use it: + F_inv = invert(F); + } + + // We offer an interface to retrieve + // certain data. Here are the + // kinematic variables: + double get_J_tilde() const + { + return material->get_J_tilde(); } -}; -// @sect3{Quadrature point history} -// As seen in step-18, the PointHistory class offers -// a method for storing data at the quadrature points. -// Here each quadrature point holds a pointer to a Material. -// Thus, different material models can be used in different regions -// of the domain. -// Among other data, we choose to store the -// Kirchhoff stress $\boldsymbol{\tau}$ and -// the tangent $J\mathfrak{c}$ for the quadrature points. - -template -class PointHistory { -public: - PointHistory(void) : - material(NULL), - F_inv(AdditionalTools::StandardTensors::I), - tau(SymmetricTensor<2, dim>()), - d2Psi_vol_dJ2(0.0), - dPsi_vol_dJ(0.0), - Jc(SymmetricTensor<4, dim>()) { - } - virtual ~PointHistory(void) { - delete material; - material = NULL; + double get_det_F() const + { + return material->get_det_F(); } - // We first create a material object. - void setup_lqp(Parameters::AllParameters & parameters) { - - // Create an instance of a three field - // compressible neo-Hookean material - material = new Material_Compressible_Neo_Hook_Three_Field( - parameters.mu, parameters.nu); - - // Initialise all tensors correctly - update_values(Tensor<2, dim>(), 0.0, 1.0); + Tensor<2, dim> get_F_inv() const + { + return F_inv; } - // Update the stored values and stresses based on the current - // deformation measure $\textrm{Grad}\mathbf{u}_{\textrm{n}}$, - // pressure $\widetilde{p}$ and - // dilation $\widetilde{J}$ field values. - void update_values(const Tensor<2, dim> & Grad_u_n, - const double p_tilde, - const double J_tilde) { - - // Calculate the 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); - - // Now we update the material model with the new deformation measure, - // pressure and dilatation - material->update_material_data(F, p_tilde, J_tilde); - - // 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(); - + // ...and the kinetic variables. These + // are used in the material and global + // tangent matrix and residual assembly + // operations: + double get_p_tilde() const + { + return material->get_p_tilde(); } - // We offer an interface to retrieve certain data. - // Here are the kinematic variables - double get_J_tilde(void) const { - return material->get_J_tilde(); - } - double get_det_F(void) const { - return material->get_det_F(); - } - Tensor<2, dim> get_F_inv(void) const { - return F_inv; + SymmetricTensor<2, dim> get_tau() const + { + return tau; } - // and the kinetic variables. - // These are used in the material and global - // tangent matrix and residual assembly operations. - double get_p_tilde(void) const { - return material->get_p_tilde(); - } - SymmetricTensor<2, dim> get_tau(void) const { - return tau; + double get_dPsi_vol_dJ() const + { + return dPsi_vol_dJ; } - double get_dPsi_vol_dJ(void) const { - return dPsi_vol_dJ; + double get_d2Psi_vol_dJ2() const + { + return d2Psi_vol_dJ2; } - double get_d2Psi_vol_dJ2(void) const { - return d2Psi_vol_dJ2; + // and finally the tangent + SymmetricTensor<4, dim> get_Jc() const + { + return Jc; } - // and finally the tangent - SymmetricTensor<4, dim> get_Jc(void) const { - return Jc; - } + // In terms of member functions, this + // class stores for the quadrature + // point it represents a copy of a + // material type in case different + // materials are used in different + // regions of the domain, as well as + // the inverse of the deformation + // gradient... + private: + Material_Compressible_Neo_Hook_Three_Field* material; -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. - Material_Compressible_Neo_Hook_Three_Field* material; + Tensor<2, dim> F_inv; - // The inverse of the deformation gradient - Tensor<2, dim> F_inv; + // ... and stress-type variables along + // with the tangent $J\mathfrak{c}$: + SymmetricTensor<2, dim> tau; + double d2Psi_vol_dJ2; + double dPsi_vol_dJ; - // the stress-type variables - SymmetricTensor<2, dim> tau; - double d2Psi_vol_dJ2; - double dPsi_vol_dJ; + SymmetricTensor<4, dim> Jc; + }; - // and the tangent - SymmetricTensor<4, dim> Jc; -}; // @sect3{Quasi-static quasi-incompressible finite-strain solid} -// The Solid class is the central class in that it represents -// the problem at hand. -template -class Solid { -public: - Solid(const std::string & input_file); - virtual - ~Solid(void); - void - run(void); - -private: - - // Threaded building-blocks data structures - // for the tangent matrix. - // (see the module on @ref distributed - // for a definition, as well as the discussion in step-40) - 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_tangent(void); - void - assemble_system_tangent_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 conditions on the displacement field - void - make_constraints(const int & it_nr, ConstraintMatrix & constraints); - - // Create and update the quadrature points - 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-Raphson 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) const; - - // Post-processing and writing data to file - void - output_results(void) const; - - // A collection of the parameters used to describe the problem setup - Parameters::AllParameters parameters; - - // The volume of the reference and current configurations - double vol_reference; - double vol_current; - - // 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. - // See step-18 for more on this - 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, - // $\widetilde{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. - // We make use of a sparsity pattern designed for a block system. - 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); - std::pair - get_error_dil(void); - - // Print information to screen - void - print_conv_header(void); - void - print_conv_footer(void); -}; +// The Solid class is the central class in that it represents the problem at +// hand. It follows the usual scheme in that all it really has is a +// constructor, destructor and a run() function that dispatches +// all the work to private functions of this class: + template + class Solid + { + public: + Solid(const std::string & input_file); + + virtual + ~Solid(); + + void + run(); + + private: + + // In the private section of this + // class, we first forward declare a + // number of objects that are used in + // parallelizing work using the + // WorkStream object (see the @ref + // threads module for more information + // on this.) + // + // We declare such structures for the + // computation of tangent (stiffness) + // matrix, right hand side, static + // condensation, and for updating + // quadrature points: + struct PerTaskData_K; + struct ScratchData_K; + + struct PerTaskData_RHS; + struct ScratchData_RHS; + + struct PerTaskData_SC; + struct ScratchData_SC; + + struct PerTaskData_UQPH; + struct ScratchData_UQPH; + + // We start the collection of member + // functions with one that builds the + // grid: + void + make_grid(); + + // Set up the finite element system to + // be solved: + void + system_setup(); + + void + determine_component_extractors(); + + // Several functions to assemble the + // system and right hand side matrices + // using multi-threading. Each of them + // comes as a wrapper function, one + // that is executed to do the work in + // the WorkStream model on one cell, + // and one that copies the work done on + // this one cell into the global object + // that represents it: + void + assemble_system_tangent(); + + void + assemble_system_tangent_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 + 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 + 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 conditions on + // the displacement field + void + make_constraints(const int & it_nr, + ConstraintMatrix & constraints); + + // Create and update the quadrature + // points. Here, no data needs to be + // copied into a global object, so the + // copy_local_to_global function is + // empty: + void + setup_qph(); + + 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-Raphson method. We break this + // function into the nonlinear loop and + // the function that solves the + // linearized Newton-Raphson step: + void + solve_nonlinear_timestep(BlockVector & solution_delta); + + std::pair + solve_linear_system(BlockVector & newton_update); + + // Solution retrieval as well as + // post-processing and writing data to + // file: + BlockVector + get_solution_total(const BlockVector & solution_delta) const; + + void + output_results() const; + + // Finally, some member variables that + // describe the current state: A + // collection of the parameters used to + // describe the problem setup... + Parameters::AllParameters parameters; + + // ...the volume of the reference and + // current configurations... + double vol_reference; + double vol_current; + + // ...and description of the geometry on which + // the problem is solved: + Triangulation triangulation; + + // Also, keep track of the current time and the + // time spent evaluating certain + // functions + Time time; + TimerOutput timer; + + // A storage object for quadrature point + // information. See step-18 for more on + // this: + 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; + const 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, $\widetilde{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 = 1, + J_dof = 2 + }; + + 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. + const QGauss qf_cell; + const QGauss qf_face; + const unsigned int n_q_points; + const 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. We make + // use of a sparsity pattern designed for + // a block system. + 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() + : + norm(1.0), u(1.0), p(1.0), J(1.0) + {} + + void reset() + { + norm = 1.0; + u = 1.0; + p = 1.0; + J = 1.0; + } + void 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; + } + + double norm, u, p, J; + }; + + Errors 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); + + std::pair + get_error_dil(); + + // Print information to screen + static + void + print_conv_header(); + + void + print_conv_footer(); + }; // @sect3{Implementation of the Solid class} // @sect4{Public interface} // 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 - // 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(); -} + 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 \times + // DGP_{n-1} \times DGP_{n-1}$ + // system. $Q_2 \times DGP_1 + // \times DGP_1$ elements + // satisfy this condition, + // while $Q_1 \times DGP_0 + // \times DGP_0$ 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), + dofs_per_cell (fe.dofs_per_cell), + u_fe(first_u_component), + p_fe(p_component), + J_fe(J_component), + dofs_per_block(n_blocks), + qf_cell(parameters.quad_order), + qf_face(parameters.quad_order), + n_q_points (qf_cell.size()), + n_q_points_f (qf_face.size()) + { + determine_component_extractors(); + } // The class destructor simply clears the data held by the DOFHandler -template -Solid::~Solid(void) { - dof_handler_ref.clear(); -} + template + Solid::~Solid() + { + dof_handler_ref.clear(); + } + -// In solving the quasi-static problem, the time -// becomes a loading parameter. We choose to 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(); +// In solving the quasi-static problem, the time becomes a loading parameter, +// i.e. we increasing the loading linearly with time, making the two concepts +// interchangeable. We choose to increment time linearly using a constant time +// step size. +// +// We start the function with preprocessing, and then output the initial grid +// before starting the simulation proper with the first time (and loading) +// increment: + template + void Solid::run() + { + make_grid(); + system_setup(); + output_results(); + time.increment(); + + // Here we define the incremental solution + // update $\varDelta \mathbf{\Xi}:= + // \{\varDelta \mathbf{u},\varDelta + // \widetilde{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; + // and plot the results output_results(); + // we then move on happily to the next time step. time.increment(); - - // Here we define the incremental solution update - // $\varDelta \mathbf{\Xi}:= \{\varDelta \mathbf{u},\varDelta \widetilde{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; - // and plot the results - output_results(); - // we then move on happily to the next time step. - time.increment(); - } -} + } + } // @sect3{Private interface} @@ -1140,769 +1416,893 @@ void Solid::run(void) { // 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) { +// 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() + { + cell_matrix = 0.0; } + }; + - 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 // gradient and symmetric gradient vector which we will compute later. -template -struct Solid::ScratchData_K { - FEValues fe_values_ref; - - // interpolation function - std::vector > Nx; - // their gradients - std::vector > > grad_Nx; - // and their symmetric gradients. - 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) { - const unsigned int n_q_points = Nx.size(); - const unsigned int n_dofs_per_cell = Nx[0].size(); - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { - Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); - Assert( grad_Nx[q_point].size() == n_dofs_per_cell, - ExcInternalError()); - Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell, - ExcInternalError()); - for (unsigned int k = 0; k < n_dofs_per_cell; ++k) { - Nx[q_point][k] = 0.0; - grad_Nx[q_point][k] = 0.0; - symm_grad_Nx[q_point][k] = 0.0; - } + template + struct Solid::ScratchData_K + { + FEValues fe_values_ref; + + // interpolation function + std::vector > Nx; + // their gradients + std::vector > > grad_Nx; + // and their symmetric gradients. + 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() + { + const unsigned int n_q_points = Nx.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); + Assert( grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) + { + Nx[q_point][k] = 0.0; + grad_Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; } + } } -}; + }; // Next are the same approach is used for the // right-hand side assembly. // The PerTaskData object again stores local contributions -template -struct Solid::PerTaskData_RHS { - Vector cell_rhs; - std::vector local_dof_indices; - - 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; + template + struct Solid::PerTaskData_RHS + { + Vector cell_rhs; + std::vector local_dof_indices; + + PerTaskData_RHS(const unsigned int dofs_per_cell) + : + cell_rhs(dofs_per_cell), + local_dof_indices(dofs_per_cell) + {} + + void reset() + { + cell_rhs = 0.0; } -}; + }; // and the ScratchData object the shape function object // and precomputed values vector -template -struct Solid::ScratchData_RHS { - FEValues fe_values_ref; - FEFaceValues fe_face_values_ref; - - std::vector > Nx; - std::vector > > symm_grad_Nx; - - 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) { - const unsigned int n_q_points = Nx.size(); - const unsigned int n_dofs_per_cell = Nx[0].size(); - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { - Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); - Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell, - ExcInternalError()); - for (unsigned int k = 0; k < n_dofs_per_cell; ++k) { - Nx[q_point][k] = 0.0; - symm_grad_Nx[q_point][k] = 0.0; - } + template + struct Solid::ScratchData_RHS + { + FEValues fe_values_ref; + FEFaceValues fe_face_values_ref; + + std::vector > Nx; + std::vector > > symm_grad_Nx; + + 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() + { + const unsigned int n_q_points = Nx.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); + Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) + { + Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; } + } } -}; + }; // Here we define structures to assemble the statically -// condensed tangent matrix. Recall that we wish to solve -// for a displacement-based formulation. +// condensed tangent matrix. Recall that we wish to solve +// for a displacement-based formulation. // We do the condensation at the element // level as the $\widetilde{p}$ and $\widetilde{J}$ -// fields are element-wise discontinuous. -// As these operations are matrix-based, +// 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 +// 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; - - FullMatrix k_orig; - FullMatrix k_pu; - FullMatrix k_pJ; - FullMatrix k_JJ; - 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_orig(dofs_per_cell, dofs_per_cell), - k_pu(n_p, n_u), - k_pJ(n_p, n_J), - k_JJ(n_J, n_J), - k_pJ_inv(n_p, n_J), - k_bbar(n_u, n_u), - A(n_J,n_u), - B(n_J, n_u), - C(n_p, n_u) { - } + template + struct Solid::PerTaskData_SC + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + + FullMatrix k_orig; + FullMatrix k_pu; + FullMatrix k_pJ; + FullMatrix k_JJ; + 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_orig(dofs_per_cell, dofs_per_cell), + k_pu(n_p, n_u), + k_pJ(n_p, n_J), + k_JJ(n_J, n_J), + k_pJ_inv(n_p, n_J), + k_bbar(n_u, n_u), + A(n_J,n_u), + B(n_J, n_u), + C(n_p, n_u) + {} + + // We choose not to reset any data as the + // matrix extraction and replacement + // tools will take care of this + void reset() + {} + }; + + +// 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 create a dummy struct for this purpose. + template + struct Solid::ScratchData_SC + { + ScratchData_SC() + {} + + ScratchData_SC(const ScratchData_SC & rhs) {} + + void reset() + {} + }; - // We 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 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) { - } -}; -// 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 > 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()) { +// the PerTaskData object to store any information but must define one +// nonetheless. + template + struct Solid::PerTaskData_UQPH + { + PerTaskData_UQPH() + {} + + void reset() + {} + }; + + +// 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 > 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() + { + const unsigned int n_q_points = solution_grads_u_total.size(); + for (unsigned int q = 0; q < n_q_points; ++q) + { + solution_grads_u_total[q] = 0.0; + solution_values_p_total[q] = 0.0; + solution_values_J_total[q] = 0.0; + } } + }; - void reset(void) { - const unsigned int n_q_points = solution_grads_u_total.size(); - for (unsigned int q = 0; q < n_q_points; ++q) { - solution_grads_u_total[q] = 0.0; - solution_values_p_total[q] = 0.0; - solution_values_J_total[q] = 0.0; - } - } -}; // @sect4{Solid::make_grid} // 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); - - // determine the volume of the reference configuration - vol_reference = GridTools::volume(triangulation); - vol_current = vol_reference; - std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl; - - // 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 - } - } - } + template + void Solid::make_grid() + { + // 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); + + // determine the volume of the reference + // configuration + vol_reference = GridTools::volume(triangulation); + vol_current = vol_reference; + std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl; + + // 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[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(); + template + void Solid::system_setup() + { + 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) { - 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 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); + 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 { - Assert(k_group <= J_dof, ExcInternalError()); - } - } -} + 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 sub-blocks from the global matrix. + template + void + Solid::determine_component_extractors() + { + element_indices_u.clear(); + element_indices_p.clear(); + element_indices_J.clear(); + + for (unsigned int k = 0; k < fe.dofs_per_cell; ++k) + { + // 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 // 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. + template + void Solid::setup_qph() + { + 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. + { + triangulation.clear_user_data(); + { + std::vector > tmp; + tmp.swap(quadrature_point_history); + } + + quadrature_point_history.resize( + triangulation.n_active_cells() * n_q_points); + + unsigned int history_index = 0; + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); cell != triangulation.end(); + ++cell) { - triangulation.clear_user_data(); - { - std::vector > tmp; - tmp.swap(quadrature_point_history); - } - - quadrature_point_history.resize( - triangulation.n_active_cells() * n_q_points); - - unsigned int history_index = 0; - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); cell != triangulation.end(); - ++cell) { - cell->set_user_pointer(&quadrature_point_history[history_index]); - history_index += n_q_points; - } - - Assert(history_index == quadrature_point_history.size(), - ExcInternalError()); + cell->set_user_pointer(&quadrature_point_history[history_index]); + history_index += n_q_points; } - // Next we setup the initial QP data - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); cell != triangulation.end(); ++cell) { - PointHistory* lqph = - reinterpret_cast*>(cell->user_pointer()); + Assert(history_index == quadrature_point_history.size(), + ExcInternalError()); + } - Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); - Assert(lqph <= &quadrature_point_history.back(), ExcInternalError()); + // Next we setup the initial QP data + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); cell != triangulation.end(); ++cell) + { + PointHistory* lqph = + reinterpret_cast*>(cell->user_pointer()); - // Setup any initial information at Gauss points - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { - lqph[q_point].setup_lqp(parameters); - } - } -} + Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); + Assert(lqph <= &quadrature_point_history.back(), ExcInternalError()); + + // Setup any initial information at 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 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(); -} + 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_J_total.size() == n_q_points, - ExcInternalError()); - - scratch.reset(); - - // 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]); - } -} + 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()); + + scratch.reset(); + + // 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} // The driver method for the Newton-Raphson scheme -template -void Solid::solve_nonlinear_timestep( - BlockVector & solution_delta) { - 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(); - 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; - } - - // Now we assemble the tangent - assemble_system_tangent(); - // and make and impose the Dirichlet constraints - make_constraints(it_nr, constraints); - constraints.condense(tangent_matrix, system_rhs); - - // Now we actually solve the linearised problem - 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; - } + template + void + Solid::solve_nonlinear_timestep(BlockVector & solution_delta) + { + 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; - throw(ExcMessage("No convergence in nonlinear solver!")); -} + // 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(); + 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; + } + + // Now we assemble the tangent + assemble_system_tangent(); + // and make and impose the Dirichlet + // constraints + make_constraints(it_nr, constraints); + constraints.condense(tangent_matrix, system_rhs); + + // Now we actually solve the linearised + // problem + 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_J " << " NU_NORM " - << " NU_U " << " NU_P " << " NU_J " << std::endl; - - for (unsigned int i = 0; i < l_width; ++i) - std::cout << "_"; - std::cout << std::endl; -} -// and here the footer -template -void Solid::print_conv_footer(void) { - static const unsigned int l_width = 155; + template + void Solid::print_conv_header() + { + static const unsigned int l_width = 155; - for (unsigned int i = 0; i < l_width; ++i) - std::cout << "_"; - std::cout << std::endl; + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; - const std::pair error_dil = get_error_dil(); + std::cout << " SOLVER STEP " + << " | LIN_IT LIN_RES RES_NORM " + << " RES_U RES_P RES_J NU_NORM " + << " NU_U NU_P NU_J " << std::endl; - std::cout << "Relative errors:" << std::endl - << "Displacement:\t" << error_update.u / error_update_0.u << std::endl - << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl - << "Dilatation:\t" << error_dil.first << std::endl - << "v / V_0:\t" << vol_current << " / " << vol_reference << " = " << error_dil.second << std::endl; + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + } -} -// 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$. -// We also return the ratio of the current volume of the domain -// to the reference volume. This is of interest for incompressible media -// where we want to check how well the isochoric constraint has been +// and here the footer + template + void Solid::print_conv_footer() + { + static const unsigned int l_width = 155; + + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + + const std::pair error_dil = get_error_dil(); + + std::cout << "Relative errors:" << std::endl + << "Displacement:\t" << error_update.u / error_update_0.u << std::endl + << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl + << "Dilatation:\t" << error_dil.first << std::endl + << "v / V_0:\t" << vol_current << " / " << vol_reference + << " = " << error_dil.second << std::endl; + } + + +// 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$. We also return the ratio of the current volume of the +// domain to the reference volume. This is of interest for incompressible +// media where we want to check how well the isochoric constraint has been // enforced. -template -std::pair Solid::get_error_dil(void) { + template + std::pair Solid::get_error_dil() + { + double dil_L2_error = 0.0; + vol_current = 0.0; - double dil_L2_error = 0.0; - vol_current = 0.0; + FEValues fe_values_ref(fe, qf_cell, update_JxW_values); - FEValues fe_values_ref(fe, qf_cell, update_JxW_values); + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + { + fe_values_ref.reinit(cell); - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); cell != triangulation.end(); ++cell) { - fe_values_ref.reinit(cell); - - PointHistory* lqph = - reinterpret_cast*>(cell->user_pointer()); + PointHistory* lqph = + reinterpret_cast*>(cell->user_pointer()); - Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); - Assert(lqph <= &quadrature_point_history.back(), ExcInternalError()); + 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) { + 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); - 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_current += det_F_qp * JxW; + } + Assert(vol_current > 0, ExcInternalError()); + } - dil_L2_error += the_error_qp_squared * JxW; - vol_current += det_F_qp * JxW; - }Assert(vol_current > 0, ExcInternalError()); - } + std::pair error_dil; + error_dil.first = std::sqrt(dil_L2_error); + error_dil.second = vol_current / vol_reference; - std::pair error_dil; - error_dil.first = std::sqrt(dil_L2_error); - error_dil.second = vol_current / vol_reference; - return error_dil; -} + return error_dil; + } -// Determine the true residual error for the problem. +// 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(); -} + 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.J = error_ud.block(J_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) const { - BlockVector solution_total(solution_n); - solution_total += solution_delta; - return solution_total; - -} + template + BlockVector + Solid::get_solution_total(const BlockVector & solution_delta) const + { + BlockVector solution_total(solution_n); + solution_total += solution_delta; + return solution_total; + } // @sect4{Solid::assemble_system_tangent} // Since we use TBB for assembly, we simply setup a copy of the @@ -1910,157 +2310,198 @@ BlockVector Solid::get_solution_total( // 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_tangent(void) { - timer.enter_subsection("Assemble tangent matrix"); - std::cout << " ASM_K " << std::flush; + template + void Solid::assemble_system_tangent() + { + 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_tangent_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_tangent_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 separate functions. -template -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_tangent_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()); - } - } - } + template + 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_tangent_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::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()); + } + } + } + + // 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 = 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); - // 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 + 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()); - } - } - } - - // 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); - } - } -} + } + // 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()); + } + } + } + + // Here we copy the lower half of the local + // matrix in the upper half of the local + // matrix + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + { + data.cell_matrix(i, j) = data.cell_matrix(j, i); + } + } + } // @sect4{Solid::assemble_system_rhs} // The assembly of the right-hand side process is similar to the @@ -2068,161 +2509,214 @@ void Solid::assemble_system_tangent_one_cell( // Note that since we are describing a problem with Neumann BCs, // we will need the face normals and so must specify this in the // update flags. -template -void Solid::assemble_system_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_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_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_tilde = lqph[q_point].get_p_tilde(); - 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) { - 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_tilde) * 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; - } - } - } - } - } - } -} + template + void Solid::assemble_system_rhs() + { + timer.enter_subsection("Assemble system right-hand side"); + std::cout << " ASM_R " << std::flush; + + system_rhs = 0.0; + + const UpdateFlags uf_cell(update_values | + update_gradients | + update_JxW_values); + const UpdateFlags uf_face(update_values | + update_normal_vectors | + update_JxW_values); + + PerTaskData_RHS per_task_data(dofs_per_cell); + ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face); + + WorkStream::run(dof_handler_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_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_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_tilde = lqph[q_point].get_p_tilde(); + 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) + { + 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_tilde) * 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} // The constraints for this problem are simple to describe. @@ -2234,108 +2728,142 @@ void Solid::assemble_system_rhs_one_cell( // 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. + 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) { - 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); - } - } + 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) { - 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); - } - } + 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) { - 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); - } - } + 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) { - 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); - } - } + 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) { - 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(); -} + 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 @@ -2344,141 +2872,158 @@ void Solid::make_constraints(const int & it_nr, // 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 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 - // the (hopefully converged) residual - unsigned int lin_it = 0; - double lin_res = 0.0; - - // | K_con | K_up | 0 | | du | | F_u | - // K_store = | K_pu | 0 | K_pJ^-1 | , dXi = | dp | , R = | F_p | - // | 0 | K_Jp | K_JJ | | dJ | | F_J | - - // Solve for the incremental displacement du + template + std::pair + Solid::solve_linear_system(BlockVector & newton_update) + { + // Need 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 the (hopefully converged) + // residual + unsigned int lin_it = 0; + double lin_res = 0.0; + + // | K_con | K_up | 0 | | du | | F_u | + // K_store = | K_pu | 0 | K_pJ^-1 | , dXi = | dp | , R = | F_p | + // | 0 | K_Jp | K_JJ | | dJ | | F_J | + + // Solve for the incremental displacement du + { + // Perform static condensation to make + // K_con = K_uu + K_bbar, 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") { - // Perform static condensation to make - // K_con = K_uu + K_bbar, - // 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 by default a SSOR preconditioner as it appears to provide - // the fastest solver convergence characteristics for this problem. - // However, for multicore computing, the Jacobi preconditioner - // which is multithreaded may converge quicker for larger linear systems. - PreconditionSelector, Vector > preconditioner ( - parameters.preconditioner_type, - parameters.preconditioner_relaxation); - preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof)); - - solver_CG.solve(tangent_matrix.block(u_dof, u_dof), - newton_update.block(u_dof), system_rhs.block(u_dof), - preconditioner); - - lin_it = solver_control.last_step(); - lin_res = solver_control.last_value(); - } else if (parameters.type_lin == "Direct") { - // 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(); - } - - // distribute the constrained dof back to the Newton update - constraints.distribute(newton_update); - - 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 - // dJ = KpJ^{-1} (F_p - K_pu du ) + 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 by default a SSOR + // preconditioner as it appears to + // provide the fastest solver + // convergence characteristics for this + // problem. However, for multicore + // computing, the Jacobi preconditioner + // which is multithreaded may converge + // quicker for larger linear systems. + PreconditionSelector, Vector > + preconditioner (parameters.preconditioner_type, + parameters.preconditioner_relaxation); + preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof)); + + solver_CG.solve(tangent_matrix.block(u_dof, u_dof), + newton_update.block(u_dof), + system_rhs.block(u_dof), + preconditioner); + + lin_it = solver_control.last_step(); + lin_res = solver_control.last_value(); + } + else if (parameters.type_lin == "Direct") { - // 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)); - } - - constraints.distribute(newton_update); - - // and finally we solve for the pressure update with the substitution - // dp = KJp^{-1} [ R_J - K_JJ dJ ] - { - // 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)); - } - - // distribute the constrained dof back to the Newton update - constraints.distribute(newton_update); - - timer.leave_subsection(); - - return std::make_pair(lin_it, lin_res); -} + // 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(); + } + + // distribute the constrained dof back to + // the Newton update + constraints.distribute(newton_update); + + 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 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)); + } + + constraints.distribute(newton_update); + + // and finally we solve for the pressure + // update with the substitution dp = + // KJp^{-1} [ R_J - K_JJ dJ ] + { + // 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)); + } + + // distribute the constrained dof back to + // the Newton update + constraints.distribute(newton_update); + + timer.leave_subsection(); + + return std::make_pair(lin_it, lin_res); + } // @sect4{Solid::assemble_system_SC} // The static condensation process could be performed at a global level @@ -2488,208 +3033,275 @@ std::pair Solid::solve_linear_system( // 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_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(); -} + template + void Solid::assemble_sc() + { + timer.enter_subsection("Perform static condensation"); + std::cout << " ASM_SC " << std::flush; + + PerTaskData_SC per_task_data(dofs_per_cell, element_indices_u.size(), + element_indices_p.size(), element_indices_J.size()); // 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 contributions to K to form K_store + template + 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)); + } -// We need to describe how to add the local contributions -// to K to form K_store -template -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); - - // We now extract the contribution of - // the dof associated with the current cell - // to the global stiffness matrix. - // The discontinuous nature of the $\widetilde{p}$ - // and $\widetilde{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. - // - // This is the strategy we will employ to get the subblocks we want: - // k_store: Since we don't have access to k_{uu}, but we know its contribution is added to the global - // K_{uu} matrix, we just want to add the element wise static-condensation k_bbar. - // k_{pJ}^-1: Similarly, k_pJ exists in the subblock. Since the copy operation is a += operation, we need - // to subtract the existing k_pJ submatrix in addition to "adding" that which we wish to - // replace it with. - // k_{Jp}^-1: Since the global matrix is symmetric, this block is the same as the one above - // and we can simply use k_pJ^-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); -} + 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 $\widetilde{p}$ + // and $\widetilde{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. + // + // This is the strategy we will employ to + // get the subblocks we want: k_store: + // Since we don't have access to k_{uu}, + // but we know its contribution is added to + // the global K_{uu} matrix, we just want + // to add the element wise + // static-condensation k_bbar. + // + // - $k_{pJ}^-1$: Similarly, k_pJ exists in + // the subblock. Since the copy + // operation is a += operation, we + // need to subtract the existing + // k_pJ submatrix in addition to + // "adding" that which we wish to + // replace it with. + // + // - $k_{Jp}^-1$: Since the global matrix + // is symmetric, this block is the + // same as the one above and we + // can simply use k_pJ^-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, then + // 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 // tutorials so will not be discussed in detail. -template -void Solid::output_results(void) const { - 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 Eulerian mapping. We also specify the polynomial degree - // to the DataOut object in order to produce a more refined output data set - // when higher order polynomials are used. - 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() const + { + 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 Eulerian mapping. We + // also specify the polynomial degree to + // the DataOut object in order to produce a + // more refined output data set when higher + // order polynomials are used. + 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 (int argc, char *argv[]) +{ + using namespace dealii; + using namespace Step44; + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + + 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; }