From 3fe8f6cb9301d91d952ea679101c80ab58b85d67 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Tue, 26 Apr 2022 20:55:39 +0200 Subject: [PATCH] step-44: Change the FE space used for p,J --- .../changes/minor/20220426Jean-PaulPelteret | 6 + examples/step-44/step-44.cc | 24 +- .../step-44-helper_res_lin_01.h | 6 +- .../step-44-helper_scalar_01.h | 6 +- .../step-44-helper_var_form_01.h | 6 +- .../step-44-helper_vector_01.h | 6 +- tests/adolc/step-44.cc | 6 +- tests/optimization/step-44.h | 6 +- tests/physics/step-44-dgp_monomial.cc | 2132 +++++++++++++++++ tests/physics/step-44-dgp_monomial.output | 11 + ...-standard_tensors-material_push_forward.cc | 6 +- .../step-44-standard_tensors-spatial.cc | 6 +- tests/physics/step-44.cc | 10 +- .../step-44-sd-quadrature_level.h | 6 +- 14 files changed, 2195 insertions(+), 42 deletions(-) create mode 100644 doc/news/changes/minor/20220426Jean-PaulPelteret create mode 100644 tests/physics/step-44-dgp_monomial.cc create mode 100644 tests/physics/step-44-dgp_monomial.output diff --git a/doc/news/changes/minor/20220426Jean-PaulPelteret b/doc/news/changes/minor/20220426Jean-PaulPelteret new file mode 100644 index 0000000000..3ac6fa63d0 --- /dev/null +++ b/doc/news/changes/minor/20220426Jean-PaulPelteret @@ -0,0 +1,6 @@ +Improved: The step-44 tutorial now uses FE_DGP elements to discretize the +pressure and dilatation fields, rather than FE_DGPMonomial. This improves +the condition number of the linear system, necessitating fewer iterations +to compute its iterative inverse. +
+(Jean-Paul Pelteret, 2022/04/26) diff --git a/examples/step-44/step-44.cc b/examples/step-44/step-44.cc index 40dd349253..26454322b8 100644 --- a/examples/step-44/step-44.cc +++ b/examples/step-44/step-44.cc @@ -43,7 +43,7 @@ #include #include -#include +#include #include #include #include @@ -1037,19 +1037,18 @@ namespace Step44 // The Finite Element System is composed of dim continuous displacement // DOFs, and discontinuous pressure and dilatation DOFs. In an attempt to // satisfy the Babuska-Brezzi or LBB stability conditions (see Hughes - // (2000)), we setup a $Q_n \times DGPM_{n-1} \times DGPM_{n-1}$ - // system. $Q_2 \times DGPM_1 \times DGPM_1$ elements satisfy this - // condition, while $Q_1 \times DGPM_0 \times DGPM_0$ elements do + // (2000)), we setup a $Q_n \times DGP_{n-1} \times DGP_{n-1}$ + // system. $Q_2 \times DGP_1 \times DGP_1$ elements satisfy this + // condition, while $Q_1 \times DGP_0 \times DGP_0$ elements do // not. However, it has been shown that the latter demonstrate good // convergence characteristics nonetheless. fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), - 1) - , // dilatation - dof_handler(triangulation) + FE_DGP(parameters.poly_degree - 1), + 1) // dilatation + , dof_handler(triangulation) , dofs_per_cell(fe.n_dofs_per_cell()) , u_fe(first_u_component) , p_fe(p_component) @@ -1080,10 +1079,11 @@ namespace Step44 // constraint $\widetilde{J}=1$ on the initial solution field. The constraint // corresponds to the determinant of the deformation gradient in the // undeformed configuration, which is the identity tensor. We use - // FE_DGPMonomial bases to interpolate the dilatation field, thus we can't + // FE_DGP bases to interpolate the dilatation field, thus we can't // simply set the corresponding dof to unity as they correspond to the - // monomial coefficients. Thus we use the VectorTools::project function to do - // the work for us. The VectorTools::project function requires an argument + // coefficients of a truncated Legendre polynomial. + // Thus we use the VectorTools::project function to do the work for us. + // The VectorTools::project function requires an argument // indicating the hanging node constraints. We have none in this program // So we have to create a constraint object. In its original state, constraint // objects are unsorted, and have to be sorted (using the diff --git a/tests/ad_common_tests/step-44-helper_res_lin_01.h b/tests/ad_common_tests/step-44-helper_res_lin_01.h index 52ed155120..c8d7ffa961 100644 --- a/tests/ad_common_tests/step-44-helper_res_lin_01.h +++ b/tests/ad_common_tests/step-44-helper_res_lin_01.h @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include #include @@ -654,9 +654,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/ad_common_tests/step-44-helper_scalar_01.h b/tests/ad_common_tests/step-44-helper_scalar_01.h index eb36eed6b9..5590d0e0a2 100644 --- a/tests/ad_common_tests/step-44-helper_scalar_01.h +++ b/tests/ad_common_tests/step-44-helper_scalar_01.h @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include #include @@ -934,9 +934,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/ad_common_tests/step-44-helper_var_form_01.h b/tests/ad_common_tests/step-44-helper_var_form_01.h index 7dbc207505..e471c8e3fc 100644 --- a/tests/ad_common_tests/step-44-helper_var_form_01.h +++ b/tests/ad_common_tests/step-44-helper_var_form_01.h @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include #include @@ -620,9 +620,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/ad_common_tests/step-44-helper_vector_01.h b/tests/ad_common_tests/step-44-helper_vector_01.h index d0a735e26c..40a21c0d66 100644 --- a/tests/ad_common_tests/step-44-helper_vector_01.h +++ b/tests/ad_common_tests/step-44-helper_vector_01.h @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include #include @@ -655,9 +655,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/adolc/step-44.cc b/tests/adolc/step-44.cc index 3790385964..42a2c4b926 100644 --- a/tests/adolc/step-44.cc +++ b/tests/adolc/step-44.cc @@ -29,7 +29,7 @@ #include #include -#include +#include #include #include #include @@ -769,9 +769,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/optimization/step-44.h b/tests/optimization/step-44.h index 1e329be61f..ed401f55bc 100644 --- a/tests/optimization/step-44.h +++ b/tests/optimization/step-44.h @@ -29,7 +29,7 @@ #include #include -#include +#include #include #include #include @@ -795,9 +795,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/physics/step-44-dgp_monomial.cc b/tests/physics/step-44-dgp_monomial.cc new file mode 100644 index 0000000000..b589c8fbf4 --- /dev/null +++ b/tests/physics/step-44-dgp_monomial.cc @@ -0,0 +1,2132 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// This is a copy of step-44 (git rev 3f7e617) to use as a base-line for +// results produced via different approaches to be compared to. +// +// This variant of the test uses the FE_DGPMonomial class to discretize the +// pressure and dilatation fields, as was originally done in commit 3f7e617. + +#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 "../tests.h" +namespace Step44 +{ + 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; + 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(); + } + 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(); + 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(); + 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); + 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; + ConditionalOStream pcout; + 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; + 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()) + , pcout(std::cout) + { + 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 (const unsigned int v : GeometryInfo::vertex_indices()) + 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; + 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; + } + }; + 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; + 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; + } + }; + 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); + pcout << "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 (const unsigned int face : GeometryInfo::face_indices()) + { + 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); + dofs_per_block = + DoFTools::count_dofs_per_fe_block(dof_handler_ref, block_component); + pcout << "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() + { + pcout << " 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) + { + timer.enter_subsection("Update QPH data"); + pcout << " 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 (const unsigned int q_point : + scratch.fe_values_ref.quadrature_point_indices()) + 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) + { + pcout << 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) + { + pcout << ' ' << std::setw(2) << newton_iteration << ' ' << std::flush; + tangent_matrix = 0.0; + system_rhs = 0.0; + assemble_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) + { + pcout << " CONVERGED! " << std::endl; + print_conv_footer(); + break; + } + assemble_system_tangent(); + make_constraints(newton_iteration); + constraints.condense(tangent_matrix, system_rhs); + const std::pair lin_solver_output = + solve_linear_system(newton_update); + 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); + solution_delta += newton_update; + update_qph_incremental(solution_delta); + pcout << " | " << 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; + } + 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 = 155; + for (unsigned int i = 0; i < l_width; ++i) + pcout << '_'; + pcout << std::endl; + pcout << " 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) + pcout << '_'; + pcout << std::endl; + } + template + void + Solid::print_conv_footer() + { + static const unsigned int l_width = 155; + for (unsigned int i = 0; i < l_width; ++i) + pcout << '_'; + pcout << std::endl; + const std::pair error_dil = get_error_dilation(); + pcout << "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 (const unsigned int q_point : + fe_values_ref.quadrature_point_indices()) + { + 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 (const unsigned int q_point : + fe_values_ref.quadrature_point_indices()) + { + 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() + { + timer.enter_subsection("Assemble tangent matrix"); + pcout << " 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); + 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) + { + 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 (const unsigned int q_point : + scratch.fe_values_ref.quadrature_point_indices()) + { + const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv(); + for (const unsigned int k : scratch.fe_values_ref.dof_indices()) + { + 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 (const unsigned int q_point : + scratch.fe_values_ref.quadrature_point_indices()) + { + 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 (const unsigned int i : scratch.fe_values_ref.dof_indices()) + { + 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 (const unsigned int j : + scratch.fe_values_ref.dof_indices_ending_at(i)) + { + 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 (const unsigned int i : scratch.fe_values_ref.dof_indices()) + for (const unsigned int j : + scratch.fe_values_ref.dof_indices_starting_at(i + 1)) + data.cell_matrix(i, j) = data.cell_matrix(j, i); + } + template + void + Solid::assemble_system_rhs() + { + timer.enter_subsection("Assemble system right-hand side"); + pcout << " 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(), + 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) + { + 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 (const unsigned int q_point : + scratch.fe_values_ref.quadrature_point_indices()) + { + const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv(); + for (const unsigned int k : scratch.fe_values_ref.dof_indices()) + { + 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 (const unsigned int q_point : + scratch.fe_values_ref.quadrature_point_indices()) + { + 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 (const unsigned int i : scratch.fe_values_ref.dof_indices()) + { + 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 (const unsigned int face : GeometryInfo::face_indices()) + if (cell->face(face)->at_boundary() == true && + cell->face(face)->boundary_id() == 6) + { + scratch.fe_face_values_ref.reinit(cell, face); + for (const unsigned int f_q_point : + scratch.fe_face_values_ref.quadrature_point_indices()) + { + 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 (const unsigned int i : scratch.fe_values_ref.dof_indices()) + { + 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) + { + pcout << " 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, + Functions::ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + Functions::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, + Functions::ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + Functions::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, + Functions::ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + Functions::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, + Functions::ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + Functions::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, + Functions::ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement) | + fe.component_mask(z_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + Functions::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, + Functions::ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + Functions::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, + Functions::ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + else + VectorTools::interpolate_boundary_values( + dof_handler_ref, + boundary_id, + Functions::ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + } + } + constraints.close(); + } + template + void + Solid::assemble_sc() + { + timer.enter_subsection("Perform static condensation"); + pcout << " 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"); + pcout << " SLV " << std::flush; + if (parameters.type_lin == "CG") + { + const auto solver_its = static_cast( + 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"); + pcout << " 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 + { + pcout << " ------ " << std::flush; + timer.enter_subsection("Linear solver"); + pcout << " 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( + static_cast(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( + static_cast(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"); + pcout << " 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; + pcout << " -- " << 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 +int +main(int argc, char **argv) +{ + initlog(); + deallog << std::setprecision(3); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + using namespace Step44; + try + { + const unsigned int dim = 3; + Solid solid(SOURCE_DIR "/prm/parameters-step-44.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/physics/step-44-dgp_monomial.output b/tests/physics/step-44-dgp_monomial.output new file mode 100644 index 0000000000..7ed68e5e88 --- /dev/null +++ b/tests/physics/step-44-dgp_monomial.output @@ -0,0 +1,11 @@ + +DEAL::Timestep 1: 0.00 -0.000110 0.00 +DEAL::Timestep 2: 0.00 -0.000220 0.00 +DEAL::Timestep 3: 0.00 -0.000325 0.00 +DEAL::Timestep 4: 0.00 -0.000422 0.00 +DEAL::Timestep 5: 0.00 -0.000509 0.00 +DEAL::Timestep 6: 0.00 -0.000585 0.00 +DEAL::Timestep 7: 0.00 -0.000650 0.00 +DEAL::Timestep 8: 0.00 -0.000705 0.00 +DEAL::Timestep 9: 0.00 -0.000752 0.00 +DEAL::Timestep 10: 0.00 -0.000790 0.00 diff --git a/tests/physics/step-44-standard_tensors-material_push_forward.cc b/tests/physics/step-44-standard_tensors-material_push_forward.cc index fe16b083ac..07fdc739b5 100644 --- a/tests/physics/step-44-standard_tensors-material_push_forward.cc +++ b/tests/physics/step-44-standard_tensors-material_push_forward.cc @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include #include @@ -964,9 +964,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/physics/step-44-standard_tensors-spatial.cc b/tests/physics/step-44-standard_tensors-spatial.cc index 72f7f83b10..e1ef3b2b52 100644 --- a/tests/physics/step-44-standard_tensors-spatial.cc +++ b/tests/physics/step-44-standard_tensors-spatial.cc @@ -30,7 +30,7 @@ #include #include -#include +#include #include #include #include @@ -764,9 +764,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/physics/step-44.cc b/tests/physics/step-44.cc index 3df381bfcd..32867428f1 100644 --- a/tests/physics/step-44.cc +++ b/tests/physics/step-44.cc @@ -15,6 +15,10 @@ // This is a copy of step-44 (git rev 3f7e617) to use as a base-line for // results produced via different approaches to be compared to. +// +// This variant of the test uses the FE_DGP class to discretize the +// pressure and dilatation fields, contrary to what is done in commit 3f7e617 +// but in line with the current implementation of the tutorial. #include #include @@ -30,7 +34,7 @@ #include #include -#include +#include #include #include #include @@ -777,9 +781,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) diff --git a/tests/symengine/sd_common_tests/step-44-sd-quadrature_level.h b/tests/symengine/sd_common_tests/step-44-sd-quadrature_level.h index 962ce30367..5561a9f350 100644 --- a/tests/symengine/sd_common_tests/step-44-sd-quadrature_level.h +++ b/tests/symengine/sd_common_tests/step-44-sd-quadrature_level.h @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include #include @@ -721,9 +721,9 @@ namespace Step44 , degree(parameters.poly_degree) , fe(FE_Q(parameters.poly_degree), dim, // displacement - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1, // pressure - FE_DGPMonomial(parameters.poly_degree - 1), + FE_DGP(parameters.poly_degree - 1), 1) , // dilatation dof_handler_ref(triangulation) -- 2.39.5