From 536a725326fa5574ecb7b9d57676506c398cb478 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Fri, 1 Dec 2017 07:38:04 +0200 Subject: [PATCH] Add assembly for cook membrane problem using automatic differentiation --- .../CMakeLists.txt | 2 +- .../cook_membrane.cc | 1560 ++++++++--------- .../doc/builds-on | 2 +- .../parameters.prm | 8 + 4 files changed, 789 insertions(+), 783 deletions(-) diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt b/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt index 724337c..b6b6434 100644 --- a/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt @@ -20,7 +20,7 @@ SET(CLEAN_UP_FILES CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) -FIND_PACKAGE(deal.II 8.0 QUIET +FIND_PACKAGE(deal.II 8.5 QUIET HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} ) IF(NOT ${deal.II_FOUND}) diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc index 55a75a6..521f093 100644 --- a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc @@ -16,7 +16,7 @@ /* * Authors: Jean-Paul Pelteret, University of Erlangen-Nuremberg, - * Andrew McBride, University of Cape Town, 2015 + * Andrew McBride, University of Cape Town, 2015, 2017 */ @@ -31,6 +31,7 @@ #include #include #include +#include #include #include @@ -52,7 +53,7 @@ #include #include -#include +#include #include #include #include @@ -62,6 +63,18 @@ #include #include +#include +#if DEAL_II_VERSION_MAJOR >= 9 && defined(DEAL_II_WITH_TRILINOS) + #include + #define ENABLE_SACADO_FORMULATION +#endif + +// These must be included below the AD headers so that +// their math functions are available for use in the +// definition of tensors and kinematic quantities +#include +#include + #include #include @@ -79,6 +92,46 @@ namespace Cook_Membrane // ParameterHandler object to read in the choices at run-time. namespace Parameters { +// @sect4{Assembly method} + +// Here we specify whether automatic differentiation is to be used to assemble +// the linear system, and if so then what order of differentiation is to be +// employed. + struct AssemblyMethod + { + unsigned int automatic_differentiation_order; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + void AssemblyMethod::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Assembly method"); + { + prm.declare_entry("Automatic differentiation order", "0", + Patterns::Integer(0,2), + "The automatic differentiation order to be used in the assembly of the linear system.\n" + "# Order = 0: Both the residual and linearisation are computed manually.\n" + "# Order = 1: The residual is computed manually but the linearisation is performed using AD.\n" + "# Order = 2: Both the residual and linearisation are computed using AD."); + } + prm.leave_subsection(); + } + + void AssemblyMethod::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Assembly method"); + { + automatic_differentiation_order = prm.get_integer("Automatic differentiation order"); + } + prm.leave_subsection(); + } + // @sect4{Finite Element system} // Here we specify the polynomial order used to approximate the solution. @@ -181,11 +234,11 @@ namespace Cook_Membrane { prm.enter_subsection("Material properties"); { - prm.declare_entry("Poisson's ratio", "0.4999", + prm.declare_entry("Poisson's ratio", "0.3", Patterns::Double(-1.0,0.5), "Poisson's ratio"); - prm.declare_entry("Shear modulus", "80.194e6", + prm.declare_entry("Shear modulus", "0.4225e6", Patterns::Double(), "Shear modulus"); } @@ -354,7 +407,9 @@ namespace Cook_Membrane // Finally we consolidate all of the above structures into a single container // that holds all of our run-time selections. - struct AllParameters : public FESystem, + struct AllParameters : + public AssemblyMethod, + public FESystem, public Geometry, public Materials, public LinearSolver, @@ -375,12 +430,13 @@ namespace Cook_Membrane { ParameterHandler prm; declare_parameters(prm); - prm.read_input(input_file); + prm.parse_input(input_file); parse_parameters(prm); } void AllParameters::declare_parameters(ParameterHandler &prm) { + AssemblyMethod::declare_parameters(prm); FESystem::declare_parameters(prm); Geometry::declare_parameters(prm); Materials::declare_parameters(prm); @@ -391,6 +447,7 @@ namespace Cook_Membrane void AllParameters::parse_parameters(ParameterHandler &prm) { + AssemblyMethod::parse_parameters(prm); FESystem::parse_parameters(prm); Geometry::parse_parameters(prm); Materials::parse_parameters(prm); @@ -517,7 +574,7 @@ namespace Cook_Membrane // store the current state (characterized by the values or measures of the // displacement field) so that we can compute the elastic coefficients // linearized around the current state. - template + template class Material_Compressible_Neo_Hook_One_Field { public: @@ -525,9 +582,7 @@ namespace Cook_Membrane 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), - b_bar(StandardTensors::I) + c_1(mu / 2.0) { Assert(kappa > 0, ExcInternalError()); } @@ -535,24 +590,23 @@ namespace Cook_Membrane ~Material_Compressible_Neo_Hook_One_Field() {} - // We update the material model with various deformation dependent data - // based on the deformation gradient $\mathbf{F}$ and the volumetric Jacobian - // $J = \text{det} \mathbf{F}$, and at the end of the function include a - // physical check for internal consistency: - void update_material_data(const Tensor<2, dim> &F) + // The first function is the total energy + // $\Psi = \Psi_{\textrm{iso}} + \Psi_{\textrm{vol}}$. + NumberType + get_Psi(const NumberType &det_F, + const SymmetricTensor<2,dim,NumberType> &b_bar) const { - det_F = determinant(F); - b_bar = std::pow(det_F, -2.0 / dim) * symmetrize(F * transpose(F)); - - Assert(det_F > 0, ExcInternalError()); + return get_Psi_vol(det_F) + get_Psi_iso(b_bar); } // The second function determines the Kirchhoff stress $\boldsymbol{\tau} // = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}$ - SymmetricTensor<2, dim> get_tau() + SymmetricTensor<2,dim,NumberType> + get_tau(const NumberType &det_F, + const SymmetricTensor<2,dim,NumberType> &b_bar) { - // See Holzapfel p231 eq6.98 onwards - return get_tau_iso() + get_tau_vol(); + // See Holzapfel p231 eq6.98 onwards + return get_tau_vol(det_F) + get_tau_iso(b_bar); } // The fourth-order elasticity tensor in the spatial setting @@ -560,89 +614,101 @@ namespace Cook_Membrane // \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}$ // where $ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial // \mathbf{C} \partial \mathbf{C}}$ - SymmetricTensor<4, dim> get_Jc() const + SymmetricTensor<4,dim,NumberType> + get_Jc(const NumberType &det_F, + const SymmetricTensor<2,dim,NumberType> &b_bar) const { - return get_Jc_vol() + get_Jc_iso(); + return get_Jc_vol(det_F) + get_Jc_iso(b_bar); } - // Derivative of the volumetric free energy with respect to - // $J$ return $\frac{\partial - // \Psi_{\text{vol}}(J)}{\partial J}$ - double get_dPsi_vol_dJ() const + private: + // Define constitutive model parameters $\kappa$ (bulk modulus) and the + // neo-Hookean model parameter $c_1$: + const double kappa; + const double c_1; + + // Value of the volumetric free energy + NumberType + get_Psi_vol(const NumberType &det_F) const { - return (kappa / 2.0) * (det_F - 1.0 / det_F); + return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F)); } - // Second derivative of the volumetric free energy wrt $J$. We - // need the following computation explicitly in the tangent so we make it - // public. We calculate $\frac{\partial^2 - // \Psi_{\textrm{vol}}(J)}{\partial J \partial - // J}$ - double get_d2Psi_vol_dJ2() const + // Value of the isochoric free energy + NumberType + get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const { - return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F))); + return c_1 * (trace(b_bar) - dim); } - // The next few functions return various data that we choose to store with - // the material: - double get_det_F() const + // Derivative of the volumetric free energy with respect to + // $J$ return $\frac{\partial + // \Psi_{\text{vol}}(J)}{\partial J}$ + NumberType + get_dPsi_vol_dJ(const NumberType &det_F) const { - return det_F; + return (kappa / 2.0) * (det_F - 1.0 / det_F); } - private: - // Define constitutive model parameters $\kappa$ (bulk modulus) and the - // neo-Hookean model parameter $c_1$: - const double kappa; - const double c_1; - - // Model specific data that is convenient to store with the material: - double det_F; - SymmetricTensor<2, dim> b_bar; - // The following functions are used internally in determining the result // of some of the public functions above. The first one determines the // volumetric Kirchhoff stress $\boldsymbol{\tau}_{\textrm{vol}}$. // Note the difference in its definition when compared to step-44. - SymmetricTensor<2, dim> get_tau_vol() const + SymmetricTensor<2,dim,NumberType> + get_tau_vol(const NumberType &det_F) const { - return get_dPsi_vol_dJ() * det_F * StandardTensors::I; + return NumberType(get_dPsi_vol_dJ(det_F) * det_F) * StandardTensors::I; } // Next, determine the isochoric Kirchhoff stress // $\boldsymbol{\tau}_{\textrm{iso}} = // \mathcal{P}:\overline{\boldsymbol{\tau}}$: - SymmetricTensor<2, dim> get_tau_iso() const + SymmetricTensor<2,dim,NumberType> + get_tau_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const { - return StandardTensors::dev_P * get_tau_bar(); + return StandardTensors::dev_P * get_tau_bar(b_bar); } // Then, determine the fictitious Kirchhoff stress // $\overline{\boldsymbol{\tau}}$: - SymmetricTensor<2, dim> get_tau_bar() const + SymmetricTensor<2,dim,NumberType> + get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const { return 2.0 * c_1 * b_bar; } + // Second derivative of the volumetric free energy wrt $J$. We + // need the following computation explicitly in the tangent so we make it + // public. We calculate $\frac{\partial^2 + // \Psi_{\textrm{vol}}(J)}{\partial J \partial + // J}$ + NumberType + get_d2Psi_vol_dJ2(const NumberType &det_F) const + { + return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F))); + } + // Calculate the volumetric part of the tangent $J // \mathfrak{c}_\textrm{vol}$. Again, note the difference in its // definition when compared to step-44. The extra terms result from two // quantities in $\boldsymbol{\tau}_{\textrm{vol}}$ being dependent on // $\boldsymbol{F}$. - SymmetricTensor<4, dim> get_Jc_vol() const + SymmetricTensor<4,dim,NumberType> + get_Jc_vol(const NumberType &det_F) const { // See Holzapfel p265 return det_F - * ( (get_dPsi_vol_dJ() + det_F * get_d2Psi_vol_dJ2())*StandardTensors::IxI - - (2.0 * get_dPsi_vol_dJ())*StandardTensors::II ); + * ( (get_dPsi_vol_dJ(det_F) + det_F * get_d2Psi_vol_dJ2(det_F))*StandardTensors::IxI + - (2.0 * get_dPsi_vol_dJ(det_F))*StandardTensors::II ); } // Calculate the isochoric part of the tangent $J // \mathfrak{c}_\textrm{iso}$: - SymmetricTensor<4, dim> get_Jc_iso() const + SymmetricTensor<4,dim,NumberType> + get_Jc_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const { - const SymmetricTensor<2, dim> tau_bar = get_tau_bar(); - const SymmetricTensor<2, dim> tau_iso = get_tau_iso(); + const SymmetricTensor<2, dim> tau_bar = get_tau_bar(b_bar); + const SymmetricTensor<2, dim> tau_iso = get_tau_iso(b_bar); const SymmetricTensor<4, dim> tau_iso_x_I = outer_product(tau_iso, StandardTensors::I); @@ -660,7 +726,8 @@ namespace Cook_Membrane // Calculate the fictitious elasticity tensor $\overline{\mathfrak{c}}$. // For the material model chosen this is simply zero: - SymmetricTensor<4, dim> get_c_bar() const + SymmetricTensor<4,dim,double> + get_c_bar() const { return SymmetricTensor<4, dim>(); } @@ -674,15 +741,11 @@ namespace Cook_Membrane // can be used in different regions of the domain. Among other data, we // choose to store the Kirchhoff stress $\boldsymbol{\tau}$ and the tangent // $J\mathfrak{c}$ for the quadrature points. - template + template class PointHistory { public: PointHistory() - : - F_inv(StandardTensors::I), - tau(SymmetricTensor<2, dim>()), - Jc(SymmetricTensor<4, dim>()) {} virtual ~PointHistory() @@ -694,68 +757,35 @@ namespace Cook_Membrane // $\textrm{Grad}\mathbf{u}_{\textrm{n}}$. void setup_lqp (const Parameters::AllParameters ¶meters) { - material.reset(new Material_Compressible_Neo_Hook_One_Field(parameters.mu, + material.reset(new Material_Compressible_Neo_Hook_One_Field(parameters.mu, parameters.nu)); - update_values(Tensor<2, dim>()); } - // To this end, we calculate the deformation gradient $\mathbf{F}$ from - // the displacement gradient $\textrm{Grad}\ \mathbf{u}$, i.e. - // $\mathbf{F}(\mathbf{u}) = \mathbf{I} + \textrm{Grad}\ \mathbf{u}$ and - // then let the material model associated with this quadrature point - // update itself. When computing the deformation gradient, we have to take - // care with which data types we compare the sum $\mathbf{I} + - // \textrm{Grad}\ \mathbf{u}$: Since $I$ has data type SymmetricTensor, - // just writing I + Grad_u_n would convert the second - // argument to a symmetric tensor, perform the sum, and then cast the - // result to a Tensor (i.e., the type of a possibly nonsymmetric - // tensor). However, since Grad_u_n is nonsymmetric in - // general, the conversion to SymmetricTensor will fail. We can avoid this - // back and forth by converting $I$ to Tensor first, and then performing - // the addition as between nonsymmetric tensors: - void update_values (const Tensor<2, dim> &Grad_u_n) - { - const Tensor<2, dim> F - = (Tensor<2, dim>(StandardTensors::I) + - Grad_u_n); - material->update_material_data(F); - - // The material has been updated so we now calculate the Kirchhoff - // stress $\mathbf{\tau}$, the tangent $J\mathfrak{c}$ and the first and - // second derivatives of the volumetric free energy. - // - // We also store the inverse of the deformation gradient since we - // frequently use it: - F_inv = invert(F); - tau = material->get_tau(); - Jc = material->get_Jc(); - } - - // We offer an interface to retrieve certain data. Here are the kinematic - // variables: - - double get_det_F() const + // We offer an interface to retrieve certain data. + // This is the strain energy: + NumberType + get_Psi(const NumberType &det_F, + const SymmetricTensor<2,dim,NumberType> &b_bar) const { - return material->get_det_F(); + return material->get_Psi(det_F,b_bar); } - const Tensor<2, dim> &get_F_inv() const - { - return F_inv; - } - - // ...and the kinetic variables. These are used in the material and + // Here are the kinetic variables. These are used in the material and // global tangent matrix and residual assembly operations: - - const SymmetricTensor<2, dim> &get_tau() const + // First is the Kirchhoff stress: + SymmetricTensor<2,dim,NumberType> + get_tau(const NumberType &det_F, + const SymmetricTensor<2,dim,NumberType> &b_bar) const { - return tau; + return material->get_tau(det_F,b_bar); } - // And finally the tangent: - const SymmetricTensor<4, dim> &get_Jc() const + // And the tangent: + SymmetricTensor<4,dim,NumberType> + get_Jc(const NumberType &det_F, + const SymmetricTensor<2,dim,NumberType> &b_bar) const { - return Jc; + return material->get_Jc(det_F,b_bar); } // In terms of member functions, this class stores for the quadrature @@ -763,27 +793,28 @@ namespace Cook_Membrane // materials are used in different regions of the domain, as well as the // inverse of the deformation gradient... private: - std_cxx11::shared_ptr< Material_Compressible_Neo_Hook_One_Field > material; - - Tensor<2, dim> F_inv; - - // ... and stress-type variables along with the tangent $J\mathfrak{c}$: - SymmetricTensor<2, dim> tau; - SymmetricTensor<4, dim> Jc; + std_cxx11::shared_ptr< Material_Compressible_Neo_Hook_One_Field > material; }; // @sect3{Quasi-static compressible finite-strain solid} + // Forward declarations for classes that will + // perform assembly of the linear system. + template + struct Assembler_Base; + template + struct Assembler; + // The Solid class is the central class in that it represents the problem at // hand. It follows the usual scheme in that all it really has is a // constructor, destructor and a run() function that dispatches // all the work to private functions of this class: - template + template class Solid { public: - Solid(const std::string &input_file); + Solid(const Parameters::AllParameters ¶meters); virtual ~Solid(); @@ -793,21 +824,6 @@ namespace Cook_Membrane private: - // In the private section of this class, we first forward declare a number - // of objects that are used in parallelizing work using the WorkStream - // object (see the @ref threads module for more information on this). - // - // We declare such structures for the computation of tangent (stiffness) - // matrix, right hand side, and for updating quadrature points: - struct PerTaskData_K; - struct ScratchData_K; - - struct PerTaskData_RHS; - struct ScratchData_RHS; - - struct PerTaskData_UQPH; - struct ScratchData_UQPH; - // We start the collection of member functions with one that builds the // grid: void @@ -823,26 +839,13 @@ namespace Cook_Membrane // and one that copies the work done on this one cell into the global // object that represents it: void - assemble_system_tangent(); - - void - assemble_system_tangent_one_cell(const typename DoFHandler::active_cell_iterator &cell, - ScratchData_K &scratch, - PerTaskData_K &data); - - void - copy_local_to_global_K(const PerTaskData_K &data); - - void - assemble_system_rhs(); - - void - assemble_system_rhs_one_cell(const typename DoFHandler::active_cell_iterator &cell, - ScratchData_RHS &scratch, - PerTaskData_RHS &data); + assemble_system(const BlockVector &solution_delta); - void - copy_local_to_global_rhs(const PerTaskData_RHS &data); + // We use a separate data structure to perform the assembly. It needs access + // to some low-level data, so we simply befriend the class instead of + // creating a complex interface to provide access as necessary. + friend struct Assembler_Base; + friend struct Assembler; // Apply Dirichlet boundary conditions on the displacement field void @@ -854,18 +857,6 @@ namespace Cook_Membrane void setup_qph(); - void - update_qph_incremental(const BlockVector &solution_delta); - - void - update_qph_incremental_one_cell(const typename DoFHandler::active_cell_iterator &cell, - ScratchData_UQPH &scratch, - PerTaskData_UQPH &data); - - void - copy_local_to_global_UQPH(const PerTaskData_UQPH &/*data*/) - {} - // Solve for the displacement using a Newton-Raphson method. We break this // function into the nonlinear loop and the function that solves the // linearized Newton-Raphson step: @@ -884,7 +875,7 @@ namespace Cook_Membrane // Finally, some member variables that describe the current state: A // collection of the parameters used to describe the problem setup... - Parameters::AllParameters parameters; + const Parameters::AllParameters ¶meters; // ...the volume of the reference and current configurations... double vol_reference; @@ -898,9 +889,10 @@ namespace Cook_Membrane Time time; TimerOutput timer; - // A storage object for quadrature point information. See step-18 for - // more on this: - std::vector > quadrature_point_history; + // A storage object for quadrature point information. As opposed to + // step-18, deal.II's native quadrature point data manager is employed here. + CellDataStorage::cell_iterator, + PointHistory > quadrature_point_history; // A description of the finite-element system including the displacement // polynomial degree, the degree-of-freedom handler, number of DoFs per @@ -1000,10 +992,12 @@ namespace Cook_Membrane // @sect4{Public interface} // We initialise the Solid class using data extracted from the parameter file. - template - Solid::Solid(const std::string &input_file) + template + Solid::Solid(const Parameters::AllParameters ¶meters) : - parameters(input_file), + parameters(parameters), + vol_reference (0.0), + vol_current (0.0), triangulation(Triangulation::maximum_smoothing), time(parameters.end_time, parameters.delta_t), timer(std::cout, @@ -1026,8 +1020,8 @@ namespace Cook_Membrane } // The class destructor simply clears the data held by the DOFHandler - template - Solid::~Solid() + template + Solid::~Solid() { dof_handler_ref.clear(); } @@ -1042,8 +1036,8 @@ namespace Cook_Membrane // before starting the simulation proper with the first time (and loading) // increment. // - template - void Solid::run() + template + void Solid::run() { make_grid(); system_setup(); @@ -1080,237 +1074,6 @@ namespace Cook_Membrane // @sect3{Private interface} -// @sect4{Threading-building-blocks structures} - -// The first group of private member functions is related to parallization. -// We use the Threading Building Blocks library (TBB) to perform as many -// computationally intensive distributed tasks as possible. In particular, we -// assemble the tangent matrix and right hand side vector, and update data -// stored at the quadrature points using TBB. Our main tool for this is the -// WorkStream class (see the @ref threads module for more information). - -// Firstly we deal with the tangent matrix assembly structures. The -// PerTaskData object stores local contributions. - template - struct Solid::PerTaskData_K - { - FullMatrix cell_matrix; - std::vector local_dof_indices; - - PerTaskData_K(const unsigned int dofs_per_cell) - : - cell_matrix(dofs_per_cell, dofs_per_cell), - local_dof_indices(dofs_per_cell) - {} - - void reset() - { - cell_matrix = 0.0; - } - }; - - -// On the other hand, the ScratchData object stores the larger objects such as -// the shape-function values array (Nx) and a shape function -// gradient and symmetric gradient vector which we will use during the -// assembly. - 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; - } - } - } - - }; - -// Next, the same approach is used for the right-hand side assembly. The -// PerTaskData object again stores local contributions and the ScratchData -// object the shape function object and precomputed values vector: - 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; - } - } - } - - }; - - -// And finally we define the structures to assist with updating the quadrature -// point information. We do not need the PerTaskData object (since there is -// nothing to store here) but must define one nonetheless. -// Note that this is because for the operation that we have here -- updating -// the data on quadrature points -- the operation is purely local: -// the things we do on every cell get consumed on every cell, without -// any global aggregation operation as is usually the case when using the -// WorkStream class. The fact that we still have to define a per-task data -// structure points to the fact that the WorkStream class may be ill-suited to -// this operation (we could, in principle simply create a new task using -// Threads::new_task for each cell) but there is not much harm done to doing -// it this way anyway. -// Furthermore, should there be different material models associated with a -// quadrature point, requiring varying levels of computational expense, then -// the method used here could be advantageous. - template - struct Solid::PerTaskData_UQPH - { - void reset() - {} - }; - - -// The ScratchData object will be used to store an alias for the solution -// vector so that we don't have to copy this large data structure. We then -// define a number of vectors to extract the solution values and gradients at -// the quadrature points. - template - struct Solid::ScratchData_UQPH - { - const BlockVector &solution_total; - - std::vector > solution_grads_u_total; - - 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()), - 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), - 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; - } - } - }; - - // @sect4{Solid::make_grid} // On to the first of the private member functions. Here we create the @@ -1339,8 +1102,8 @@ Point grid_y_transform (const Point &pt_in) return pt_out; } - template - void Solid::make_grid() + template + void Solid::make_grid() { // Divide the beam, but only along the x- and y-coordinate directions std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge); @@ -1388,7 +1151,6 @@ Point grid_y_transform (const Point &pt_in) vol_reference = GridTools::volume(triangulation); vol_current = vol_reference; std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl; - } @@ -1397,8 +1159,8 @@ Point grid_y_transform (const Point &pt_in) // Next we describe how the FE system is setup. We first determine the number // of components per block. Since the displacement is a vector component, the // first dim components belong to it. - template - void Solid::system_setup() + template + void Solid::system_setup() { timer.enter_subsection("Setup system"); @@ -1464,123 +1226,39 @@ Point grid_y_transform (const Point &pt_in) // // Firstly the actual QPH data objects are created. This must be done only // once the grid is refined to its finest level. - template - void Solid::setup_qph() + template + void Solid::setup_qph() { std::cout << " Setting up quadrature point data..." << std::endl; - { - triangulation.clear_user_data(); - { - std::vector > tmp; - tmp.swap(quadrature_point_history); - } - - quadrature_point_history - .resize(triangulation.n_active_cells() * n_q_points); - - unsigned int history_index = 0; - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); cell != triangulation.end(); - ++cell) - { - cell->set_user_pointer(&quadrature_point_history[history_index]); - history_index += n_q_points; - } - - Assert(history_index == quadrature_point_history.size(), - ExcInternalError()); - } + quadrature_point_history.initialize(triangulation.begin_active(), + triangulation.end(), + n_q_points); - // Next we setup the initial quadrature - // point data: + // Next we setup the initial quadrature point data. Note that when + // the quadrature point data is retrieved, it is returned as a vector + // of smart pointers. for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) { - PointHistory *lqph = - reinterpret_cast*>(cell->user_pointer()); + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); - Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); - Assert(lqph <= &quadrature_point_history.back(), ExcInternalError()); - - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - lqph[q_point].setup_lqp(parameters); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point]->setup_lqp(parameters); } } -// @sect4{Solid::update_qph_incremental} -// As the update of QP information occurs frequently and involves a number of -// expensive operations, we define a multithreaded approach to distributing -// the task across a number of CPU cores. -// -// To start this, we first we need to obtain the total solution as it stands -// at this Newton increment and then create the initial copy of the scratch and -// copy data objects: - template - void Solid::update_qph_incremental(const BlockVector &solution_delta) - { - timer.enter_subsection("Update QPH data"); - std::cout << " UQPH " << std::flush; - - const BlockVector solution_total(get_total_solution(solution_delta)); - - const UpdateFlags uf_UQPH(update_gradients); - PerTaskData_UQPH per_task_data_UQPH; - ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total); - - // We then pass them and the one-cell update function to the WorkStream to - // be processed: - WorkStream::run(dof_handler_ref.begin_active(), - dof_handler_ref.end(), - *this, - &Solid::update_qph_incremental_one_cell, - &Solid::copy_local_to_global_UQPH, - scratch_data_UQPH, - per_task_data_UQPH); - - timer.leave_subsection(); - } - - -// Now we describe how we extract data from the solution vector and pass it -// along to each QP storage object for processing. - template - void - Solid::update_qph_incremental_one_cell(const typename DoFHandler::active_cell_iterator &cell, - ScratchData_UQPH &scratch, - PerTaskData_UQPH &/*data*/) - { - PointHistory *lqph = - reinterpret_cast*>(cell->user_pointer()); - - Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); - Assert(lqph <= &quadrature_point_history.back(), ExcInternalError()); - - Assert(scratch.solution_grads_u_total.size() == n_q_points, - ExcInternalError()); - - scratch.reset(); - - // We first need to find the solution gradients at quadrature points - // inside the current cell and then we update each local QP using the - // displacement gradient: - scratch.fe_values_ref.reinit(cell); - scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total, - scratch.solution_grads_u_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]); - } - // @sect4{Solid::solve_nonlinear_timestep} // The next function is the driver method for the Newton-Raphson scheme. At // its top we create a new vector to store the current Newton update step, // reset the error storage objects and print solver header. - template + template void - Solid::solve_nonlinear_timestep(BlockVector &solution_delta) + Solid::solve_nonlinear_timestep(BlockVector &solution_delta) { std::cout << std::endl << "Timestep " << time.get_timestep() << " @ " << time.current() << "s" << std::endl; @@ -1613,10 +1291,12 @@ Point grid_y_transform (const Point &pt_in) { std::cout << " " << std::setw(2) << newton_iteration << " " << std::flush; - tangent_matrix = 0.0; - system_rhs = 0.0; + // If we have decided that we want to continue with the iteration, we + // assemble the tangent, make and impose the Dirichlet constraints, + // and do the solve of the linearized system: + make_constraints(newton_iteration); + assemble_system(solution_delta); - assemble_system_rhs(); get_error_residual(error_residual); if (newton_iteration == 0) @@ -1636,13 +1316,6 @@ Point grid_y_transform (const Point &pt_in) break; } - // If we have decided that we want to continue with the iteration, we - // assemble the tangent, make and impose the Dirichlet constraints, - // and do the solve of the linearized system: - assemble_system_tangent(); - make_constraints(newton_iteration); - constraints.condense(tangent_matrix, system_rhs); - const std::pair lin_solver_output = solve_linear_system(newton_update); @@ -1658,7 +1331,6 @@ Point grid_y_transform (const Point &pt_in) error_update_norm.normalise(error_update_0); 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 << " " @@ -1686,16 +1358,16 @@ Point grid_y_transform (const Point &pt_in) // This program prints out data in a nice table that is updated // on a per-iteration basis. The next two functions set up the table // header and footer: - template - void Solid::print_conv_header() + template + void Solid::print_conv_header() { - static const unsigned int l_width = 102; + static const unsigned int l_width = 87; for (unsigned int i = 0; i < l_width; ++i) std::cout << "_"; std::cout << std::endl; - std::cout << " SOLVER STEP " + std::cout << " SOLVER STEP " << " | LIN_IT LIN_RES RES_NORM " << " RES_U NU_NORM " << " NU_U " << std::endl; @@ -1707,10 +1379,10 @@ Point grid_y_transform (const Point &pt_in) - template - void Solid::print_conv_footer() + template + void Solid::print_conv_footer() { - static const unsigned int l_width = 102; + static const unsigned int l_width = 87; for (unsigned int i = 0; i < l_width; ++i) std::cout << "_"; @@ -1726,10 +1398,10 @@ Point grid_y_transform (const Point &pt_in) // At the end we also output the result that can be compared to that found in // the literature, namely the displacement at the upper right corner of the // beam. - template - void Solid::print_vertical_tip_displacement() + template + void Solid::print_vertical_tip_displacement() { - static const unsigned int l_width = 102; + static const unsigned int l_width = 87; for (unsigned int i = 0; i < l_width; ++i) std::cout << "_"; @@ -1788,8 +1460,8 @@ Point grid_y_transform (const Point &pt_in) // error in the residual for the unconstrained degrees of freedom. Note that to // do so, we need to ignore constrained DOFs by setting the residual in these // vector components to zero. - template - void Solid::get_error_residual(Errors &error_residual) + template + void Solid::get_error_residual(Errors &error_residual) { BlockVector error_res(dofs_per_block); @@ -1805,8 +1477,8 @@ Point grid_y_transform (const Point &pt_in) // @sect4{Solid::get_error_udpate} // Determine the true Newton update error for the problem - template - void Solid::get_error_update(const BlockVector &newton_update, + template + void Solid::get_error_update(const BlockVector &newton_update, Errors &error_update) { BlockVector error_ud(dofs_per_block); @@ -1825,9 +1497,9 @@ Point grid_y_transform (const Point &pt_in) // This function provides the total solution, which is valid at any Newton step. // This is required as, to reduce computational error, the total solution is // only updated at the end of the timestep. - template + template BlockVector - Solid::get_total_solution(const BlockVector &solution_delta) const + Solid::get_total_solution(const BlockVector &solution_delta) const { BlockVector solution_total(solution_n); solution_total += solution_delta; @@ -1835,302 +1507,585 @@ Point grid_y_transform (const Point &pt_in) } -// @sect4{Solid::assemble_system_tangent} +// @sect4{Solid::assemble_system} -// Since we use TBB for assembly, we simply setup a copy of the -// data structures required for the process and pass them, along -// with the memory addresses of the assembly functions to the -// WorkStream object for processing. Note that we must ensure that -// the matrix is reset before any assembly operations can occur. - template - void Solid::assemble_system_tangent() + template + struct Assembler_Base { - timer.enter_subsection("Assemble tangent matrix"); - std::cout << " ASM_K " << std::flush; + virtual ~Assembler_Base() {} - tangent_matrix = 0.0; + // Here we deal with the tangent matrix assembly structures. The + // PerTaskData object stores local contributions. + struct PerTaskData_ASM + { + Solid &solid; + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; - const UpdateFlags uf_cell(update_gradients | - update_JxW_values); + PerTaskData_ASM(Solid &solid) + : + solid (solid), + cell_matrix(solid.dofs_per_cell, solid.dofs_per_cell), + cell_rhs(solid.dofs_per_cell), + local_dof_indices(solid.dofs_per_cell) + {} - PerTaskData_K per_task_data(dofs_per_cell); - ScratchData_K scratch_data(fe, qf_cell, uf_cell); + void reset() + { + cell_matrix = 0.0; + cell_rhs = 0.0; + } + }; - WorkStream::run(dof_handler_ref.begin_active(), - dof_handler_ref.end(), - *this, - &Solid::assemble_system_tangent_one_cell, - &Solid::copy_local_to_global_K, - scratch_data, - per_task_data); + // On the other hand, the ScratchData object stores the larger objects such as + // the shape-function values array (Nx) and a shape function + // gradient and symmetric gradient vector which we will use during the + // assembly. + struct ScratchData_ASM + { + const BlockVector &solution_total; + std::vector > solution_grads_u_total; - timer.leave_subsection(); - } + FEValues fe_values_ref; + FEFaceValues fe_face_values_ref; -// This function adds the local contribution to the system matrix. -// Note that we choose not to use the constraint matrix to do the -// job for us because the tangent matrix and residual processes have -// been split up into two separate functions. - template - void Solid::copy_local_to_global_K(const PerTaskData_K &data) - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = 0; j < dofs_per_cell; ++j) - tangent_matrix.add(data.local_dof_indices[i], - data.local_dof_indices[j], - data.cell_matrix(i, j)); - } + std::vector > > grad_Nx; + std::vector > > symm_grad_Nx; -// Of course, we still have to define how we assemble the tangent matrix -// contribution for a single cell. We first need to reset and initialise some -// of the scratch data structures and retrieve some basic information -// regarding the DOF numbering on this cell. We can precalculate the cell -// shape function gradients. Note that the shape function gradients -// are defined with regard to the current configuration. That is -// $\textrm{grad}\ \boldsymbol{\varphi} = \textrm{Grad}\ \boldsymbol{\varphi} -// \ \mathbf{F}^{-1}$. - template - void - Solid::assemble_system_tangent_one_cell(const typename DoFHandler::active_cell_iterator &cell, - ScratchData_K &scratch, - PerTaskData_K &data) - { - data.reset(); - scratch.reset(); - scratch.fe_values_ref.reinit(cell); - cell->get_dof_indices(data.local_dof_indices); - PointHistory *lqph = - reinterpret_cast*>(cell->user_pointer()); - - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + ScratchData_ASM(const FiniteElement &fe_cell, + const QGauss &qf_cell, + const UpdateFlags uf_cell, + const QGauss & qf_face, + const UpdateFlags uf_face, + const BlockVector &solution_total) + : + solution_total(solution_total), + solution_grads_u_total(qf_cell.size()), + fe_values_ref(fe_cell, qf_cell, uf_cell), + fe_face_values_ref(fe_cell, qf_face, uf_face), + 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_ASM(const ScratchData_ASM &rhs) + : + solution_total (rhs.solution_total), + solution_grads_u_total(rhs.solution_grads_u_total), + 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()), + grad_Nx(rhs.grad_Nx), + symm_grad_Nx(rhs.symm_grad_Nx) + {} + + void reset() { - 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 n_q_points = fe_values_ref.get_quadrature().size(); + const unsigned int n_dofs_per_cell = fe_values_ref.dofs_per_cell; + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { - const unsigned int k_group = fe.system_to_base_index(k).first.first; + Assert( grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); - if (k_group == u_dof) + solution_grads_u_total[q_point] = Tensor<2,dim,NumberType>(); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) { - 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]); + grad_Nx[q_point][k] = Tensor<2,dim,NumberType>(); + symm_grad_Nx[q_point][k] = SymmetricTensor<2,dim,NumberType>(); } - else - Assert(k_group <= u_dof, ExcInternalError()); } } - // Now we build the local cell stiffness matrix. Since the global and - // local system matrices are symmetric, we can exploit this property by - // building only the lower half of the local matrix and copying the values - // to the upper half. - // - // In doing so, we first extract some configuration dependent variables - // from our QPH history objects for the current quadrature point. - 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 det_F = lqph[q_point].get_det_F(); - - // Next we define some aliases to make the assembly process easier to - // follow - const std::vector - &N = scratch.Nx[q_point]; - const std::vector > - &symm_grad_Nx = scratch.symm_grad_Nx[q_point]; - const std::vector > - &grad_Nx = scratch.grad_Nx[q_point]; - const double JxW = scratch.fe_values_ref.JxW(q_point); + }; - for (unsigned int i = 0; i < dofs_per_cell; ++i) + // Of course, we still have to define how we assemble the tangent matrix + // contribution for a single cell. + void + assemble_system_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_ASM &scratch, + PerTaskData_ASM &data) + { + // Due to the C++ specialization rules, we need one more + // level of indirection in order to define the assembly + // routine for all different number. The next function call + // is specialized for each NumberType, but to prevent having + // to specialize the whole class along with it we have inlined + // the definition of the other functions that are common to + // all implementations. + assemble_system_tangent_residual_one_cell(cell, scratch, data); + assemble_neumann_contribution_one_cell(cell, scratch, data); + } + + // This function adds the local contribution to the system matrix. + void + copy_local_to_global_ASM(const PerTaskData_ASM &data) + { + const ConstraintMatrix &constraints = data.solid.constraints; + BlockSparseMatrix &tangent_matrix = data.solid.tangent_matrix; + BlockVector &system_rhs = data.solid.system_rhs; + + constraints.distribute_local_to_global( + data.cell_matrix, data.cell_rhs, + data.local_dof_indices, + tangent_matrix, system_rhs); + } + + protected: + + // This function needs to exist in the base class for + // Workstream to work with a reference to the base class. + virtual void + assemble_system_tangent_residual_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_ASM &scratch, + PerTaskData_ASM &data) + { + AssertThrow(false, ExcPureFunctionCalled()); + } + + void + assemble_neumann_contribution_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_ASM &scratch, + PerTaskData_ASM &data) + { + // Aliases for data referenced from the Solid class + const unsigned int &n_q_points_f = data.solid.n_q_points_f; + const unsigned int &dofs_per_cell = data.solid.dofs_per_cell; + const Parameters::AllParameters ¶meters = data.solid.parameters; + const Time &time = data.solid.time; + const FESystem &fe = data.solid.fe; + const unsigned int &u_dof = data.solid.u_dof; + + // Next we assemble the Neumann contribution. We first check to see it the + // cell face exists on a boundary on which a traction is applied and add + // the contribution if this is the case. + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + if (cell->face(face)->at_boundary() == true + && cell->face(face)->boundary_id() == 11) { - 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; + scratch.fe_face_values_ref.reinit(cell, face); - for (unsigned int j = 0; j <= i; ++j) + for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; + ++f_q_point) { - const unsigned int component_j = fe.system_to_component_index(j).first; - const unsigned int j_group = fe.system_to_base_index(j).first.first; - - // This is the $\mathsf{\mathbf{k}}_{\mathbf{u} \mathbf{u}}$ - // contribution. It comprises a material contribution, and a - // geometrical stress contribution which is only added along - // the local matrix diagonals: - if ((i_group == j_group) && (i_group == u_dof)) + // We specify the traction in reference configuration. + // For this problem, a defined total vertical force is applied + // in the reference configuration. + // The direction of the applied traction is assumed not to + // evolve with the deformation of the domain. + + // Note that the contributions to the right hand side vector we + // compute here only exist in the displacement components of the + // vector. + const double time_ramp = (time.current() / time.end()); + const double magnitude = (1.0/(16.0*parameters.scale*1.0*parameters.scale))*time_ramp; // (Total force) / (RHS surface area) + Tensor<1,dim> dir; + dir[1] = 1.0; + const Tensor<1, dim> traction = magnitude*dir; + + for (unsigned int i = 0; i < dofs_per_cell; ++i) { - 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; + 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; + } } - else - Assert((i_group <= u_dof) && (j_group <= u_dof), - ExcInternalError()); } } - } + } - // Finally, we need to copy the lower half of the local matrix into the - // upper half: - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = i + 1; j < dofs_per_cell; ++j) - data.cell_matrix(i, j) = data.cell_matrix(j, i); - } + }; -// @sect4{Solid::assemble_system_rhs} -// The assembly of the right-hand side process is similar to the -// tangent matrix, so we will not describe it in too much detail. -// Note that since we are describing a problem with Neumann BCs, -// we will need the face normals and so must specify this in the -// update flags. template - void Solid::assemble_system_rhs() + struct Assembler : Assembler_Base { - timer.enter_subsection("Assemble system right-hand side"); - std::cout << " ASM_R " << std::flush; + typedef double NumberType; + using typename Assembler_Base::ScratchData_ASM; + using typename Assembler_Base::PerTaskData_ASM; - system_rhs = 0.0; + virtual ~Assembler() {} - const UpdateFlags uf_cell(update_values | - update_gradients | - update_JxW_values); - const UpdateFlags uf_face(update_values | - update_JxW_values); + virtual void + assemble_system_tangent_residual_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_ASM &scratch, + PerTaskData_ASM &data) + { + // Aliases for data referenced from the Solid class + const unsigned int &n_q_points = data.solid.n_q_points; + const unsigned int &dofs_per_cell = data.solid.dofs_per_cell; + const Parameters::AllParameters ¶meters = data.solid.parameters; + const FESystem &fe = data.solid.fe; + const unsigned int &u_dof = data.solid.u_dof; + const FEValuesExtractors::Vector &u_fe = data.solid.u_fe; + + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + + const std::vector > > lqph = + const_cast &>(data.solid).quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + + // We first need to find the solution gradients at quadrature points + // inside the current cell and then we update each local QP using the + // displacement gradient: + scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total, + scratch.solution_grads_u_total); + + // Now we build the local cell stiffness matrix. Since the global and + // local system matrices are symmetric, we can exploit this property by + // building only the lower half of the local matrix and copying the values + // to the upper half. + // + // In doing so, we first extract some configuration dependent variables + // from our QPH history objects for the current quadrature point. + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2,dim,NumberType> &grad_u = scratch.solution_grads_u_total[q_point]; + const Tensor<2,dim,NumberType> F = Physics::Elasticity::Kinematics::F(grad_u); + const NumberType det_F = determinant(F); + const Tensor<2,dim,NumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F); + const SymmetricTensor<2,dim,NumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar); + const Tensor<2,dim,NumberType> F_inv = invert(F); + Assert(det_F > NumberType(0.0), ExcInternalError()); + + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + const unsigned int k_group = fe.system_to_base_index(k).first.first; - PerTaskData_RHS per_task_data(dofs_per_cell); - ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face); + 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 + Assert(k_group <= u_dof, ExcInternalError()); + } - WorkStream::run(dof_handler_ref.begin_active(), - dof_handler_ref.end(), - *this, - &Solid::assemble_system_rhs_one_cell, - &Solid::copy_local_to_global_rhs, - scratch_data, - per_task_data); + const SymmetricTensor<2,dim,NumberType> tau = lqph[q_point]->get_tau(det_F,b_bar); + const SymmetricTensor<4,dim,NumberType> Jc = lqph[q_point]->get_Jc(det_F,b_bar); + const Tensor<2,dim,NumberType> tau_ns (tau); - timer.leave_subsection(); - } + // Next we define some aliases to make the assembly process easier to + // follow + 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; + if (i_group == u_dof) + data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW; + else + Assert(i_group <= u_dof, ExcInternalError()); - 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); - } + for (unsigned int j = 0; j <= i; ++j) + { + const unsigned int component_j = fe.system_to_component_index(j).first; + const unsigned int j_group = fe.system_to_base_index(j).first.first; + + // This is the $\mathsf{\mathbf{k}}_{\mathbf{u} \mathbf{u}}$ + // contribution. It comprises a material contribution, and a + // geometrical stress contribution which is only added along + // the local matrix diagonals: + if ((i_group == j_group) && (i_group == u_dof)) + { + 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_ns + * grad_Nx[j][component_j] * JxW; + } + else + Assert((i_group <= u_dof) && (j_group <= u_dof), + ExcInternalError()); + } + } + } + // Finally, we need to copy the lower half of the local matrix into the + // upper half: + 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); + } + + }; + +#ifdef ENABLE_SACADO_FORMULATION + template - void - Solid::assemble_system_rhs_one_cell(const typename DoFHandler::active_cell_iterator &cell, - ScratchData_RHS &scratch, - PerTaskData_RHS &data) + struct Assembler > : Assembler_Base > { - data.reset(); - scratch.reset(); - scratch.fe_values_ref.reinit(cell); - cell->get_dof_indices(data.local_dof_indices); - PointHistory *lqph = - reinterpret_cast*>(cell->user_pointer()); - - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + typedef Sacado::Fad::DFad ADNumberType; + using typename Assembler_Base::ScratchData_ASM; + using typename Assembler_Base::PerTaskData_ASM; + + virtual ~Assembler() {} + + virtual void + assemble_system_tangent_residual_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_ASM &scratch, + PerTaskData_ASM &data) + { + // Aliases for data referenced from the Solid class + const unsigned int &n_q_points = data.solid.n_q_points; + const unsigned int &dofs_per_cell = data.solid.dofs_per_cell; + const Parameters::AllParameters ¶meters = data.solid.parameters; + const FESystem &fe = data.solid.fe; + const unsigned int &u_dof = data.solid.u_dof; + const FEValuesExtractors::Vector &u_fe = data.solid.u_fe; + + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + + const std::vector > > lqph = + const_cast &>(data.solid).quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + + const unsigned int n_independent_variables = data.local_dof_indices.size(); + std::vector local_dof_values(n_independent_variables); + cell->get_dof_values(scratch.solution_total, + local_dof_values.begin(), + local_dof_values.end()); + + // We now retrieve a set of degree-of-freedom values that + // have the operations that are performed with them tracked. + std::vector local_dof_values_ad (n_independent_variables); + for (unsigned int i=0; i residual_ad (dofs_per_cell, ADNumberType(0.0)); + 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(); + const Tensor<2,dim,ADNumberType> &grad_u = scratch.solution_grads_u_total[q_point]; + const Tensor<2,dim,ADNumberType> F = Physics::Elasticity::Kinematics::F(grad_u); + const ADNumberType det_F = determinant(F); +// const Tensor<2,dim,ADNumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F); + const Tensor<2,dim,ADNumberType> F_bar = ADNumberType(std::pow(determinant(F),-1.0/dim))*F; + const SymmetricTensor<2,dim,ADNumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar); + const Tensor<2,dim,ADNumberType> F_inv = invert(F); + Assert(det_F > ADNumberType(0.0), ExcInternalError()); 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); + { + 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 Assert(k_group <= u_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 SymmetricTensor<2,dim,ADNumberType> tau = lqph[q_point]->get_tau(det_F,b_bar); - const std::vector - &N = scratch.Nx[q_point]; - const std::vector > - &symm_grad_Nx = scratch.symm_grad_Nx[q_point]; + // Next we define some position-dependent aliases, again to + // make the assembly process easier to follow. + const std::vector > &symm_grad_Nx = scratch.symm_grad_Nx[q_point]; const double JxW = scratch.fe_values_ref.JxW(q_point); - // We first compute the contributions - // from the internal forces. Note, by - // definition of the rhs as the negative - // of the residual, these contributions - // are subtracted. for (unsigned int i = 0; i < dofs_per_cell; ++i) { - const unsigned int i_group = fe.system_to_base_index(i).first.first; + 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; if (i_group == u_dof) - data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW; + residual_ad[i] += (symm_grad_Nx[i] * tau) * JxW; else Assert(i_group <= u_dof, ExcInternalError()); } } - // Next we assemble the Neumann contribution. We first check to see it the - // cell face exists on a boundary on which a traction is applied and add - // the contribution if this is the case. - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; - ++face) - if (cell->face(face)->at_boundary() == true - && cell->face(face)->boundary_id() == 11) + for (unsigned int I=0; I dir; - dir[1] = 1.0; - const Tensor<1, dim> traction = magnitude*dir; - - 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 + struct Assembler > > : Assembler_Base > > + { + typedef Sacado::Fad::DFad ADDerivType; + typedef Sacado::Rad::ADvar ADNumberType; + using typename Assembler_Base::ScratchData_ASM; + using typename Assembler_Base::PerTaskData_ASM; + + virtual ~Assembler() {} + + virtual void + assemble_system_tangent_residual_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_ASM &scratch, + PerTaskData_ASM &data) + { + // Aliases for data referenced from the Solid class + const unsigned int &n_q_points = data.solid.n_q_points; + const unsigned int &dofs_per_cell = data.solid.dofs_per_cell; + const Parameters::AllParameters ¶meters = data.solid.parameters; + const FESystem &fe = data.solid.fe; + const unsigned int &u_dof = data.solid.u_dof; + const FEValuesExtractors::Vector &u_fe = data.solid.u_fe; + + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + + const std::vector > > lqph = + const_cast &>(data.solid).quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + + const unsigned int n_independent_variables = data.local_dof_indices.size(); + std::vector local_dof_values(n_independent_variables); + cell->get_dof_values(scratch.solution_total, + local_dof_values.begin(), + local_dof_values.end()); + + // We now retrieve a set of degree-of-freedom values that + // have the operations that are performed with them tracked. + std::vector local_dof_values_ad (n_independent_variables); + for (unsigned int i=0; i &grad_u = scratch.solution_grads_u_total[q_point]; + const Tensor<2,dim,ADNumberType> F = Physics::Elasticity::Kinematics::F(grad_u); + const ADNumberType det_F = determinant(F); + const Tensor<2,dim,ADNumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F); + const SymmetricTensor<2,dim,ADNumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar); + Assert(det_F > ADNumberType(0.0), ExcInternalError()); + + // Next we define some position-dependent aliases, again to + // make the assembly process easier to follow. + const double JxW = scratch.fe_values_ref.JxW(q_point); + + const ADNumberType Psi = lqph[q_point]->get_Psi(det_F,b_bar); + + // We extract the configuration-dependent material energy + // from our QPH history objects for the current quadrature point + // and integrate its contribution to increment the total + // cell energy. + cell_energy_ad += Psi * JxW; + } + + // Compute derivatives of reverse-mode AD variables + ADNumberType::Gradcomp(); + + for (unsigned int I=0; I + void Solid::assemble_system(const BlockVector &solution_delta) + { + timer.enter_subsection("Assemble linear system"); + std::cout << " ASM " << std::flush; + + tangent_matrix = 0.0; + system_rhs = 0.0; + + const UpdateFlags uf_cell(update_gradients | + update_JxW_values); + const UpdateFlags uf_face(update_values | + update_JxW_values); + + const BlockVector solution_total(get_total_solution(solution_delta)); + typename Assembler_Base::PerTaskData_ASM per_task_data(*this); + typename Assembler_Base::ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total); + Assembler assembler; + + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + static_cast&>(assembler), + &Assembler_Base::assemble_system_one_cell, + &Assembler_Base::copy_local_to_global_ASM, + scratch_data, + per_task_data); + + timer.leave_subsection(); } + // @sect4{Solid::make_constraints} // The constraints for this problem are simple to describe. // However, since we are dealing with an iterative Newton method, @@ -2138,8 +2093,8 @@ Point grid_y_transform (const Point &pt_in) // be specified at the zeroth iteration and subsequently no // additional contributions are to be made since the constraints // are already exactly satisfied. - template - void Solid::make_constraints(const int &it_nr) + template + void Solid::make_constraints(const int &it_nr) { std::cout << " CST " << std::flush; @@ -2214,9 +2169,9 @@ Point grid_y_transform (const Point &pt_in) // @sect4{Solid::solve_linear_system} // As the system is composed of a single block, defining a solution scheme // for the linear problem is straight-forward. - template + template std::pair - Solid::solve_linear_system(BlockVector &newton_update) + Solid::solve_linear_system(BlockVector &newton_update) { BlockVector A(dofs_per_block); BlockVector B(dofs_per_block); @@ -2287,8 +2242,8 @@ Point grid_y_transform (const Point &pt_in) // Here we present how the results are written to file to be viewed // using ParaView or Visit. The method is similar to that shown in the // tutorials so will not be discussed in detail. - template - void Solid::output_results() const + template + void Solid::output_results() const { DataOut data_out; std::vector @@ -2335,14 +2290,57 @@ int main (int argc, char *argv[]) using namespace dealii; using namespace Cook_Membrane; - Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, dealii::numbers::invalid_unsigned_int); - + const unsigned int dim = 2; try { deallog.depth_console(0); + Parameters::AllParameters parameters("parameters.prm"); + if (parameters.automatic_differentiation_order == 0) + { + std::cout << "Assembly method: Residual and linearisation are computed manually." << std::endl; - Solid<3> solid_3d("parameters.prm"); - solid_3d.run(); + // Allow multi-threading + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, + dealii::numbers::invalid_unsigned_int); + + typedef double NumberType; + Solid solid_3d(parameters); + solid_3d.run(); + } + #ifdef ENABLE_SACADO_FORMULATION + else if (parameters.automatic_differentiation_order == 1) + { + std::cout << "Assembly method: Residual computed manually; linearisation performed using AD." << std::endl; + + // Allow multi-threading + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, + dealii::numbers::invalid_unsigned_int); + + typedef Sacado::Fad::DFad NumberType; + Solid solid_3d(parameters); + solid_3d.run(); + } + else if (parameters.automatic_differentiation_order == 2) + { + std::cout << "Assembly method: Residual and linearisation computed using AD." << std::endl; + + // Sacado Rad-Fad is not thread-safe, so disable threading. + // Parallisation using MPI would be possible though. + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, + 1); + + typedef Sacado::Rad::ADvar< Sacado::Fad::DFad > NumberType; + Solid solid_3d(parameters); + solid_3d.run(); + } + #endif + else + { + AssertThrow(false, + ExcMessage("The selected assembly method is not supported. " + "You need deal.II 9.0 and Trilinos with the Sacado package " + "to enable assembly using automatic differentiation.")); + } } catch (std::exception &exc) { diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/doc/builds-on b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/builds-on index 6f6a518..6a43f4c 100644 --- a/Quasi_static_Finite_strain_Compressible_Elasticity/doc/builds-on +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/builds-on @@ -1 +1 @@ -step-18 step-44 +step-18 step-33 step-44 diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/parameters.prm b/Quasi_static_Finite_strain_Compressible_Elasticity/parameters.prm index 595f06a..c21fe8f 100644 --- a/Quasi_static_Finite_strain_Compressible_Elasticity/parameters.prm +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/parameters.prm @@ -1,5 +1,13 @@ # Listing of Parameters # --------------------- +subsection Assembly method + # The automatic differentiation order to be used in the assembly of the linear system. + # Order = 0: Both the residual and linearisation are computed manually. + # Order = 1: The residual is computed manually but the linearisation is performed using AD. + # Order = 2: Both the residual and linearisation are computed using AD. + set Automatic differentiation order = 0 +end + subsection Finite element system # Displacement system polynomial order set Polynomial degree = 1 -- 2.39.5