From 867e6169f8e4dafd349830909b8996a87f476662 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Fri, 14 Sep 2018 17:31:39 +0200 Subject: [PATCH] Add initial version of test for linesearch in FE problem (step-44) --- .../parameters-step-44-with_linesearch.prm | 81 + ...arameters-step-44-with_linesearch_ptrb.prm | 81 + .../parameters-step-44-without_linesearch.prm | 81 + .../optimization/step-44-with_line_search.cc | 67 + .../step-44-with_line_search.output | 3 + .../step-44-with_line_search_ptrb.cc | 68 + .../step-44-with_line_search_ptrb.output | 3 + .../step-44-without_line_search.cc | 66 + .../step-44-without_line_search.output | 3 + tests/optimization/step-44.h | 2247 +++++++++++++++++ 10 files changed, 2700 insertions(+) create mode 100644 tests/optimization/prm/parameters-step-44-with_linesearch.prm create mode 100644 tests/optimization/prm/parameters-step-44-with_linesearch_ptrb.prm create mode 100644 tests/optimization/prm/parameters-step-44-without_linesearch.prm create mode 100644 tests/optimization/step-44-with_line_search.cc create mode 100644 tests/optimization/step-44-with_line_search.output create mode 100644 tests/optimization/step-44-with_line_search_ptrb.cc create mode 100644 tests/optimization/step-44-with_line_search_ptrb.output create mode 100644 tests/optimization/step-44-without_line_search.cc create mode 100644 tests/optimization/step-44-without_line_search.output create mode 100644 tests/optimization/step-44.h diff --git a/tests/optimization/prm/parameters-step-44-with_linesearch.prm b/tests/optimization/prm/parameters-step-44-with_linesearch.prm new file mode 100644 index 0000000000..79a29e9f44 --- /dev/null +++ b/tests/optimization/prm/parameters-step-44-with_linesearch.prm @@ -0,0 +1,81 @@ +# Listing of Parameters +# --------------------- +subsection Finite element system + # Displacement system polynomial order + set Polynomial degree = 1 + + # Gauss quadrature order + set Quadrature order = 2 +end + + +subsection Geometry + # Global refinement level + set Global refinement = 2 + + # Global grid scaling factor + set Grid scale = 1e-3 + + # Ratio of applied pressure to reference pressure + set Pressure ratio p/p0 = 100 +end + + +subsection Linear solver + # Linear solver iterations (multiples of the system matrix size) + # In 2-d, this value is best set at 2. In 3-d, a value of 1 work fine. + set Max iteration multiplier = 2 + + # Linear solver residual (scaled by residual norm) + set Residual = 1e-6 + + # Use static condensation and solve a 1-block system, or solve + # the full 3-block system using Linear Operators and the Schur + # complement + set Use static condensation = true + + # Preconditioner type + set Preconditioner type = ssor + + # Preconditioner relaxation value + set Preconditioner relaxation = 0.65 + + # Type of solver used to solve the linear system + set Solver type = CG +end + + +subsection Material properties + # Poisson's ratio + set Poisson's ratio = 0.4999 + + # Shear modulus + set Shear modulus = 80.194e6 +end + + +subsection Nonlinear solver + # Number of Newton-Raphson iterations allowed + set Max iterations Newton-Raphson = 10 + + # Displacement error tolerance + set Tolerance displacement = 1.0e-6 + + # Force residual tolerance + set Tolerance force = 1.0e-9 + + # Use line-search minimization algorithm + set Use line search = true + + # Use an approximation for the gradient during line search + set Use line search approximate gradient = false +end + + +subsection Time + # End time + set End time = 1 + + # Time step size + set Time step size = 0.5 +end diff --git a/tests/optimization/prm/parameters-step-44-with_linesearch_ptrb.prm b/tests/optimization/prm/parameters-step-44-with_linesearch_ptrb.prm new file mode 100644 index 0000000000..042f1a8722 --- /dev/null +++ b/tests/optimization/prm/parameters-step-44-with_linesearch_ptrb.prm @@ -0,0 +1,81 @@ +# Listing of Parameters +# --------------------- +subsection Finite element system + # Displacement system polynomial order + set Polynomial degree = 1 + + # Gauss quadrature order + set Quadrature order = 2 +end + + +subsection Geometry + # Global refinement level + set Global refinement = 2 + + # Global grid scaling factor + set Grid scale = 1e-3 + + # Ratio of applied pressure to reference pressure + set Pressure ratio p/p0 = 100 +end + + +subsection Linear solver + # Linear solver iterations (multiples of the system matrix size) + # In 2-d, this value is best set at 2. In 3-d, a value of 1 work fine. + set Max iteration multiplier = 2 + + # Linear solver residual (scaled by residual norm) + set Residual = 1e-6 + + # Use static condensation and solve a 1-block system, or solve + # the full 3-block system using Linear Operators and the Schur + # complement + set Use static condensation = true + + # Preconditioner type + set Preconditioner type = ssor + + # Preconditioner relaxation value + set Preconditioner relaxation = 0.65 + + # Type of solver used to solve the linear system + set Solver type = CG +end + + +subsection Material properties + # Poisson's ratio + set Poisson's ratio = 0.4999 + + # Shear modulus + set Shear modulus = 80.194e6 +end + + +subsection Nonlinear solver + # Number of Newton-Raphson iterations allowed + set Max iterations Newton-Raphson = 10 + + # Displacement error tolerance + set Tolerance displacement = 1.0e-6 + + # Force residual tolerance + set Tolerance force = 1.0e-9 + + # Use line-search minimization algorithm + set Use line search = true + + # Use an approximation for the gradient during line search + set Use line search approximate gradient = true +end + + +subsection Time + # End time + set End time = 1 + + # Time step size + set Time step size = 0.5 +end diff --git a/tests/optimization/prm/parameters-step-44-without_linesearch.prm b/tests/optimization/prm/parameters-step-44-without_linesearch.prm new file mode 100644 index 0000000000..72f908c3c0 --- /dev/null +++ b/tests/optimization/prm/parameters-step-44-without_linesearch.prm @@ -0,0 +1,81 @@ +# Listing of Parameters +# --------------------- +subsection Finite element system + # Displacement system polynomial order + set Polynomial degree = 1 + + # Gauss quadrature order + set Quadrature order = 2 +end + + +subsection Geometry + # Global refinement level + set Global refinement = 2 + + # Global grid scaling factor + set Grid scale = 1e-3 + + # Ratio of applied pressure to reference pressure + set Pressure ratio p/p0 = 100 +end + + +subsection Linear solver + # Linear solver iterations (multiples of the system matrix size) + # In 2-d, this value is best set at 2. In 3-d, a value of 1 work fine. + set Max iteration multiplier = 2 + + # Linear solver residual (scaled by residual norm) + set Residual = 1e-6 + + # Use static condensation and solve a 1-block system, or solve + # the full 3-block system using Linear Operators and the Schur + # complement + set Use static condensation = true + + # Preconditioner type + set Preconditioner type = ssor + + # Preconditioner relaxation value + set Preconditioner relaxation = 0.65 + + # Type of solver used to solve the linear system + set Solver type = CG +end + + +subsection Material properties + # Poisson's ratio + set Poisson's ratio = 0.4999 + + # Shear modulus + set Shear modulus = 80.194e6 +end + + +subsection Nonlinear solver + # Number of Newton-Raphson iterations allowed + set Max iterations Newton-Raphson = 10 + + # Displacement error tolerance + set Tolerance displacement = 1.0e-6 + + # Force residual tolerance + set Tolerance force = 1.0e-9 + + # Use line-search minimization algorithm + set Use line search = false + + # Use an approximation for the gradient during line search + set Use line search approximate gradient = false +end + + +subsection Time + # End time + set End time = 1 + + # Time step size + set Time step size = 0.5 +end diff --git a/tests/optimization/step-44-with_line_search.cc b/tests/optimization/step-44-with_line_search.cc new file mode 100644 index 0000000000..f164b72e36 --- /dev/null +++ b/tests/optimization/step-44-with_line_search.cc @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// A reduced version of step-44 (based off of git rev 3f7e617) that utilizes +// line-search to select the optimal step-size at each Newton iteration. +// The chosen parameters dictate that the exact value of the gradient +// for the function to be minimized is computed. + +#include "../tests.h" +#include "step-44.h" + +int +main(int argc, char **argv) +{ + initlog(true); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + using namespace dealii; + using namespace Step44; + try + { + const unsigned int dim = 3; + Solid solid(SOURCE_DIR + "/prm/parameters-step-44-with_linesearch.prm"); + solid.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; +} diff --git a/tests/optimization/step-44-with_line_search.output b/tests/optimization/step-44-with_line_search.output new file mode 100644 index 0000000000..4c17c7fac4 --- /dev/null +++ b/tests/optimization/step-44-with_line_search.output @@ -0,0 +1,3 @@ + +DEAL::Timestep 1: 0.00000 -0.000509022 0.00000 +DEAL::Timestep 2: 0.00000 -0.000789713 0.00000 diff --git a/tests/optimization/step-44-with_line_search_ptrb.cc b/tests/optimization/step-44-with_line_search_ptrb.cc new file mode 100644 index 0000000000..884c7e4a83 --- /dev/null +++ b/tests/optimization/step-44-with_line_search_ptrb.cc @@ -0,0 +1,68 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// A reduced version of step-44 (based off of git rev 3f7e617) that utilizes +// line-search to select the optimal step-size at each Newton iteration. +// This is the same as optimization/step-44-with_line_search.cc, except that +// a perturbation method is used to compute the gradient of the function to +// be minimized. + +#include "../tests.h" +#include "step-44.h" + +int +main(int argc, char **argv) +{ + initlog(true); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + using namespace dealii; + using namespace Step44; + try + { + const unsigned int dim = 3; + Solid solid(SOURCE_DIR + "/prm/parameters-step-44-with_linesearch_ptrb.prm"); + solid.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; +} diff --git a/tests/optimization/step-44-with_line_search_ptrb.output b/tests/optimization/step-44-with_line_search_ptrb.output new file mode 100644 index 0000000000..4c17c7fac4 --- /dev/null +++ b/tests/optimization/step-44-with_line_search_ptrb.output @@ -0,0 +1,3 @@ + +DEAL::Timestep 1: 0.00000 -0.000509022 0.00000 +DEAL::Timestep 2: 0.00000 -0.000789713 0.00000 diff --git a/tests/optimization/step-44-without_line_search.cc b/tests/optimization/step-44-without_line_search.cc new file mode 100644 index 0000000000..f5437b7054 --- /dev/null +++ b/tests/optimization/step-44-without_line_search.cc @@ -0,0 +1,66 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// A reduced version of step-44 (based off of git rev 3f7e617) that does not +// employ line-search. This serves as a compartive result for +// optimization/step-44-with_line_search.cc . + +#include "../tests.h" +#include "step-44.h" + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + using namespace dealii; + using namespace Step44; + try + { + const unsigned int dim = 3; + Solid solid(SOURCE_DIR + "/prm/parameters-step-44-without_linesearch.prm"); + solid.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; +} diff --git a/tests/optimization/step-44-without_line_search.output b/tests/optimization/step-44-without_line_search.output new file mode 100644 index 0000000000..4c17c7fac4 --- /dev/null +++ b/tests/optimization/step-44-without_line_search.output @@ -0,0 +1,3 @@ + +DEAL::Timestep 1: 0.00000 -0.000509022 0.00000 +DEAL::Timestep 2: 0.00000 -0.000789713 0.00000 diff --git a/tests/optimization/step-44.h b/tests/optimization/step-44.h new file mode 100644 index 0000000000..649d203a1e --- /dev/null +++ b/tests/optimization/step-44.h @@ -0,0 +1,2247 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// This is a copy of step-44 (git rev 3f7e617), with modifications made +// to test the implementation of the line-search algorithm. + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include + +namespace Step44 +{ + using namespace dealii; + namespace Parameters + { + 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(); + } + struct Geometry + { + unsigned 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(); + } + 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(); + } + struct LinearSolver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + bool use_static_condensation; + std::string preconditioner_type; + 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("Use static condensation", + "true", + Patterns::Bool(), + "Solve the full block system or a reduced problem"); + 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"); + use_static_condensation = prm.get_bool("Use static condensation"); + preconditioner_type = prm.get("Preconditioner type"); + preconditioner_relaxation = prm.get_double("Preconditioner relaxation"); + } + prm.leave_subsection(); + } + struct NonlinearSolver + { + unsigned int max_iterations_NR; + double tol_f; + double tol_u; + bool use_line_search; + bool use_line_search_ptrb; + 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.declare_entry("Use line search", + "false", + Patterns::Bool(), + "Use line-search minimization algorithm"); + prm.declare_entry( + "Use line search approximate gradient", + "false", + Patterns::Bool(), + "Use an approximation for the gradient during line search"); + } + 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"); + use_line_search = prm.get_bool("Use line search"); + use_line_search_ptrb = + prm.get_bool("Use line search approximate gradient"); + } + prm.leave_subsection(); + } + 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(); + } + 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.parse_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); + } + } // namespace Parameters + template + class StandardTensors + { + public: + static const SymmetricTensor<2, dim> I; + static const SymmetricTensor<4, dim> IxI; + static const SymmetricTensor<4, dim> II; + static const SymmetricTensor<4, dim> dev_P; + }; + template + const SymmetricTensor<2, dim> + StandardTensors::I = unit_symmetric_tensor(); + template + const SymmetricTensor<4, dim> StandardTensors::IxI = outer_product(I, I); + template + const SymmetricTensor<4, dim> + StandardTensors::II = identity_tensor(); + template + const SymmetricTensor<4, dim> + StandardTensors::dev_P = deviator_tensor(); + 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() const + { + return time_end; + } + double + get_delta_t() const + { + return delta_t; + } + unsigned int + get_timestep() const + { + return timestep; + } + void + increment() + { + time_current += delta_t; + ++timestep; + } + + private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; + }; + template + class Material_Compressible_Neo_Hook_Three_Field + { + public: + Material_Compressible_Neo_Hook_Three_Field(const double mu, const double nu) + : kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))) + , c_1(mu / 2.0) + , det_F(1.0) + , p_tilde(0.0) + , J_tilde(1.0) + , b_bar(StandardTensors::I) + { + Assert(kappa > 0, ExcInternalError()); + } + ~Material_Compressible_Neo_Hook_Three_Field() + {} + 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 / dim) * symmetrize(F * transpose(F)); + p_tilde = p_tilde_in; + J_tilde = J_tilde_in; + Assert(det_F > 0, ExcInternalError()); + } + SymmetricTensor<2, dim> + get_tau() + { + return get_tau_iso() + get_tau_vol(); + } + SymmetricTensor<4, dim> + get_Jc() const + { + return get_Jc_vol() + get_Jc_iso(); + } + double + get_dPsi_vol_dJ() const + { + return (kappa / 2.0) * (J_tilde - 1.0 / J_tilde); + } + double + get_d2Psi_vol_dJ2() const + { + return ((kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde))); + } + double + get_det_F() const + { + return det_F; + } + double + get_p_tilde() const + { + return p_tilde; + } + double + get_J_tilde() const + { + return J_tilde; + } + + protected: + const double kappa; + const double c_1; + double det_F; + double p_tilde; + double J_tilde; + SymmetricTensor<2, dim> b_bar; + SymmetricTensor<2, dim> + get_tau_vol() const + { + return p_tilde * det_F * StandardTensors::I; + } + SymmetricTensor<2, dim> + get_tau_iso() const + { + return StandardTensors::dev_P * get_tau_bar(); + } + SymmetricTensor<2, dim> + get_tau_bar() const + { + return 2.0 * c_1 * b_bar; + } + SymmetricTensor<4, dim> + get_Jc_vol() const + { + return p_tilde * det_F * + (StandardTensors::IxI - (2.0 * StandardTensors::II)); + } + 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, StandardTensors::I); + const SymmetricTensor<4, dim> I_x_tau_iso = + outer_product(StandardTensors::I, tau_iso); + const SymmetricTensor<4, dim> c_bar = get_c_bar(); + return (2.0 / dim) * trace(tau_bar) * StandardTensors::dev_P - + (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso) + + StandardTensors::dev_P * c_bar * StandardTensors::dev_P; + } + SymmetricTensor<4, dim> + get_c_bar() const + { + return SymmetricTensor<4, dim>(); + } + }; + template + class PointHistory + { + public: + PointHistory() + : F_inv(StandardTensors::I) + , tau(SymmetricTensor<2, dim>()) + , d2Psi_vol_dJ2(0.0) + , dPsi_vol_dJ(0.0) + , Jc(SymmetricTensor<4, dim>()) + {} + virtual ~PointHistory() + {} + void + setup_lqp(const Parameters::AllParameters ¶meters) + { + material.reset( + new Material_Compressible_Neo_Hook_Three_Field(parameters.mu, + parameters.nu)); + update_values(Tensor<2, dim>(), 0.0, 1.0); + } + void + update_values(const Tensor<2, dim> &Grad_u_n, + const double p_tilde, + const double J_tilde) + { + const Tensor<2, dim> F = + (Tensor<2, dim>(StandardTensors::I) + Grad_u_n); + material->update_material_data(F, p_tilde, J_tilde); + F_inv = invert(F); + tau = material->get_tau(); + Jc = material->get_Jc(); + dPsi_vol_dJ = material->get_dPsi_vol_dJ(); + d2Psi_vol_dJ2 = material->get_d2Psi_vol_dJ2(); + } + double + get_J_tilde() const + { + return material->get_J_tilde(); + } + double + get_det_F() const + { + return material->get_det_F(); + } + const Tensor<2, dim> & + get_F_inv() const + { + return F_inv; + } + double + get_p_tilde() const + { + return material->get_p_tilde(); + } + const SymmetricTensor<2, dim> & + get_tau() const + { + return tau; + } + double + get_dPsi_vol_dJ() const + { + return dPsi_vol_dJ; + } + double + get_d2Psi_vol_dJ2() const + { + return d2Psi_vol_dJ2; + } + const SymmetricTensor<4, dim> & + get_Jc() const + { + return Jc; + } + + private: + std::shared_ptr> material; + Tensor<2, dim> F_inv; + SymmetricTensor<2, dim> tau; + double d2Psi_vol_dJ2; + double dPsi_vol_dJ; + SymmetricTensor<4, dim> Jc; + }; + template + class Solid + { + public: + Solid(const std::string &input_file); + virtual ~Solid(); + void + run(); + + private: + struct PerTaskData_K; + struct ScratchData_K; + struct PerTaskData_RHS; + struct ScratchData_RHS; + struct PerTaskData_SC; + struct ScratchData_SC; + struct PerTaskData_UQPH; + struct ScratchData_UQPH; + void + make_grid(); + void + system_setup(); + void + determine_component_extractors(); + void + assemble_system_tangent(BlockSparseMatrix &tangent_matrix, + const bool print = true); + void + assemble_system_tangent_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_K & scratch, + PerTaskData_K & data) const; + void + copy_local_to_global_K(const PerTaskData_K &data); + void + assemble_system_rhs(BlockVector &system_rhs, + const bool print = true); + void + assemble_system_rhs_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_RHS & scratch, + PerTaskData_RHS & data) const; + void + copy_local_to_global_rhs(const PerTaskData_RHS &data); + void + assemble_sc(); + void + assemble_sc_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_SC & scratch, + PerTaskData_SC & data); + void + copy_local_to_global_sc(const PerTaskData_SC &data); + void + make_constraints(const int &it_nr); + void + setup_qph(); + void + update_qph_incremental(const BlockVector &solution_delta, + const bool print = true); + void + update_qph_incremental_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_UQPH & scratch, + PerTaskData_UQPH & data); + void + copy_local_to_global_UQPH(const PerTaskData_UQPH & /*data*/) + {} + void + solve_nonlinear_timestep(BlockVector &solution_delta); + std::pair + solve_linear_system(BlockVector &newton_update); + BlockVector + get_total_solution(const BlockVector &solution_delta) const; + void + output_results() const; + Parameters::AllParameters parameters; + double vol_reference; + Triangulation triangulation; + Time time; + mutable TimerOutput timer; + CellDataStorage::cell_iterator, + PointHistory> + quadrature_point_history; + const unsigned int degree; + const FESystem fe; + DoFHandler dof_handler_ref; + const unsigned int dofs_per_cell; + const FEValuesExtractors::Vector u_fe; + const FEValuesExtractors::Scalar p_fe; + const FEValuesExtractors::Scalar J_fe; + static const unsigned int n_blocks = 3; + static const unsigned int n_components = dim + 2; + static const unsigned int first_u_component = 0; + static const unsigned int p_component = dim; + static const unsigned int J_component = dim + 1; + enum + { + u_dof = 0, + p_dof = 1, + J_dof = 2 + }; + std::vector dofs_per_block; + std::vector element_indices_u; + std::vector element_indices_p; + std::vector element_indices_J; + const QGauss qf_cell; + const QGauss qf_face; + const unsigned int n_q_points; + const unsigned int n_q_points_f; + AffineConstraints constraints; + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix tangent_matrix; + BlockVector system_rhs; + BlockVector solution_n; + 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; + void + get_error_residual(Errors &error_residual); + void + get_error_update(const BlockVector &newton_update, + Errors & error_update); + std::pair + get_error_dilation() const; + double + compute_vol_current() const; + static void + print_conv_header(); + void + print_conv_footer(); + }; + template + Solid::Solid(const std::string &input_file) + : parameters(input_file) + , triangulation(Triangulation::maximum_smoothing) + , time(parameters.end_time, parameters.delta_t) + , timer(std::cout, TimerOutput::never, TimerOutput::wall_times) + , degree(parameters.poly_degree) + , fe(FE_Q(parameters.poly_degree), + dim, // displacement + FE_DGPMonomial(parameters.poly_degree - 1), + 1, // pressure + FE_DGPMonomial(parameters.poly_degree - 1), + 1) + , // dilatation + 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()) + { + Assert(dim == 2 || dim == 3, + ExcMessage("This problem only works in 2 or 3 space dimensions.")); + determine_component_extractors(); + } + template + Solid::~Solid() + { + dof_handler_ref.clear(); + } + template + void + Solid::run() + { + make_grid(); + system_setup(); + { + AffineConstraints constraints; + constraints.close(); + const ComponentSelectFunction J_mask(J_component, n_components); + VectorTools::project(dof_handler_ref, + constraints, + QGauss(degree + 2), + J_mask, + solution_n); + } + // output_results(); + time.increment(); + BlockVector solution_delta(dofs_per_block); + while (time.current() <= time.end()) + { + solution_delta = 0.0; + solve_nonlinear_timestep(solution_delta); + solution_n += solution_delta; + // output_results(); + // Output displacement at centre of traction surface + { + const Point soln_pt( + dim == 3 ? Point(0.0, 1.0, 0.0) * parameters.scale : + Point(0.0, 1.0) * parameters.scale); + typename DoFHandler::active_cell_iterator + cell = dof_handler_ref.begin_active(), + endc = dof_handler_ref.end(); + for (; cell != endc; ++cell) + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; + ++v) + if (cell->vertex(v).distance(soln_pt) < 1e-6 * parameters.scale) + { + Tensor<1, dim> soln; + for (unsigned int d = 0; d < dim; ++d) + soln[d] = solution_n(cell->vertex_dof_index(v, u_dof + d)); + deallog << "Timestep " << time.get_timestep() << ": " << soln + << std::endl; + } + } + time.increment(); + } + } + template + struct Solid::PerTaskData_K + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + BlockSparseMatrix * tangent_matrix; + PerTaskData_K(const unsigned int dofs_per_cell, + BlockSparseMatrix &tangent_matrix) + : cell_matrix(dofs_per_cell, dofs_per_cell) + , local_dof_indices(dofs_per_cell) + , tangent_matrix(&tangent_matrix) + {} + void + reset() + { + cell_matrix = 0.0; + } + }; + template + struct Solid::ScratchData_K + { + FEValues fe_values_ref; + std::vector> Nx; + std::vector>> grad_Nx; + std::vector>> symm_grad_Nx; + ScratchData_K(const FiniteElement &fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell) + : fe_values_ref(fe_cell, qf_cell, uf_cell) + , Nx(qf_cell.size(), std::vector(fe_cell.dofs_per_cell)) + , grad_Nx(qf_cell.size(), + std::vector>(fe_cell.dofs_per_cell)) + , symm_grad_Nx(qf_cell.size(), + std::vector>( + fe_cell.dofs_per_cell)) + {} + ScratchData_K(const ScratchData_K &rhs) + : fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()) + , Nx(rhs.Nx) + , grad_Nx(rhs.grad_Nx) + , symm_grad_Nx(rhs.symm_grad_Nx) + {} + void + reset() + { + const unsigned int n_q_points = Nx.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); + Assert(grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) + { + Nx[q_point][k] = 0.0; + grad_Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + } + } + } + }; + template + struct Solid::PerTaskData_RHS + { + Vector cell_rhs; + std::vector local_dof_indices; + BlockVector * system_rhs; + PerTaskData_RHS(const unsigned int dofs_per_cell, + BlockVector &system_rhs) + : cell_rhs(dofs_per_cell) + , local_dof_indices(dofs_per_cell) + , system_rhs(&system_rhs) + {} + void + reset() + { + cell_rhs = 0.0; + } + }; + template + struct Solid::ScratchData_RHS + { + FEValues fe_values_ref; + FEFaceValues fe_face_values_ref; + std::vector> Nx; + std::vector>> symm_grad_Nx; + ScratchData_RHS(const FiniteElement &fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell, + const QGauss & qf_face, + const UpdateFlags uf_face) + : fe_values_ref(fe_cell, qf_cell, uf_cell) + , fe_face_values_ref(fe_cell, qf_face, uf_face) + , Nx(qf_cell.size(), std::vector(fe_cell.dofs_per_cell)) + , symm_grad_Nx(qf_cell.size(), + std::vector>( + fe_cell.dofs_per_cell)) + {} + ScratchData_RHS(const ScratchData_RHS &rhs) + : fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()) + , fe_face_values_ref(rhs.fe_face_values_ref.get_fe(), + rhs.fe_face_values_ref.get_quadrature(), + rhs.fe_face_values_ref.get_update_flags()) + , Nx(rhs.Nx) + , symm_grad_Nx(rhs.symm_grad_Nx) + {} + void + reset() + { + const unsigned int n_q_points = Nx.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); + Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) + { + Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + } + } + } + }; + template + struct Solid::PerTaskData_SC + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + FullMatrix k_orig; + FullMatrix k_pu; + FullMatrix k_pJ; + FullMatrix k_JJ; + FullMatrix k_pJ_inv; + FullMatrix k_bbar; + FullMatrix A; + FullMatrix B; + FullMatrix C; + PerTaskData_SC(const unsigned int dofs_per_cell, + const unsigned int n_u, + const unsigned int n_p, + const unsigned int n_J) + : cell_matrix(dofs_per_cell, dofs_per_cell) + , local_dof_indices(dofs_per_cell) + , k_orig(dofs_per_cell, dofs_per_cell) + , k_pu(n_p, n_u) + , k_pJ(n_p, n_J) + , k_JJ(n_J, n_J) + , k_pJ_inv(n_p, n_J) + , k_bbar(n_u, n_u) + , A(n_J, n_u) + , B(n_J, n_u) + , C(n_p, n_u) + {} + void + reset() + {} + }; + template + struct Solid::ScratchData_SC + { + void + reset() + {} + }; + template + struct Solid::PerTaskData_UQPH + { + void + reset() + {} + }; + template + struct Solid::ScratchData_UQPH + { + const BlockVector & solution_total; + std::vector> solution_grads_u_total; + std::vector solution_values_p_total; + std::vector solution_values_J_total; + FEValues fe_values_ref; + ScratchData_UQPH(const FiniteElement & fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell, + const BlockVector &solution_total) + : solution_total(solution_total) + , solution_grads_u_total(qf_cell.size()) + , solution_values_p_total(qf_cell.size()) + , solution_values_J_total(qf_cell.size()) + , fe_values_ref(fe_cell, qf_cell, uf_cell) + {} + ScratchData_UQPH(const ScratchData_UQPH &rhs) + : solution_total(rhs.solution_total) + , solution_grads_u_total(rhs.solution_grads_u_total) + , solution_values_p_total(rhs.solution_values_p_total) + , solution_values_J_total(rhs.solution_values_J_total) + , fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()) + {} + void + reset() + { + const unsigned int n_q_points = solution_grads_u_total.size(); + for (unsigned int q = 0; q < n_q_points; ++q) + { + solution_grads_u_total[q] = 0.0; + solution_values_p_total[q] = 0.0; + solution_values_J_total[q] = 0.0; + } + } + }; + template + void + Solid::make_grid() + { + GridGenerator::hyper_rectangle( + triangulation, + (dim == 3 ? Point(0.0, 0.0, 0.0) : Point(0.0, 0.0)), + (dim == 3 ? Point(1.0, 1.0, 1.0) : Point(1.0, 1.0)), + true); + GridTools::scale(parameters.scale, triangulation); + triangulation.refine_global(std::max(1U, parameters.global_refinement)); + vol_reference = GridTools::volume(triangulation); + std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl; + typename Triangulation::active_cell_iterator cell = triangulation + .begin_active(), + endc = + triangulation.end(); + for (; cell != endc; ++cell) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + { + if (cell->face(face)->at_boundary() == true && + cell->face(face)->center()[1] == 1.0 * parameters.scale) + { + if (dim == 3) + { + if (cell->face(face)->center()[0] < 0.5 * parameters.scale && + cell->face(face)->center()[2] < 0.5 * parameters.scale) + cell->face(face)->set_boundary_id(6); + } + else + { + if (cell->face(face)->center()[0] < 0.5 * parameters.scale) + cell->face(face)->set_boundary_id(6); + } + } + } + } + template + void + Solid::system_setup() + { + timer.enter_subsection("Setup system"); + std::vector block_component(n_components, + u_dof); // Displacement + block_component[p_component] = p_dof; // Pressure + block_component[J_component] = J_dof; // Dilatation + dof_handler_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; + tangent_matrix.clear(); + { + const types::global_dof_index n_dofs_u = dofs_per_block[u_dof]; + const types::global_dof_index n_dofs_p = dofs_per_block[p_dof]; + const types::global_dof_index n_dofs_J = dofs_per_block[J_dof]; + BlockDynamicSparsityPattern dsp(n_blocks, n_blocks); + dsp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u); + dsp.block(u_dof, p_dof).reinit(n_dofs_u, n_dofs_p); + dsp.block(u_dof, J_dof).reinit(n_dofs_u, n_dofs_J); + dsp.block(p_dof, u_dof).reinit(n_dofs_p, n_dofs_u); + dsp.block(p_dof, p_dof).reinit(n_dofs_p, n_dofs_p); + dsp.block(p_dof, J_dof).reinit(n_dofs_p, n_dofs_J); + dsp.block(J_dof, u_dof).reinit(n_dofs_J, n_dofs_u); + dsp.block(J_dof, p_dof).reinit(n_dofs_J, n_dofs_p); + dsp.block(J_dof, J_dof).reinit(n_dofs_J, n_dofs_J); + dsp.collect_sizes(); + 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, dsp, constraints, false); + sparsity_pattern.copy_from(dsp); + } + tangent_matrix.reinit(sparsity_pattern); + system_rhs.reinit(dofs_per_block); + system_rhs.collect_sizes(); + solution_n.reinit(dofs_per_block); + solution_n.collect_sizes(); + setup_qph(); + timer.leave_subsection(); + } + template + void + Solid::determine_component_extractors() + { + element_indices_u.clear(); + element_indices_p.clear(); + element_indices_J.clear(); + for (unsigned int k = 0; k < fe.dofs_per_cell; ++k) + { + 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()); + } + } + } + template + void + Solid::setup_qph() + { + std::cout << " Setting up quadrature point data..." << std::endl; + quadrature_point_history.initialize(triangulation.begin_active(), + triangulation.end(), + n_q_points); + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); + cell != triangulation.end(); + ++cell) + { + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point]->setup_lqp(parameters); + } + } + template + void + Solid::update_qph_incremental(const BlockVector &solution_delta, + const bool print) + { + timer.enter_subsection("Update QPH data"); + if (print) + std::cout << " UQPH " << std::flush; + const BlockVector solution_total( + get_total_solution(solution_delta)); + 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); + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::update_qph_incremental_one_cell, + &Solid::copy_local_to_global_UQPH, + scratch_data_UQPH, + per_task_data_UQPH); + timer.leave_subsection(); + } + template + void + Solid::update_qph_incremental_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_UQPH & scratch, + PerTaskData_UQPH & /*data*/) + { + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, 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(); + 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); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point]->update_values(scratch.solution_grads_u_total[q_point], + scratch.solution_values_p_total[q_point], + scratch.solution_values_J_total[q_point]); + } + template + void + Solid::solve_nonlinear_timestep(BlockVector &solution_delta) + { + std::cout << std::endl + << "Timestep " << time.get_timestep() << " @ " << time.current() + << "s" << std::endl; + BlockVector newton_update(dofs_per_block); + error_residual.reset(); + error_residual_0.reset(); + error_residual_norm.reset(); + error_update.reset(); + error_update_0.reset(); + error_update_norm.reset(); + print_conv_header(); + unsigned int newton_iteration = 0; + for (; newton_iteration < parameters.max_iterations_NR; ++newton_iteration) + { + std::cout << " " << std::setw(2) << newton_iteration << " " + << std::flush; + tangent_matrix = 0.0; + system_rhs = 0.0; + assemble_system_rhs(system_rhs); + get_error_residual(error_residual); + if (newton_iteration == 0) + error_residual_0 = error_residual; + error_residual_norm = error_residual; + error_residual_norm.normalise(error_residual_0); + if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u && + error_residual_norm.u <= parameters.tol_f) + { + std::cout << " CONVERGED! " << std::endl; + print_conv_footer(); + break; + } + assemble_system_tangent(this->tangent_matrix); + make_constraints(newton_iteration); + constraints.condense(tangent_matrix, system_rhs); + const std::pair lin_solver_output = + solve_linear_system(newton_update); + + // Here (in a most inefficient way) we check the effect of + // taking the proposed Newton step, and adjust it using + // the line-search algorith. + + // Previous step residual + BlockVector residual_0 = this->system_rhs; + residual_0 *= -1.0; // Residual = -RHS + constraints.set_zero(residual_0); + const double tang_mtrx_norm_old = + tangent_matrix.block(0, 0).frobenius_norm(); + + // Minimization function (exact) + auto ls_minimization_function = + [this, + &residual_0, + &solution_delta, + &newton_update, + &tang_mtrx_norm_old](const double ss /*step size*/) { + // Ensure that the constraints for the Dirichlet BC's are correct, + // irrespective of the chosen step size. + BlockVector solution_delta_trial(newton_update); + solution_delta_trial *= ss; + this->constraints.distribute(solution_delta_trial); + // Now add the constribution from the previously accepted solution + // history. + solution_delta_trial += solution_delta; + + // Compute residual for proposed step + this->update_qph_incremental(solution_delta_trial, false); + this->assemble_system_rhs(this->system_rhs, false); + this->assemble_system_tangent(this->tangent_matrix, false); + constraints.condense(this->tangent_matrix, this->system_rhs); + BlockVector residual_trial = this->system_rhs; + residual_trial *= -1.0; // Residual = -RHS + + + const double tang_mtrx_norm_new = + tangent_matrix.block(0, 0).frobenius_norm(); + + if (ss != 0.0) + Assert(tang_mtrx_norm_new != tang_mtrx_norm_old, + ExcInternalError()); + + // Negelect the constrained entries in the consideration + this->constraints.set_zero(residual_trial); + + // Wriggers p159 eq 5.11 + const double f = 0.5 * (residual_trial * residual_trial); // Value + // Wriggers p159 eq 5.14 is wrong... + // const double g = -(residual_0*residual_trial); // Gradient + // This should be g = G(V + alpha*delta)*[ K(V + alpha*delta)*delta + // ] + BlockVector tmp(newton_update); + tangent_matrix.vmult(tmp, newton_update); + const double g = tmp * residual_trial; // Gradient + + return std::make_pair(f, g); + }; + + // Minimization function (perturbation method, to confirm result + // of exact minimization function's gradient -- this verifies + // that Wriggers 5.14 is incorrect) + auto ls_minimization_function_ptrb = + [&ls_minimization_function](const double ss) { + const double ptrb = 0.005 * 1e-2; + const auto f_m = ls_minimization_function(ss); + const auto f_minus = ls_minimization_function(ss - ptrb); + const auto f_plus = ls_minimization_function(ss + ptrb); + return std::make_pair(f_m.first, + (f_plus.first - f_minus.first) / (2. * ptrb)); + }; + + // Invoke line search + const bool use_ptrb = parameters.use_line_search_ptrb; + auto perform_linesearch = [&ls_minimization_function, + &ls_minimization_function_ptrb, + &use_ptrb]() { + const auto res_0 = (use_ptrb ? ls_minimization_function_ptrb(0.0) : + ls_minimization_function(0.0)); + Assert(res_0.second < 0.0, + ExcMessage("Gradient should be negative. Current value: " + + std::to_string(res_0.second))); + const auto res_1 = (use_ptrb ? ls_minimization_function_ptrb(1.0) : + ls_minimization_function(1.0)); + + // Wriggers discussion after 5.14 + if (res_0.second * res_1.second > 0.0) + return 1.0; + + // The values for eta, mu are chosen such that + // more strict convergence conditions are enforced. + const double a1 = 1.0; + const double eta = 0.5; + const double mu = 0.49; + const double a_max = 1.25; + const double max_evals = 20; + const bool debug_linesearch = false; + const auto res = (use_ptrb ? LineMinimization::line_search( + ls_minimization_function_ptrb, + res_0.first, + res_0.second, + LineMinimization::poly_fit, + a1, + eta, + mu, + a_max, + max_evals, + debug_linesearch) : + LineMinimization::line_search( + ls_minimization_function, + res_0.first, + res_0.second, + LineMinimization::poly_fit, + a1, + eta, + mu, + a_max, + max_evals, + debug_linesearch)); + + return res.first; // Final stepsize + }; + const double linesearch_step_size = + (parameters.use_line_search ? perform_linesearch() : 1.0); + + get_error_update(newton_update, error_update); + if (newton_iteration == 0) + error_update_0 = error_update; + error_update_norm = error_update; + error_update_norm.normalise(error_update_0); + if (linesearch_step_size != 1.0) + newton_update *= linesearch_step_size; + 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 << " " << linesearch_step_size + << " " << 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; + } + AssertThrow(newton_iteration <= parameters.max_iterations_NR, + ExcMessage("No convergence in nonlinear solver!")); + } + template + void + Solid::print_conv_header() + { + static const unsigned int l_width = 166; + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + std::cout << " SOLVER STEP " + << " | LIN_IT LIN_RES LS_STP_SZ 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; + } + template + void + Solid::print_conv_footer() + { + static const unsigned int l_width = 166; + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + const std::pair error_dil = get_error_dilation(); + 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" << error_dil.second * vol_reference << " / " + << vol_reference << " = " << error_dil.second << std::endl; + } + template + double + Solid::compute_vol_current() const + { + double vol_current = 0.0; + FEValues fe_values_ref(fe, qf_cell, update_JxW_values); + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); + cell != triangulation.end(); + ++cell) + { + fe_values_ref.reinit(cell); + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double det_F_qp = lqph[q_point]->get_det_F(); + const double JxW = fe_values_ref.JxW(q_point); + vol_current += det_F_qp * JxW; + } + } + Assert(vol_current > 0.0, ExcInternalError()); + return vol_current; + } + template + std::pair + Solid::get_error_dilation() const + { + double dil_L2_error = 0.0; + FEValues fe_values_ref(fe, qf_cell, update_JxW_values); + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); + cell != triangulation.end(); + ++cell) + { + fe_values_ref.reinit(cell); + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double det_F_qp = lqph[q_point]->get_det_F(); + const double J_tilde_qp = lqph[q_point]->get_J_tilde(); + const double the_error_qp_squared = + std::pow((det_F_qp - J_tilde_qp), 2); + const double JxW = fe_values_ref.JxW(q_point); + dil_L2_error += the_error_qp_squared * JxW; + } + } + return std::make_pair(std::sqrt(dil_L2_error), + compute_vol_current() / vol_reference); + } + template + void + Solid::get_error_residual(Errors &error_residual) + { + BlockVector error_res(dofs_per_block); + for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i) + if (!constraints.is_constrained(i)) + error_res(i) = system_rhs(i); + error_residual.norm = error_res.l2_norm(); + error_residual.u = error_res.block(u_dof).l2_norm(); + error_residual.p = error_res.block(p_dof).l2_norm(); + error_residual.J = error_res.block(J_dof).l2_norm(); + } + template + void + Solid::get_error_update(const BlockVector &newton_update, + Errors & error_update) + { + BlockVector error_ud(dofs_per_block); + 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 + BlockVector + Solid::get_total_solution( + const BlockVector &solution_delta) const + { + BlockVector solution_total(solution_n); + solution_total += solution_delta; + return solution_total; + } + template + void + Solid::assemble_system_tangent(BlockSparseMatrix &tangent_matrix, + const bool print) + { + timer.enter_subsection("Assemble tangent matrix"); + if (print) + std::cout << " ASM_K " << std::flush; + tangent_matrix = 0.0; + const UpdateFlags uf_cell(update_values | update_gradients | + update_JxW_values); + PerTaskData_K per_task_data(dofs_per_cell, tangent_matrix); + ScratchData_K scratch_data(fe, qf_cell, uf_cell); + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + std::bind(&Solid::assemble_system_tangent_one_cell, + this, + std::placeholders::_1, + std::placeholders::_2, + std::placeholders::_3), + std::bind(&Solid::copy_local_to_global_K, + this, + std::placeholders::_1), + scratch_data, + per_task_data); + timer.leave_subsection(); + } + template + void + Solid::copy_local_to_global_K(const PerTaskData_K &data) + { + BlockSparseMatrix &tangent_matrix = *data.tangent_matrix; + 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)); + } + template + void + Solid::assemble_system_tangent_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_K & scratch, + PerTaskData_K & data) const + { + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + 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()); + } + } + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + 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(); + const std::vector & N = scratch.Nx[q_point]; + const std::vector> &symm_grad_Nx = + scratch.symm_grad_Nx[q_point]; + const std::vector> &grad_Nx = scratch.grad_Nx[q_point]; + const double JxW = scratch.fe_values_ref.JxW(q_point); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int component_i = + fe.system_to_component_index(i).first; + 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; + if ((i_group == j_group) && (i_group == u_dof)) + { + data.cell_matrix(i, j) += symm_grad_Nx[i] * + Jc // The material contribution: + * 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; + } + else if ((i_group == p_dof) && (j_group == u_dof)) + { + data.cell_matrix(i, j) += + N[i] * det_F * + (symm_grad_Nx[j] * StandardTensors::I) * JxW; + } + else if ((i_group == J_dof) && (j_group == p_dof)) + data.cell_matrix(i, j) -= N[i] * N[j] * JxW; + 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()); + } + } + } + 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); + } + template + void + Solid::assemble_system_rhs(BlockVector &system_rhs, + const bool print) + { + timer.enter_subsection("Assemble system right-hand side"); + if (print) + 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, system_rhs); + ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face); + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + std::bind(&Solid::assemble_system_rhs_one_cell, + this, + std::placeholders::_1, + std::placeholders::_2, + std::placeholders::_3), + std::bind(&Solid::copy_local_to_global_rhs, + this, + std::placeholders::_1), + scratch_data, + per_task_data); + timer.leave_subsection(); + } + template + void + Solid::copy_local_to_global_rhs(const PerTaskData_RHS &data) + { + BlockVector &system_rhs = *data.system_rhs; + for (unsigned int i = 0; i < dofs_per_cell; ++i) + system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i); + } + template + void + Solid::assemble_system_rhs_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_RHS & scratch, + PerTaskData_RHS & data) const + { + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + const std::vector>> lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + 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()); + } + } + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + 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(); + const std::vector & N = scratch.Nx[q_point]; + const std::vector> &symm_grad_Nx = + scratch.symm_grad_Nx[q_point]; + const double JxW = scratch.fe_values_ref.JxW(q_point); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int i_group = fe.system_to_base_index(i).first.first; + if (i_group == u_dof) + data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW; + else if (i_group == p_dof) + data.cell_rhs(i) -= N[i] * (det_F - J_tilde) * JxW; + 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()); + } + } + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + if (cell->face(face)->at_boundary() == true && + cell->face(face)->boundary_id() == 6) + { + scratch.fe_face_values_ref.reinit(cell, face); + for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; + ++f_q_point) + { + const Tensor<1, dim> &N = + scratch.fe_face_values_ref.normal_vector(f_q_point); + 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) + { + 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); + data.cell_rhs(i) += (Ni * traction[component_i]) * JxW; + } + } + } + } + } + template + void + Solid::make_constraints(const int &it_nr) + { + std::cout << " CST " << std::flush; + if (it_nr > 1) + return; + constraints.clear(); + const bool apply_dirichlet_bc = (it_nr == 0); + const FEValuesExtractors::Scalar x_displacement(0); + const FEValuesExtractors::Scalar y_displacement(1); + { + const int boundary_id = 0; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement)); + } + { + const int boundary_id = 2; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + } + if (dim == 3) + { + const FEValuesExtractors::Scalar z_displacement(2); + { + const int boundary_id = 3; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + } + { + const int boundary_id = 4; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + } + { + const int boundary_id = 6; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + } + } + else + { + { + const int boundary_id = 3; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + } + { + const int boundary_id = 6; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + } + } + constraints.close(); + } + template + void + Solid::assemble_sc() + { + timer.enter_subsection("Perform static condensation"); + std::cout << " ASM_SC " << std::flush; + PerTaskData_SC per_task_data(dofs_per_cell, + element_indices_u.size(), + element_indices_p.size(), + element_indices_J.size()); + ScratchData_SC scratch_data; + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::assemble_sc_one_cell, + &Solid::copy_local_to_global_sc, + scratch_data, + per_task_data); + timer.leave_subsection(); + } + template + void + Solid::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)); + } + template + void + Solid::assemble_sc_one_cell( + const typename DoFHandler::active_cell_iterator &cell, + ScratchData_SC & scratch, + PerTaskData_SC & data) + { + data.reset(); + scratch.reset(); + cell->get_dof_indices(data.local_dof_indices); + data.k_orig.extract_submatrix_from(tangent_matrix, + data.local_dof_indices, + data.local_dof_indices); + data.k_pu.extract_submatrix_from(data.k_orig, + element_indices_p, + element_indices_u); + data.k_pJ.extract_submatrix_from(data.k_orig, + element_indices_p, + element_indices_J); + data.k_JJ.extract_submatrix_from(data.k_orig, + element_indices_J, + element_indices_J); + data.k_pJ_inv.invert(data.k_pJ); + data.k_pJ_inv.mmult(data.A, data.k_pu); + data.k_JJ.mmult(data.B, data.A); + data.k_pJ_inv.Tmmult(data.C, data.B); + data.k_pu.Tmmult(data.k_bbar, data.C); + data.k_bbar.scatter_matrix_to(element_indices_u, + element_indices_u, + data.cell_matrix); + data.k_pJ_inv.add(-1.0, data.k_pJ); + data.k_pJ_inv.scatter_matrix_to(element_indices_p, + element_indices_J, + data.cell_matrix); + } + template + std::pair + Solid::solve_linear_system(BlockVector &newton_update) + { + unsigned int lin_it = 0; + double lin_res = 0.0; + if (parameters.use_static_condensation == true) + { + BlockVector A(dofs_per_block); + BlockVector B(dofs_per_block); + { + assemble_sc(); + tangent_matrix.block(p_dof, J_dof) + .vmult(A.block(J_dof), system_rhs.block(p_dof)); + tangent_matrix.block(J_dof, J_dof) + .vmult(B.block(J_dof), A.block(J_dof)); + A.block(J_dof) = system_rhs.block(J_dof); + A.block(J_dof) -= B.block(J_dof); + tangent_matrix.block(p_dof, J_dof) + .Tvmult(A.block(p_dof), A.block(J_dof)); + tangent_matrix.block(u_dof, p_dof) + .vmult(A.block(u_dof), A.block(p_dof)); + 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, false, false); + GrowingVectorMemory> GVM; + SolverCG> solver_CG(solver_control, GVM); + PreconditionSelector, Vector> + preconditioner(parameters.preconditioner_type, + parameters.preconditioner_relaxation); + preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof)); + solver_CG.solve(tangent_matrix.block(u_dof, u_dof), + newton_update.block(u_dof), + system_rhs.block(u_dof), + preconditioner); + lin_it = solver_control.last_step(); + lin_res = solver_control.last_value(); + } + else if (parameters.type_lin == "Direct") + { + 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 + Assert(false, ExcMessage("Linear solver type not implemented")); + timer.leave_subsection(); + } + constraints.distribute(newton_update); + timer.enter_subsection("Linear solver postprocessing"); + std::cout << " PP " << std::flush; + { + tangent_matrix.block(p_dof, u_dof) + .vmult(A.block(p_dof), newton_update.block(u_dof)); + A.block(p_dof) *= -1.0; + A.block(p_dof) += system_rhs.block(p_dof); + tangent_matrix.block(p_dof, J_dof) + .vmult(newton_update.block(J_dof), A.block(p_dof)); + } + constraints.distribute(newton_update); + { + tangent_matrix.block(J_dof, J_dof) + .vmult(A.block(J_dof), newton_update.block(J_dof)); + A.block(J_dof) *= -1.0; + A.block(J_dof) += system_rhs.block(J_dof); + tangent_matrix.block(p_dof, J_dof) + .Tvmult(newton_update.block(p_dof), A.block(J_dof)); + } + constraints.distribute(newton_update); + timer.leave_subsection(); + } + else + { + std::cout << " ------ " << std::flush; + timer.enter_subsection("Linear solver"); + std::cout << " SLV " << std::flush; + if (parameters.type_lin == "CG") + { + const Vector &f_u = system_rhs.block(u_dof); + const Vector &f_p = system_rhs.block(p_dof); + const Vector &f_J = system_rhs.block(J_dof); + Vector & d_u = newton_update.block(u_dof); + Vector & d_p = newton_update.block(p_dof); + Vector & d_J = newton_update.block(J_dof); + const auto K_uu = + linear_operator(tangent_matrix.block(u_dof, u_dof)); + const auto K_up = + linear_operator(tangent_matrix.block(u_dof, p_dof)); + const auto K_pu = + linear_operator(tangent_matrix.block(p_dof, u_dof)); + const auto K_Jp = + linear_operator(tangent_matrix.block(J_dof, p_dof)); + const auto K_JJ = + linear_operator(tangent_matrix.block(J_dof, J_dof)); + PreconditionSelector, Vector> + preconditioner_K_Jp_inv("jacobi"); + preconditioner_K_Jp_inv.use_matrix( + tangent_matrix.block(J_dof, p_dof)); + ReductionControl solver_control_K_Jp_inv( + tangent_matrix.block(J_dof, p_dof).m() * + parameters.max_iterations_lin, + 1.0e-30, + parameters.tol_lin); + SolverSelector> solver_K_Jp_inv; + solver_K_Jp_inv.select("cg"); + solver_K_Jp_inv.set_control(solver_control_K_Jp_inv); + const auto K_Jp_inv = + inverse_operator(K_Jp, solver_K_Jp_inv, preconditioner_K_Jp_inv); + const auto K_pJ_inv = transpose_operator(K_Jp_inv); + const auto K_pp_bar = K_Jp_inv * K_JJ * K_pJ_inv; + const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu; + const auto K_uu_con = K_uu + K_uu_bar_bar; + PreconditionSelector, Vector> + preconditioner_K_con_inv(parameters.preconditioner_type, + parameters.preconditioner_relaxation); + preconditioner_K_con_inv.use_matrix( + tangent_matrix.block(u_dof, u_dof)); + ReductionControl solver_control_K_con_inv( + tangent_matrix.block(u_dof, u_dof).m() * + parameters.max_iterations_lin, + 1.0e-30, + parameters.tol_lin); + SolverSelector> solver_K_con_inv; + solver_K_con_inv.select("cg"); + solver_K_con_inv.set_control(solver_control_K_con_inv); + const auto K_uu_con_inv = + inverse_operator(K_uu_con, + solver_K_con_inv, + preconditioner_K_con_inv); + d_u = + K_uu_con_inv * (f_u - K_up * (K_Jp_inv * f_J - K_pp_bar * f_p)); + timer.leave_subsection(); + timer.enter_subsection("Linear solver postprocessing"); + std::cout << " PP " << std::flush; + d_J = K_pJ_inv * (f_p - K_pu * d_u); + d_p = K_Jp_inv * (f_J - K_JJ * d_J); + lin_it = solver_control_K_con_inv.last_step(); + lin_res = solver_control_K_con_inv.last_value(); + } + else if (parameters.type_lin == "Direct") + { + SparseDirectUMFPACK A_direct; + A_direct.initialize(tangent_matrix); + A_direct.vmult(newton_update, system_rhs); + lin_it = 1; + lin_res = 0.0; + std::cout << " -- " << std::flush; + } + else + Assert(false, ExcMessage("Linear solver type not implemented")); + timer.leave_subsection(); + constraints.distribute(newton_update); + } + return std::make_pair(lin_it, lin_res); + } + template + void + Solid::output_results() const + { + DataOut data_out; + std::vector + data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + std::vector solution_name(dim, "displacement"); + solution_name.push_back("pressure"); + solution_name.push_back("dilatation"); + data_out.attach_dof_handler(dof_handler_ref); + data_out.add_data_vector(solution_n, + solution_name, + DataOut::type_dof_data, + data_component_interpretation); + Vector soln(solution_n.size()); + for (unsigned int i = 0; i < soln.size(); ++i) + soln(i) = solution_n(i); + MappingQEulerian q_mapping(degree, dof_handler_ref, soln); + data_out.build_patches(q_mapping, degree); + std::ostringstream filename; + filename << "solution-" << dim << "d-" << time.get_timestep() << ".vtk"; + std::ofstream output(filename.str().c_str()); + data_out.write_vtk(output); + } +} // namespace Step44 -- 2.39.5