/* 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 */
/* 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 <deal.II/base/function.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/point.h>
#include <iostream>
#include <fstream>
-// 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
// 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
// We place these common operations
// in a separate namespace for convenience.
// We also include some widely used operators
-namespace AdditionalTools {
-// Define an operation that takes two
-// symmetric second-order tensors
-// $\mathbf{A}$ and $\mathbf{B}$
-// such that their outer-product
-// $ \mathbf{A} \overline{\otimes} \mathbf{B} \Rightarrow C_{ijkl} = A_{ik} B_{jl} $
-template<int dim>
-SymmetricTensor<4, dim> outer_product_T23(const SymmetricTensor<2, dim> & A,
- const SymmetricTensor<2, dim> & B) {
- SymmetricTensor<4, dim> A_ik_B_jl;
-
- for (unsigned int i = 0; i < dim; ++i) {
- for (unsigned int j = i; j < dim; ++j) {
- for (unsigned int k = 0; k < dim; ++k) {
- for (unsigned int l = k; k < dim; ++k) {
- A_ik_B_jl[i][j][k][l] += A[i][k] * B[j][l];
- }
- }
- }
- }
-
- return A_ik_B_jl;
-}
+ namespace AdditionalTools
+ {
// The extract_submatrix function
// takes specific entries from a matrix,
// first two parameters which hold the
// row and columns to be extracted.
// The matrix is automatically resized
-// to size $ r \times c $.
-template<typename MatrixType>
-void extract_submatrix(const std::vector<unsigned int> &row_index_set,
- const std::vector<unsigned int> &column_index_set,
- const MatrixType &matrix, FullMatrix<double> &sub_matrix) {
-
- const unsigned int n_rows_submatrix = row_index_set.size();
- const unsigned int n_cols_submatrix = column_index_set.size();
- // check the size of the input vectors
- Assert(n_rows_submatrix > 0, ExcInternalError());
- Assert(n_cols_submatrix > 0, ExcInternalError());
-
- sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix);
-
- for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) {
- const unsigned int row = row_index_set[sub_row];
- Assert(row<=matrix.m(), ExcInternalError());
+// to size $ r \times c $. At the beginning we
+// check the size of the input vectors
+ template <typename MatrixType>
+ void extract_submatrix (const std::vector<unsigned int> &row_index_set,
+ const std::vector<unsigned int> &column_index_set,
+ const MatrixType &matrix,
+ FullMatrix<double> &sub_matrix)
+ {
+
+ const unsigned int n_rows_submatrix = row_index_set.size();
+ const unsigned int n_cols_submatrix = column_index_set.size();
+ 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 <code> BlockSparseMatrix </code>.
-template<>
-void extract_submatrix<dealii::BlockSparseMatrix<double> >(
- const std::vector<unsigned int> &row_index_set,
- const std::vector<unsigned int> &column_index_set,
- const dealii::BlockSparseMatrix<double> &matrix,
- FullMatrix<double> &sub_matrix) {
+ template <>
+ void
+ extract_submatrix<dealii::BlockSparseMatrix<double> >
+ (const std::vector<unsigned int> &row_index_set,
+ const std::vector<unsigned int> &column_index_set,
+ const dealii::BlockSparseMatrix<double> &matrix,
+ FullMatrix<double> &sub_matrix)
+ {
- const unsigned int n_rows_submatrix = row_index_set.size();
- const unsigned int n_cols_submatrix = column_index_set.size();
+ const unsigned int n_rows_submatrix = row_index_set.size();
+ const unsigned int n_cols_submatrix = column_index_set.size();
- // 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,
// first two parameters which hold the
// row and column entries to be replaced.
// The matrix expected to be of the correct size.
-template<typename MatrixType>
-void replace_submatrix(const std::vector<unsigned int> &row_index_set,
- const std::vector<unsigned int> &column_index_set,
- const MatrixType &sub_matrix, FullMatrix<double> &matrix) {
- const unsigned int n_rows_submatrix = row_index_set.size();
- Assert(n_rows_submatrix<=sub_matrix.m(), ExcInternalError());
- const unsigned int n_cols_submatrix = column_index_set.size();
- Assert(n_cols_submatrix<=sub_matrix.n(), ExcInternalError());
-
- for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) {
- const unsigned int row = row_index_set[sub_row];
- Assert(row<=matrix.m(), ExcInternalError());
+ template <typename MatrixType>
+ void
+ replace_submatrix(const std::vector<unsigned int> &row_index_set,
+ const std::vector<unsigned int> &column_index_set,
+ const MatrixType &sub_matrix,
+ FullMatrix<double> &matrix)
+ {
+ const unsigned int n_rows_submatrix = row_index_set.size();
+ Assert(n_rows_submatrix<=sub_matrix.m(), ExcInternalError());
+ const unsigned int n_cols_submatrix = column_index_set.size();
+ Assert(n_cols_submatrix<=sub_matrix.n(), ExcInternalError());
+
+ for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
+ {
+ const unsigned int row = row_index_set[sub_row];
+ Assert(row<=matrix.m(), ExcInternalError());
- for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) {
- const unsigned int col = column_index_set[sub_col];
- Assert(col<=matrix.n(), ExcInternalError());
+ 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<int dim>
-class StandardTensors {
-public:
-
- // $\mathbf{I}$
- static SymmetricTensor<2, dim> const I;
- // $\mathbf{I} \otimes \mathbf{I}$
- static SymmetricTensor<4, dim> const IxI;
- // $\mathcal{S}$, note that as we only use
- // this fourth-order unit tensor to operate
- // on symmetric second-order tensors.
- // To maintain notation consistent with Holzapfel (2001)
- // we name the tensor $\mathcal{I}$
- static SymmetricTensor<4, dim> const II;
- // Fourth-order deviatoric such that
- // $\textrm{dev}(\bullet) = (\bullet) - (1/3)[(\bullet):\mathbf{I}]\mathbf{I}$
- static SymmetricTensor<4, dim> const dev_P;
-};
-
-template<int dim>
-SymmetricTensor<2, dim> const StandardTensors<dim>::I = SymmetricTensor<2, dim>(
- unit_symmetric_tensor<dim>());
-template<int dim>
-SymmetricTensor<4, dim> const StandardTensors<dim>::IxI =
- SymmetricTensor<4, dim>(outer_product(I, I));
-template<int dim>
-SymmetricTensor<4, dim> const StandardTensors<dim>::II =
- SymmetricTensor<4, dim>(identity_tensor<dim>());
-template<int dim>
-SymmetricTensor<4, dim> const StandardTensors<dim>::dev_P = (II
- - (1.0 / dim) * IxI);
-}
+ template <int dim>
+ 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 <int dim>
+ const SymmetricTensor<2, dim>
+ StandardTensors<dim>::I = unit_symmetric_tensor<dim>();
+
+ template <int dim>
+ const SymmetricTensor<4, dim>
+ StandardTensors<dim>::IxI = outer_product(I, I);
+
+ template <int dim>
+ const SymmetricTensor<4, dim>
+ StandardTensors<dim>::II = identity_tensor<dim>();
+
+ template <int dim>
+ const SymmetricTensor<4, dim>
+ StandardTensors<dim>::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 <a href="http://en.wikipedia.org/wiki/Bulk_modulus">bulk modulus</a> and
-// $\lambda$ is <a href="http://en.wikipedia.org/wiki/Lam%C3%A9_parameters">Lame's first parameter</a>.
-template<int dim>
-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<dim>::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 <a href="http://en.wikipedia.org/wiki/Bulk_modulus">bulk modulus</a>
+// and $\lambda$ is <a
+// href="http://en.wikipedia.org/wiki/Lam%C3%A9_parameters">Lame's first
+// parameter</a>.
+//
+// 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 <int dim>
+ 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<dim>::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<dim>::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<dim>::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<dim>::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<dim>::dev_P * get_tau_bar();
+ return p_tilde * det_F
+ * ( AdditionalTools::StandardTensors<dim>::IxI
+ - (2.0 * AdditionalTools::StandardTensors<dim>::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<dim>::I);
+ const SymmetricTensor<4, dim> I_x_tau_iso
+ = outer_product(AdditionalTools::StandardTensors<dim>::I,
+ tau_iso);
+ const SymmetricTensor<4, dim> c_bar = get_c_bar();
+
+ return (2.0 / 3.0) * trace(tau_bar)
+ * AdditionalTools::StandardTensors<dim>::dev_P
+ - (2.0 / 3.0) * (tau_iso_x_I + I_x_tau_iso)
+ + AdditionalTools::StandardTensors<dim>::dev_P * c_bar
+ * AdditionalTools::StandardTensors<dim>::dev_P;
+ }
+
+ // 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<dim>::IxI
- - (2.0 * AdditionalTools::StandardTensors<dim>::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<dim>::I);
- const SymmetricTensor<4, dim> I_x_tau_iso = outer_product(
- AdditionalTools::StandardTensors<dim>::I, tau_iso);
- const SymmetricTensor<4, dim> c_bar = get_c_bar();
-
- return (2.0 / 3.0) * trace(tau_bar)
- * AdditionalTools::StandardTensors<dim>::dev_P
- - (2.0 / 3.0) * (tau_iso_x_I + I_x_tau_iso)
- + AdditionalTools::StandardTensors<dim>::dev_P * c_bar
- * AdditionalTools::StandardTensors<dim>::dev_P;
+// As seen in step-18, the <code>
+// PointHistory </code> 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 <int dim>
+ class PointHistory
+ {
+ public:
+ PointHistory()
+ :
+ material(NULL),
+ F_inv(AdditionalTools::StandardTensors<dim>::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<dim>(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<dim>::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 <code> PointHistory </code> 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<int dim>
-class PointHistory {
-public:
- PointHistory(void) :
- material(NULL),
- F_inv(AdditionalTools::StandardTensors<dim>::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<dim>(
- 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<Tensor<2, dim> >(AdditionalTools::StandardTensors<
- dim>::I);
- const Tensor<2, dim> F = I + Grad_u_n;
-
- // We use the inverse of $\mathbf{F}$ frequently so we store it
- F_inv = invert(F);
-
- // 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<dim>* 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<dim>* 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<int dim>
-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<dim>::active_cell_iterator & cell,
- ScratchData_K & scratch, PerTaskData_K & data);
- void
- copy_local_to_global_K(const PerTaskData_K & data);
- void
- assemble_system_rhs(void);
- void
- assemble_system_rhs_one_cell(
- const typename DoFHandler<dim>::active_cell_iterator & cell,
- ScratchData_RHS & scratch, PerTaskData_RHS & data);
- void
- copy_local_to_global_rhs(const PerTaskData_RHS & data);
- void
- assemble_sc(void);
- void
- assemble_sc_one_cell(
- const typename DoFHandler<dim>::active_cell_iterator & cell,
- ScratchData_SC & scratch, PerTaskData_SC & data);
- void
- copy_local_to_global_sc(const PerTaskData_SC & data);
- // Apply Dirichlet boundary 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<double> & solution_delta);
- void
- update_qph_incremental_one_cell(
- const typename DoFHandler<dim>::active_cell_iterator & cell,
- ScratchData_UQPH & scratch, PerTaskData_UQPH & data);
- void copy_local_to_global_UQPH(const PerTaskData_UQPH & data) {
- }
-
- // Solve for the displacement using a Newton-Raphson method
- void
- solve_nonlinear_timestep(BlockVector<double> & solution_delta);
- std::pair<unsigned int, double>
- solve_linear_system(BlockVector<double> & newton_update);
-
- // Solution retrieval
- BlockVector<double>
- get_solution_total(const BlockVector<double> & solution_delta) 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<dim> triangulation;
-
- // Keep track of the current time and the time spent evaluating certain functions
- Time time;
- TimerOutput timer;
-
- // A storage object for quadrature point information.
- // See step-18 for more on this
- std::vector<PointHistory<dim> > quadrature_point_history;
-
- // A description of the finite-element system including the displacement polynomial degree,
- // the degree-of-freedom handler, number of dof's per cell and the extractor objects used
- // to retrieve information from the solution vectors
- const unsigned int degree;
- const FESystem<dim> fe;
- DoFHandler<dim> dof_handler_ref;
- unsigned int dofs_per_cell;
- const FEValuesExtractors::Vector u_fe;
- const FEValuesExtractors::Scalar p_fe;
- const FEValuesExtractors::Scalar J_fe;
-
- // Description of how the block-system is arranged
- // There are 3 blocks, the first contains a vector DOF $\mathbf{u}$
- // while the other two describe scalar DOFs,
- // $\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<unsigned int> dofs_per_block;
- std::vector<unsigned int> element_indices_u;
- std::vector<unsigned int> element_indices_p;
- std::vector<unsigned int> element_indices_J;
-
- // Rules for Gauss-quadrature on both the cell and faces. The
- // number of quadrature points on both cells and faces is
- // recorded.
- QGauss<dim> qf_cell;
- QGauss<dim - 1> qf_face;
- unsigned int n_q_points;
- unsigned int n_q_points_f;
-
- // Objects that store the converged solution and right-hand side vectors,
- // as well as the tangent matrix. There is a ConstraintMatrix object
- // used to keep track of constraints.
- // We make use of a sparsity pattern designed for a block system.
- ConstraintMatrix constraints;
- BlockSparsityPattern sparsity_pattern;
- BlockSparseMatrix<double> tangent_matrix;
- BlockVector<double> system_rhs;
- BlockVector<double> solution_n;
-
- // Then define a number of variables to store norms and update
- // norms and normalisation factors.
- struct Errors {
- Errors(void) :
- norm(1.0), u(1.0), p(1.0), J(1.0) {
- }
- double norm, u, p, J;
- void reset(void) {
- norm = 1.0;
- u = 1.0;
- p = 1.0;
- J = 1.0;
- }
- void normalise(const Errors & rhs) {
- if (rhs.norm != 0.0)
- norm /= rhs.norm;
- if (rhs.u != 0.0)
- u /= rhs.u;
- if (rhs.p != 0.0)
- p /= rhs.p;
- if (rhs.J != 0.0)
- J /= rhs.J;
- }
- } error_residual, error_residual_0, error_residual_norm, error_update,
- error_update_0, error_update_norm;
-
- // Methods to calculate error measures
- void
- get_error_residual(Errors & error_residual);
- void
- get_error_update(const BlockVector<double> & newton_update,
- Errors & error_update);
- std::pair<double, double>
- 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 <code>run()</code> function that dispatches
+// all the work to private functions of this class:
+ template <int dim>
+ 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<dim>::active_cell_iterator & cell,
+ ScratchData_K & scratch,
+ PerTaskData_K & data);
+
+ void
+ copy_local_to_global_K(const PerTaskData_K & data);
+
+ void
+ assemble_system_rhs();
+
+ void
+ assemble_system_rhs_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,
+ ScratchData_RHS & scratch,
+ PerTaskData_RHS & data);
+
+ void
+ copy_local_to_global_rhs(const PerTaskData_RHS & data);
+
+ void
+ assemble_sc();
+
+ void
+ assemble_sc_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,
+ ScratchData_SC & scratch,
+ PerTaskData_SC & data);
+
+ void
+ copy_local_to_global_sc(const PerTaskData_SC & data);
+
+ // Apply Dirichlet boundary 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<double> & solution_delta);
+
+ void
+ update_qph_incremental_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,
+ ScratchData_UQPH & scratch,
+ PerTaskData_UQPH & data);
+
+ void
+ copy_local_to_global_UQPH(const PerTaskData_UQPH & data)
+ {}
+
+ // Solve for the displacement using a
+ // Newton-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<double> & solution_delta);
+
+ std::pair<unsigned int, double>
+ solve_linear_system(BlockVector<double> & newton_update);
+
+ // Solution retrieval as well as
+ // post-processing and writing data to
+ // file:
+ BlockVector<double>
+ get_solution_total(const BlockVector<double> & 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<dim> 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<PointHistory<dim> > quadrature_point_history;
+
+ // A description of the finite-element
+ // system including the displacement
+ // polynomial degree, the
+ // degree-of-freedom handler, number of
+ // dof's per cell and the extractor
+ // objects used to retrieve information
+ // from the solution vectors:
+ const unsigned int degree;
+ const FESystem<dim> fe;
+ DoFHandler<dim> dof_handler_ref;
+ 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<unsigned int> dofs_per_block;
+ std::vector<unsigned int> element_indices_u;
+ std::vector<unsigned int> element_indices_p;
+ std::vector<unsigned int> element_indices_J;
+
+ // Rules for Gauss-quadrature on both the
+ // cell and faces. The number of
+ // quadrature points on both cells and
+ // faces is recorded.
+ const QGauss<dim> qf_cell;
+ const QGauss<dim - 1> 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<double> tangent_matrix;
+ BlockVector<double> system_rhs;
+ BlockVector<double> solution_n;
+
+ // Then define a number of variables to
+ // store norms and update norms and
+ // normalisation factors.
+ struct Errors
+ {
+ Errors()
+ :
+ 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<double> & newton_update,
+ Errors & error_update);
+
+ std::pair<double, double>
+ get_error_dil();
+
+ // Print information to screen
+ static
+ void
+ print_conv_header();
+
+ void
+ print_conv_footer();
+ };
// @sect3{Implementation of the <code>Solid</code> class}
// @sect4{Public interface}
// We initialise the Solid class using data extracted
// from the parameter file.
-template<int dim>
-Solid<dim>::Solid(const std::string & input_file) :
-parameters(input_file), triangulation(
- Triangulation<dim>::maximum_smoothing), time(
- parameters.end_time, parameters.delta_t), timer(std::cout,
- TimerOutput::summary, TimerOutput::wall_times), degree(
- parameters.poly_degree),
- // The Finite Element System is composed of dim continuous
- // displacement DOFs, and discontinuous pressure and
- // dilatation DOFs. In an attempt to satisfy the LBB conditions,
- // we setup a Q(n)-P(n-1)-P(n-1) system. Q2-P1-P1 elements satisfy
- // this condition, while Q1-P0-P0 elements do not. However, it
- // has been shown that the latter demonstrate good convergence
- // characteristics nonetheless.
- fe(FE_Q<dim>(parameters.poly_degree), dim, // displacement
- FE_DGPMonomial<dim>(parameters.poly_degree - 1), 1, // pressure
- FE_DGPMonomial<dim>(parameters.poly_degree - 1), 1), // dilatation
- dof_handler_ref(triangulation), u_fe(first_u_component), p_fe(
- p_component), J_fe(J_component), dofs_per_block(n_blocks), qf_cell(
- parameters.quad_order), qf_face(parameters.quad_order) {
- n_q_points = qf_cell.size();
- n_q_points_f = qf_face.size();
- dofs_per_cell = fe.dofs_per_cell;
- determine_component_extractors();
-}
+ template <int dim>
+ Solid<dim>::Solid(const std::string & input_file)
+ :
+ parameters(input_file),
+ triangulation(Triangulation<dim>::maximum_smoothing),
+ time(parameters.end_time, parameters.delta_t),
+ timer(std::cout,
+ TimerOutput::summary,
+ TimerOutput::wall_times),
+ degree(parameters.poly_degree),
+ // The Finite Element
+ // System is composed of
+ // dim continuous
+ // displacement DOFs, and
+ // discontinuous pressure
+ // and dilatation DOFs. In
+ // an attempt to satisfy
+ // the LBB conditions, we
+ // setup a $Q_n \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<dim>(parameters.poly_degree), dim, // displacement
+ FE_DGPMonomial<dim>(parameters.poly_degree - 1), 1, // pressure
+ FE_DGPMonomial<dim>(parameters.poly_degree - 1), 1), // dilatation
+ dof_handler_ref(triangulation),
+ 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<int dim>
-Solid<dim>::~Solid(void) {
- dof_handler_ref.clear();
-}
+ template <int dim>
+ Solid<dim>::~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<int dim>
-void Solid<dim>::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 <int dim>
+ void Solid<dim>::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<double> solution_delta(dofs_per_block);
+ solution_delta.collect_sizes();
+
+ // Now we loop over the time domain
+ while (time.current() < time.end())
+ {
+ // We need to reset the solution update
+ // for this time step
+ solution_delta = 0.0;
+
+ // Solve the current time step and update total
+ // solution vector
+ solve_nonlinear_timestep(solution_delta);
+ // $\varDelta \mathbf{\Xi}_{\textrm{n}} =
+ // \varDelta \mathbf{\Xi}_{\textrm{n-1}}
+ // + \varDelta \mathbf{\Xi}$
+ solution_n += solution_delta;
+ // 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<double> solution_delta(dofs_per_block);
- solution_delta.collect_sizes();
-
- // Now we loop over the time domain
- while (time.current() < time.end()) {
- // We need to reset the solution update
- // for this time step
- solution_delta = 0.0;
-
- // Solve the current time step and update total
- // solution vector
- solve_nonlinear_timestep(solution_delta);
- // $\varDelta \mathbf{\Xi}_{\textrm{n}} = \varDelta \mathbf{\Xi}_{\textrm{n-1}} + \varDelta \mathbf{\Xi}$
- solution_n += solution_delta;
- // and plot the results
- output_results();
- // we then move on happily to the next time step.
- time.increment();
- }
-}
+ }
+ }
// @sect3{Private interface}
// at the quadrature points using TBB.
// Firstly we deal with the tangent matrix assembly structures.
-// The PerTaskData object stores local contributions.
-template<int dim>
-struct Solid<dim>::PerTaskData_K {
- FullMatrix<double> cell_matrix;
- std::vector<unsigned int> local_dof_indices;
-
- PerTaskData_K(const unsigned int dofs_per_cell) :
- cell_matrix(dofs_per_cell, dofs_per_cell), local_dof_indices(
- dofs_per_cell) {
+// The PerTaskData object stores local contributions.
+ template <int dim>
+ struct Solid<dim>::PerTaskData_K
+ {
+ FullMatrix<double> cell_matrix;
+ std::vector<unsigned int> local_dof_indices;
+
+ PerTaskData_K(const unsigned int dofs_per_cell)
+ :
+ cell_matrix(dofs_per_cell, dofs_per_cell),
+ local_dof_indices(dofs_per_cell)
+ {}
+
+ void reset()
+ {
+ 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<int dim>
-struct Solid<dim>::ScratchData_K {
- FEValues<dim> fe_values_ref;
-
- // interpolation function
- std::vector<std::vector<double> > Nx;
- // their gradients
- std::vector<std::vector<Tensor<2, dim> > > grad_Nx;
- // and their symmetric gradients.
- std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
-
- ScratchData_K(const FiniteElement<dim> & fe_cell,
- const QGauss<dim> & qf_cell, const UpdateFlags uf_cell) :
- fe_values_ref(fe_cell, qf_cell, uf_cell), Nx(qf_cell.size(),
- std::vector<double>(fe_cell.dofs_per_cell)), grad_Nx(
- qf_cell.size(),
- std::vector<Tensor<2, dim> >(fe_cell.dofs_per_cell)), symm_grad_Nx(
- qf_cell.size(),
- std::vector<SymmetricTensor<2, dim> >(
- fe_cell.dofs_per_cell)) {
- }
-
- ScratchData_K(const ScratchData_K & rhs) :
- fe_values_ref(rhs.fe_values_ref.get_fe(),
- rhs.fe_values_ref.get_quadrature(),
- rhs.fe_values_ref.get_update_flags()), Nx(rhs.Nx), grad_Nx(
- rhs.grad_Nx), symm_grad_Nx(rhs.symm_grad_Nx) {
- }
-
- void reset(void) {
- 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 <int dim>
+ struct Solid<dim>::ScratchData_K
+ {
+ FEValues<dim> fe_values_ref;
+
+ // interpolation function
+ std::vector<std::vector<double> > Nx;
+ // their gradients
+ std::vector<std::vector<Tensor<2, dim> > > grad_Nx;
+ // and their symmetric gradients.
+ std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
+
+ ScratchData_K(const FiniteElement<dim> & fe_cell,
+ const QGauss<dim> & qf_cell,
+ const UpdateFlags uf_cell)
+ :
+ fe_values_ref(fe_cell, qf_cell, uf_cell),
+ Nx(qf_cell.size(),
+ std::vector<double>(fe_cell.dofs_per_cell)),
+ grad_Nx(qf_cell.size(),
+ std::vector<Tensor<2, dim> >(fe_cell.dofs_per_cell)),
+ symm_grad_Nx(qf_cell.size(),
+ std::vector<SymmetricTensor<2, dim> >
+ (fe_cell.dofs_per_cell))
+ {}
+
+ ScratchData_K(const ScratchData_K & rhs)
+ :
+ fe_values_ref(rhs.fe_values_ref.get_fe(),
+ rhs.fe_values_ref.get_quadrature(),
+ rhs.fe_values_ref.get_update_flags()),
+ Nx(rhs.Nx),
+ grad_Nx(rhs.grad_Nx),
+ symm_grad_Nx(rhs.symm_grad_Nx)
+ {}
+
+ void reset()
+ {
+ 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<int dim>
-struct Solid<dim>::PerTaskData_RHS {
- Vector<double> cell_rhs;
- std::vector<unsigned int> 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 <int dim>
+ struct Solid<dim>::PerTaskData_RHS
+ {
+ Vector<double> cell_rhs;
+ std::vector<unsigned int> 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<int dim>
-struct Solid<dim>::ScratchData_RHS {
- FEValues<dim> fe_values_ref;
- FEFaceValues<dim> fe_face_values_ref;
-
- std::vector<std::vector<double> > Nx;
- std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
-
- ScratchData_RHS(const FiniteElement<dim> & fe_cell,
- const QGauss<dim> & qf_cell, const UpdateFlags uf_cell,
- const QGauss<dim - 1> & qf_face, const UpdateFlags uf_face) :
- fe_values_ref(fe_cell, qf_cell, uf_cell), fe_face_values_ref(
- fe_cell, qf_face, uf_face), Nx(qf_cell.size(),
- std::vector<double>(fe_cell.dofs_per_cell)), symm_grad_Nx(
- qf_cell.size(),
- std::vector<SymmetricTensor<2, dim> >(
- fe_cell.dofs_per_cell)) {
- }
-
- ScratchData_RHS(const ScratchData_RHS & rhs) :
- fe_values_ref(rhs.fe_values_ref.get_fe(),
- rhs.fe_values_ref.get_quadrature(),
- rhs.fe_values_ref.get_update_flags()), fe_face_values_ref(
- rhs.fe_face_values_ref.get_fe(),
- rhs.fe_face_values_ref.get_quadrature(),
- rhs.fe_face_values_ref.get_update_flags()), Nx(rhs.Nx), symm_grad_Nx(
- rhs.symm_grad_Nx) {
- }
-
- void reset(void) {
- 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 <int dim>
+ struct Solid<dim>::ScratchData_RHS
+ {
+ FEValues<dim> fe_values_ref;
+ FEFaceValues<dim> fe_face_values_ref;
+
+ std::vector<std::vector<double> > Nx;
+ std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
+
+ ScratchData_RHS(const FiniteElement<dim> & fe_cell,
+ const QGauss<dim> & qf_cell, const UpdateFlags uf_cell,
+ const QGauss<dim - 1> & qf_face, const UpdateFlags uf_face)
+ :
+ fe_values_ref(fe_cell, qf_cell, uf_cell),
+ fe_face_values_ref(fe_cell, qf_face, uf_face),
+ Nx(qf_cell.size(),
+ std::vector<double>(fe_cell.dofs_per_cell)),
+ symm_grad_Nx(qf_cell.size(),
+ std::vector<SymmetricTensor<2, dim> >
+ (fe_cell.dofs_per_cell))
+ {}
+
+ ScratchData_RHS(const ScratchData_RHS & rhs)
+ :
+ fe_values_ref(rhs.fe_values_ref.get_fe(),
+ rhs.fe_values_ref.get_quadrature(),
+ rhs.fe_values_ref.get_update_flags()),
+ fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
+ rhs.fe_face_values_ref.get_quadrature(),
+ rhs.fe_face_values_ref.get_update_flags()),
+ Nx(rhs.Nx),
+ symm_grad_Nx(rhs.symm_grad_Nx)
+ {}
+
+ void reset()
+ {
+ 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<int dim>
-struct Solid<dim>::PerTaskData_SC {
- FullMatrix<double> cell_matrix;
- std::vector<unsigned int> local_dof_indices;
-
- FullMatrix<double> k_orig;
- FullMatrix<double> k_pu;
- FullMatrix<double> k_pJ;
- FullMatrix<double> k_JJ;
- FullMatrix<double> k_pJ_inv;
- FullMatrix<double> k_bbar;
- FullMatrix<double> A;
- FullMatrix<double> B;
- FullMatrix<double> C;
-
- PerTaskData_SC(const unsigned int dofs_per_cell, const unsigned int n_u,
- const unsigned int n_p, const unsigned int n_J) :
- cell_matrix(dofs_per_cell, dofs_per_cell),
- local_dof_indices(dofs_per_cell),
- k_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 <int dim>
+ struct Solid<dim>::PerTaskData_SC
+ {
+ FullMatrix<double> cell_matrix;
+ std::vector<unsigned int> local_dof_indices;
+
+ FullMatrix<double> k_orig;
+ FullMatrix<double> k_pu;
+ FullMatrix<double> k_pJ;
+ FullMatrix<double> k_JJ;
+ FullMatrix<double> k_pJ_inv;
+ FullMatrix<double> k_bbar;
+ FullMatrix<double> A;
+ FullMatrix<double> B;
+ FullMatrix<double> C;
+
+ PerTaskData_SC(const unsigned int dofs_per_cell,
+ const unsigned int n_u,
+ const unsigned int n_p,
+ const unsigned int n_J)
+ :
+ cell_matrix(dofs_per_cell, dofs_per_cell),
+ local_dof_indices(dofs_per_cell),
+ k_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 <int dim>
+ struct Solid<dim>::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<int dim>
-struct Solid<dim>::ScratchData_SC {
- ScratchData_SC(void) {
- }
- ScratchData_SC(const ScratchData_SC & rhs) {
- }
- void reset(void) {
- }
-};
// And finally we define the structures to assist with updating the quadrature
// point information. Similar to the SC assembly process, we choose not to use
-// the PerTaskData object to store any information but must define one nonetheless.
-template<int dim>
-struct Solid<dim>::PerTaskData_UQPH {
- PerTaskData_UQPH(void) {
- }
- void reset(void) {
- }
-};
-// The ScratchData object will be used to store an alias for the solution vector
-// so that we don't have to copy this large data structure. We then define
-// a number of vectors to extract the solution values and gradients at the
-// quadrature points.
-template<int dim>
-struct Solid<dim>::ScratchData_UQPH {
-
- const BlockVector<double> & solution_total;
-
- std::vector<Tensor<2, dim> > solution_grads_u_total;
- std::vector<double> solution_values_p_total;
- std::vector<double> solution_values_J_total;
-
- FEValues<dim> fe_values_ref;
-
- ScratchData_UQPH(const FiniteElement<dim> & fe_cell,
- const QGauss<dim> & qf_cell,
- const UpdateFlags uf_cell,
- const BlockVector<double> & solution_total) :
- solution_total(solution_total),
- solution_grads_u_total(qf_cell.size()),
- solution_values_p_total(qf_cell.size()),
- solution_values_J_total(qf_cell.size()),
- fe_values_ref(fe_cell, qf_cell, uf_cell) {
- }
-
- ScratchData_UQPH(const ScratchData_UQPH & rhs) :
- solution_total(rhs.solution_total), solution_grads_u_total(
- rhs.solution_grads_u_total), solution_values_p_total(
- rhs.solution_values_p_total), solution_values_J_total(
- rhs.solution_values_J_total), fe_values_ref(
- rhs.fe_values_ref.get_fe(),
- rhs.fe_values_ref.get_quadrature(),
- rhs.fe_values_ref.get_update_flags()) {
+// the PerTaskData object to store any information but must define one
+// nonetheless.
+ template <int dim>
+ struct Solid<dim>::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 <int dim>
+ struct Solid<dim>::ScratchData_UQPH
+ {
+ const BlockVector<double> &solution_total;
+
+ std::vector<Tensor<2, dim> > solution_grads_u_total;
+ std::vector<double> solution_values_p_total;
+ std::vector<double> solution_values_J_total;
+
+ FEValues<dim> fe_values_ref;
+
+ ScratchData_UQPH(const FiniteElement<dim> & fe_cell,
+ const QGauss<dim> & qf_cell,
+ const UpdateFlags uf_cell,
+ const BlockVector<double> & solution_total)
+ :
+ solution_total(solution_total),
+ solution_grads_u_total(qf_cell.size()),
+ solution_values_p_total(qf_cell.size()),
+ solution_values_J_total(qf_cell.size()),
+ fe_values_ref(fe_cell, qf_cell, uf_cell)
+ {}
+
+ ScratchData_UQPH(const ScratchData_UQPH & rhs)
+ :
+ solution_total(rhs.solution_total),
+ solution_grads_u_total(rhs.solution_grads_u_total),
+ solution_values_p_total(rhs.solution_values_p_total),
+ solution_values_J_total(rhs.solution_values_J_total),
+ fe_values_ref(rhs.fe_values_ref.get_fe(),
+ rhs.fe_values_ref.get_quadrature(),
+ rhs.fe_values_ref.get_update_flags())
+ {}
+
+ void reset()
+ {
+ 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<int dim>
-void Solid<dim>::make_grid(void) {
- // Create a unit cube with each face given a boundary ID number
- GridGenerator::hyper_rectangle(triangulation, Point<dim>(0.0, 0.0, 0.0),
- Point<dim>(1.0, 1.0, 1.0), true);
- GridTools::scale(parameters.scale, triangulation);
-
- // The grid must be refined at least once for the indentation problem
- if (parameters.global_refinement == 0)
- triangulation.refine_global(1);
- else
- triangulation.refine_global(parameters.global_refinement);
-
- // 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<dim>::active_cell_iterator cell =
- triangulation.begin_active(), endc = triangulation.end();
- for (; cell != endc; ++cell) {
- if (cell->at_boundary() == true) {
- for (unsigned int face = 0;
- face < GeometryInfo<dim>::faces_per_cell; ++face) {
- // Find faces on the +y surface
- if (cell->face(face)->at_boundary() == true
- && cell->face(face)->center()[2]
- == 1.0 * parameters.scale) {
- if (cell->face(face)->center()[0] < 0.5 * parameters.scale
- && cell->face(face)->center()[1]
- < 0.5 * parameters.scale) {
- cell->face(face)->set_boundary_indicator(6); // Set a new boundary id on a patch
- }
- }
- }
+ template <int dim>
+ void Solid<dim>::make_grid()
+ {
+ // Create a unit cube with each face given
+ // a boundary ID number
+ GridGenerator::hyper_rectangle(triangulation, Point<dim>(0.0, 0.0, 0.0),
+ Point<dim>(1.0, 1.0, 1.0), true);
+ GridTools::scale(parameters.scale, triangulation);
+
+ // The grid must be refined at least once
+ // for the indentation problem
+ if (parameters.global_refinement == 0)
+ triangulation.refine_global(1);
+ else
+ triangulation.refine_global(parameters.global_refinement);
+
+ // 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<dim>::active_cell_iterator cell =
+ triangulation.begin_active(), endc = triangulation.end();
+ for (; cell != endc; ++cell)
+ {
+ if (cell->at_boundary() == true)
+ {
+ for (unsigned int face = 0;
+ face < GeometryInfo<dim>::faces_per_cell; ++face)
+ {
+ // Find faces on the +y surface
+ if (cell->face(face)->at_boundary() == true
+ && cell->face(face)->center()[2]
+ == 1.0 * parameters.scale) {
+ if (cell->face(face)->center()[0] < 0.5 * parameters.scale
+ && cell->face(face)->center()[1]
+ < 0.5 * parameters.scale) {
+ cell->face(face)->set_boundary_indicator(6); // Set a new boundary id on a patch
+ }
}
- }
-}
+ }
+ }
+ }
+ }
+
// @sect4{Solid::system_setup}
// Next we describe how the FE system is setup.
-template<int dim>
-void Solid<dim>::system_setup(void) {
- timer.enter_subsection("Setup system");
-
- // We first describe the number of components per block. Since the
- // displacement is a vector component, the first dim components
- // belong to it, while the next two describe scalar pressure and
- // dilatation DOFs.
- std::vector<unsigned int> block_component(n_components, u_dof); // Displacement
- block_component[p_component] = p_dof; // Pressure
- block_component[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 <int dim>
+ void Solid<dim>::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<unsigned int> block_component(n_components, u_dof); // Displacement
+ block_component[p_component] = p_dof; // Pressure
+ block_component[J_component] = J_dof; // Dilatation
+
+ // DOF handler is then initialised and we
+ // renumber the grid in an efficient
+ // manner. We also record the number of
+ // DOF's per block.
+ dof_handler_ref.distribute_dofs(fe);
+ DoFRenumbering::Cuthill_McKee(dof_handler_ref);
+ DoFRenumbering::component_wise(dof_handler_ref, block_component);
+ DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block,
+ block_component);
+
+ std::cout << "Triangulation:"
+ << "\n\t Number of active cells: " << triangulation.n_active_cells()
+ << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs()
+ << std::endl;
+
+ // Setup the sparsity pattern and tangent matrix
+ tangent_matrix.clear();
+ {
+ const unsigned int n_dofs_u = dofs_per_block[u_dof];
+ const unsigned int n_dofs_p = dofs_per_block[p_dof];
+ const unsigned int n_dofs_J = dofs_per_block[J_dof];
+
+ BlockCompressedSimpleSparsityPattern csp(n_blocks, n_blocks);
+
+ csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
+ csp.block(u_dof, p_dof).reinit(n_dofs_u, n_dofs_p);
+ csp.block(u_dof, J_dof).reinit(n_dofs_u, n_dofs_J);
+
+ csp.block(p_dof, u_dof).reinit(n_dofs_p, n_dofs_u);
+ csp.block(p_dof, p_dof).reinit(n_dofs_p, n_dofs_p);
+ csp.block(p_dof, J_dof).reinit(n_dofs_p, n_dofs_J);
+
+ csp.block(J_dof, u_dof).reinit(n_dofs_J, n_dofs_u);
+ csp.block(J_dof, p_dof).reinit(n_dofs_J, n_dofs_p);
+ csp.block(J_dof, J_dof).reinit(n_dofs_J, n_dofs_J);
+ csp.collect_sizes();
+
+ // In order to perform the static condensation efficiently,
+ // we choose to exploit the symmetry of the the system matrix.
+ // The global system matrix has the following structure
+ // | K_con | K_up | 0 | | dU_u | | R_u |
+ // K = | K_pu | 0 | K_pJ^-1 | , dU = | dU_p | , R = | R_p |
+ // | 0 | K_Jp | K_JJ | | dU_J | | R_J |
+ // We optimise the sparsity pattern to reflect this structure
+ // and prevent unnecessary data creation for the right-diagonal
+ // block components.
+ Table<2, DoFTools::Coupling> coupling(n_components, n_components);
+ for (unsigned int ii = 0; ii < n_components; ++ii)
{
- 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<int dim>
-void Solid<dim>::determine_component_extractors(void) {
- element_indices_u.clear();
- element_indices_p.clear();
- element_indices_J.clear();
-
- for (unsigned int k = 0; k < fe.dofs_per_cell; ++k) {
- // The next call has the FE System indicate to which block component
- // the current DOF is attached to.
- // Currently, the interpolation fields are setup such that
- // 0 indicates a displacement DOF, 1 a pressure DOF and 2 a dilatation DOF.
- const unsigned int k_group = fe.system_to_base_index(k).first.first;
- if (k_group == u_dof) {
- element_indices_u.push_back(k);
- } else if (k_group == p_dof) {
- element_indices_p.push_back(k);
- } else if (k_group == J_dof) {
- element_indices_J.push_back(k);
+ 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 <int dim>
+ void
+ Solid<dim>::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<int dim>
-void Solid<dim>::setup_qph(void) {
- std::cout << " Setting up quadrature point data..." << std::endl;
-
- // Firstly the actual QPH data objects are created. This must be done
- // only once the grid is refined to its finest level.
+ template <int dim>
+ void Solid<dim>::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<PointHistory<dim> > 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<dim>::active_cell_iterator cell =
+ triangulation.begin_active(); cell != triangulation.end();
+ ++cell)
{
- triangulation.clear_user_data();
- {
- std::vector<PointHistory<dim> > 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<dim>::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<dim>::active_cell_iterator cell =
- triangulation.begin_active(); cell != triangulation.end(); ++cell) {
- PointHistory<dim>* lqph =
- reinterpret_cast<PointHistory<dim>*>(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<dim>::active_cell_iterator cell =
+ triangulation.begin_active(); cell != triangulation.end(); ++cell)
+ {
+ PointHistory<dim>* lqph =
+ reinterpret_cast<PointHistory<dim>*>(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<int dim>
-void Solid<dim>::update_qph_incremental(
- const BlockVector<double> & solution_delta) {
- timer.enter_subsection("Update QPH data");
- std::cout << " UQPH " << std::flush;
-
- // Firstly we need to obtain the total solution as it stands
- // at this Newton increment
- const BlockVector<double> solution_total(
- get_solution_total(solution_delta));
-
- // Next we create the initial copy of TBB objects
- const UpdateFlags uf_UQPH(update_values | update_gradients);
- PerTaskData_UQPH per_task_data_UQPH;
- ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total);
-
- // and pass them and the one-cell update function to the WorkStream to be processed
- WorkStream::run(dof_handler_ref.begin_active(), dof_handler_ref.end(),
- *this, &Solid::update_qph_incremental_one_cell,
- &Solid::copy_local_to_global_UQPH, scratch_data_UQPH,
- per_task_data_UQPH);
-
- timer.leave_subsection();
-}
+ template <int dim>
+ void Solid<dim>::update_qph_incremental(const BlockVector<double> & solution_delta)
+ {
+ timer.enter_subsection("Update QPH data");
+ std::cout << " UQPH " << std::flush;
+
+ // Firstly we need to obtain the total
+ // solution as it stands at this Newton
+ // increment
+ const BlockVector<double> solution_total(
+ get_solution_total(solution_delta));
+
+ // Next we create the initial copy of TBB
+ // objects
+ const UpdateFlags uf_UQPH(update_values | update_gradients);
+ PerTaskData_UQPH per_task_data_UQPH;
+ ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total);
+
+ // and pass them and the one-cell update
+ // function to the WorkStream to be
+ // processed
+ WorkStream::run(dof_handler_ref.begin_active(),
+ dof_handler_ref.end(),
+ *this,
+ &Solid::update_qph_incremental_one_cell,
+ &Solid::copy_local_to_global_UQPH,
+ scratch_data_UQPH,
+ per_task_data_UQPH);
+
+ timer.leave_subsection();
+ }
// Now we describe how we extract data from the solution vector and pass it
// along to each QP storage object for processing.
-template<int dim>
-void Solid<dim>::update_qph_incremental_one_cell(
- const typename DoFHandler<dim>::active_cell_iterator & cell,
- ScratchData_UQPH & scratch, PerTaskData_UQPH & data) {
- PointHistory<dim>* lqph =
- reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
-
- Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
- Assert(lqph <= &quadrature_point_history.back(), ExcInternalError());
-
- Assert(scratch.solution_grads_u_total.size() == n_q_points,
- ExcInternalError());
- Assert(scratch.solution_values_p_total.size() == n_q_points,
- ExcInternalError());
- Assert(scratch.solution_values_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 <int dim>
+ void
+ Solid<dim>::update_qph_incremental_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,
+ ScratchData_UQPH & scratch,
+ PerTaskData_UQPH & data)
+ {
+ PointHistory<dim>* lqph =
+ reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+
+ Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
+ Assert(lqph <= &quadrature_point_history.back(), ExcInternalError());
+
+ Assert(scratch.solution_grads_u_total.size() == n_q_points,
+ ExcInternalError());
+ Assert(scratch.solution_values_p_total.size() == n_q_points,
+ ExcInternalError());
+ Assert(scratch.solution_values_J_total.size() == n_q_points,
+ ExcInternalError());
+
+ 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<int dim>
-void Solid<dim>::solve_nonlinear_timestep(
- BlockVector<double> & 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<double> newton_update(dofs_per_block);
- newton_update.collect_sizes();
-
- // Reset the error storage objects
- error_residual.reset();
- error_residual_0.reset();
- error_residual_norm.reset();
- error_update.reset();
- error_update_0.reset();
- error_update_norm.reset();
-
- // Print solver header
- print_conv_header();
-
- // We now perform a number of Newton iterations to iteratively solve
- // the nonlinear problem.
- for (unsigned int it_nr = 0; it_nr < parameters.max_iterations_NR;
- ++it_nr) {
- // Print Newton iteration
- std::cout << " " << std::setw(2) << it_nr << " " << std::flush;
-
- // Since the problem is fully nonlinear and we are using a
- // full Newton method, the data stored in the tangent matrix
- // and right-hand side vector is not reusable and must be cleared
- // at each Newton step.
- tangent_matrix = 0.0;
- system_rhs = 0.0;
-
- // We initially build the right-hand side vector to check for convergence.
- // The unconstrained DOF's of the rhs vector hold the out-of-balance
- // forces. The building is done before assembling the system matrix as the latter
- // is an expensive operation and we can potentially avoid an extra
- // assembly process by not assembling the tangent matrix when convergence
- // is attained.
- assemble_system_rhs();
- 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<unsigned int, double> lin_solver_output =
- solve_linear_system(newton_update);
-
- get_error_update(newton_update, error_update);
- if (it_nr == 0)
- error_update_0 = error_update;
-
- // We can now determine the normalised Newton update error
- error_update_norm = error_update;
- error_update_norm.normalise(error_update_0);
-
- // The current solution state is unacceptable, so we need to update
- // the solution increment for this time step, update all quadrature
- // point information pertaining to this new displacement and stress state
- // and continue iterating.
- solution_delta += newton_update;
- update_qph_incremental(solution_delta);
-
- std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
- << std::scientific << lin_solver_output.first << " "
- << lin_solver_output.second << " " << error_residual_norm.norm
- << " " << error_residual_norm.u << " "
- << error_residual_norm.p << " " << error_residual_norm.J
- << " " << error_update_norm.norm << " " << error_update_norm.u
- << " " << error_update_norm.p << " " << error_update_norm.J
- << " " << std::endl;
- }
+ template <int dim>
+ void
+ Solid<dim>::solve_nonlinear_timestep(BlockVector<double> & 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<double> newton_update(dofs_per_block);
+ newton_update.collect_sizes();
+
+ // Reset the error storage objects
+ error_residual.reset();
+ error_residual_0.reset();
+ error_residual_norm.reset();
+ error_update.reset();
+ error_update_0.reset();
+ error_update_norm.reset();
+
+ // Print solver header
+ print_conv_header();
+
+ // We now perform a number of Newton
+ // iterations to iteratively solve the
+ // nonlinear problem.
+ for (unsigned int it_nr = 0; it_nr < parameters.max_iterations_NR;
+ ++it_nr)
+ {
+ // Print Newton iteration
+ std::cout << " " << std::setw(2) << it_nr << " " << std::flush;
+
+ // Since the problem is fully nonlinear
+ // and we are using a full Newton method,
+ // the data stored in the tangent matrix
+ // and right-hand side vector is not
+ // reusable and must be cleared at each
+ // Newton step.
+ tangent_matrix = 0.0;
+ system_rhs = 0.0;
- 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<unsigned int, double>
+ lin_solver_output = solve_linear_system(newton_update);
+
+ get_error_update(newton_update, error_update);
+ if (it_nr == 0)
+ error_update_0 = error_update;
+
+ // We can now determine the normalised
+ // Newton update error
+ error_update_norm = error_update;
+ error_update_norm.normalise(error_update_0);
+
+ // The current solution state is
+ // unacceptable, so we need to update the
+ // solution increment for this time step,
+ // update all quadrature point
+ // information pertaining to this new
+ // displacement and stress state and
+ // continue iterating.
+ solution_delta += newton_update;
+ update_qph_incremental(solution_delta);
+
+ std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
+ << std::scientific << lin_solver_output.first << " "
+ << lin_solver_output.second << " " << error_residual_norm.norm
+ << " " << error_residual_norm.u << " "
+ << error_residual_norm.p << " " << error_residual_norm.J
+ << " " << error_update_norm.norm << " " << error_update_norm.u
+ << " " << error_update_norm.p << " " << error_update_norm.J
+ << " " << std::endl;
+ }
+
+ throw ExcMessage("No convergence in nonlinear solver!");
+ }
// We print out data in a nice table that is updated
// on a per-iteration basis. Here we set up the table
// header
-template<int dim>
-void Solid<dim>::print_conv_header(void) {
- static const unsigned int l_width = 155;
-
- for (unsigned int i = 0; i < l_width; ++i)
- std::cout << "_";
- std::cout << std::endl;
-
- std::cout << " " << "SOLVER STEP" << " "
- << " | " << " LIN_IT " << " LIN_RES " << " RES_NORM "
- << " RES_U " << " RES_P " << " RES_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<int dim>
-void Solid<dim>::print_conv_footer(void) {
- static const unsigned int l_width = 155;
+ template <int dim>
+ void Solid<dim>::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 <double,double> 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 <int dim>
+ void Solid<dim>::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 <double,double> 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<int dim>
-std::pair<double, double> Solid<dim>::get_error_dil(void) {
+ template <int dim>
+ std::pair<double, double> Solid<dim>::get_error_dil()
+ {
+ double dil_L2_error = 0.0;
+ vol_current = 0.0;
- double dil_L2_error = 0.0;
- vol_current = 0.0;
+ FEValues<dim> fe_values_ref(fe, qf_cell, update_JxW_values);
- FEValues<dim> fe_values_ref(fe, qf_cell, update_JxW_values);
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ {
+ fe_values_ref.reinit(cell);
- for (typename Triangulation<dim>::active_cell_iterator cell =
- triangulation.begin_active(); cell != triangulation.end(); ++cell) {
- fe_values_ref.reinit(cell);
-
- PointHistory<dim>* lqph =
- reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+ PointHistory<dim>* lqph =
+ reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
- Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
- Assert(lqph <= &quadrature_point_history.back(), ExcInternalError());
+ Assert(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<double, double> error_dil;
+ error_dil.first = std::sqrt(dil_L2_error);
+ error_dil.second = vol_current / vol_reference;
- std::pair<double, double> 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<int dim>
-void Solid<dim>::get_error_residual(Errors & error_residual) {
- BlockVector<double> error_res(dofs_per_block);
- error_res.collect_sizes();
-
- // Need to ignore constrained DOFs
- for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
- if (!constraints.is_constrained(i))
- error_res(i) = system_rhs(i);
-
- error_residual.norm = error_res.l2_norm();
- error_residual.u = error_res.block(u_dof).l2_norm();
- error_residual.p = error_res.block(p_dof).l2_norm();
- error_residual.J = error_res.block(J_dof).l2_norm();
-}
+ template <int dim>
+ void Solid<dim>::get_error_residual(Errors & error_residual)
+ {
+ BlockVector<double> error_res(dofs_per_block);
+ error_res.collect_sizes();
+
+ // Need to ignore constrained DOFs
+ for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
+ if (!constraints.is_constrained(i))
+ error_res(i) = system_rhs(i);
+
+ error_residual.norm = error_res.l2_norm();
+ error_residual.u = error_res.block(u_dof).l2_norm();
+ error_residual.p = error_res.block(p_dof).l2_norm();
+ error_residual.J = error_res.block(J_dof).l2_norm();
+ }
+
// Determine the true Newton update error for the problem
-template<int dim>
-void Solid<dim>::get_error_update(const BlockVector<double> & newton_update,
- Errors & error_update) {
- BlockVector<double> error_ud(dofs_per_block);
- error_ud.collect_sizes();
-
- // Need to ignore constrained DOFs as they have a prescribed
- // value
- for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
- if (!constraints.is_constrained(i))
- error_ud(i) = newton_update(i);
-
- error_update.norm = error_ud.l2_norm();
- error_update.u = error_ud.block(u_dof).l2_norm();
- error_update.p = error_ud.block(p_dof).l2_norm();
- error_update.J = error_ud.block(J_dof).l2_norm();
-}
+ template <int dim>
+ void Solid<dim>::get_error_update(const BlockVector<double> & newton_update,
+ Errors & error_update)
+ {
+ BlockVector<double> error_ud(dofs_per_block);
+ error_ud.collect_sizes();
+
+ // Need to ignore constrained DOFs as they
+ // have a prescribed value
+ for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
+ if (!constraints.is_constrained(i))
+ error_ud(i) = newton_update(i);
+
+ error_update.norm = error_ud.l2_norm();
+ error_update.u = error_ud.block(u_dof).l2_norm();
+ error_update.p = error_ud.block(p_dof).l2_norm();
+ error_update.J = error_ud.block(J_dof).l2_norm();
+ }
// This function provides the total solution, which is valid at any Newton step.
// This is required as, to reduce computational error, the total solution is
// only updated at the end of the timestep.
-template<int dim>
-BlockVector<double> Solid<dim>::get_solution_total(
- const BlockVector<double> & solution_delta) const {
- BlockVector<double> solution_total(solution_n);
- solution_total += solution_delta;
- return solution_total;
-
-}
+ template <int dim>
+ BlockVector<double>
+ Solid<dim>::get_solution_total(const BlockVector<double> & solution_delta) const
+ {
+ BlockVector<double> 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
// with the memory addresses of the assembly functions to the
// WorkStream object for processing. Note that we must ensure that
// the matrix is reset before any assembly operations can occur.
-template<int dim>
-void Solid<dim>::assemble_system_tangent(void) {
- timer.enter_subsection("Assemble tangent matrix");
- std::cout << " ASM_K " << std::flush;
+ template <int dim>
+ void Solid<dim>::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<int dim>
-void Solid<dim>::copy_local_to_global_K(const PerTaskData_K & data) {
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
- tangent_matrix.add(data.local_dof_indices[i],
- data.local_dof_indices[j], data.cell_matrix(i, j));
-}
-
-// Here we define how we assemble the tangent matrix contribution for a
-// single cell.
-template<int dim>
-void Solid<dim>::assemble_system_tangent_one_cell(
- const typename DoFHandler<dim>::active_cell_iterator & cell,
- ScratchData_K & scratch, PerTaskData_K & data) {
- // We first need to reset and initialise some
- // of the data structures and retrieve some
- // basic information regarding the DOF numbering on this cell
- data.reset();
- scratch.reset();
- scratch.fe_values_ref.reinit(cell);
- cell->get_dof_indices(data.local_dof_indices);
- PointHistory<dim> *lqph =
- reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
-
- // We can precalculate the cell shape function values and gradients. Note that the
- // shape function gradients are defined wrt the current configuration.
- // That is
- // $\textrm{grad}\boldsymbol{\varphi} = \textrm{Grad}\boldsymbol{\varphi} \mathbf{F}^{-1}$
- static const SymmetricTensor<2, dim> I = AdditionalTools::StandardTensors<
- dim>::I;
- for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
- const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
- for (unsigned int k = 0; k < dofs_per_cell; ++k) {
- const unsigned int k_group = fe.system_to_base_index(k).first.first;
-
- if (k_group == u_dof) {
- scratch.grad_Nx[q_point][k] =
- scratch.fe_values_ref[u_fe].gradient(k, q_point)
- * F_inv;
- scratch.symm_grad_Nx[q_point][k] = symmetrize(
- scratch.grad_Nx[q_point][k]);
- } else if (k_group == p_dof) {
- scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
- q_point);
- } else if (k_group == J_dof) {
- scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
- q_point);
- } else {
- Assert(k_group <= J_dof, ExcInternalError());
- }
- }
- }
+ template <int dim>
+ void Solid<dim>::copy_local_to_global_K(const PerTaskData_K & data)
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ tangent_matrix.add(data.local_dof_indices[i],
+ data.local_dof_indices[j],
+ data.cell_matrix(i, j));
+ }
+
+// Here we define how we assemble the tangent matrix contribution for a single
+// cell.
+ template <int dim>
+ void
+ Solid<dim>::assemble_system_tangent_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,
+ ScratchData_K & scratch,
+ PerTaskData_K & data)
+ {
+ // We first need to reset and initialise
+ // some of the data structures and retrieve
+ // some basic information regarding the DOF
+ // numbering on this cell
+ data.reset();
+ scratch.reset();
+ scratch.fe_values_ref.reinit(cell);
+ cell->get_dof_indices(data.local_dof_indices);
+ PointHistory<dim> *lqph =
+ reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+
+ // We can precalculate the cell shape
+ // function values and gradients. Note that
+ // the shape function gradients are defined
+ // wrt the current configuration. That is
+ // $\textrm{grad}\boldsymbol{\varphi} =
+ // \textrm{Grad}\boldsymbol{\varphi}
+ // \mathbf{F}^{-1}$
+ static const SymmetricTensor<2, dim> I = AdditionalTools::StandardTensors<dim>::I;
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
+ for (unsigned int k = 0; k < dofs_per_cell; ++k)
+ {
+ const unsigned int k_group = fe.system_to_base_index(k).first.first;
+
+ if (k_group == u_dof)
+ {
+ scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point)
+ * F_inv;
+ scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
+ }
+ else if (k_group == p_dof)
+ {
+ scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
+ q_point);
+ }
+ else if (k_group == J_dof)
+ {
+ scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
+ q_point);
+ }
+ else
+ {
+ Assert(k_group <= J_dof, ExcInternalError());
+ }
+ }
+ }
+
+ // 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<double>
+ & N = scratch.Nx[q_point];
+ const std::vector<SymmetricTensor<2, dim> >
+ & symm_grad_Nx = scratch.symm_grad_Nx[q_point];
+ const std::vector<Tensor<2, dim> >
+ & grad_Nx = scratch.grad_Nx[q_point];
+ const double JxW = scratch.fe_values_ref.JxW(q_point);
- // Now we build the local cell stiffness matrix. Since the global and local system
- // matrices are symmetric, we can exploit this property by building only the lower
- // half of the local matrix and copying the values to the upper half.
- // So we only assemble half of the K_uu, K_pp (= 0), K_JJ blocks, while the whole
- // K_pJ, K_uJ (=0), K_up blocks are built.
- for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
- // We first extract some configuration dependent variables from our
- // QPH history objects that for the current q_point.
- // Get the current stress state $\boldsymbol{\tau}$
- const Tensor<2, dim> tau =
- static_cast<Tensor<2, dim> >(lqph[q_point].get_tau());
- const SymmetricTensor<4, dim> Jc = lqph[q_point].get_Jc();
- const double d2Psi_vol_dJ2 = lqph[q_point].get_d2Psi_vol_dJ2();
- const double det_F = lqph[q_point].get_det_F();
-
- // Next we define some aliases to make the assembly process easier to follow
- const std::vector<double> & N = scratch.Nx[q_point];
- const std::vector<SymmetricTensor<2, dim> > & symm_grad_Nx =
- scratch.symm_grad_Nx[q_point];
- const std::vector<Tensor<2, dim> > & grad_Nx = scratch.grad_Nx[q_point];
- const double JxW = scratch.fe_values_ref.JxW(q_point);
-
- for (unsigned int i = 0; i < dofs_per_cell; ++i) {
- const unsigned int component_i =
- fe.system_to_component_index(i).first;
- // Determine the dimensional component that matches the dof component (i.e. i % dim)
- const unsigned int i_group = fe.system_to_base_index(i).first.first;
-
- for (unsigned int j = 0; j <= i; ++j) {
- const unsigned int component_j =
- fe.system_to_component_index(j).first;
- const unsigned int j_group =
- fe.system_to_base_index(j).first.first;
-
- // This is the K_{uu} contribution. It comprises of a material
- // contribution and a geometrical stress contribution which is only
- // added along the local matrix diagonals
- if ((i_group == j_group) && (i_group == u_dof)) {
- // The material contribution:
- data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc
- * symm_grad_Nx[j] * JxW;
- if (component_i == component_j) // geometrical stress contribution
- data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ const unsigned int component_i = fe.system_to_component_index(i).first;
+ // Determine the dimensional component
+ // that matches the dof component
+ // (i.e. i % dim)
+ const unsigned int i_group = fe.system_to_base_index(i).first.first;
+
+ for (unsigned int j = 0; j <= i; ++j)
+ {
+ const unsigned int component_j = fe.system_to_component_index(j).first;
+ const unsigned int j_group = fe.system_to_base_index(j).first.first;
+
+ // This is the K_{uu}
+ // contribution. It comprises of a
+ // material contribution and a
+ // geometrical stress contribution
+ // which is only added along the
+ // local matrix diagonals
+ if ((i_group == j_group) && (i_group == u_dof))
+ {
+ // The material contribution:
+ data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc
+ * symm_grad_Nx[j] * JxW;
+ if (component_i == component_j) // geometrical stress contribution
+ data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau
* grad_Nx[j][component_j] * JxW;
- }
- // Next is the K_{pu} contribution
- else if ((i_group == p_dof) && (j_group == u_dof)) {
- data.cell_matrix(i, j) += N[i] * det_F
- * (symm_grad_Nx[j]
- * AdditionalTools::StandardTensors<dim>::I)
- * JxW;
- }
- // and the K_{Jp} contribution
- else if ((i_group == J_dof) && (j_group == p_dof)) {
- data.cell_matrix(i, j) -= N[i] * N[j] * JxW;
- }
- // and lastly the K_{JJ} contribution
- else if ((i_group == j_group) && (i_group == J_dof)) {
- data.cell_matrix(i, j) += N[i] * d2Psi_vol_dJ2 * N[j] * JxW;
- } else
- Assert((i_group <= J_dof) && (j_group <= J_dof),
- ExcInternalError());
- }
- }
- }
-
- // 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<dim>::I)
+ * JxW;
+ }
+ // and the K_{Jp} contribution
+ else if ((i_group == J_dof) && (j_group == p_dof))
+ {
+ data.cell_matrix(i, j) -= N[i] * N[j] * JxW;
+ }
+ // and lastly the K_{JJ} contribution
+ else if ((i_group == j_group) && (i_group == J_dof))
+ {
+ data.cell_matrix(i, j) += N[i] * d2Psi_vol_dJ2 * N[j] * JxW;
+ }
+ else
+ Assert((i_group <= J_dof) && (j_group <= J_dof),
+ ExcInternalError());
+ }
+ }
+ }
+
+ // 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
// Note that since we are describing a problem with Neumann BCs,
// we will need the face normals and so must specify this in the
// update flags.
-template<int dim>
-void Solid<dim>::assemble_system_rhs(void) {
- timer.enter_subsection("Assemble system right-hand side");
- std::cout << " ASM_R " << std::flush;
-
- system_rhs = 0.0;
-
- const UpdateFlags uf_cell(
- update_values | update_gradients | update_JxW_values);
- const UpdateFlags uf_face(
- update_values | update_normal_vectors | update_JxW_values);
-
- PerTaskData_RHS per_task_data(dofs_per_cell);
- ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face);
-
- WorkStream::run(dof_handler_ref.begin_active(), dof_handler_ref.end(),
- *this, &Solid::assemble_system_rhs_one_cell,
- &Solid::copy_local_to_global_rhs, scratch_data, per_task_data);
-
- timer.leave_subsection();
-}
-
-template<int dim>
-void Solid<dim>::copy_local_to_global_rhs(const PerTaskData_RHS & data) {
- for (unsigned int i = 0; i < dofs_per_cell; ++i) {
- system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);
- }
-}
-
-template<int dim>
-void Solid<dim>::assemble_system_rhs_one_cell(
- const typename DoFHandler<dim>::active_cell_iterator & cell,
- ScratchData_RHS & scratch, PerTaskData_RHS & data) {
- // Again we reset the data structures
- data.reset();
- scratch.reset();
- scratch.fe_values_ref.reinit(cell);
- cell->get_dof_indices(data.local_dof_indices);
- PointHistory<dim> *lqph =
- reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
-
- // and then precompute some shape function data
- for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
- const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
-
- for (unsigned int k = 0; k < dofs_per_cell; ++k) {
- const unsigned int k_group = fe.system_to_base_index(k).first.first;
-
- if (k_group == u_dof) {
- scratch.symm_grad_Nx[q_point][k] = symmetrize(
- scratch.fe_values_ref[u_fe].gradient(k, q_point)
- * F_inv);
- } else if (k_group == p_dof) {
- scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
- q_point);
- } else if (k_group == J_dof) {
- scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
- q_point);
- } else
- Assert(k_group <= J_dof, ExcInternalError());
- }
- }
-
- // and can now assemble the right-hand side contribution
- for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
- // We fist retrieve data stored at the qp
- const SymmetricTensor<2, dim> tau = lqph[q_point].get_tau();
- const double det_F = lqph[q_point].get_det_F();
- const double J_tilde = lqph[q_point].get_J_tilde();
- const double p_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<double> & N = scratch.Nx[q_point];
- const std::vector<SymmetricTensor<2, dim> > & symm_grad_Nx =
- scratch.symm_grad_Nx[q_point];
- const double JxW = scratch.fe_values_ref.JxW(q_point);
-
- // We first compute the contributions from the internal forces.
- // Note, by definition of the rhs as the negative of the residual,
- // these contributions are subtracted.
- for (unsigned int i = 0; i < dofs_per_cell; ++i) {
- const unsigned int i_group = fe.system_to_base_index(i).first.first;
- // Add the contribution to the F_u block
- if (i_group == u_dof) {
- 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<dim>::faces_per_cell;
- ++face) {
- if (cell->face(face)->at_boundary() == true
- && cell->face(face)->boundary_indicator() == 6) {
- scratch.fe_face_values_ref.reinit(cell, face);
-
- for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
- ++f_q_point) {
- // We retrieve the face normal at this QP
- const Tensor<1, dim> & N =
- scratch.fe_face_values_ref.normal_vector(f_q_point);
-
- // and specify the traction in reference configuration. For this problem,
- // a defined pressure is applied in the reference configuration.
- // The direction of the applied traction is assumed
- // not to evolve with the deformation of the domain. The
- // traction is defined using the first Piola-Kirchhoff stress is simply
- // t_0 = P*N = (pI)*N = p*N
- // We choose to use the time variable to linearly ramp up the pressure
- // load.
- static const double p0 = -4.0
- / (parameters.scale * parameters.scale);
- const double time_ramp = (time.current() / time.end());
- const double pressure = p0 * parameters.p_p0 * time_ramp;
- const Tensor<1, dim> traction = pressure * N;
-
- for (unsigned int i = 0; i < dofs_per_cell; ++i) {
- const unsigned int i_group =
- fe.system_to_base_index(i).first.first;
-
- if (i_group == u_dof) {
- // More shortcuts being assigned
- const unsigned int component_i =
- fe.system_to_component_index(i).first;
- const double Ni =
- scratch.fe_face_values_ref.shape_value(i,
- f_q_point);
- const double JxW = scratch.fe_face_values_ref.JxW(
- f_q_point);
-
- // And finally we can add the traction vector contribution to
- // the local RHS vector. Note that this contribution is present
- // on displacement DOFs only.
- data.cell_rhs(i) += (Ni * traction[component_i])
- * JxW;
- }
- }
- }
- }
- }
- }
-}
+ template <int dim>
+ void Solid<dim>::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 <int dim>
+ void Solid<dim>::copy_local_to_global_rhs(const PerTaskData_RHS & data)
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);
+ }
+ }
+
+
+
+ template <int dim>
+ void
+ Solid<dim>::assemble_system_rhs_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,
+ ScratchData_RHS & scratch,
+ PerTaskData_RHS & data)
+ {
+ // Again we reset the data structures
+ data.reset();
+ scratch.reset();
+ scratch.fe_values_ref.reinit(cell);
+ cell->get_dof_indices(data.local_dof_indices);
+ PointHistory<dim> *lqph =
+ reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+
+ // and then precompute some shape function
+ // data
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
+
+ for (unsigned int k = 0; k < dofs_per_cell; ++k) {
+ const unsigned int k_group = fe.system_to_base_index(k).first.first;
+
+ if (k_group == u_dof)
+ {
+ scratch.symm_grad_Nx[q_point][k]
+ = symmetrize(scratch.fe_values_ref[u_fe].gradient(k, q_point)
+ * F_inv);
+ }
+ else if (k_group == p_dof)
+ {
+ scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
+ q_point);
+ }
+ else if (k_group == J_dof)
+ {
+ scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
+ q_point);
+ }
+ else
+ Assert(k_group <= J_dof, ExcInternalError());
+ }
+ }
+
+ // and can now assemble the right-hand side
+ // contribution
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ // We fist retrieve data stored at the qp
+ const SymmetricTensor<2, dim> tau = lqph[q_point].get_tau();
+ const double det_F = lqph[q_point].get_det_F();
+ const double J_tilde = lqph[q_point].get_J_tilde();
+ const double p_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<double>
+ & N = scratch.Nx[q_point];
+ const std::vector<SymmetricTensor<2, dim> >
+ & symm_grad_Nx = scratch.symm_grad_Nx[q_point];
+ const double JxW = scratch.fe_values_ref.JxW(q_point);
+
+ // We first compute the contributions
+ // from the internal forces. Note, by
+ // definition of the rhs as the negative
+ // of the residual, these contributions
+ // are subtracted.
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ const unsigned int i_group = fe.system_to_base_index(i).first.first;
+ // Add the contribution to the F_u
+ // block
+ if (i_group == u_dof)
+ {
+ 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<dim>::faces_per_cell;
+ ++face)
+ {
+ if (cell->face(face)->at_boundary() == true
+ && cell->face(face)->boundary_indicator() == 6)
+ {
+ scratch.fe_face_values_ref.reinit(cell, face);
+
+ for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
+ ++f_q_point)
+ {
+ // We retrieve the face normal at
+ // this QP
+ const Tensor<1, dim> & N =
+ scratch.fe_face_values_ref.normal_vector(f_q_point);
+
+ // and specify the traction in
+ // reference configuration. For
+ // this problem, a defined pressure
+ // is applied in the reference
+ // configuration. The direction of
+ // the applied traction is assumed
+ // not to evolve with the
+ // deformation of the domain. The
+ // traction is defined using the
+ // first Piola-Kirchhoff stress is
+ // simply t_0 = P*N = (pI)*N = p*N
+ // We choose to use the time
+ // variable to linearly ramp up the
+ // pressure load.
+ static const double p0 = -4.0
+ /
+ (parameters.scale * parameters.scale);
+ const double time_ramp = (time.current() / time.end());
+ const double pressure = p0 * parameters.p_p0 * time_ramp;
+ const Tensor<1, dim> traction = pressure * N;
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ const unsigned int i_group =
+ fe.system_to_base_index(i).first.first;
+
+ if (i_group == u_dof)
+ {
+ // More shortcuts being assigned
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+ const double Ni =
+ scratch.fe_face_values_ref.shape_value(i,
+ f_q_point);
+ const double JxW = scratch.fe_face_values_ref.JxW(
+ f_q_point);
+
+ // And finally we can add the
+ // traction vector contribution
+ // to the local RHS
+ // vector. Note that this
+ // contribution is present on
+ // displacement DOFs only.
+ data.cell_rhs(i) += (Ni * traction[component_i])
+ * JxW;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
// @sect4{Solid::make_constraints}
// The constraints for this problem are simple to describe.
// completeness although for this problem the constraints are
// trivial and it would not have made a difference if this had
// not been accounted for in this problem.
-template<int dim>
-void Solid<dim>::make_constraints(const int & it_nr,
- ConstraintMatrix & constraints) {
- std::cout << " CST " << std::flush;
-
- // Since the constraints are different at Newton iterations,
- // we need to clear the constraints matrix and completely
- // rebuild it. However, after the first iteration, the
- // constraints remain the same and we can simply skip the
- // rebuilding step if we do not clear it.
- if (it_nr > 1)
- return;
- constraints.clear();
- const bool apply_dirichlet_bc = (it_nr == 0);
-
- // The boundary conditions for the indentation problem are as follows:
- // On the -x, -y and -z faces (ID's 0,2,4) we set up a symmetry condition
- // to allow only planar movement while the +x and +y faces (ID's 1,3) are
- // traction free. In this contrived problem, part of the +z face (ID 5) is
- // set to have no motion in the x- and y-component. Finally, as described
- // earlier, the other part of the +z face has an the applied pressure but
- // is also constrained in the x- and y-directions.
+ template <int dim>
+ void Solid<dim>::make_constraints(const int & it_nr,
+ ConstraintMatrix & constraints)
+ {
+ std::cout << " CST " << std::flush;
+
+ // Since the constraints are different at
+ // Newton iterations, we need to clear the
+ // constraints matrix and completely
+ // rebuild it. However, after the first
+ // iteration, the constraints remain the
+ // same and we can simply skip the
+ // rebuilding step if we do not clear it.
+ if (it_nr > 1)
+ return;
+ constraints.clear();
+ const bool apply_dirichlet_bc = (it_nr == 0);
+
+ // The boundary conditions for the
+ // indentation problem are as follows: On
+ // the -x, -y and -z faces (ID's 0,2,4) we
+ // set up a symmetry condition to allow
+ // only planar movement while the +x and +y
+ // faces (ID's 1,3) are traction free. In
+ // this contrived problem, part of the +z
+ // face (ID 5) is set to have no motion in
+ // the x- and y-component. Finally, as
+ // described earlier, the other part of the
+ // +z face has an the applied pressure but
+ // is also constrained in the x- and
+ // y-directions.
+ {
+ const int boundary_id = 0;
+
+ std::vector<bool> components(n_components, false);
+ components[0] = true;
+
+ if (apply_dirichlet_bc == true)
{
- const int boundary_id = 0;
-
- std::vector<bool> components(n_components, false);
- components[0] = true;
-
- if (apply_dirichlet_bc == true) {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- } else {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- }
- }
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ } else {
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ }
+ }
+ {
+ const int boundary_id = 2;
+
+ std::vector<bool> components(n_components, false);
+ components[1] = true;
+
+ if (apply_dirichlet_bc == true)
{
- const int boundary_id = 2;
-
- std::vector<bool> components(n_components, false);
- components[1] = true;
-
- if (apply_dirichlet_bc == true) {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- } else {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- }
- }
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ } else {
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ }
+ }
+ {
+ const int boundary_id = 4;
+ std::vector<bool> components(n_components, false);
+ components[2] = true;
+
+ if (apply_dirichlet_bc == true)
{
- const int boundary_id = 4;
- std::vector<bool> components(n_components, false);
- components[2] = true;
-
- if (apply_dirichlet_bc == true) {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- } else {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- }
- }
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ } else {
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ }
+ }
+ {
+ const int boundary_id = 5;
+ std::vector<bool> components(n_components, true);
+ components[2] = false;
+
+ if (apply_dirichlet_bc == true)
{
- const int boundary_id = 5;
- std::vector<bool> components(n_components, true);
- components[2] = false;
-
- if (apply_dirichlet_bc == true) {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- } else {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- }
- }
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ } else {
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ }
+ }
+ {
+ const int boundary_id = 6;
+ std::vector<bool> components(n_components, true);
+ components[2] = false;
+
+ if (apply_dirichlet_bc == true)
{
- const int boundary_id = 6;
- std::vector<bool> components(n_components, true);
- components[2] = false;
-
- if (apply_dirichlet_bc == true) {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- } else {
- VectorTools::interpolate_boundary_values(dof_handler_ref,
- boundary_id, ZeroFunction<dim>(n_components), constraints,
- components);
- }
- }
-
- constraints.close();
-}
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ } else {
+ VectorTools::interpolate_boundary_values(dof_handler_ref,
+ boundary_id,
+ ZeroFunction<dim>(n_components),
+ constraints,
+ components);
+ }
+ }
+
+ constraints.close();
+ }
// @sect4{Solid::solve_linear_system}
// Solving the entire block system is a bit problematic as there are no
// condense them out to form a smaller displacement-only system which
// we will then solve and subsequently post-process to retrieve the
// pressure and dilatation solutions.
-template<int dim>
-std::pair<unsigned int, double> Solid<dim>::solve_linear_system(
- BlockVector<double> & newton_update) {
- // Need two temporary vectors to help
- // with the static condensation.
- BlockVector<double> A(dofs_per_block);
- BlockVector<double> B(dofs_per_block);
- A.collect_sizes();
- B.collect_sizes();
-
- // Store the number of linear solver iterations
- // 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 <int dim>
+ std::pair<unsigned int, double>
+ Solid<dim>::solve_linear_system(BlockVector<double> & newton_update)
+ {
+ // Need two temporary vectors to help with
+ // the static condensation.
+ BlockVector<double> A(dofs_per_block);
+ BlockVector<double> B(dofs_per_block);
+ A.collect_sizes();
+ B.collect_sizes();
+
+ // Store the number of linear solver
+ // iterations 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<Vector<double> > GVM;
- SolverCG<Vector<double> > 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<SparseMatrix<double>, Vector<double> > 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<Vector<double> > GVM;
+ SolverCG<Vector<double> > 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<SparseMatrix<double>, Vector<double> >
+ 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
// block-diagonal K_{pt} block by inverting the local blocks. We can
// again use TBB to do this since each operation will be independent of
// one another.
-template<int dim>
-void Solid<dim>::assemble_sc(void) {
- timer.enter_subsection("Perform static condensation");
- std::cout << " ASM_SC " << std::flush;
-
- PerTaskData_SC per_task_data(dofs_per_cell, element_indices_u.size(),
- element_indices_p.size(), element_indices_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 <int dim>
+ void Solid<dim>::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 <int dim>
+ void Solid<dim>::copy_local_to_global_sc(const PerTaskData_SC & data)
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ tangent_matrix.add(data.local_dof_indices[i],
+ data.local_dof_indices[j],
+ data.cell_matrix(i, j));
+ }
-// We need to describe how to add the local contributions
-// to K to form K_store
-template<int dim>
-void Solid<dim>::copy_local_to_global_sc(const PerTaskData_SC & data) {
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
- tangent_matrix.add(data.local_dof_indices[i],
- data.local_dof_indices[j], data.cell_matrix(i, j));
-}
// Now we describe the static condensation process.
-template<int dim>
-void Solid<dim>::assemble_sc_one_cell(
- const typename DoFHandler<dim>::active_cell_iterator & cell,
- ScratchData_SC & scratch, PerTaskData_SC & data) {
- // As per usual, we must first find out which global numbers the
- // degrees of freedom on this cell have and reset some data structures
- data.reset();
- scratch.reset();
- cell->get_dof_indices(data.local_dof_indices);
-
- // 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 <int dim>
+ void
+ Solid<dim>::assemble_sc_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,
+ ScratchData_SC & scratch,
+ PerTaskData_SC & data)
+ {
+ // As per usual, we must first find out
+ // which global numbers the degrees of
+ // freedom on this cell have and reset some
+ // data structures
+ data.reset();
+ scratch.reset();
+ cell->get_dof_indices(data.local_dof_indices);
+
+ // We now extract the contribution of
+ // the dof associated with the current cell
+ // to the global stiffness matrix.
+ // The discontinuous nature of the $\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<int dim>
-void Solid<dim>::output_results(void) const {
- DataOut<dim> data_out;
- std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(
- dim, DataComponentInterpretation::component_is_part_of_vector);
- data_component_interpretation.push_back(
- DataComponentInterpretation::component_is_scalar);
- data_component_interpretation.push_back(
- DataComponentInterpretation::component_is_scalar);
-
- std::vector<std::string> solution_name(dim, "displacement");
- solution_name.push_back("pressure");
- solution_name.push_back("dilatation");
-
- data_out.attach_dof_handler(dof_handler_ref);
- data_out.add_data_vector(solution_n, solution_name,
- DataOut<dim>::type_dof_data, data_component_interpretation);
-
- // Since we are dealing with a large deformation problem, it would be nice
- // to display the result on a displaced grid! The MappingQEulerian class
- // linked with the DataOut class provides an interface through which this
- // can be achieved without physically moving the grid points ourselves.
- // We first need to copy the solution to a temporary vector and then
- // create the 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<double> soln(solution_n.size());
- for (unsigned int i = 0; i < soln.size(); ++i)
- soln(i) = solution_n(i);
- MappingQEulerian<dim> q_mapping(degree, soln, dof_handler_ref);
- data_out.build_patches(q_mapping, degree);
-
- std::ostringstream filename;
- filename << "solution-" << time.get_timestep() << ".vtk";
-
- std::ofstream output(filename.str().c_str());
- data_out.write_vtk(output);
+ template <int dim>
+ void Solid<dim>::output_results() const
+ {
+ DataOut<dim> data_out;
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ data_component_interpretation(dim,
+ DataComponentInterpretation::component_is_part_of_vector);
+ data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
+ data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
+
+ std::vector<std::string> solution_name(dim, "displacement");
+ solution_name.push_back("pressure");
+ solution_name.push_back("dilatation");
+
+ data_out.attach_dof_handler(dof_handler_ref);
+ data_out.add_data_vector(solution_n,
+ solution_name,
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
+
+ // Since we are dealing with a large
+ // deformation problem, it would be nice to
+ // display the result on a displaced grid!
+ // The MappingQEulerian class linked with
+ // the DataOut class provides an interface
+ // through which this can be achieved
+ // without physically moving the grid
+ // points ourselves. We first need to copy
+ // the solution to a temporary vector and
+ // then create the 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<double> soln(solution_n.size());
+ for (unsigned int i = 0; i < soln.size(); ++i)
+ soln(i) = solution_n(i);
+ MappingQEulerian<dim> q_mapping(degree, soln, dof_handler_ref);
+ data_out.build_patches(q_mapping, degree);
+
+ std::ostringstream filename;
+ filename << "solution-" << time.get_timestep() << ".vtk";
+
+ std::ofstream output(filename.str().c_str());
+ data_out.write_vtk(output);
+ }
+
}
+
// @sect3{Main function}
// Lastly we provide the main driver function which appears
// no different to the other tutorials.
-int main(void) {
- try {
- deallog.depth_console(0);
-
- Solid<3> solid_3d("parameters.prm");
- solid_3d.run();
- } catch (std::exception &exc) {
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Exception on processing: " << std::endl << exc.what()
- << std::endl << "Aborting!" << std::endl
- << "----------------------------------------------------"
- << std::endl;
-
- return 1;
- } catch (...) {
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Unknown exception!" << std::endl << "Aborting!"
- << std::endl
- << "----------------------------------------------------"
- << std::endl;
- return 1;
- }
-
- return 0;
+int main (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;
}