207 * We start by including all the necessary deal.II header files and some
C++
208 * related ones. They have been discussed in detail in previous tutorial
209 * programs, so you need only refer to past tutorials
for details.
212 * #include <deal.II/base/
function.h>
213 * #include <deal.II/base/parameter_handler.h>
214 * #include <deal.II/base/
point.h>
215 * #include <deal.II/base/quadrature_lib.h>
216 * #include <deal.II/base/symmetric_tensor.h>
217 * #include <deal.II/base/tensor.h>
218 * #include <deal.II/base/timer.h>
219 * #include <deal.II/base/work_stream.h>
220 * #include <deal.II/base/quadrature_point_data.h>
221 * #include <deal.II/base/
std_cxx11/shared_ptr.h>
223 * #include <deal.II/dofs/dof_renumbering.h>
224 * #include <deal.II/dofs/dof_tools.h>
226 * #include <deal.II/grid/grid_generator.h>
227 * #include <deal.II/grid/grid_tools.h>
228 * #include <deal.II/grid/grid_in.h>
229 * #include <deal.II/grid/tria.h>
230 * #include <deal.II/grid/tria_boundary_lib.h>
232 * #include <deal.II/fe/fe_dgp_monomial.h>
233 * #include <deal.II/fe/fe_q.h>
234 * #include <deal.II/fe/fe_system.h>
235 * #include <deal.II/fe/fe_tools.h>
236 * #include <deal.II/fe/fe_values.h>
237 * #include <deal.II/fe/mapping_q_eulerian.h>
238 * #include <deal.II/fe/mapping_q.h>
240 * #include <deal.II/lac/block_sparse_matrix.h>
241 * #include <deal.II/lac/block_vector.h>
242 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
243 * #include <deal.II/lac/full_matrix.h>
244 * #include <deal.II/lac/precondition_selector.h>
245 * #include <deal.II/lac/solver_cg.h>
246 * #include <deal.II/lac/sparse_direct.h>
247 * #include <deal.II/lac/constraint_matrix.h>
249 * #include <deal.II/numerics/data_out.h>
250 * #include <deal.II/numerics/vector_tools.h>
252 * #include <deal.II/base/config.h>
254 * #include <deal.II/differentiation/ad.h>
255 * #define ENABLE_SACADO_FORMULATION
260 * These must be included below the AD headers so that
261 * their math
functions are available
for use in the
262 * definition of tensors and kinematic quantities
265 * #include <deal.II/physics/elasticity/kinematics.h>
266 * #include <deal.II/physics/elasticity/standard_tensors.h>
268 * #include <iostream>
274 * We then stick everything that relates to
this tutorial program into a
275 *
namespace of its own, and import all the deal.II function and class names
279 *
namespace Cook_Membrane
286 * <a name=
"Runtimeparameters"></a>
287 * <h3>Run-time parameters</h3>
291 * There are several parameters that can be
set in the code so we
set up a
295 *
namespace Parameters
300 * <a name=
"Assemblymethod"></a>
301 * <h4>Assembly method</h4>
305 * Here we specify whether automatic differentiation is to be used to
assemble
306 * the linear system, and
if so then what order of differentiation is to be
310 *
struct AssemblyMethod
312 *
unsigned int automatic_differentiation_order;
324 * prm.enter_subsection(
"Assembly method");
326 * prm.declare_entry(
"Automatic differentiation order",
"0",
328 *
"The automatic differentiation order to be used in the assembly of the linear system.\n"
329 *
"# Order = 0: Both the residual and linearisation are computed manually.\n"
330 *
"# Order = 1: The residual is computed manually but the linearisation is performed using AD.\n"
331 *
"# Order = 2: Both the residual and linearisation are computed using AD.");
333 * prm.leave_subsection();
338 * prm.enter_subsection(
"Assembly method");
340 * automatic_differentiation_order = prm.get_integer(
"Automatic differentiation order");
342 * prm.leave_subsection();
348 * <a name=
"FiniteElementsystem"></a>
349 * <h4>Finite Element system</h4>
353 * Here we specify the polynomial order used to
approximate the solution.
354 * The quadrature order should be adjusted accordingly.
359 *
unsigned int poly_degree;
360 *
unsigned int quad_order;
372 * prm.enter_subsection(
"Finite element system");
374 * prm.declare_entry(
"Polynomial degree",
"2",
376 *
"Displacement system polynomial order");
378 * prm.declare_entry(
"Quadrature order",
"3",
380 *
"Gauss quadrature order");
382 * prm.leave_subsection();
387 * prm.enter_subsection(
"Finite element system");
389 * poly_degree = prm.get_integer(
"Polynomial degree");
390 * quad_order = prm.get_integer(
"Quadrature order");
392 * prm.leave_subsection();
398 * <a name=
"Geometry"></a>
403 * Make adjustments to the problem geometry and its discretisation.
408 *
unsigned int elements_per_edge;
420 * prm.enter_subsection(
"Geometry");
422 * prm.declare_entry(
"Elements per edge",
"32",
424 *
"Number of elements per long edge of the beam");
426 * prm.declare_entry(
"Grid scale",
"1e-3",
428 *
"Global grid scaling factor");
430 * prm.leave_subsection();
435 * prm.enter_subsection(
"Geometry");
437 * elements_per_edge = prm.get_integer(
"Elements per edge");
438 *
scale = prm.get_double(
"Grid scale");
440 * prm.leave_subsection();
446 * <a name=
"Materials"></a>
451 * We also need the shear modulus @f$ \mu @f$ and Poisson ration @f$ \nu @f$
for the
452 * neo-Hookean material.
469 * prm.enter_subsection(
"Material properties");
471 * prm.declare_entry(
"Poisson's ratio",
"0.3",
473 *
"Poisson's ratio");
475 * prm.declare_entry(
"Shear modulus",
"0.4225e6",
479 * prm.leave_subsection();
484 * prm.enter_subsection(
"Material properties");
486 * nu = prm.get_double(
"Poisson's ratio");
487 *
mu = prm.get_double(
"Shear modulus");
489 * prm.leave_subsection();
495 * <a name=
"Linearsolver"></a>
496 * <h4>Linear solver</h4>
500 * Next, we choose both solver and preconditioner settings. The use of an
501 * effective preconditioner is critical to ensure convergence when a large
502 * nonlinear motion occurs within a Newton increment.
505 *
struct LinearSolver
507 * std::string type_lin;
509 *
double max_iterations_lin;
510 * std::string preconditioner_type;
511 *
double preconditioner_relaxation;
522 * prm.enter_subsection(
"Linear solver");
524 * prm.declare_entry(
"Solver type",
"CG",
526 *
"Type of solver used to solve the linear system");
528 * prm.declare_entry(
"Residual",
"1e-6",
530 *
"Linear solver residual (scaled by residual norm)");
532 * prm.declare_entry(
"Max iteration multiplier",
"1",
534 *
"Linear solver iterations (multiples of the system matrix size)");
536 * prm.declare_entry(
"Preconditioner type",
"ssor",
538 *
"Type of preconditioner");
540 * prm.declare_entry(
"Preconditioner relaxation",
"0.65",
542 *
"Preconditioner relaxation value");
544 * prm.leave_subsection();
549 * prm.enter_subsection(
"Linear solver");
551 * type_lin = prm.get(
"Solver type");
552 * tol_lin = prm.get_double(
"Residual");
553 * max_iterations_lin = prm.get_double(
"Max iteration multiplier");
554 * preconditioner_type = prm.get(
"Preconditioner type");
555 * preconditioner_relaxation = prm.get_double(
"Preconditioner relaxation");
557 * prm.leave_subsection();
563 * <a name=
"Nonlinearsolver"></a>
564 * <h4>Nonlinear solver</h4>
568 *
A Newton-Raphson scheme is used to solve the nonlinear system of governing
569 * equations. We now define the tolerances and the maximum number of
570 * iterations
for the Newton-Raphson nonlinear solver.
573 *
struct NonlinearSolver
575 *
unsigned int max_iterations_NR;
588 * prm.enter_subsection(
"Nonlinear solver");
590 * prm.declare_entry(
"Max iterations Newton-Raphson",
"10",
592 *
"Number of Newton-Raphson iterations allowed");
594 * prm.declare_entry(
"Tolerance force",
"1.0e-9",
596 *
"Force residual tolerance");
598 * prm.declare_entry(
"Tolerance displacement",
"1.0e-6",
600 *
"Displacement error tolerance");
602 * prm.leave_subsection();
607 * prm.enter_subsection(
"Nonlinear solver");
609 * max_iterations_NR = prm.get_integer(
"Max iterations Newton-Raphson");
610 * tol_f = prm.get_double(
"Tolerance force");
611 * tol_u = prm.get_double(
"Tolerance displacement");
613 * prm.leave_subsection();
619 * <a name=
"Time"></a>
624 * Set the timestep size @f$ \varDelta t @f$ and the simulation
end-time.
641 * prm.enter_subsection(
"Time");
643 * prm.declare_entry(
"End time",
"1",
647 * prm.declare_entry(
"Time step size",
"0.1",
651 * prm.leave_subsection();
656 * prm.enter_subsection(
"Time");
658 * end_time = prm.get_double(
"End time");
659 * delta_t = prm.get_double(
"Time step size");
661 * prm.leave_subsection();
667 * <a name=
"Allparameters"></a>
668 * <h4>All parameters</h4>
672 * Finally we consolidate all of the above structures into a single container
673 * that holds all of our
run-time selections.
676 *
struct AllParameters :
677 *
public AssemblyMethod,
681 *
public LinearSolver,
682 *
public NonlinearSolver,
686 * AllParameters(
const std::string &input_file);
695 * AllParameters::AllParameters(
const std::string &input_file)
698 * declare_parameters(prm);
700 * parse_parameters(prm);
705 * AssemblyMethod::declare_parameters(prm);
706 * FESystem::declare_parameters(prm);
707 * Geometry::declare_parameters(prm);
708 * Materials::declare_parameters(prm);
709 * LinearSolver::declare_parameters(prm);
710 * NonlinearSolver::declare_parameters(prm);
711 * Time::declare_parameters(prm);
716 * AssemblyMethod::parse_parameters(prm);
717 * FESystem::parse_parameters(prm);
718 * Geometry::parse_parameters(prm);
719 * Materials::parse_parameters(prm);
720 * LinearSolver::parse_parameters(prm);
721 * NonlinearSolver::parse_parameters(prm);
722 * Time::parse_parameters(prm);
730 * <a name=
"Timeclass"></a>
731 * <h3>Time
class</h3>
735 *
A simple
class to store time data. Its functioning is transparent so no
736 * discussion is necessary. For simplicity we assume a constant time step
743 * Time (
const double time_end,
744 *
const double delta_t)
748 * time_end(time_end),
755 *
double current() const
757 *
return time_current;
763 *
double get_delta_t() const
767 *
unsigned int get_timestep() const
773 * time_current += delta_t;
778 *
unsigned int timestep;
779 *
double time_current;
780 *
const double time_end;
781 *
const double delta_t;
787 * <a name=
"CompressibleneoHookeanmaterialwithinaonefieldformulation"></a>
788 * <h3>Compressible neo-Hookean material within a
one-field formulation</h3>
792 * As discussed in the literature and @ref step_44
"step-44", Neo-Hookean materials are a type
793 * of hyperelastic materials. The entire domain is assumed to be composed of a
794 * compressible neo-Hookean material. This
class defines the behaviour of
795 *
this material within a
one-field formulation. Compressible neo-Hookean
796 * materials can be described by a strain-energy
function (SEF) @f$ \Psi =
797 * \Psi_{\text{iso}}(\overline{\mathbf{
b}}) + \Psi_{\text{vol}}(J)
802 * The isochoric response is given by @f$
803 * \Psi_{\text{iso}}(\overline{\mathbf{
b}}) = c_{1} [\overline{I}_{1} - 3] @f$
804 * where @f$ c_{1} = \frac{\mu}{2} @f$ and @f$\overline{I}_{1}@f$ is the
first
805 * invariant of the left- or right-isochoric Cauchy-Green deformation tensors.
806 * That is @f$\overline{I}_1 :=\textrm{tr}(\overline{\mathbf{
b}})@f$. In
this
807 * example the SEF that governs the volumetric response is defined as @f$
808 * \Psi_{\text{vol}}(J) = \kappa \frac{1}{4} [ J^2 - 1
809 * - 2\textrm{ln}\; J ]@f$, where @f$\kappa:= \lambda + 2/3 \mu@f$ is
810 * the <a href=
"http://en.wikipedia.org/wiki/Bulk_modulus">bulk modulus</a>
811 * and @f$\lambda@f$ is <a
812 * href=
"http://en.wikipedia.org/wiki/Lam%C3%A9_parameters">Lame
's first
817 * The following class will be used to characterize the material we work with,
818 * and provides a central point that one would need to modify if one were to
819 * implement a different material model. For it to work, we will store one
820 * object of this type per quadrature point, and in each of these objects
821 * store the current state (characterized by the values or measures of the
822 * displacement field) so that we can compute the elastic coefficients
823 * linearized around the current state.
826 * template <int dim,typename NumberType>
827 * class Material_Compressible_Neo_Hook_One_Field
830 * Material_Compressible_Neo_Hook_One_Field(const double mu,
833 * kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
836 * Assert(kappa > 0, ExcInternalError());
839 * ~Material_Compressible_Neo_Hook_One_Field()
844 * The first function is the total energy
845 * @f$\Psi = \Psi_{\textrm{iso}} + \Psi_{\textrm{vol}}@f$.
849 * get_Psi(const NumberType &det_F,
850 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
852 * return get_Psi_vol(det_F) + get_Psi_iso(b_bar);
857 * The second function determines the Kirchhoff stress @f$\boldsymbol{\tau}
858 * = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}@f$
861 * SymmetricTensor<2,dim,NumberType>
862 * get_tau(const NumberType &det_F,
863 * const SymmetricTensor<2,dim,NumberType> &b_bar)
867 * See Holzapfel p231 eq6.98 onwards
870 * return get_tau_vol(det_F) + get_tau_iso(b_bar);
875 * The fourth-order elasticity tensor in the spatial setting
876 * @f$\mathfrak{c}@f$ is calculated from the SEF @f$\Psi@f$ as @f$ J
877 * \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}@f$
878 * where @f$ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial
879 * \mathbf{C} \partial \mathbf{C}}@f$
882 * SymmetricTensor<4,dim,NumberType>
883 * get_Jc(const NumberType &det_F,
884 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
886 * return get_Jc_vol(det_F) + get_Jc_iso(b_bar);
892 * Define constitutive model parameters @f$\kappa@f$ (bulk modulus) and the
893 * neo-Hookean model parameter @f$c_1@f$:
896 * const double kappa;
901 * Value of the volumetric free energy
905 * get_Psi_vol(const NumberType &det_F) const
907 * return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
912 * Value of the isochoric free energy
916 * get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
918 * return c_1 * (trace(b_bar) - dim);
923 * Derivative of the volumetric free energy with respect to
924 * @f$J@f$ return @f$\frac{\partial
925 * \Psi_{\text{vol}}(J)}{\partial J}@f$
929 * get_dPsi_vol_dJ(const NumberType &det_F) const
931 * return (kappa / 2.0) * (det_F - 1.0 / det_F);
936 * The following functions are used internally in determining the result
937 * of some of the public functions above. The first one determines the
938 * volumetric Kirchhoff stress @f$\boldsymbol{\tau}_{\textrm{vol}}@f$.
939 * Note the difference in its definition when compared to @ref step_44 "step-44".
942 * SymmetricTensor<2,dim,NumberType>
943 * get_tau_vol(const NumberType &det_F) const
945 * return NumberType(get_dPsi_vol_dJ(det_F) * det_F) * Physics::Elasticity::StandardTensors<dim>::I;
950 * Next, determine the isochoric Kirchhoff stress
951 * @f$\boldsymbol{\tau}_{\textrm{iso}} =
952 * \mathcal{P}:\overline{\boldsymbol{\tau}}@f$:
955 * SymmetricTensor<2,dim,NumberType>
956 * get_tau_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
958 * return Physics::Elasticity::StandardTensors<dim>::dev_P * get_tau_bar(b_bar);
963 * Then, determine the fictitious Kirchhoff stress
964 * @f$\overline{\boldsymbol{\tau}}@f$:
967 * SymmetricTensor<2,dim,NumberType>
968 * get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const
970 * return 2.0 * c_1 * b_bar;
975 * Second derivative of the volumetric free energy wrt @f$J@f$. We
976 * need the following computation explicitly in the tangent so we make it
977 * public. We calculate @f$\frac{\partial^2
978 * \Psi_{\textrm{vol}}(J)}{\partial J \partial
983 * get_d2Psi_vol_dJ2(const NumberType &det_F) const
985 * return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
990 * Calculate the volumetric part of the tangent @f$J
991 * \mathfrak{c}_\textrm{vol}@f$. Again, note the difference in its
992 * definition when compared to @ref step_44 "step-44". The extra terms result from two
993 * quantities in @f$\boldsymbol{\tau}_{\textrm{vol}}@f$ being dependent on
994 * @f$\boldsymbol{F}@f$.
997 * SymmetricTensor<4,dim,NumberType>
998 * get_Jc_vol(const NumberType &det_F) const
1002 * See Holzapfel p265
1006 * * ( (get_dPsi_vol_dJ(det_F) + det_F * get_d2Psi_vol_dJ2(det_F))*Physics::Elasticity::StandardTensors<dim>::IxI
1007 * - (2.0 * get_dPsi_vol_dJ(det_F))*Physics::Elasticity::StandardTensors<dim>::S );
1012 * Calculate the isochoric part of the tangent @f$J
1013 * \mathfrak{c}_\textrm{iso}@f$:
1016 * SymmetricTensor<4,dim,NumberType>
1017 * get_Jc_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
1019 * const SymmetricTensor<2, dim> tau_bar = get_tau_bar(b_bar);
1020 * const SymmetricTensor<2, dim> tau_iso = get_tau_iso(b_bar);
1021 * const SymmetricTensor<4, dim> tau_iso_x_I
1022 * = outer_product(tau_iso,
1023 * Physics::Elasticity::StandardTensors<dim>::I);
1024 * const SymmetricTensor<4, dim> I_x_tau_iso
1025 * = outer_product(Physics::Elasticity::StandardTensors<dim>::I,
1027 * const SymmetricTensor<4, dim> c_bar = get_c_bar();
1029 * return (2.0 / dim) * trace(tau_bar)
1030 * * Physics::Elasticity::StandardTensors<dim>::dev_P
1031 * - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso)
1032 * + Physics::Elasticity::StandardTensors<dim>::dev_P * c_bar
1033 * * Physics::Elasticity::StandardTensors<dim>::dev_P;
1038 * Calculate the fictitious elasticity tensor @f$\overline{\mathfrak{c}}@f$.
1039 * For the material model chosen this is simply zero:
1042 * SymmetricTensor<4,dim,double>
1045 * return SymmetricTensor<4, dim>();
1052 * <a name="Quadraturepointhistory"></a>
1053 * <h3>Quadrature point history</h3>
1057 * As seen in @ref step_18 "step-18", the <code> PointHistory </code> class offers a method
1058 * for storing data at the quadrature points. Here each quadrature point
1059 * holds a pointer to a material description. Thus, different material models
1060 * can be used in different regions of the domain. Among other data, we
1061 * choose to store the Kirchhoff stress @f$\boldsymbol{\tau}@f$ and the tangent
1062 * @f$J\mathfrak{c}@f$ for the quadrature points.
1065 * template <int dim,typename NumberType>
1066 * class PointHistory
1072 * virtual ~PointHistory()
1077 * The first function is used to create a material object and to
1078 * initialize all tensors correctly: The second one updates the stored
1079 * values and stresses based on the current deformation measure
1080 * @f$\textrm{Grad}\mathbf{u}_{\textrm{n}}@f$.
1083 * void setup_lqp (const Parameters::AllParameters ¶meters)
1085 * material.reset(new Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu,
1091 * We offer an interface to retrieve certain data.
1092 * This is the strain energy:
1096 * get_Psi(const NumberType &det_F,
1097 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1099 * return material->get_Psi(det_F,b_bar);
1104 * Here are the kinetic variables. These are used in the material and
1105 * global tangent matrix and residual assembly operations:
1106 * First is the Kirchhoff stress:
1109 * SymmetricTensor<2,dim,NumberType>
1110 * get_tau(const NumberType &det_F,
1111 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1113 * return material->get_tau(det_F,b_bar);
1121 * SymmetricTensor<4,dim,NumberType>
1122 * get_Jc(const NumberType &det_F,
1123 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1125 * return material->get_Jc(det_F,b_bar);
1130 * In terms of member functions, this class stores for the quadrature
1131 * point it represents a copy of a material type in case different
1132 * materials are used in different regions of the domain, as well as the
1133 * inverse of the deformation gradient...
1137 * std::shared_ptr< Material_Compressible_Neo_Hook_One_Field<dim,NumberType> > material;
1144 * <a name="Quasistaticcompressiblefinitestrainsolid"></a>
1145 * <h3>Quasi-static compressible finite-strain solid</h3>
1149 * Forward declarations for classes that will
1150 * perform assembly of the linear system.
1153 * template <int dim,typename NumberType>
1154 * struct Assembler_Base;
1155 * template <int dim,typename NumberType>
1160 * The Solid class is the central class in that it represents the problem at
1161 * hand. It follows the usual scheme in that all it really has is a
1162 * constructor, destructor and a <code>run()</code> function that dispatches
1163 * all the work to private functions of this class:
1166 * template <int dim,typename NumberType>
1170 * Solid(const Parameters::AllParameters ¶meters);
1182 * We start the collection of member functions with one that builds the
1191 * Set up the finite element system to be solved:
1199 * Several functions to assemble the system and right hand side matrices
1200 * using multithreading. Each of them comes as a wrapper function, one
1201 * that is executed to do the work in the WorkStream model on one cell,
1202 * and one that copies the work done on this one cell into the global
1203 * object that represents it:
1207 * assemble_system(const BlockVector<double> &solution_delta);
1211 * We use a separate data structure to perform the assembly. It needs access
1212 * to some low-level data, so we simply befriend the class instead of
1213 * creating a complex interface to provide access as necessary.
1216 * friend struct Assembler_Base<dim,NumberType>;
1217 * friend struct Assembler<dim,NumberType>;
1221 * Apply Dirichlet boundary conditions on the displacement field
1225 * make_constraints(const int &it_nr);
1229 * Create and update the quadrature points. Here, no data needs to be
1230 * copied into a global object, so the copy_local_to_global function is
1239 * Solve for the displacement using a Newton-Raphson method. We break this
1240 * function into the nonlinear loop and the function that solves the
1241 * linearized Newton-Raphson step:
1245 * solve_nonlinear_timestep(BlockVector<double> &solution_delta);
1247 * std::pair<unsigned int, double>
1248 * solve_linear_system(BlockVector<double> &newton_update);
1252 * Solution retrieval as well as post-processing and writing data to file :
1255 * BlockVector<double>
1256 * get_total_solution(const BlockVector<double> &solution_delta) const;
1259 * output_results() const;
1263 * Finally, some member variables that describe the current state: A
1264 * collection of the parameters used to describe the problem setup...
1267 * const Parameters::AllParameters ¶meters;
1271 * ...the volume of the reference and current configurations...
1274 * double vol_reference;
1275 * double vol_current;
1279 * ...and description of the geometry on which the problem is solved:
1282 * Triangulation<dim> triangulation;
1286 * Also, keep track of the current time and the time spent evaluating
1291 * TimerOutput timer;
1295 * A storage object for quadrature point information. As opposed to
1296 * @ref step_18 "step-18", deal.II's native quadrature
point data manager is employed here.
1300 * PointHistory<dim,NumberType> > quadrature_point_history;
1304 *
A description of the finite-element system including the displacement
1305 * polynomial degree, the degree-of-freedom handler, number of DoFs per
1306 * cell and the extractor objects used to retrieve information from the
1310 *
const unsigned int degree;
1313 *
const unsigned int dofs_per_cell;
1318 * Description of how the block-system is arranged. There is just 1 block,
1319 * that contains a vector DOF @f$\mathbf{u}@f$.
1320 * There are two reasons that we retain the block system in
this problem.
1321 * The
first is pure laziness to perform further modifications to the
1322 * code from which
this work originated. The
second is that a block system
1323 * would typically necessary when extending
this code to multiphysics
1327 *
static const unsigned int n_blocks = 1;
1329 *
static const unsigned int first_u_component = 0;
1336 * std::vector<types::global_dof_index> dofs_per_block;
1340 * Rules
for Gauss-quadrature on both the cell and faces. The number of
1341 * quadrature points on both cells and faces is recorded.
1345 *
const QGauss<dim - 1> qf_face;
1346 *
const unsigned int n_q_points;
1347 *
const unsigned int n_q_points_f;
1351 * Objects that store the converged solution and right-hand side vectors,
1353 * to keep track of constraints. We make use of a sparsity pattern
1354 * designed
for a block system.
1365 * Then define a number of variables to store norms and update norms and
1366 * normalisation factors.
1381 *
void normalise(
const Errors &rhs)
1383 *
if (rhs.norm != 0.0)
1392 * Errors error_residual, error_residual_0, error_residual_norm, error_update,
1393 * error_update_0, error_update_norm;
1397 * Methods to calculate error measures
1401 * get_error_residual(Errors &error_residual);
1405 * Errors &error_update);
1409 * Print information to screen in a pleasing way...
1414 * print_conv_header();
1417 * print_conv_footer();
1420 * print_vertical_tip_displacement();
1426 * <a name=
"ImplementationofthecodeSolidcodeclass"></a>
1427 * <h3>Implementation of the <code>Solid</code>
class</h3>
1432 * <a name=
"Publicinterface"></a>
1433 * <h4>Public interface</h4>
1437 * We initialise the Solid
class using data extracted from the parameter file.
1440 *
template <
int dim,
typename NumberType>
1441 * Solid<dim,NumberType>::Solid(
const Parameters::AllParameters ¶meters)
1443 * parameters(parameters),
1444 * vol_reference (0.0),
1445 * vol_current (0.0),
1447 * time(parameters.end_time, parameters.delta_t),
1451 * degree(parameters.poly_degree),
1454 * The Finite Element System is composed of dim continuous displacement
1458 * fe(
FE_Q<dim>(parameters.poly_degree), dim),
1460 * dofs_per_cell (fe.dofs_per_cell),
1461 * u_fe(first_u_component),
1463 * qf_cell(parameters.quad_order),
1464 * qf_face(parameters.quad_order),
1465 * n_q_points (qf_cell.size()),
1466 * n_q_points_f (qf_face.size())
1473 * The
class destructor simply clears the data held by the DOFHandler
1476 *
template <
int dim,
typename NumberType>
1477 * Solid<dim,NumberType>::~Solid()
1479 * dof_handler_ref.
clear();
1485 * In solving the quasi-
static problem, the time becomes a loading parameter,
1486 * i.e. we increasing the loading linearly with time, making the two concepts
1487 * interchangeable. We choose to increment time linearly
using a constant time
1492 * We start the
function with preprocessing, and then output the
initial grid
1493 * before starting the simulation proper with the
first time (and loading)
1500 *
template <
int dim,
typename NumberType>
1510 * We then declare the incremental solution update @f$\varDelta
1511 * \mathbf{\Xi}:= \{\varDelta \mathbf{u}\}@f$ and start the
loop over the
1516 * At the beginning, we reset the solution update
for this time step...
1520 *
while (time.current() <= time.end())
1522 * solution_delta = 0.0;
1526 * ...solve the current time step and update total solution vector
1527 * @f$\mathbf{\Xi}_{\textrm{n}} = \mathbf{\Xi}_{\textrm{n-1}} +
1528 * \varDelta \mathbf{\Xi}@f$...
1531 * solve_nonlinear_timestep(solution_delta);
1532 * solution_n += solution_delta;
1536 * ...and plot the results before moving on happily to the next time
1546 * Lastly, we print the vertical tip displacement of the Cook cantilever
1547 * after the full load is applied
1550 * print_vertical_tip_displacement();
1557 * <a name=
"Privateinterface"></a>
1558 * <h3>Private interface</h3>
1563 * <a name=
"Solidmake_grid"></a>
1564 * <h4>Solid::make_grid</h4>
1568 * On to the
first of the
private member
functions. Here we create the
1569 *
triangulation of the domain,
for which we choose a scaled an anisotripically
1570 * discretised rectangle which is subsequently transformed into the correct
1571 * of the Cook cantilever. Each relevant boundary face is then given a boundary
1576 * We then determine the
volume of the reference configuration and print it
1583 *
template <
int dim>
1586 *
const double &x = pt_in[0];
1587 *
const double &y = pt_in[1];
1589 *
const double y_upper = 44.0 + (16.0/48.0)*x;
1590 *
const double y_lower = 0.0 + (44.0/48.0)*x;
1591 *
const double theta = y/44.0;
1592 *
const double y_transform = (1-theta)*y_lower + theta*y_upper;
1595 * pt_out[1] = y_transform;
1600 *
template <
int dim,
typename NumberType>
1601 *
void Solid<dim,NumberType>::make_grid()
1605 * Divide the beam, but only along the x- and y-coordinate directions
1608 * std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge);
1611 * Only allow
one element through the thickness
1612 * (modelling a plane strain condition)
1616 * repetitions[dim-1] = 1;
1628 * Since we wish to
apply a Neumann BC to the right-hand surface, we
1629 * must find the cell faces in
this part of the domain and mark them with
1630 * a distinct boundary ID number. The faces we are looking
for are on the
1631 * +x surface and will get boundary ID 11.
1632 * Dirichlet boundaries exist on the left-hand face of the beam (
this fixed
1633 * boundary will get ID 1) and on the +Z and -Z faces (which correspond to
1634 * ID 2 and we will use to impose the plane strain condition)
1637 *
const double tol_boundary = 1
e-6;
1640 *
for (; cell != endc; ++cell)
1641 *
for (
unsigned int face = 0;
1642 * face < GeometryInfo<dim>::faces_per_cell; ++face)
1643 *
if (cell->face(face)->at_boundary() ==
true)
1645 *
if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary)
1646 * cell->face(face)->set_boundary_id(1);
1647 *
else if (std::abs(cell->face(face)->center()[0] - 48.0) < tol_boundary)
1648 * cell->face(face)->set_boundary_id(11);
1649 *
else if (std::abs(std::abs(cell->face(face)->center()[0]) - 0.5) < tol_boundary)
1650 * cell->face(face)->set_boundary_id(2);
1655 * Transform the hyper-rectangle into the beam shape
1663 * vol_current = vol_reference;
1664 * std::cout <<
"Grid:\n\t Reference volume: " << vol_reference << std::endl;
1671 * <a name=
"Solidsystem_setup"></a>
1672 * <h4>Solid::system_setup</h4>
1676 * Next we describe how the FE system is setup. We
first determine the number
1677 * of components per block. Since the displacement is a vector component, the
1678 *
first dim components belong to it.
1681 *
template <
int dim,
typename NumberType>
1682 *
void Solid<dim,NumberType>::system_setup()
1684 * timer.enter_subsection(
"Setup system");
1686 * std::vector<unsigned int> block_component(
n_components, u_dof);
1690 * The DOF handler is then initialised and we renumber the grid in an
1691 * efficient manner. We also record the number of DOFs per block.
1700 * std::cout <<
"Triangulation:"
1701 * <<
"\n\t Number of active cells: " <<
triangulation.n_active_cells()
1702 * <<
"\n\t Number of degrees of freedom: " << dof_handler_ref.
n_dofs()
1707 * Setup the sparsity pattern and tangent
matrix
1710 * tangent_matrix.
clear();
1716 * csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
1717 * csp.collect_sizes();
1721 * Naturally,
for a
one-field vector-valued problem, all of the
1722 * components of the system are coupled.
1737 * tangent_matrix.
reinit(sparsity_pattern);
1741 * We then
set up storage vectors
1744 * system_rhs.
reinit(dofs_per_block);
1747 * solution_n.
reinit(dofs_per_block);
1752 * ...and
finally set up the quadrature
1758 * timer.leave_subsection();
1765 * <a name=
"Solidsetup_qph"></a>
1766 * <h4>Solid::setup_qph</h4>
1767 * The method used to store quadrature information is already described in
1768 * @ref step_18
"step-18" and @ref step_44
"step-44". Here we implement a similar setup
for a SMP machine.
1772 * Firstly the actual QPH data objects are created. This must be done only
1773 * once the grid is refined to its finest
level.
1776 *
template <
int dim,
typename NumberType>
1777 *
void Solid<dim,NumberType>::setup_qph()
1779 * std::cout <<
" Setting up quadrature point data..." << std::endl;
1781 * quadrature_point_history.initialize(
triangulation.begin_active(),
1787 * Next we setup the
initial quadrature
point data. Note that when
1788 * the quadrature
point data is retrieved, it is returned as a vector
1789 * of smart pointers.
1795 *
const std::vector<std::shared_ptr<PointHistory<dim,NumberType> > > lqph =
1796 * quadrature_point_history.get_data(cell);
1799 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1800 * lqph[q_point]->setup_lqp(parameters);
1808 * <a name=
"Solidsolve_nonlinear_timestep"></a>
1809 * <h4>Solid::solve_nonlinear_timestep</h4>
1813 * The next
function is the driver method
for the Newton-Raphson scheme. At
1814 * its top we create a
new vector to store the current Newton update step,
1815 * reset the error storage objects and print solver header.
1818 *
template <
int dim,
typename NumberType>
1822 * std::cout << std::endl <<
"Timestep " << time.get_timestep() <<
" @ "
1823 * << time.current() <<
"s" << std::endl;
1827 * error_residual.reset();
1828 * error_residual_0.reset();
1829 * error_residual_norm.reset();
1830 * error_update.reset();
1831 * error_update_0.reset();
1832 * error_update_norm.reset();
1834 * print_conv_header();
1838 * We now perform a number of Newton iterations to iteratively solve the
1839 * nonlinear problem. Since the problem is fully nonlinear and we are
1840 *
using a full Newton method, the data stored in the tangent
matrix and
1841 * right-hand side vector is not reusable and must be cleared at each
1842 * Newton step. We then initially build the right-hand side vector to
1843 * check
for convergence (and store
this value in the
first iteration).
1844 * The unconstrained DOFs of the rhs vector hold the out-of-balance
1845 * forces. The building is done before assembling the system
matrix as the
1846 * latter is an expensive operation and we can potentially avoid an extra
1847 * assembly process by not assembling the tangent
matrix when convergence
1851 *
unsigned int newton_iteration = 0;
1852 *
for (; newton_iteration < parameters.max_iterations_NR;
1853 * ++newton_iteration)
1855 * std::cout <<
" " << std::setw(2) << newton_iteration <<
" " << std::flush;
1859 * If we have decided that we want to
continue with the iteration, we
1860 *
assemble the tangent, make and impose the Dirichlet constraints,
1861 * and
do the solve of the linearized system:
1864 * make_constraints(newton_iteration);
1865 * assemble_system(solution_delta);
1867 * get_error_residual(error_residual);
1869 *
if (newton_iteration == 0)
1870 * error_residual_0 = error_residual;
1874 * We can now determine the normalised residual error and check
for
1875 * solution convergence:
1878 * error_residual_norm = error_residual;
1879 * error_residual_norm.normalise(error_residual_0);
1881 *
if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u
1882 * && error_residual_norm.u <= parameters.tol_f)
1884 * std::cout <<
" CONVERGED! " << std::endl;
1885 * print_conv_footer();
1890 *
const std::pair<unsigned int, double>
1891 * lin_solver_output = solve_linear_system(newton_update);
1893 * get_error_update(newton_update, error_update);
1894 *
if (newton_iteration == 0)
1895 * error_update_0 = error_update;
1899 * We can now determine the normalised Newton update error, and
1900 * perform the actual update of the solution increment
for the current
1901 * time step, update all quadrature
point information pertaining to
1902 *
this new displacement and stress state and
continue iterating:
1905 * error_update_norm = error_update;
1906 * error_update_norm.normalise(error_update_0);
1908 * solution_delta += newton_update;
1910 * std::cout <<
" | " << std::fixed << std::setprecision(3) << std::setw(7)
1911 * << std::scientific << lin_solver_output.first <<
" "
1912 * << lin_solver_output.second <<
" " << error_residual_norm.
norm
1913 * <<
" " << error_residual_norm.u <<
" "
1914 * <<
" " << error_update_norm.norm <<
" " << error_update_norm.u
1915 * <<
" " << std::endl;
1920 * At the
end,
if it turns out that we have in fact done more iterations
1921 * than the parameter file allowed, we
raise an exception that can be
1922 * caught in the main() function. The
call <code>
AssertThrow(condition,
1923 * exc_object)</code> is in essence equivalent to <code>if (!cond) throw
1924 * exc_object;</code> but the former form fills certain fields in the
1925 * exception
object that identify the location (filename and line number)
1926 * where the exception was raised to make it simpler to identify where the
1930 *
AssertThrow (newton_iteration <= parameters.max_iterations_NR,
1931 *
ExcMessage("No convergence in nonlinear solver!"));
1938 * <a name="Solidprint_conv_headerSolidprint_conv_footerandSolidprint_vertical_tip_displacement"></a>
1939 * <h4>Solid::print_conv_header, Solid::print_conv_footer and Solid::print_vertical_tip_displacement</h4>
1943 * This program prints out data in a nice table that is updated
1944 * on a per-iteration basis. The next two
functions set up the table
1945 * header and footer:
1948 * template <
int dim,typename NumberType>
1949 *
void Solid<dim,NumberType>::print_conv_header()
1951 *
static const unsigned int l_width = 87;
1953 *
for (
unsigned int i = 0; i < l_width; ++i)
1955 * std::cout << std::endl;
1957 * std::cout <<
" SOLVER STEP "
1958 * <<
" | LIN_IT LIN_RES RES_NORM "
1959 * <<
" RES_U NU_NORM "
1960 * <<
" NU_U " << std::endl;
1962 *
for (
unsigned int i = 0; i < l_width; ++i)
1964 * std::cout << std::endl;
1969 *
template <
int dim,
typename NumberType>
1970 *
void Solid<dim,NumberType>::print_conv_footer()
1972 *
static const unsigned int l_width = 87;
1974 *
for (
unsigned int i = 0; i < l_width; ++i)
1976 * std::cout << std::endl;
1978 * std::cout <<
"Relative errors:" << std::endl
1979 * <<
"Displacement:\t" << error_update.u / error_update_0.u << std::endl
1980 * <<
"Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
1981 * <<
"v / V_0:\t" << vol_current <<
" / " << vol_reference
1987 * At the
end we also output the result that can be compared to that found in
1988 * the literature, namely the displacement at the upper right corner of the
1992 *
template <
int dim,
typename NumberType>
1993 *
void Solid<dim,NumberType>::print_vertical_tip_displacement()
1995 *
static const unsigned int l_width = 87;
1997 *
for (
unsigned int i = 0; i < l_width; ++i)
1999 * std::cout << std::endl;
2001 *
Point<dim> soln_pt (48.0*parameters.scale,60.0*parameters.scale);
2003 * soln_pt[2] = 0.5*parameters.scale;
2004 *
double vertical_tip_displacement = 0.0;
2005 *
double vertical_tip_displacement_check = 0.0;
2009 *
for (; cell != endc; ++cell)
2013 *
if (cell->point_inside(soln_pt) ==
true)
2016 *
for (
unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
2017 *
if (cell->vertex(v).distance(soln_pt) < 1
e-6)
2021 * Extract y-component of solution at the given
point
2022 * This
point is coindicent with a vertex, so we can
2023 *
extract it directly as we
're using FE_Q finite elements
2024 * that have support at the vertices
2027 * vertical_tip_displacement = solution_n(cell->vertex_dof_index(v,u_dof+1));
2031 * Sanity check using alternate method to extract the solution
2032 * at the given point. To do this, we must create an FEValues instance
2033 * to help us extract the solution value at the desired point
2036 * const MappingQ<dim> mapping (parameters.poly_degree);
2037 * const Point<dim> qp_unit = mapping.transform_real_to_unit_cell(cell,soln_pt);
2038 * const Quadrature<dim> soln_qrule (qp_unit);
2039 * AssertThrow(soln_qrule.size() == 1, ExcInternalError());
2040 * FEValues<dim> fe_values_soln (fe, soln_qrule, update_values);
2041 * fe_values_soln.reinit(cell);
2045 * Extract y-component of solution at given point
2048 * std::vector< Tensor<1,dim> > soln_values (soln_qrule.size());
2049 * fe_values_soln[u_fe].get_function_values(solution_n,
2051 * vertical_tip_displacement_check = soln_values[0][u_dof+1];
2056 * AssertThrow(vertical_tip_displacement > 0.0, ExcMessage("Found no cell with point inside!"))
2058 * std::cout << "Vertical tip displacement: " << vertical_tip_displacement
2059 * << "\t Check: " << vertical_tip_displacement_check
2067 * <a name="Solidget_error_residual"></a>
2068 * <h4>Solid::get_error_residual</h4>
2072 * Determine the true residual error for the problem. That is, determine the
2073 * error in the residual for the unconstrained degrees of freedom. Note that to
2074 * do so, we need to ignore constrained DOFs by setting the residual in these
2075 * vector components to zero.
2078 * template <int dim,typename NumberType>
2079 * void Solid<dim,NumberType>::get_error_residual(Errors &error_residual)
2081 * BlockVector<double> error_res(dofs_per_block);
2083 * for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
2084 * if (!constraints.is_constrained(i))
2085 * error_res(i) = system_rhs(i);
2087 * error_residual.norm = error_res.l2_norm();
2088 * error_residual.u = error_res.block(u_dof).l2_norm();
2095 * <a name="Solidget_error_udpate"></a>
2096 * <h4>Solid::get_error_udpate</h4>
2100 * Determine the true Newton update error for the problem
2103 * template <int dim,typename NumberType>
2104 * void Solid<dim,NumberType>::get_error_update(const BlockVector<double> &newton_update,
2105 * Errors &error_update)
2107 * BlockVector<double> error_ud(dofs_per_block);
2108 * for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
2109 * if (!constraints.is_constrained(i))
2110 * error_ud(i) = newton_update(i);
2112 * error_update.norm = error_ud.l2_norm();
2113 * error_update.u = error_ud.block(u_dof).l2_norm();
2121 * <a name="Solidget_total_solution"></a>
2122 * <h4>Solid::get_total_solution</h4>
2126 * This function provides the total solution, which is valid at any Newton step.
2127 * This is required as, to reduce computational error, the total solution is
2128 * only updated at the end of the timestep.
2131 * template <int dim,typename NumberType>
2132 * BlockVector<double>
2133 * Solid<dim,NumberType>::get_total_solution(const BlockVector<double> &solution_delta) const
2135 * BlockVector<double> solution_total(solution_n);
2136 * solution_total += solution_delta;
2137 * return solution_total;
2144 * <a name="Solidassemble_system"></a>
2145 * <h4>Solid::assemble_system</h4>
2151 * template <int dim,typename NumberType>
2152 * struct Assembler_Base
2154 * virtual ~Assembler_Base() {}
2158 * Here we deal with the tangent matrix assembly structures. The
2159 * PerTaskData object stores local contributions.
2162 * struct PerTaskData_ASM
2164 * const Solid<dim,NumberType> *solid;
2165 * FullMatrix<double> cell_matrix;
2166 * Vector<double> cell_rhs;
2167 * std::vector<types::global_dof_index> local_dof_indices;
2169 * PerTaskData_ASM(const Solid<dim,NumberType> *solid)
2172 * cell_matrix(solid->dofs_per_cell, solid->dofs_per_cell),
2173 * cell_rhs(solid->dofs_per_cell),
2174 * local_dof_indices(solid->dofs_per_cell)
2179 * cell_matrix = 0.0;
2186 * On the other hand, the ScratchData object stores the larger objects such as
2187 * the shape-function values array (<code>Nx</code>) and a shape function
2188 * gradient and symmetric gradient vector which we will use during the
2192 * struct ScratchData_ASM
2194 * const BlockVector<double> &solution_total;
2195 * std::vector<Tensor<2, dim,NumberType> > solution_grads_u_total;
2197 * FEValues<dim> fe_values_ref;
2198 * FEFaceValues<dim> fe_face_values_ref;
2200 * std::vector<std::vector<Tensor<2, dim,NumberType> > > grad_Nx;
2201 * std::vector<std::vector<SymmetricTensor<2,dim,NumberType> > > symm_grad_Nx;
2203 * ScratchData_ASM(const FiniteElement<dim> &fe_cell,
2204 * const QGauss<dim> &qf_cell,
2205 * const UpdateFlags uf_cell,
2206 * const QGauss<dim-1> & qf_face,
2207 * const UpdateFlags uf_face,
2208 * const BlockVector<double> &solution_total)
2210 * solution_total(solution_total),
2211 * solution_grads_u_total(qf_cell.size()),
2212 * fe_values_ref(fe_cell, qf_cell, uf_cell),
2213 * fe_face_values_ref(fe_cell, qf_face, uf_face),
2214 * grad_Nx(qf_cell.size(),
2215 * std::vector<Tensor<2,dim,NumberType> >(fe_cell.dofs_per_cell)),
2216 * symm_grad_Nx(qf_cell.size(),
2217 * std::vector<SymmetricTensor<2,dim,NumberType> >
2218 * (fe_cell.dofs_per_cell))
2221 * ScratchData_ASM(const ScratchData_ASM &rhs)
2223 * solution_total (rhs.solution_total),
2224 * solution_grads_u_total(rhs.solution_grads_u_total),
2225 * fe_values_ref(rhs.fe_values_ref.get_fe(),
2226 * rhs.fe_values_ref.get_quadrature(),
2227 * rhs.fe_values_ref.get_update_flags()),
2228 * fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
2229 * rhs.fe_face_values_ref.get_quadrature(),
2230 * rhs.fe_face_values_ref.get_update_flags()),
2231 * grad_Nx(rhs.grad_Nx),
2232 * symm_grad_Nx(rhs.symm_grad_Nx)
2237 * const unsigned int n_q_points = fe_values_ref.get_quadrature().size();
2238 * const unsigned int n_dofs_per_cell = fe_values_ref.dofs_per_cell;
2239 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2241 * Assert( grad_Nx[q_point].size() == n_dofs_per_cell,
2242 * ExcInternalError());
2243 * Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell,
2244 * ExcInternalError());
2246 * solution_grads_u_total[q_point] = Tensor<2,dim,NumberType>();
2247 * for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
2249 * grad_Nx[q_point][k] = Tensor<2,dim,NumberType>();
2250 * symm_grad_Nx[q_point][k] = SymmetricTensor<2,dim,NumberType>();
2259 * Of course, we still have to define how we assemble the tangent matrix
2260 * contribution for a single cell.
2264 * assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2265 * ScratchData_ASM &scratch,
2266 * PerTaskData_ASM &data)
2270 * Due to the C++ specialization rules, we need one more
2271 * level of indirection in order to define the assembly
2272 * routine for all different number. The next function call
2273 * is specialized for each NumberType, but to prevent having
2274 * to specialize the whole class along with it we have inlined
2275 * the definition of the other functions that are common to
2276 * all implementations.
2279 * assemble_system_tangent_residual_one_cell(cell, scratch, data);
2280 * assemble_neumann_contribution_one_cell(cell, scratch, data);
2285 * This function adds the local contribution to the system matrix.
2289 * copy_local_to_global_ASM(const PerTaskData_ASM &data)
2291 * const ConstraintMatrix &constraints = data.solid->constraints;
2292 * BlockSparseMatrix<double> &tangent_matrix = const_cast<Solid<dim,NumberType> *>(data.solid)->tangent_matrix;
2293 * BlockVector<double> &system_rhs = const_cast<Solid<dim,NumberType> *>(data.solid)->system_rhs;
2295 * constraints.distribute_local_to_global(
2296 * data.cell_matrix, data.cell_rhs,
2297 * data.local_dof_indices,
2298 * tangent_matrix, system_rhs);
2305 * This function needs to exist in the base class for
2306 * Workstream to work with a reference to the base class.
2310 * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &/*cell*/,
2311 * ScratchData_ASM &/*scratch*/,
2312 * PerTaskData_ASM &/*data*/)
2314 * AssertThrow(false, ExcPureFunctionCalled());
2318 * assemble_neumann_contribution_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2319 * ScratchData_ASM &scratch,
2320 * PerTaskData_ASM &data)
2324 * Aliases for data referenced from the Solid class
2327 * const unsigned int &n_q_points_f = data.solid->n_q_points_f;
2328 * const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2329 * const Parameters::AllParameters ¶meters = data.solid->parameters;
2330 * const Time &time = data.solid->time;
2331 * const FESystem<dim> &fe = data.solid->fe;
2332 * const unsigned int &u_dof = data.solid->u_dof;
2336 * Next we assemble the Neumann contribution. We first check to see it the
2337 * cell face exists on a boundary on which a traction is applied and add
2338 * the contribution if this is the case.
2341 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2343 * if (cell->face(face)->at_boundary() == true
2344 * && cell->face(face)->boundary_id() == 11)
2346 * scratch.fe_face_values_ref.reinit(cell, face);
2348 * for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
2353 * We specify the traction in reference configuration.
2354 * For this problem, a defined total vertical force is applied
2355 * in the reference configuration.
2356 * The direction of the applied traction is assumed not to
2357 * evolve with the deformation of the domain.
2361 * Note that the contributions to the right hand side vector we
2362 * compute here only exist in the displacement components of the
2366 * const double time_ramp = (time.current() / time.end());
2367 * const double magnitude = (1.0/(16.0*parameters.scale*1.0*parameters.scale))*time_ramp; // (Total force) / (RHS surface area)
2368 * Tensor<1,dim> dir;
2370 * const Tensor<1, dim> traction = magnitude*dir;
2372 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2374 * const unsigned int i_group =
2375 * fe.system_to_base_index(i).first.first;
2377 * if (i_group == u_dof)
2379 * const unsigned int component_i =
2380 * fe.system_to_component_index(i).first;
2382 * scratch.fe_face_values_ref.shape_value(i,
2384 * const double JxW = scratch.fe_face_values_ref.JxW(
2387 * data.cell_rhs(i) += (Ni * traction[component_i])
2397 * template <int dim>
2398 * struct Assembler<dim,double> : Assembler_Base<dim,double>
2400 * typedef double NumberType;
2401 * using typename Assembler_Base<dim,NumberType>::ScratchData_ASM;
2402 * using typename Assembler_Base<dim,NumberType>::PerTaskData_ASM;
2404 * virtual ~Assembler() {}
2407 * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2408 * ScratchData_ASM &scratch,
2409 * PerTaskData_ASM &data)
2413 * Aliases for data referenced from the Solid class
2416 * const unsigned int &n_q_points = data.solid->n_q_points;
2417 * const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2418 * const FESystem<dim> &fe = data.solid->fe;
2419 * const unsigned int &u_dof = data.solid->u_dof;
2420 * const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
2424 * scratch.fe_values_ref.reinit(cell);
2425 * cell->get_dof_indices(data.local_dof_indices);
2427 * const std::vector<std::shared_ptr<const PointHistory<dim,NumberType> > > lqph =
2428 * const_cast<const Solid<dim,NumberType> *>(data.solid)->quadrature_point_history.get_data(cell);
2429 * Assert(lqph.size() == n_q_points, ExcInternalError());
2433 * We first need to find the solution gradients at quadrature points
2434 * inside the current cell and then we update each local QP using the
2435 * displacement gradient:
2438 * scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
2439 * scratch.solution_grads_u_total);
2443 * Now we build the local cell stiffness matrix. Since the global and
2444 * local system matrices are symmetric, we can exploit this property by
2445 * building only the lower half of the local matrix and copying the values
2446 * to the upper half.
2450 * In doing so, we first extract some configuration dependent variables
2451 * from our QPH history objects for the current quadrature point.
2454 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2456 * const Tensor<2,dim,NumberType> &grad_u = scratch.solution_grads_u_total[q_point];
2457 * const Tensor<2,dim,NumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
2458 * const NumberType det_F = determinant(F);
2459 * const Tensor<2,dim,NumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
2460 * const SymmetricTensor<2,dim,NumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
2461 * const Tensor<2,dim,NumberType> F_inv = invert(F);
2462 * Assert(det_F > NumberType(0.0), ExcInternalError());
2464 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2466 * const unsigned int k_group = fe.system_to_base_index(k).first.first;
2468 * if (k_group == u_dof)
2470 * scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
2471 * scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
2474 * Assert(k_group <= u_dof, ExcInternalError());
2477 * const SymmetricTensor<2,dim,NumberType> tau = lqph[q_point]->get_tau(det_F,b_bar);
2478 * const SymmetricTensor<4,dim,NumberType> Jc = lqph[q_point]->get_Jc(det_F,b_bar);
2479 * const Tensor<2,dim,NumberType> tau_ns (tau);
2483 * Next we define some aliases to make the assembly process easier to
2487 * const std::vector<SymmetricTensor<2, dim> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
2488 * const std::vector<Tensor<2, dim> > &grad_Nx = scratch.grad_Nx[q_point];
2489 * const double JxW = scratch.fe_values_ref.JxW(q_point);
2491 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2493 * const unsigned int component_i = fe.system_to_component_index(i).first;
2494 * const unsigned int i_group = fe.system_to_base_index(i).first.first;
2496 * if (i_group == u_dof)
2497 * data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
2499 * Assert(i_group <= u_dof, ExcInternalError());
2501 * for (unsigned int j = 0; j <= i; ++j)
2503 * const unsigned int component_j = fe.system_to_component_index(j).first;
2504 * const unsigned int j_group = fe.system_to_base_index(j).first.first;
2508 * This is the @f$\mathsf{\mathbf{k}}_{\mathbf{u} \mathbf{u}}@f$
2509 * contribution. It comprises a material contribution, and a
2510 * geometrical stress contribution which is only added along
2511 * the local matrix diagonals:
2514 * if ((i_group == j_group) && (i_group == u_dof))
2516 * data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
2517 * * symm_grad_Nx[j] * JxW;
2518 * if (component_i == component_j) // geometrical stress contribution
2519 * data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau_ns
2520 * * grad_Nx[j][component_j] * JxW;
2523 * Assert((i_group <= u_dof) && (j_group <= u_dof),
2524 * ExcInternalError());
2532 * Finally, we need to copy the lower half of the local matrix into the
2536 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2537 * for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
2538 * data.cell_matrix(i, j) = data.cell_matrix(j, i);
2543 * #ifdef ENABLE_SACADO_FORMULATION
2546 * template <int dim>
2547 * struct Assembler<dim,Sacado::Fad::DFad<double> > : Assembler_Base<dim,Sacado::Fad::DFad<double> >
2549 * typedef Sacado::Fad::DFad<double> ADNumberType;
2550 * using typename Assembler_Base<dim,ADNumberType>::ScratchData_ASM;
2551 * using typename Assembler_Base<dim,ADNumberType>::PerTaskData_ASM;
2553 * virtual ~Assembler() {}
2556 * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2557 * ScratchData_ASM &scratch,
2558 * PerTaskData_ASM &data)
2562 * Aliases for data referenced from the Solid class
2565 * const unsigned int &n_q_points = data.solid->n_q_points;
2566 * const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2567 * const FESystem<dim> &fe = data.solid->fe;
2568 * const unsigned int &u_dof = data.solid->u_dof;
2569 * const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
2573 * scratch.fe_values_ref.reinit(cell);
2574 * cell->get_dof_indices(data.local_dof_indices);
2576 * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
2577 * const_cast<const Solid<dim,ADNumberType> *>(data.solid)->quadrature_point_history.get_data(cell);
2578 * Assert(lqph.size() == n_q_points, ExcInternalError());
2580 * const unsigned int n_independent_variables = data.local_dof_indices.size();
2581 * std::vector<double> local_dof_values(n_independent_variables);
2582 * cell->get_dof_values(scratch.solution_total,
2583 * local_dof_values.begin(),
2584 * local_dof_values.end());
2588 * We now retrieve a set of degree-of-freedom values that
2589 * have the operations that are performed with them tracked.
2592 * std::vector<ADNumberType> local_dof_values_ad (n_independent_variables);
2593 * for (unsigned int i=0; i<n_independent_variables; ++i)
2594 * local_dof_values_ad[i] = ADNumberType(n_independent_variables, i, local_dof_values[i]);
2598 * Compute all values, gradients etc. based on sensitive
2599 * AD degree-of-freedom values.
2602 * scratch.fe_values_ref[u_fe].get_function_gradients_from_local_dof_values(
2603 * local_dof_values_ad,
2604 * scratch.solution_grads_u_total);
2608 * Accumulate the residual value for each degree of freedom.
2609 * Note: Its important that the vectors is initialised (zero'd) correctly.
2612 * std::vector<ADNumberType> residual_ad (dofs_per_cell, ADNumberType(0.0));
2613 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2623 *
for (
unsigned int k = 0; k < dofs_per_cell; ++k)
2625 *
const unsigned int k_group = fe.system_to_base_index(k).first.first;
2627 *
if (k_group == u_dof)
2629 * scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
2630 * scratch.symm_grad_Nx[q_point][k] =
symmetrize(scratch.grad_Nx[q_point][k]);
2640 * Next we define some position-dependent aliases, again to
2641 * make the assembly process easier to follow.
2644 *
const std::vector<SymmetricTensor<2, dim,ADNumberType> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
2645 *
const double JxW = scratch.fe_values_ref.JxW(q_point);
2647 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
2649 *
const unsigned int i_group = fe.system_to_base_index(i).first.first;
2651 *
if (i_group == u_dof)
2652 * residual_ad[i] += (symm_grad_Nx[i] * tau) * JxW;
2658 *
for (
unsigned int I=0; I<n_independent_variables; ++I)
2660 *
const ADNumberType &residual_I = residual_ad[I];
2661 * data.cell_rhs(I) = -residual_I.val();
2662 *
for (
unsigned int J=0; J<n_independent_variables; ++J)
2666 * Compute the gradients of the residual entry [forward-mode]
2669 * data.cell_matrix(I,J) = residual_I.dx(J);
2677 *
template <
int dim>
2678 *
struct Assembler<dim,Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > > : Assembler_Base<dim,Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > >
2680 *
typedef Sacado::Fad::DFad<double> ADDerivType;
2681 *
typedef Sacado::Rad::ADvar<ADDerivType> ADNumberType;
2682 *
using typename Assembler_Base<dim,ADNumberType>::ScratchData_ASM;
2683 *
using typename Assembler_Base<dim,ADNumberType>::PerTaskData_ASM;
2685 *
virtual ~Assembler() {}
2689 * ScratchData_ASM &scratch,
2690 * PerTaskData_ASM &data)
2694 * Aliases
for data referenced from the Solid
class
2697 *
const unsigned int &n_q_points = data.solid->n_q_points;
2702 * scratch.fe_values_ref.reinit(cell);
2703 * cell->get_dof_indices(data.local_dof_indices);
2705 *
const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
2706 * data.solid->quadrature_point_history.get_data(cell);
2709 *
const unsigned int n_independent_variables = data.local_dof_indices.size();
2710 * std::vector<double> local_dof_values(n_independent_variables);
2711 * cell->get_dof_values(scratch.solution_total,
2712 * local_dof_values.begin(),
2713 * local_dof_values.end());
2717 * We now retrieve a
set of degree-of-freedom values that
2718 * have the operations that are performed with them tracked.
2721 * std::vector<ADNumberType> local_dof_values_ad (n_independent_variables);
2722 *
for (
unsigned int i=0; i<n_independent_variables; ++i)
2723 * local_dof_values_ad[i] = ADDerivType(n_independent_variables, i, local_dof_values[i]);
2727 * Compute all values, gradients etc. based on sensitive
2728 * AD degree-of-freedom values.
2731 * scratch.fe_values_ref[u_fe].get_function_gradients_from_local_dof_values(
2732 * local_dof_values_ad,
2733 * scratch.solution_grads_u_total);
2737 * Next we compute the total potential energy of the element.
2738 * This is defined as follows:
2739 * Total energy = (
internal - external) energies
2740 * Note: Its important that
this value is initialised (
zero'd) correctly.
2743 * ADNumberType cell_energy_ad = ADNumberType(0.0);
2744 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2746 * const Tensor<2,dim,ADNumberType> &grad_u = scratch.solution_grads_u_total[q_point];
2747 * const Tensor<2,dim,ADNumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
2748 * const ADNumberType det_F = determinant(F);
2749 * const Tensor<2,dim,ADNumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
2750 * const SymmetricTensor<2,dim,ADNumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
2751 * Assert(det_F > ADNumberType(0.0), ExcInternalError());
2755 * Next we define some position-dependent aliases, again to
2756 * make the assembly process easier to follow.
2759 * const double JxW = scratch.fe_values_ref.JxW(q_point);
2761 * const ADNumberType Psi = lqph[q_point]->get_Psi(det_F,b_bar);
2765 * We extract the configuration-dependent material energy
2766 * from our QPH history objects for the current quadrature point
2767 * and integrate its contribution to increment the total
2771 * cell_energy_ad += Psi * JxW;
2776 * Compute derivatives of reverse-mode AD variables
2779 * ADNumberType::Gradcomp();
2781 * for (unsigned int I=0; I<n_independent_variables; ++I)
2785 * This computes the adjoint df/dX_{i} [reverse-mode]
2788 * const ADDerivType residual_I = local_dof_values_ad[I].adj();
2789 * data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
2790 * for (unsigned int J=0; J<n_independent_variables; ++J)
2794 * Compute the gradients of the residual entry [forward-mode]
2797 * data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
2810 * Since we use TBB for assembly, we simply setup a copy of the
2811 * data structures required for the process and pass them, along
2812 * with the memory addresses of the assembly functions to the
2813 * WorkStream object for processing. Note that we must ensure that
2814 * the matrix is reset before any assembly operations can occur.
2817 * template <int dim,typename NumberType>
2818 * void Solid<dim,NumberType>::assemble_system(const BlockVector<double> &solution_delta)
2820 * timer.enter_subsection("Assemble linear system");
2821 * std::cout << " ASM " << std::flush;
2823 * tangent_matrix = 0.0;
2826 * const UpdateFlags uf_cell(update_gradients |
2827 * update_JxW_values);
2828 * const UpdateFlags uf_face(update_values |
2829 * update_JxW_values);
2831 * const BlockVector<double> solution_total(get_total_solution(solution_delta));
2832 * typename Assembler_Base<dim,NumberType>::PerTaskData_ASM per_task_data(this);
2833 * typename Assembler_Base<dim,NumberType>::ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total);
2834 * Assembler<dim,NumberType> assembler;
2836 * WorkStream::run(dof_handler_ref.begin_active(),
2837 * dof_handler_ref.end(),
2838 * static_cast<Assembler_Base<dim,NumberType>&>(assembler),
2839 * &Assembler_Base<dim,NumberType>::assemble_system_one_cell,
2840 * &Assembler_Base<dim,NumberType>::copy_local_to_global_ASM,
2844 * timer.leave_subsection();
2851 * <a name="Solidmake_constraints"></a>
2852 * <h4>Solid::make_constraints</h4>
2853 * The constraints for this problem are simple to describe.
2854 * However, since we are dealing with an iterative Newton method,
2855 * it should be noted that any displacement constraints should only
2856 * be specified at the zeroth iteration and subsequently no
2857 * additional contributions are to be made since the constraints
2858 * are already exactly satisfied.
2861 * template <int dim,typename NumberType>
2862 * void Solid<dim,NumberType>::make_constraints(const int &it_nr)
2864 * std::cout << " CST " << std::flush;
2868 * Since the constraints are different at different Newton iterations, we
2869 * need to clear the constraints matrix and completely rebuild
2870 * it. However, after the first iteration, the constraints remain the same
2871 * and we can simply skip the rebuilding step if we do not clear it.
2876 * constraints.clear();
2877 * const bool apply_dirichlet_bc = (it_nr == 0);
2881 * The boundary conditions for the indentation problem are as follows: On
2882 * the -x, -y and -z faces (ID's 0,2,4) we
set up a symmetry condition to
2883 * allow only planar movement
while the +x and +y faces (ID
's 1,3) are
2884 * traction free. In this contrived problem, part of the +z face (ID 5) is
2885 * set to have no motion in the x- and y-component. Finally, as described
2886 * earlier, the other part of the +z face has an the applied pressure but
2887 * is also constrained in the x- and y-directions.
2891 * In the following, we will have to tell the function interpolation
2892 * boundary values which components of the solution vector should be
2893 * constrained (i.e., whether it's the x-, y-, z-displacements or
2894 * combinations thereof). This is done
using ComponentMask objects (see
2895 * @ref GlossComponentMask) which we can get from the finite element
if we
2896 * provide it with an extractor
object for the component we wish to
2897 * select. To
this end we
first set up such extractor objects and later
2898 * use it when generating the relevant component masks:
2902 * Fixed left hand side of the beam
2908 *
if (apply_dirichlet_bc ==
true)
2913 * fe.component_mask(u_fe));
2919 * fe.component_mask(u_fe));
2924 * Zero Z-displacement through thickness direction
2925 * This corresponds to a plane strain condition being imposed on the beam
2933 *
if (apply_dirichlet_bc ==
true)
2938 * fe.component_mask(z_displacement));
2944 * fe.component_mask(z_displacement));
2947 * constraints.close();
2953 * <a name=
"Solidsolve_linear_system"></a>
2954 * <h4>Solid::solve_linear_system</h4>
2955 * As the system is composed of a single block, defining a solution scheme
2956 *
for the linear problem is straight-forward.
2959 *
template <
int dim,
typename NumberType>
2960 * std::pair<unsigned int, double>
2966 *
unsigned int lin_it = 0;
2967 *
double lin_res = 0.0;
2971 * We solve
for the incremental displacement @f$d\mathbf{u}@f$.
2975 * timer.enter_subsection(
"Linear solver");
2976 * std::cout <<
" SLV " << std::flush;
2977 *
if (parameters.type_lin ==
"CG")
2979 *
const int solver_its = tangent_matrix.
block(u_dof, u_dof).
m()
2980 * * parameters.max_iterations_lin;
2981 *
const double tol_sol = parameters.tol_lin
2991 * We
've chosen by default a SSOR preconditioner as it appears to
2992 * provide the fastest solver convergence characteristics for this
2993 * problem on a single-thread machine. However, for multicore
2994 * computing, the Jacobi preconditioner which is multithreaded may
2995 * converge quicker for larger linear systems.
2998 * PreconditionSelector<SparseMatrix<double>, Vector<double> >
2999 * preconditioner (parameters.preconditioner_type,
3000 * parameters.preconditioner_relaxation);
3001 * preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
3003 * solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
3004 * newton_update.block(u_dof),
3005 * system_rhs.block(u_dof),
3008 * lin_it = solver_control.last_step();
3009 * lin_res = solver_control.last_value();
3011 * else if (parameters.type_lin == "Direct")
3015 * Otherwise if the problem is small
3016 * enough, a direct solver can be
3020 * SparseDirectUMFPACK A_direct;
3021 * A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
3022 * A_direct.vmult(newton_update.block(u_dof), system_rhs.block(u_dof));
3028 * Assert (false, ExcMessage("Linear solver type not implemented"));
3030 * timer.leave_subsection();
3035 * Now that we have the displacement update, distribute the constraints
3036 * back to the Newton update:
3039 * constraints.distribute(newton_update);
3041 * return std::make_pair(lin_it, lin_res);
3047 * <a name="Solidoutput_results"></a>
3048 * <h4>Solid::output_results</h4>
3049 * Here we present how the results are written to file to be viewed
3050 * using ParaView or Visit. The method is similar to that shown in the
3051 * tutorials so will not be discussed in detail.
3054 * template <int dim,typename NumberType>
3055 * void Solid<dim,NumberType>::output_results() const
3057 * DataOut<dim> data_out;
3058 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3059 * data_component_interpretation(dim,
3060 * DataComponentInterpretation::component_is_part_of_vector);
3062 * std::vector<std::string> solution_name(dim, "displacement");
3064 * data_out.attach_dof_handler(dof_handler_ref);
3065 * data_out.add_data_vector(solution_n,
3067 * DataOut<dim>::type_dof_data,
3068 * data_component_interpretation);
3072 * Since we are dealing with a large deformation problem, it would be nice
3073 * to display the result on a displaced grid! The MappingQEulerian class
3074 * linked with the DataOut class provides an interface through which this
3075 * can be achieved without physically moving the grid points in the
3076 * Triangulation object ourselves. We first need to copy the solution to
3077 * a temporary vector and then create the Eulerian mapping. We also
3078 * specify the polynomial degree to the DataOut object in order to produce
3079 * a more refined output data set when higher order polynomials are used.
3082 * Vector<double> soln(solution_n.size());
3083 * for (unsigned int i = 0; i < soln.size(); ++i)
3084 * soln(i) = solution_n(i);
3085 * MappingQEulerian<dim> q_mapping(degree, dof_handler_ref, soln);
3086 * data_out.build_patches(q_mapping, degree);
3088 * std::ostringstream filename;
3089 * filename << "solution-" << time.get_timestep() << ".vtk";
3091 * std::ofstream output(filename.str().c_str());
3092 * data_out.write_vtk(output);
3101 * <a name="Mainfunction"></a>
3102 * <h3>Main function</h3>
3103 * Lastly we provide the main driver function which appears
3104 * no different to the other tutorials.
3107 * int main (int argc, char *argv[])
3109 * using namespace dealii;
3110 * using namespace Cook_Membrane;
3112 * const unsigned int dim = 2;
3115 * deallog.depth_console(0);
3116 * Parameters::AllParameters parameters("parameters.prm");
3117 * if (parameters.automatic_differentiation_order == 0)
3119 * std::cout << "Assembly method: Residual and linearisation are computed manually." << std::endl;
3123 * Allow multi-threading
3126 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3127 * ::numbers::invalid_unsigned_int);
3129 * typedef double NumberType;
3130 * Solid<dim,NumberType> solid_3d(parameters);
3133 * #ifdef ENABLE_SACADO_FORMULATION
3134 * else if (parameters.automatic_differentiation_order == 1)
3136 * std::cout << "Assembly method: Residual computed manually; linearisation performed using AD." << std::endl;
3140 * Allow multi-threading
3143 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3144 * ::numbers::invalid_unsigned_int);
3146 * typedef Sacado::Fad::DFad<double> NumberType;
3147 * Solid<dim,NumberType> solid_3d(parameters);
3150 * else if (parameters.automatic_differentiation_order == 2)
3152 * std::cout << "Assembly method: Residual and linearisation computed using AD." << std::endl;
3156 * Sacado Rad-Fad is not thread-safe, so disable threading.
3157 * Parallisation using MPI would be possible though.
3160 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3163 * typedef Sacado::Rad::ADvar< Sacado::Fad::DFad<double> > NumberType;
3164 * Solid<dim,NumberType> solid_3d(parameters);
3170 * AssertThrow(false,
3171 * ExcMessage("The selected assembly method is not supported. "
3172 * "You need deal.II 9.0 and Trilinos with the Sacado package "
3173 * "to enable assembly using automatic differentiation."));
3176 * catch (std::exception &exc)
3178 * std::cerr << std::endl << std::endl
3179 * << "----------------------------------------------------"
3181 * std::cerr << "Exception on processing: " << std::endl << exc.what()
3182 * << std::endl << "Aborting!" << std::endl
3183 * << "----------------------------------------------------"
3190 * std::cerr << std::endl << std::endl
3191 * << "----------------------------------------------------"
3193 * std::cerr << "Unknown exception!" << std::endl << "Aborting!"
3195 * << "----------------------------------------------------"