Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Quasi_static_Finite_strain_Compressible_Elasticity.h
Go to the documentation of this file.
1 
198  *
199  * /*
200  * * Authors: Jean-Paul Pelteret, University of Erlangen-Nuremberg,
201  * * Andrew McBride, University of Cape Town, 2015, 2017
202  * */
203  *
204  *
205  * @endcode
206  *
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.
210  *
211  * @code
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>
222  *
223  * #include <deal.II/dofs/dof_renumbering.h>
224  * #include <deal.II/dofs/dof_tools.h>
225  *
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>
231  *
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>
239  *
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>
248  *
249  * #include <deal.II/numerics/data_out.h>
250  * #include <deal.II/numerics/vector_tools.h>
251  *
252  * #include <deal.II/base/config.h>
253  * #if DEAL_II_VERSION_MAJOR >= 9 && defined(DEAL_II_WITH_TRILINOS)
254  * #include <deal.II/differentiation/ad.h>
255  * #define ENABLE_SACADO_FORMULATION
256  * #endif
257  *
258  * @endcode
259  *
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
263  *
264  * @code
265  * #include <deal.II/physics/elasticity/kinematics.h>
266  * #include <deal.II/physics/elasticity/standard_tensors.h>
267  *
268  * #include <iostream>
269  * #include <fstream>
270  *
271  *
272  * @endcode
273  *
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
276  * into it:
277  *
278  * @code
279  * namespace Cook_Membrane
280  * {
281  * using namespace dealii;
282  *
283  * @endcode
284  *
285  *
286  * <a name="Runtimeparameters"></a>
287  * <h3>Run-time parameters</h3>
288  *
289 
290  *
291  * There are several parameters that can be set in the code so we set up a
292  * ParameterHandler object to read in the choices at run-time.
293  *
294  * @code
295  * namespace Parameters
296  * {
297  * @endcode
298  *
299  *
300  * <a name="Assemblymethod"></a>
301  * <h4>Assembly method</h4>
302  *
303 
304  *
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
307  * employed.
308  *
309  * @code
310  * struct AssemblyMethod
311  * {
312  * unsigned int automatic_differentiation_order;
313  *
314  * static void
315  * declare_parameters(ParameterHandler &prm);
316  *
317  * void
318  * parse_parameters(ParameterHandler &prm);
319  * };
320  *
321  *
322  * void AssemblyMethod::declare_parameters(ParameterHandler &prm)
323  * {
324  * prm.enter_subsection("Assembly method");
325  * {
326  * prm.declare_entry("Automatic differentiation order", "0",
327  * Patterns::Integer(0,2),
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.");
332  * }
333  * prm.leave_subsection();
334  * }
335  *
336  * void AssemblyMethod::parse_parameters(ParameterHandler &prm)
337  * {
338  * prm.enter_subsection("Assembly method");
339  * {
340  * automatic_differentiation_order = prm.get_integer("Automatic differentiation order");
341  * }
342  * prm.leave_subsection();
343  * }
344  *
345  * @endcode
346  *
347  *
348  * <a name="FiniteElementsystem"></a>
349  * <h4>Finite Element system</h4>
350  *
351 
352  *
353  * Here we specify the polynomial order used to approximate the solution.
354  * The quadrature order should be adjusted accordingly.
355  *
356  * @code
357  * struct FESystem
358  * {
359  * unsigned int poly_degree;
360  * unsigned int quad_order;
361  *
362  * static void
363  * declare_parameters(ParameterHandler &prm);
364  *
365  * void
366  * parse_parameters(ParameterHandler &prm);
367  * };
368  *
369  *
370  * void FESystem::declare_parameters(ParameterHandler &prm)
371  * {
372  * prm.enter_subsection("Finite element system");
373  * {
374  * prm.declare_entry("Polynomial degree", "2",
375  * Patterns::Integer(0),
376  * "Displacement system polynomial order");
377  *
378  * prm.declare_entry("Quadrature order", "3",
379  * Patterns::Integer(0),
380  * "Gauss quadrature order");
381  * }
382  * prm.leave_subsection();
383  * }
384  *
385  * void FESystem::parse_parameters(ParameterHandler &prm)
386  * {
387  * prm.enter_subsection("Finite element system");
388  * {
389  * poly_degree = prm.get_integer("Polynomial degree");
390  * quad_order = prm.get_integer("Quadrature order");
391  * }
392  * prm.leave_subsection();
393  * }
394  *
395  * @endcode
396  *
397  *
398  * <a name="Geometry"></a>
399  * <h4>Geometry</h4>
400  *
401 
402  *
403  * Make adjustments to the problem geometry and its discretisation.
404  *
405  * @code
406  * struct Geometry
407  * {
408  * unsigned int elements_per_edge;
409  * double scale;
410  *
411  * static void
412  * declare_parameters(ParameterHandler &prm);
413  *
414  * void
415  * parse_parameters(ParameterHandler &prm);
416  * };
417  *
418  * void Geometry::declare_parameters(ParameterHandler &prm)
419  * {
420  * prm.enter_subsection("Geometry");
421  * {
422  * prm.declare_entry("Elements per edge", "32",
423  * Patterns::Integer(0),
424  * "Number of elements per long edge of the beam");
425  *
426  * prm.declare_entry("Grid scale", "1e-3",
427  * Patterns::Double(0.0),
428  * "Global grid scaling factor");
429  * }
430  * prm.leave_subsection();
431  * }
432  *
433  * void Geometry::parse_parameters(ParameterHandler &prm)
434  * {
435  * prm.enter_subsection("Geometry");
436  * {
437  * elements_per_edge = prm.get_integer("Elements per edge");
438  * scale = prm.get_double("Grid scale");
439  * }
440  * prm.leave_subsection();
441  * }
442  *
443  * @endcode
444  *
445  *
446  * <a name="Materials"></a>
447  * <h4>Materials</h4>
448  *
449 
450  *
451  * We also need the shear modulus @f$ \mu @f$ and Poisson ration @f$ \nu @f$ for the
452  * neo-Hookean material.
453  *
454  * @code
455  * struct Materials
456  * {
457  * double nu;
458  * double mu;
459  *
460  * static void
461  * declare_parameters(ParameterHandler &prm);
462  *
463  * void
464  * parse_parameters(ParameterHandler &prm);
465  * };
466  *
467  * void Materials::declare_parameters(ParameterHandler &prm)
468  * {
469  * prm.enter_subsection("Material properties");
470  * {
471  * prm.declare_entry("Poisson's ratio", "0.3",
472  * Patterns::Double(-1.0,0.5),
473  * "Poisson's ratio");
474  *
475  * prm.declare_entry("Shear modulus", "0.4225e6",
476  * Patterns::Double(),
477  * "Shear modulus");
478  * }
479  * prm.leave_subsection();
480  * }
481  *
482  * void Materials::parse_parameters(ParameterHandler &prm)
483  * {
484  * prm.enter_subsection("Material properties");
485  * {
486  * nu = prm.get_double("Poisson's ratio");
487  * mu = prm.get_double("Shear modulus");
488  * }
489  * prm.leave_subsection();
490  * }
491  *
492  * @endcode
493  *
494  *
495  * <a name="Linearsolver"></a>
496  * <h4>Linear solver</h4>
497  *
498 
499  *
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.
503  *
504  * @code
505  * struct LinearSolver
506  * {
507  * std::string type_lin;
508  * double tol_lin;
509  * double max_iterations_lin;
510  * std::string preconditioner_type;
511  * double preconditioner_relaxation;
512  *
513  * static void
514  * declare_parameters(ParameterHandler &prm);
515  *
516  * void
517  * parse_parameters(ParameterHandler &prm);
518  * };
519  *
520  * void LinearSolver::declare_parameters(ParameterHandler &prm)
521  * {
522  * prm.enter_subsection("Linear solver");
523  * {
524  * prm.declare_entry("Solver type", "CG",
525  * Patterns::Selection("CG|Direct"),
526  * "Type of solver used to solve the linear system");
527  *
528  * prm.declare_entry("Residual", "1e-6",
529  * Patterns::Double(0.0),
530  * "Linear solver residual (scaled by residual norm)");
531  *
532  * prm.declare_entry("Max iteration multiplier", "1",
533  * Patterns::Double(0.0),
534  * "Linear solver iterations (multiples of the system matrix size)");
535  *
536  * prm.declare_entry("Preconditioner type", "ssor",
537  * Patterns::Selection("jacobi|ssor"),
538  * "Type of preconditioner");
539  *
540  * prm.declare_entry("Preconditioner relaxation", "0.65",
541  * Patterns::Double(0.0),
542  * "Preconditioner relaxation value");
543  * }
544  * prm.leave_subsection();
545  * }
546  *
547  * void LinearSolver::parse_parameters(ParameterHandler &prm)
548  * {
549  * prm.enter_subsection("Linear solver");
550  * {
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");
556  * }
557  * prm.leave_subsection();
558  * }
559  *
560  * @endcode
561  *
562  *
563  * <a name="Nonlinearsolver"></a>
564  * <h4>Nonlinear solver</h4>
565  *
566 
567  *
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.
571  *
572  * @code
573  * struct NonlinearSolver
574  * {
575  * unsigned int max_iterations_NR;
576  * double tol_f;
577  * double tol_u;
578  *
579  * static void
580  * declare_parameters(ParameterHandler &prm);
581  *
582  * void
583  * parse_parameters(ParameterHandler &prm);
584  * };
585  *
586  * void NonlinearSolver::declare_parameters(ParameterHandler &prm)
587  * {
588  * prm.enter_subsection("Nonlinear solver");
589  * {
590  * prm.declare_entry("Max iterations Newton-Raphson", "10",
591  * Patterns::Integer(0),
592  * "Number of Newton-Raphson iterations allowed");
593  *
594  * prm.declare_entry("Tolerance force", "1.0e-9",
595  * Patterns::Double(0.0),
596  * "Force residual tolerance");
597  *
598  * prm.declare_entry("Tolerance displacement", "1.0e-6",
599  * Patterns::Double(0.0),
600  * "Displacement error tolerance");
601  * }
602  * prm.leave_subsection();
603  * }
604  *
605  * void NonlinearSolver::parse_parameters(ParameterHandler &prm)
606  * {
607  * prm.enter_subsection("Nonlinear solver");
608  * {
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");
612  * }
613  * prm.leave_subsection();
614  * }
615  *
616  * @endcode
617  *
618  *
619  * <a name="Time"></a>
620  * <h4>Time</h4>
621  *
622 
623  *
624  * Set the timestep size @f$ \varDelta t @f$ and the simulation end-time.
625  *
626  * @code
627  * struct Time
628  * {
629  * double delta_t;
630  * double end_time;
631  *
632  * static void
633  * declare_parameters(ParameterHandler &prm);
634  *
635  * void
636  * parse_parameters(ParameterHandler &prm);
637  * };
638  *
639  * void Time::declare_parameters(ParameterHandler &prm)
640  * {
641  * prm.enter_subsection("Time");
642  * {
643  * prm.declare_entry("End time", "1",
644  * Patterns::Double(),
645  * "End time");
646  *
647  * prm.declare_entry("Time step size", "0.1",
648  * Patterns::Double(),
649  * "Time step size");
650  * }
651  * prm.leave_subsection();
652  * }
653  *
654  * void Time::parse_parameters(ParameterHandler &prm)
655  * {
656  * prm.enter_subsection("Time");
657  * {
658  * end_time = prm.get_double("End time");
659  * delta_t = prm.get_double("Time step size");
660  * }
661  * prm.leave_subsection();
662  * }
663  *
664  * @endcode
665  *
666  *
667  * <a name="Allparameters"></a>
668  * <h4>All parameters</h4>
669  *
670 
671  *
672  * Finally we consolidate all of the above structures into a single container
673  * that holds all of our run-time selections.
674  *
675  * @code
676  * struct AllParameters :
677  * public AssemblyMethod,
678  * public FESystem,
679  * public Geometry,
680  * public Materials,
681  * public LinearSolver,
682  * public NonlinearSolver,
683  * public Time
684  *
685  * {
686  * AllParameters(const std::string &input_file);
687  *
688  * static void
689  * declare_parameters(ParameterHandler &prm);
690  *
691  * void
692  * parse_parameters(ParameterHandler &prm);
693  * };
694  *
695  * AllParameters::AllParameters(const std::string &input_file)
696  * {
697  * ParameterHandler prm;
698  * declare_parameters(prm);
699  * prm.parse_input(input_file);
700  * parse_parameters(prm);
701  * }
702  *
703  * void AllParameters::declare_parameters(ParameterHandler &prm)
704  * {
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);
712  * }
713  *
714  * void AllParameters::parse_parameters(ParameterHandler &prm)
715  * {
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);
723  * }
724  * }
725  *
726  *
727  * @endcode
728  *
729  *
730  * <a name="Timeclass"></a>
731  * <h3>Time class</h3>
732  *
733 
734  *
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
737  * size.
738  *
739  * @code
740  * class Time
741  * {
742  * public:
743  * Time (const double time_end,
744  * const double delta_t)
745  * :
746  * timestep(0),
747  * time_current(0.0),
748  * time_end(time_end),
749  * delta_t(delta_t)
750  * {}
751  *
752  * virtual ~Time()
753  * {}
754  *
755  * double current() const
756  * {
757  * return time_current;
758  * }
759  * double end() const
760  * {
761  * return time_end;
762  * }
763  * double get_delta_t() const
764  * {
765  * return delta_t;
766  * }
767  * unsigned int get_timestep() const
768  * {
769  * return timestep;
770  * }
771  * void increment()
772  * {
773  * time_current += delta_t;
774  * ++timestep;
775  * }
776  *
777  * private:
778  * unsigned int timestep;
779  * double time_current;
780  * const double time_end;
781  * const double delta_t;
782  * };
783  *
784  * @endcode
785  *
786  *
787  * <a name="CompressibleneoHookeanmaterialwithinaonefieldformulation"></a>
788  * <h3>Compressible neo-Hookean material within a one-field formulation</h3>
789  *
790 
791  *
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)
798  * @f$.
799  *
800 
801  *
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
813  * parameter</a>.
814  *
815 
816  *
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.
824  *
825  * @code
826  * template <int dim,typename NumberType>
827  * class Material_Compressible_Neo_Hook_One_Field
828  * {
829  * public:
830  * Material_Compressible_Neo_Hook_One_Field(const double mu,
831  * const double nu)
832  * :
833  * kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
834  * c_1(mu / 2.0)
835  * {
836  * Assert(kappa > 0, ExcInternalError());
837  * }
838  *
839  * ~Material_Compressible_Neo_Hook_One_Field()
840  * {}
841  *
842  * @endcode
843  *
844  * The first function is the total energy
845  * @f$\Psi = \Psi_{\textrm{iso}} + \Psi_{\textrm{vol}}@f$.
846  *
847  * @code
848  * NumberType
849  * get_Psi(const NumberType &det_F,
850  * const SymmetricTensor<2,dim,NumberType> &b_bar) const
851  * {
852  * return get_Psi_vol(det_F) + get_Psi_iso(b_bar);
853  * }
854  *
855  * @endcode
856  *
857  * The second function determines the Kirchhoff stress @f$\boldsymbol{\tau}
858  * = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}@f$
859  *
860  * @code
861  * SymmetricTensor<2,dim,NumberType>
862  * get_tau(const NumberType &det_F,
863  * const SymmetricTensor<2,dim,NumberType> &b_bar)
864  * {
865  * @endcode
866  *
867  * See Holzapfel p231 eq6.98 onwards
868  *
869  * @code
870  * return get_tau_vol(det_F) + get_tau_iso(b_bar);
871  * }
872  *
873  * @endcode
874  *
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$
880  *
881  * @code
882  * SymmetricTensor<4,dim,NumberType>
883  * get_Jc(const NumberType &det_F,
884  * const SymmetricTensor<2,dim,NumberType> &b_bar) const
885  * {
886  * return get_Jc_vol(det_F) + get_Jc_iso(b_bar);
887  * }
888  *
889  * private:
890  * @endcode
891  *
892  * Define constitutive model parameters @f$\kappa@f$ (bulk modulus) and the
893  * neo-Hookean model parameter @f$c_1@f$:
894  *
895  * @code
896  * const double kappa;
897  * const double c_1;
898  *
899  * @endcode
900  *
901  * Value of the volumetric free energy
902  *
903  * @code
904  * NumberType
905  * get_Psi_vol(const NumberType &det_F) const
906  * {
907  * return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
908  * }
909  *
910  * @endcode
911  *
912  * Value of the isochoric free energy
913  *
914  * @code
915  * NumberType
916  * get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
917  * {
918  * return c_1 * (trace(b_bar) - dim);
919  * }
920  *
921  * @endcode
922  *
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$
926  *
927  * @code
928  * NumberType
929  * get_dPsi_vol_dJ(const NumberType &det_F) const
930  * {
931  * return (kappa / 2.0) * (det_F - 1.0 / det_F);
932  * }
933  *
934  * @endcode
935  *
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".
940  *
941  * @code
942  * SymmetricTensor<2,dim,NumberType>
943  * get_tau_vol(const NumberType &det_F) const
944  * {
945  * return NumberType(get_dPsi_vol_dJ(det_F) * det_F) * Physics::Elasticity::StandardTensors<dim>::I;
946  * }
947  *
948  * @endcode
949  *
950  * Next, determine the isochoric Kirchhoff stress
951  * @f$\boldsymbol{\tau}_{\textrm{iso}} =
952  * \mathcal{P}:\overline{\boldsymbol{\tau}}@f$:
953  *
954  * @code
955  * SymmetricTensor<2,dim,NumberType>
956  * get_tau_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
957  * {
958  * return Physics::Elasticity::StandardTensors<dim>::dev_P * get_tau_bar(b_bar);
959  * }
960  *
961  * @endcode
962  *
963  * Then, determine the fictitious Kirchhoff stress
964  * @f$\overline{\boldsymbol{\tau}}@f$:
965  *
966  * @code
967  * SymmetricTensor<2,dim,NumberType>
968  * get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const
969  * {
970  * return 2.0 * c_1 * b_bar;
971  * }
972  *
973  * @endcode
974  *
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
979  * J}@f$
980  *
981  * @code
982  * NumberType
983  * get_d2Psi_vol_dJ2(const NumberType &det_F) const
984  * {
985  * return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
986  * }
987  *
988  * @endcode
989  *
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$.
995  *
996  * @code
997  * SymmetricTensor<4,dim,NumberType>
998  * get_Jc_vol(const NumberType &det_F) const
999  * {
1000  * @endcode
1001  *
1002  * See Holzapfel p265
1003  *
1004  * @code
1005  * return det_F
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 );
1008  * }
1009  *
1010  * @endcode
1011  *
1012  * Calculate the isochoric part of the tangent @f$J
1013  * \mathfrak{c}_\textrm{iso}@f$:
1014  *
1015  * @code
1016  * SymmetricTensor<4,dim,NumberType>
1017  * get_Jc_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
1018  * {
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,
1026  * tau_iso);
1027  * const SymmetricTensor<4, dim> c_bar = get_c_bar();
1028  *
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;
1034  * }
1035  *
1036  * @endcode
1037  *
1038  * Calculate the fictitious elasticity tensor @f$\overline{\mathfrak{c}}@f$.
1039  * For the material model chosen this is simply zero:
1040  *
1041  * @code
1042  * SymmetricTensor<4,dim,double>
1043  * get_c_bar() const
1044  * {
1045  * return SymmetricTensor<4, dim>();
1046  * }
1047  * };
1048  *
1049  * @endcode
1050  *
1051  *
1052  * <a name="Quadraturepointhistory"></a>
1053  * <h3>Quadrature point history</h3>
1054  *
1055 
1056  *
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.
1063  *
1064  * @code
1065  * template <int dim,typename NumberType>
1066  * class PointHistory
1067  * {
1068  * public:
1069  * PointHistory()
1070  * {}
1071  *
1072  * virtual ~PointHistory()
1073  * {}
1074  *
1075  * @endcode
1076  *
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$.
1081  *
1082  * @code
1083  * void setup_lqp (const Parameters::AllParameters &parameters)
1084  * {
1085  * material.reset(new Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu,
1086  * parameters.nu));
1087  * }
1088  *
1089  * @endcode
1090  *
1091  * We offer an interface to retrieve certain data.
1092  * This is the strain energy:
1093  *
1094  * @code
1095  * NumberType
1096  * get_Psi(const NumberType &det_F,
1097  * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1098  * {
1099  * return material->get_Psi(det_F,b_bar);
1100  * }
1101  *
1102  * @endcode
1103  *
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:
1107  *
1108  * @code
1109  * SymmetricTensor<2,dim,NumberType>
1110  * get_tau(const NumberType &det_F,
1111  * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1112  * {
1113  * return material->get_tau(det_F,b_bar);
1114  * }
1115  *
1116  * @endcode
1117  *
1118  * And the tangent:
1119  *
1120  * @code
1121  * SymmetricTensor<4,dim,NumberType>
1122  * get_Jc(const NumberType &det_F,
1123  * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1124  * {
1125  * return material->get_Jc(det_F,b_bar);
1126  * }
1127  *
1128  * @endcode
1129  *
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...
1134  *
1135  * @code
1136  * private:
1137  * std::shared_ptr< Material_Compressible_Neo_Hook_One_Field<dim,NumberType> > material;
1138  * };
1139  *
1140  *
1141  * @endcode
1142  *
1143  *
1144  * <a name="Quasistaticcompressiblefinitestrainsolid"></a>
1145  * <h3>Quasi-static compressible finite-strain solid</h3>
1146  *
1147 
1148  *
1149  * Forward declarations for classes that will
1150  * perform assembly of the linear system.
1151  *
1152  * @code
1153  * template <int dim,typename NumberType>
1154  * struct Assembler_Base;
1155  * template <int dim,typename NumberType>
1156  * struct Assembler;
1157  *
1158  * @endcode
1159  *
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:
1164  *
1165  * @code
1166  * template <int dim,typename NumberType>
1167  * class Solid
1168  * {
1169  * public:
1170  * Solid(const Parameters::AllParameters &parameters);
1171  *
1172  * virtual
1173  * ~Solid();
1174  *
1175  * void
1176  * run();
1177  *
1178  * private:
1179  *
1180  * @endcode
1181  *
1182  * We start the collection of member functions with one that builds the
1183  * grid:
1184  *
1185  * @code
1186  * void
1187  * make_grid();
1188  *
1189  * @endcode
1190  *
1191  * Set up the finite element system to be solved:
1192  *
1193  * @code
1194  * void
1195  * system_setup();
1196  *
1197  * @endcode
1198  *
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:
1204  *
1205  * @code
1206  * void
1207  * assemble_system(const BlockVector<double> &solution_delta);
1208  *
1209  * @endcode
1210  *
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.
1214  *
1215  * @code
1216  * friend struct Assembler_Base<dim,NumberType>;
1217  * friend struct Assembler<dim,NumberType>;
1218  *
1219  * @endcode
1220  *
1221  * Apply Dirichlet boundary conditions on the displacement field
1222  *
1223  * @code
1224  * void
1225  * make_constraints(const int &it_nr);
1226  *
1227  * @endcode
1228  *
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
1231  * empty:
1232  *
1233  * @code
1234  * void
1235  * setup_qph();
1236  *
1237  * @endcode
1238  *
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:
1242  *
1243  * @code
1244  * void
1245  * solve_nonlinear_timestep(BlockVector<double> &solution_delta);
1246  *
1247  * std::pair<unsigned int, double>
1248  * solve_linear_system(BlockVector<double> &newton_update);
1249  *
1250  * @endcode
1251  *
1252  * Solution retrieval as well as post-processing and writing data to file :
1253  *
1254  * @code
1255  * BlockVector<double>
1256  * get_total_solution(const BlockVector<double> &solution_delta) const;
1257  *
1258  * void
1259  * output_results() const;
1260  *
1261  * @endcode
1262  *
1263  * Finally, some member variables that describe the current state: A
1264  * collection of the parameters used to describe the problem setup...
1265  *
1266  * @code
1267  * const Parameters::AllParameters &parameters;
1268  *
1269  * @endcode
1270  *
1271  * ...the volume of the reference and current configurations...
1272  *
1273  * @code
1274  * double vol_reference;
1275  * double vol_current;
1276  *
1277  * @endcode
1278  *
1279  * ...and description of the geometry on which the problem is solved:
1280  *
1281  * @code
1282  * Triangulation<dim> triangulation;
1283  *
1284  * @endcode
1285  *
1286  * Also, keep track of the current time and the time spent evaluating
1287  * certain functions
1288  *
1289  * @code
1290  * Time time;
1291  * TimerOutput timer;
1292  *
1293  * @endcode
1294  *
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.
1297  *
1298  * @code
1300  * PointHistory<dim,NumberType> > quadrature_point_history;
1301  *
1302  * @endcode
1303  *
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
1307  * solution vectors:
1308  *
1309  * @code
1310  * const unsigned int degree;
1311  * const FESystem<dim> fe;
1312  * DoFHandler<dim> dof_handler_ref;
1313  * const unsigned int dofs_per_cell;
1314  * const FEValuesExtractors::Vector u_fe;
1315  *
1316  * @endcode
1317  *
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
1324  * problems.
1325  *
1326  * @code
1327  * static const unsigned int n_blocks = 1;
1328  * static const unsigned int n_components = dim;
1329  * static const unsigned int first_u_component = 0;
1330  *
1331  * enum
1332  * {
1333  * u_dof = 0
1334  * };
1335  *
1336  * std::vector<types::global_dof_index> dofs_per_block;
1337  *
1338  * @endcode
1339  *
1340  * Rules for Gauss-quadrature on both the cell and faces. The number of
1341  * quadrature points on both cells and faces is recorded.
1342  *
1343  * @code
1344  * const QGauss<dim> qf_cell;
1345  * const QGauss<dim - 1> qf_face;
1346  * const unsigned int n_q_points;
1347  * const unsigned int n_q_points_f;
1348  *
1349  * @endcode
1350  *
1351  * Objects that store the converged solution and right-hand side vectors,
1352  * as well as the tangent matrix. There is a ConstraintMatrix object used
1353  * to keep track of constraints. We make use of a sparsity pattern
1354  * designed for a block system.
1355  *
1356  * @code
1357  * ConstraintMatrix constraints;
1358  * BlockSparsityPattern sparsity_pattern;
1359  * BlockSparseMatrix<double> tangent_matrix;
1360  * BlockVector<double> system_rhs;
1361  * BlockVector<double> solution_n;
1362  *
1363  * @endcode
1364  *
1365  * Then define a number of variables to store norms and update norms and
1366  * normalisation factors.
1367  *
1368  * @code
1369  * struct Errors
1370  * {
1371  * Errors()
1372  * :
1373  * norm(1.0), u(1.0)
1374  * {}
1375  *
1376  * void reset()
1377  * {
1378  * norm = 1.0;
1379  * u = 1.0;
1380  * }
1381  * void normalise(const Errors &rhs)
1382  * {
1383  * if (rhs.norm != 0.0)
1384  * norm /= rhs.norm;
1385  * if (rhs.u != 0.0)
1386  * u /= rhs.u;
1387  * }
1388  *
1389  * double norm, u;
1390  * };
1391  *
1392  * Errors error_residual, error_residual_0, error_residual_norm, error_update,
1393  * error_update_0, error_update_norm;
1394  *
1395  * @endcode
1396  *
1397  * Methods to calculate error measures
1398  *
1399  * @code
1400  * void
1401  * get_error_residual(Errors &error_residual);
1402  *
1403  * void
1404  * get_error_update(const BlockVector<double> &newton_update,
1405  * Errors &error_update);
1406  *
1407  * @endcode
1408  *
1409  * Print information to screen in a pleasing way...
1410  *
1411  * @code
1412  * static
1413  * void
1414  * print_conv_header();
1415  *
1416  * void
1417  * print_conv_footer();
1418  *
1419  * void
1420  * print_vertical_tip_displacement();
1421  * };
1422  *
1423  * @endcode
1424  *
1425  *
1426  * <a name="ImplementationofthecodeSolidcodeclass"></a>
1427  * <h3>Implementation of the <code>Solid</code> class</h3>
1428  *
1429 
1430  *
1431  *
1432  * <a name="Publicinterface"></a>
1433  * <h4>Public interface</h4>
1434  *
1435 
1436  *
1437  * We initialise the Solid class using data extracted from the parameter file.
1438  *
1439  * @code
1440  * template <int dim,typename NumberType>
1441  * Solid<dim,NumberType>::Solid(const Parameters::AllParameters &parameters)
1442  * :
1443  * parameters(parameters),
1444  * vol_reference (0.0),
1445  * vol_current (0.0),
1447  * time(parameters.end_time, parameters.delta_t),
1448  * timer(std::cout,
1451  * degree(parameters.poly_degree),
1452  * @endcode
1453  *
1454  * The Finite Element System is composed of dim continuous displacement
1455  * DOFs.
1456  *
1457  * @code
1458  * fe(FE_Q<dim>(parameters.poly_degree), dim), // displacement
1459  * dof_handler_ref(triangulation),
1460  * dofs_per_cell (fe.dofs_per_cell),
1461  * u_fe(first_u_component),
1462  * dofs_per_block(n_blocks),
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())
1467  * {
1468  *
1469  * }
1470  *
1471  * @endcode
1472  *
1473  * The class destructor simply clears the data held by the DOFHandler
1474  *
1475  * @code
1476  * template <int dim,typename NumberType>
1477  * Solid<dim,NumberType>::~Solid()
1478  * {
1479  * dof_handler_ref.clear();
1480  * }
1481  *
1482  *
1483  * @endcode
1484  *
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
1488  * step size.
1489  *
1490 
1491  *
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)
1494  * increment.
1495  *
1496 
1497  *
1498  *
1499  * @code
1500  * template <int dim,typename NumberType>
1502  * {
1503  * make_grid();
1504  * system_setup();
1505  * output_results();
1506  * time.increment();
1507  *
1508  * @endcode
1509  *
1510  * We then declare the incremental solution update @f$\varDelta
1511  * \mathbf{\Xi}:= \{\varDelta \mathbf{u}\}@f$ and start the loop over the
1512  * time domain.
1513  *
1514 
1515  *
1516  * At the beginning, we reset the solution update for this time step...
1517  *
1518  * @code
1519  * BlockVector<double> solution_delta(dofs_per_block);
1520  * while (time.current() <= time.end())
1521  * {
1522  * solution_delta = 0.0;
1523  *
1524  * @endcode
1525  *
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$...
1529  *
1530  * @code
1531  * solve_nonlinear_timestep(solution_delta);
1532  * solution_n += solution_delta;
1533  *
1534  * @endcode
1535  *
1536  * ...and plot the results before moving on happily to the next time
1537  * step:
1538  *
1539  * @code
1540  * output_results();
1541  * time.increment();
1542  * }
1543  *
1544  * @endcode
1545  *
1546  * Lastly, we print the vertical tip displacement of the Cook cantilever
1547  * after the full load is applied
1548  *
1549  * @code
1550  * print_vertical_tip_displacement();
1551  * }
1552  *
1553  *
1554  * @endcode
1555  *
1556  *
1557  * <a name="Privateinterface"></a>
1558  * <h3>Private interface</h3>
1559  *
1560 
1561  *
1562  *
1563  * <a name="Solidmake_grid"></a>
1564  * <h4>Solid::make_grid</h4>
1565  *
1566 
1567  *
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
1572  * ID number.
1573  *
1574 
1575  *
1576  * We then determine the volume of the reference configuration and print it
1577  * for comparison.
1578  *
1579 
1580  *
1581  *
1582  * @code
1583  * template <int dim>
1584  * Point<dim> grid_y_transform (const Point<dim> &pt_in)
1585  * {
1586  * const double &x = pt_in[0];
1587  * const double &y = pt_in[1];
1588  *
1589  * const double y_upper = 44.0 + (16.0/48.0)*x; // Line defining upper edge of beam
1590  * const double y_lower = 0.0 + (44.0/48.0)*x; // Line defining lower edge of beam
1591  * const double theta = y/44.0; // Fraction of height along left side of beam
1592  * const double y_transform = (1-theta)*y_lower + theta*y_upper; // Final transformation
1593  *
1594  * Point<dim> pt_out = pt_in;
1595  * pt_out[1] = y_transform;
1596  *
1597  * return pt_out;
1598  * }
1599  *
1600  * template <int dim,typename NumberType>
1601  * void Solid<dim,NumberType>::make_grid()
1602  * {
1603  * @endcode
1604  *
1605  * Divide the beam, but only along the x- and y-coordinate directions
1606  *
1607  * @code
1608  * std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge);
1609  * @endcode
1610  *
1611  * Only allow one element through the thickness
1612  * (modelling a plane strain condition)
1613  *
1614  * @code
1615  * if (dim == 3)
1616  * repetitions[dim-1] = 1;
1617  *
1618  * const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, -0.5) : Point<dim>(0.0, 0.0));
1619  * const Point<dim> top_right = (dim == 3 ? Point<dim>(48.0, 44.0, 0.5) : Point<dim>(48.0, 44.0));
1620  *
1622  * repetitions,
1623  * bottom_left,
1624  * top_right);
1625  *
1626  * @endcode
1627  *
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)
1635  *
1636  * @code
1637  * const double tol_boundary = 1e-6;
1639  * triangulation.begin_active(), endc = triangulation.end();
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)
1644  * {
1645  * if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary)
1646  * cell->face(face)->set_boundary_id(1); // -X faces
1647  * else if (std::abs(cell->face(face)->center()[0] - 48.0) < tol_boundary)
1648  * cell->face(face)->set_boundary_id(11); // +X faces
1649  * else if (std::abs(std::abs(cell->face(face)->center()[0]) - 0.5) < tol_boundary)
1650  * cell->face(face)->set_boundary_id(2); // +Z and -Z faces
1651  * }
1652  *
1653  * @endcode
1654  *
1655  * Transform the hyper-rectangle into the beam shape
1656  *
1657  * @code
1658  * GridTools::transform(&grid_y_transform<dim>, triangulation);
1659  *
1660  * GridTools::scale(parameters.scale, triangulation);
1661  *
1662  * vol_reference = GridTools::volume(triangulation);
1663  * vol_current = vol_reference;
1664  * std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
1665  * }
1666  *
1667  *
1668  * @endcode
1669  *
1670  *
1671  * <a name="Solidsystem_setup"></a>
1672  * <h4>Solid::system_setup</h4>
1673  *
1674 
1675  *
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.
1679  *
1680  * @code
1681  * template <int dim,typename NumberType>
1682  * void Solid<dim,NumberType>::system_setup()
1683  * {
1684  * timer.enter_subsection("Setup system");
1685  *
1686  * std::vector<unsigned int> block_component(n_components, u_dof); // Displacement
1687  *
1688  * @endcode
1689  *
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.
1692  *
1693  * @code
1694  * dof_handler_ref.distribute_dofs(fe);
1695  * DoFRenumbering::Cuthill_McKee(dof_handler_ref);
1696  * DoFRenumbering::component_wise(dof_handler_ref, block_component);
1697  * DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block,
1698  * block_component);
1699  *
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()
1703  * << std::endl;
1704  *
1705  * @endcode
1706  *
1707  * Setup the sparsity pattern and tangent matrix
1708  *
1709  * @code
1710  * tangent_matrix.clear();
1711  * {
1712  * const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
1713  *
1715  *
1716  * csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
1717  * csp.collect_sizes();
1718  *
1719  * @endcode
1720  *
1721  * Naturally, for a one-field vector-valued problem, all of the
1722  * components of the system are coupled.
1723  *
1724  * @code
1726  * for (unsigned int ii = 0; ii < n_components; ++ii)
1727  * for (unsigned int jj = 0; jj < n_components; ++jj)
1728  * coupling[ii][jj] = DoFTools::always;
1729  * DoFTools::make_sparsity_pattern(dof_handler_ref,
1730  * coupling,
1731  * csp,
1732  * constraints,
1733  * false);
1734  * sparsity_pattern.copy_from(csp);
1735  * }
1736  *
1737  * tangent_matrix.reinit(sparsity_pattern);
1738  *
1739  * @endcode
1740  *
1741  * We then set up storage vectors
1742  *
1743  * @code
1744  * system_rhs.reinit(dofs_per_block);
1745  * system_rhs.collect_sizes();
1746  *
1747  * solution_n.reinit(dofs_per_block);
1748  * solution_n.collect_sizes();
1749  *
1750  * @endcode
1751  *
1752  * ...and finally set up the quadrature
1753  * point history:
1754  *
1755  * @code
1756  * setup_qph();
1757  *
1758  * timer.leave_subsection();
1759  * }
1760  *
1761  *
1762  * @endcode
1763  *
1764  *
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.
1769  *
1770 
1771  *
1772  * Firstly the actual QPH data objects are created. This must be done only
1773  * once the grid is refined to its finest level.
1774  *
1775  * @code
1776  * template <int dim,typename NumberType>
1777  * void Solid<dim,NumberType>::setup_qph()
1778  * {
1779  * std::cout << " Setting up quadrature point data..." << std::endl;
1780  *
1781  * quadrature_point_history.initialize(triangulation.begin_active(),
1782  * triangulation.end(),
1783  * n_q_points);
1784  *
1785  * @endcode
1786  *
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.
1790  *
1791  * @code
1792  * for (typename Triangulation<dim>::active_cell_iterator cell =
1793  * triangulation.begin_active(); cell != triangulation.end(); ++cell)
1794  * {
1795  * const std::vector<std::shared_ptr<PointHistory<dim,NumberType> > > lqph =
1796  * quadrature_point_history.get_data(cell);
1797  * Assert(lqph.size() == n_q_points, ExcInternalError());
1798  *
1799  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1800  * lqph[q_point]->setup_lqp(parameters);
1801  * }
1802  * }
1803  *
1804  *
1805  * @endcode
1806  *
1807  *
1808  * <a name="Solidsolve_nonlinear_timestep"></a>
1809  * <h4>Solid::solve_nonlinear_timestep</h4>
1810  *
1811 
1812  *
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.
1816  *
1817  * @code
1818  * template <int dim,typename NumberType>
1819  * void
1820  * Solid<dim,NumberType>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
1821  * {
1822  * std::cout << std::endl << "Timestep " << time.get_timestep() << " @ "
1823  * << time.current() << "s" << std::endl;
1824  *
1825  * BlockVector<double> newton_update(dofs_per_block);
1826  *
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();
1833  *
1834  * print_conv_header();
1835  *
1836  * @endcode
1837  *
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
1848  * is attained.
1849  *
1850  * @code
1851  * unsigned int newton_iteration = 0;
1852  * for (; newton_iteration < parameters.max_iterations_NR;
1853  * ++newton_iteration)
1854  * {
1855  * std::cout << " " << std::setw(2) << newton_iteration << " " << std::flush;
1856  *
1857  * @endcode
1858  *
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:
1862  *
1863  * @code
1864  * make_constraints(newton_iteration);
1865  * assemble_system(solution_delta);
1866  *
1867  * get_error_residual(error_residual);
1868  *
1869  * if (newton_iteration == 0)
1870  * error_residual_0 = error_residual;
1871  *
1872  * @endcode
1873  *
1874  * We can now determine the normalised residual error and check for
1875  * solution convergence:
1876  *
1877  * @code
1878  * error_residual_norm = error_residual;
1879  * error_residual_norm.normalise(error_residual_0);
1880  *
1881  * if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u
1882  * && error_residual_norm.u <= parameters.tol_f)
1883  * {
1884  * std::cout << " CONVERGED! " << std::endl;
1885  * print_conv_footer();
1886  *
1887  * break;
1888  * }
1889  *
1890  * const std::pair<unsigned int, double>
1891  * lin_solver_output = solve_linear_system(newton_update);
1892  *
1893  * get_error_update(newton_update, error_update);
1894  * if (newton_iteration == 0)
1895  * error_update_0 = error_update;
1896  *
1897  * @endcode
1898  *
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:
1903  *
1904  * @code
1905  * error_update_norm = error_update;
1906  * error_update_norm.normalise(error_update_0);
1907  *
1908  * solution_delta += newton_update;
1909  *
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;
1916  * }
1917  *
1918  * @endcode
1919  *
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
1927  * problem happened.
1928  *
1929  * @code
1930  * AssertThrow (newton_iteration <= parameters.max_iterations_NR,
1931  * ExcMessage("No convergence in nonlinear solver!"));
1932  * }
1933  *
1934  *
1935  * @endcode
1936  *
1937  *
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>
1940  *
1941 
1942  *
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:
1946  *
1947  * @code
1948  * template <int dim,typename NumberType>
1949  * void Solid<dim,NumberType>::print_conv_header()
1950  * {
1951  * static const unsigned int l_width = 87;
1952  *
1953  * for (unsigned int i = 0; i < l_width; ++i)
1954  * std::cout << "_";
1955  * std::cout << std::endl;
1956  *
1957  * std::cout << " SOLVER STEP "
1958  * << " | LIN_IT LIN_RES RES_NORM "
1959  * << " RES_U NU_NORM "
1960  * << " NU_U " << std::endl;
1961  *
1962  * for (unsigned int i = 0; i < l_width; ++i)
1963  * std::cout << "_";
1964  * std::cout << std::endl;
1965  * }
1966  *
1967  *
1968  *
1969  * template <int dim,typename NumberType>
1970  * void Solid<dim,NumberType>::print_conv_footer()
1971  * {
1972  * static const unsigned int l_width = 87;
1973  *
1974  * for (unsigned int i = 0; i < l_width; ++i)
1975  * std::cout << "_";
1976  * std::cout << std::endl;
1977  *
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
1982  * << std::endl;
1983  * }
1984  *
1985  * @endcode
1986  *
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
1989  * beam.
1990  *
1991  * @code
1992  * template <int dim,typename NumberType>
1993  * void Solid<dim,NumberType>::print_vertical_tip_displacement()
1994  * {
1995  * static const unsigned int l_width = 87;
1996  *
1997  * for (unsigned int i = 0; i < l_width; ++i)
1998  * std::cout << "_";
1999  * std::cout << std::endl;
2000  *
2001  * Point<dim> soln_pt (48.0*parameters.scale,60.0*parameters.scale);
2002  * if (dim == 3)
2003  * soln_pt[2] = 0.5*parameters.scale;
2004  * double vertical_tip_displacement = 0.0;
2005  * double vertical_tip_displacement_check = 0.0;
2006  *
2007  * typename DoFHandler<dim>::active_cell_iterator cell =
2008  * dof_handler_ref.begin_active(), endc = dof_handler_ref.end();
2009  * for (; cell != endc; ++cell)
2010  * {
2011  * @endcode
2012  *
2013  * if (cell->point_inside(soln_pt) == true)
2014  *
2015  * @code
2016  * for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
2017  * if (cell->vertex(v).distance(soln_pt) < 1e-6)
2018  * {
2019  * @endcode
2020  *
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
2025  *
2026  * @code
2027  * vertical_tip_displacement = solution_n(cell->vertex_dof_index(v,u_dof+1));
2028  *
2029  * @endcode
2030  *
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
2034  *
2035  * @code
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);
2042  *
2043  * @endcode
2044  *
2045  * Extract y-component of solution at given point
2046  *
2047  * @code
2048  * std::vector< Tensor<1,dim> > soln_values (soln_qrule.size());
2049  * fe_values_soln[u_fe].get_function_values(solution_n,
2050  * soln_values);
2051  * vertical_tip_displacement_check = soln_values[0][u_dof+1];
2052  *
2053  * break;
2054  * }
2055  * }
2056  * AssertThrow(vertical_tip_displacement > 0.0, ExcMessage("Found no cell with point inside!"))
2057  *
2058  * std::cout << "Vertical tip displacement: " << vertical_tip_displacement
2059  * << "\t Check: " << vertical_tip_displacement_check
2060  * << std::endl;
2061  * }
2062  *
2063  *
2064  * @endcode
2065  *
2066  *
2067  * <a name="Solidget_error_residual"></a>
2068  * <h4>Solid::get_error_residual</h4>
2069  *
2070 
2071  *
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.
2076  *
2077  * @code
2078  * template <int dim,typename NumberType>
2079  * void Solid<dim,NumberType>::get_error_residual(Errors &error_residual)
2080  * {
2081  * BlockVector<double> error_res(dofs_per_block);
2082  *
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);
2086  *
2087  * error_residual.norm = error_res.l2_norm();
2088  * error_residual.u = error_res.block(u_dof).l2_norm();
2089  * }
2090  *
2091  *
2092  * @endcode
2093  *
2094  *
2095  * <a name="Solidget_error_udpate"></a>
2096  * <h4>Solid::get_error_udpate</h4>
2097  *
2098 
2099  *
2100  * Determine the true Newton update error for the problem
2101  *
2102  * @code
2103  * template <int dim,typename NumberType>
2104  * void Solid<dim,NumberType>::get_error_update(const BlockVector<double> &newton_update,
2105  * Errors &error_update)
2106  * {
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);
2111  *
2112  * error_update.norm = error_ud.l2_norm();
2113  * error_update.u = error_ud.block(u_dof).l2_norm();
2114  * }
2115  *
2116  *
2117  *
2118  * @endcode
2119  *
2120  *
2121  * <a name="Solidget_total_solution"></a>
2122  * <h4>Solid::get_total_solution</h4>
2123  *
2124 
2125  *
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.
2129  *
2130  * @code
2131  * template <int dim,typename NumberType>
2132  * BlockVector<double>
2133  * Solid<dim,NumberType>::get_total_solution(const BlockVector<double> &solution_delta) const
2134  * {
2135  * BlockVector<double> solution_total(solution_n);
2136  * solution_total += solution_delta;
2137  * return solution_total;
2138  * }
2139  *
2140  *
2141  * @endcode
2142  *
2143  *
2144  * <a name="Solidassemble_system"></a>
2145  * <h4>Solid::assemble_system</h4>
2146  *
2147 
2148  *
2149  *
2150  * @code
2151  * template <int dim,typename NumberType>
2152  * struct Assembler_Base
2153  * {
2154  * virtual ~Assembler_Base() {}
2155  *
2156  * @endcode
2157  *
2158  * Here we deal with the tangent matrix assembly structures. The
2159  * PerTaskData object stores local contributions.
2160  *
2161  * @code
2162  * struct PerTaskData_ASM
2163  * {
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;
2168  *
2169  * PerTaskData_ASM(const Solid<dim,NumberType> *solid)
2170  * :
2171  * solid (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)
2175  * {}
2176  *
2177  * void reset()
2178  * {
2179  * cell_matrix = 0.0;
2180  * cell_rhs = 0.0;
2181  * }
2182  * };
2183  *
2184  * @endcode
2185  *
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
2189  * assembly.
2190  *
2191  * @code
2192  * struct ScratchData_ASM
2193  * {
2194  * const BlockVector<double> &solution_total;
2195  * std::vector<Tensor<2, dim,NumberType> > solution_grads_u_total;
2196  *
2197  * FEValues<dim> fe_values_ref;
2198  * FEFaceValues<dim> fe_face_values_ref;
2199  *
2200  * std::vector<std::vector<Tensor<2, dim,NumberType> > > grad_Nx;
2201  * std::vector<std::vector<SymmetricTensor<2,dim,NumberType> > > symm_grad_Nx;
2202  *
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)
2209  * :
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))
2219  * {}
2220  *
2221  * ScratchData_ASM(const ScratchData_ASM &rhs)
2222  * :
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)
2233  * {}
2234  *
2235  * void reset()
2236  * {
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)
2240  * {
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());
2245  *
2246  * solution_grads_u_total[q_point] = Tensor<2,dim,NumberType>();
2247  * for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
2248  * {
2249  * grad_Nx[q_point][k] = Tensor<2,dim,NumberType>();
2250  * symm_grad_Nx[q_point][k] = SymmetricTensor<2,dim,NumberType>();
2251  * }
2252  * }
2253  * }
2254  *
2255  * };
2256  *
2257  * @endcode
2258  *
2259  * Of course, we still have to define how we assemble the tangent matrix
2260  * contribution for a single cell.
2261  *
2262  * @code
2263  * void
2264  * assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2265  * ScratchData_ASM &scratch,
2266  * PerTaskData_ASM &data)
2267  * {
2268  * @endcode
2269  *
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.
2277  *
2278  * @code
2279  * assemble_system_tangent_residual_one_cell(cell, scratch, data);
2280  * assemble_neumann_contribution_one_cell(cell, scratch, data);
2281  * }
2282  *
2283  * @endcode
2284  *
2285  * This function adds the local contribution to the system matrix.
2286  *
2287  * @code
2288  * void
2289  * copy_local_to_global_ASM(const PerTaskData_ASM &data)
2290  * {
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;
2294  *
2295  * constraints.distribute_local_to_global(
2296  * data.cell_matrix, data.cell_rhs,
2297  * data.local_dof_indices,
2298  * tangent_matrix, system_rhs);
2299  * }
2300  *
2301  * protected:
2302  *
2303  * @endcode
2304  *
2305  * This function needs to exist in the base class for
2306  * Workstream to work with a reference to the base class.
2307  *
2308  * @code
2309  * virtual void
2310  * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &/*cell*/,
2311  * ScratchData_ASM &/*scratch*/,
2312  * PerTaskData_ASM &/*data*/)
2313  * {
2314  * AssertThrow(false, ExcPureFunctionCalled());
2315  * }
2316  *
2317  * void
2318  * assemble_neumann_contribution_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2319  * ScratchData_ASM &scratch,
2320  * PerTaskData_ASM &data)
2321  * {
2322  * @endcode
2323  *
2324  * Aliases for data referenced from the Solid class
2325  *
2326  * @code
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 &parameters = 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;
2333  *
2334  * @endcode
2335  *
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.
2339  *
2340  * @code
2341  * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2342  * ++face)
2343  * if (cell->face(face)->at_boundary() == true
2344  * && cell->face(face)->boundary_id() == 11)
2345  * {
2346  * scratch.fe_face_values_ref.reinit(cell, face);
2347  *
2348  * for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
2349  * ++f_q_point)
2350  * {
2351  * @endcode
2352  *
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.
2358  *
2359 
2360  *
2361  * Note that the contributions to the right hand side vector we
2362  * compute here only exist in the displacement components of the
2363  * vector.
2364  *
2365  * @code
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;
2369  * dir[1] = 1.0;
2370  * const Tensor<1, dim> traction = magnitude*dir;
2371  *
2372  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2373  * {
2374  * const unsigned int i_group =
2375  * fe.system_to_base_index(i).first.first;
2376  *
2377  * if (i_group == u_dof)
2378  * {
2379  * const unsigned int component_i =
2380  * fe.system_to_component_index(i).first;
2381  * const double Ni =
2382  * scratch.fe_face_values_ref.shape_value(i,
2383  * f_q_point);
2384  * const double JxW = scratch.fe_face_values_ref.JxW(
2385  * f_q_point);
2386  *
2387  * data.cell_rhs(i) += (Ni * traction[component_i])
2388  * * JxW;
2389  * }
2390  * }
2391  * }
2392  * }
2393  * }
2394  *
2395  * };
2396  *
2397  * template <int dim>
2398  * struct Assembler<dim,double> : Assembler_Base<dim,double>
2399  * {
2400  * typedef double NumberType;
2401  * using typename Assembler_Base<dim,NumberType>::ScratchData_ASM;
2402  * using typename Assembler_Base<dim,NumberType>::PerTaskData_ASM;
2403  *
2404  * virtual ~Assembler() {}
2405  *
2406  * virtual void
2407  * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2408  * ScratchData_ASM &scratch,
2409  * PerTaskData_ASM &data)
2410  * {
2411  * @endcode
2412  *
2413  * Aliases for data referenced from the Solid class
2414  *
2415  * @code
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;
2421  *
2422  * data.reset();
2423  * scratch.reset();
2424  * scratch.fe_values_ref.reinit(cell);
2425  * cell->get_dof_indices(data.local_dof_indices);
2426  *
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());
2430  *
2431  * @endcode
2432  *
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:
2436  *
2437  * @code
2438  * scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
2439  * scratch.solution_grads_u_total);
2440  *
2441  * @endcode
2442  *
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.
2447  *
2448 
2449  *
2450  * In doing so, we first extract some configuration dependent variables
2451  * from our QPH history objects for the current quadrature point.
2452  *
2453  * @code
2454  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2455  * {
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());
2463  *
2464  * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2465  * {
2466  * const unsigned int k_group = fe.system_to_base_index(k).first.first;
2467  *
2468  * if (k_group == u_dof)
2469  * {
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]);
2472  * }
2473  * else
2474  * Assert(k_group <= u_dof, ExcInternalError());
2475  * }
2476  *
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);
2480  *
2481  * @endcode
2482  *
2483  * Next we define some aliases to make the assembly process easier to
2484  * follow
2485  *
2486  * @code
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);
2490  *
2491  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2492  * {
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;
2495  *
2496  * if (i_group == u_dof)
2497  * data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
2498  * else
2499  * Assert(i_group <= u_dof, ExcInternalError());
2500  *
2501  * for (unsigned int j = 0; j <= i; ++j)
2502  * {
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;
2505  *
2506  * @endcode
2507  *
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:
2512  *
2513  * @code
2514  * if ((i_group == j_group) && (i_group == u_dof))
2515  * {
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;
2521  * }
2522  * else
2523  * Assert((i_group <= u_dof) && (j_group <= u_dof),
2524  * ExcInternalError());
2525  * }
2526  * }
2527  * }
2528  *
2529  *
2530  * @endcode
2531  *
2532  * Finally, we need to copy the lower half of the local matrix into the
2533  * upper half:
2534  *
2535  * @code
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);
2539  * }
2540  *
2541  * };
2542  *
2543  * #ifdef ENABLE_SACADO_FORMULATION
2544  *
2545  *
2546  * template <int dim>
2547  * struct Assembler<dim,Sacado::Fad::DFad<double> > : Assembler_Base<dim,Sacado::Fad::DFad<double> >
2548  * {
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;
2552  *
2553  * virtual ~Assembler() {}
2554  *
2555  * virtual void
2556  * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2557  * ScratchData_ASM &scratch,
2558  * PerTaskData_ASM &data)
2559  * {
2560  * @endcode
2561  *
2562  * Aliases for data referenced from the Solid class
2563  *
2564  * @code
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;
2570  *
2571  * data.reset();
2572  * scratch.reset();
2573  * scratch.fe_values_ref.reinit(cell);
2574  * cell->get_dof_indices(data.local_dof_indices);
2575  *
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());
2579  *
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());
2585  *
2586  * @endcode
2587  *
2588  * We now retrieve a set of degree-of-freedom values that
2589  * have the operations that are performed with them tracked.
2590  *
2591  * @code
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]);
2595  *
2596  * @endcode
2597  *
2598  * Compute all values, gradients etc. based on sensitive
2599  * AD degree-of-freedom values.
2600  *
2601  * @code
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);
2605  *
2606  * @endcode
2607  *
2608  * Accumulate the residual value for each degree of freedom.
2609  * Note: Its important that the vectors is initialised (zero'd) correctly.
2610  *
2611  * @code
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)
2614  * {
2615  * const Tensor<2,dim,ADNumberType> &grad_u = scratch.solution_grads_u_total[q_point];
2617  * const ADNumberType det_F = determinant(F);
2620  * const Tensor<2,dim,ADNumberType> F_inv = invert(F);
2621  * Assert(det_F > ADNumberType(0.0), ExcInternalError());
2622  *
2623  * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2624  * {
2625  * const unsigned int k_group = fe.system_to_base_index(k).first.first;
2626  *
2627  * if (k_group == u_dof)
2628  * {
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]);
2631  * }
2632  * else
2633  * Assert(k_group <= u_dof, ExcInternalError());
2634  * }
2635  *
2636  * const SymmetricTensor<2,dim,ADNumberType> tau = lqph[q_point]->get_tau(det_F,b_bar);
2637  *
2638  * @endcode
2639  *
2640  * Next we define some position-dependent aliases, again to
2641  * make the assembly process easier to follow.
2642  *
2643  * @code
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);
2646  *
2647  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2648  * {
2649  * const unsigned int i_group = fe.system_to_base_index(i).first.first;
2650  *
2651  * if (i_group == u_dof)
2652  * residual_ad[i] += (symm_grad_Nx[i] * tau) * JxW;
2653  * else
2654  * Assert(i_group <= u_dof, ExcInternalError());
2655  * }
2656  * }
2657  *
2658  * for (unsigned int I=0; I<n_independent_variables; ++I)
2659  * {
2660  * const ADNumberType &residual_I = residual_ad[I];
2661  * data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
2662  * for (unsigned int J=0; J<n_independent_variables; ++J)
2663  * {
2664  * @endcode
2665  *
2666  * Compute the gradients of the residual entry [forward-mode]
2667  *
2668  * @code
2669  * data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
2670  * }
2671  * }
2672  * }
2673  *
2674  * };
2675  *
2676  *
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> > >
2679  * {
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;
2684  *
2685  * virtual ~Assembler() {}
2686  *
2687  * virtual void
2688  * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2689  * ScratchData_ASM &scratch,
2690  * PerTaskData_ASM &data)
2691  * {
2692  * @endcode
2693  *
2694  * Aliases for data referenced from the Solid class
2695  *
2696  * @code
2697  * const unsigned int &n_q_points = data.solid->n_q_points;
2698  * const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
2699  *
2700  * data.reset();
2701  * scratch.reset();
2702  * scratch.fe_values_ref.reinit(cell);
2703  * cell->get_dof_indices(data.local_dof_indices);
2704  *
2705  * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
2706  * data.solid->quadrature_point_history.get_data(cell);
2707  * Assert(lqph.size() == n_q_points, ExcInternalError());
2708  *
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());
2714  *
2715  * @endcode
2716  *
2717  * We now retrieve a set of degree-of-freedom values that
2718  * have the operations that are performed with them tracked.
2719  *
2720  * @code
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]);
2724  *
2725  * @endcode
2726  *
2727  * Compute all values, gradients etc. based on sensitive
2728  * AD degree-of-freedom values.
2729  *
2730  * @code
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);
2734  *
2735  * @endcode
2736  *
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.
2741  *
2742  * @code
2743  * ADNumberType cell_energy_ad = ADNumberType(0.0);
2744  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2745  * {
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());
2752  *
2753  * @endcode
2754  *
2755  * Next we define some position-dependent aliases, again to
2756  * make the assembly process easier to follow.
2757  *
2758  * @code
2759  * const double JxW = scratch.fe_values_ref.JxW(q_point);
2760  *
2761  * const ADNumberType Psi = lqph[q_point]->get_Psi(det_F,b_bar);
2762  *
2763  * @endcode
2764  *
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
2768  * cell energy.
2769  *
2770  * @code
2771  * cell_energy_ad += Psi * JxW;
2772  * }
2773  *
2774  * @endcode
2775  *
2776  * Compute derivatives of reverse-mode AD variables
2777  *
2778  * @code
2779  * ADNumberType::Gradcomp();
2780  *
2781  * for (unsigned int I=0; I<n_independent_variables; ++I)
2782  * {
2783  * @endcode
2784  *
2785  * This computes the adjoint df/dX_{i} [reverse-mode]
2786  *
2787  * @code
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)
2791  * {
2792  * @endcode
2793  *
2794  * Compute the gradients of the residual entry [forward-mode]
2795  *
2796  * @code
2797  * data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
2798  * }
2799  * }
2800  * }
2801  *
2802  * };
2803  *
2804  *
2805  * #endif
2806  *
2807  *
2808  * @endcode
2809  *
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.
2815  *
2816  * @code
2817  * template <int dim,typename NumberType>
2818  * void Solid<dim,NumberType>::assemble_system(const BlockVector<double> &solution_delta)
2819  * {
2820  * timer.enter_subsection("Assemble linear system");
2821  * std::cout << " ASM " << std::flush;
2822  *
2823  * tangent_matrix = 0.0;
2824  * system_rhs = 0.0;
2825  *
2826  * const UpdateFlags uf_cell(update_gradients |
2827  * update_JxW_values);
2828  * const UpdateFlags uf_face(update_values |
2829  * update_JxW_values);
2830  *
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;
2835  *
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,
2841  * scratch_data,
2842  * per_task_data);
2843  *
2844  * timer.leave_subsection();
2845  * }
2846  *
2847  *
2848  * @endcode
2849  *
2850  *
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.
2859  *
2860  * @code
2861  * template <int dim,typename NumberType>
2862  * void Solid<dim,NumberType>::make_constraints(const int &it_nr)
2863  * {
2864  * std::cout << " CST " << std::flush;
2865  *
2866  * @endcode
2867  *
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.
2872  *
2873  * @code
2874  * if (it_nr > 1)
2875  * return;
2876  * constraints.clear();
2877  * const bool apply_dirichlet_bc = (it_nr == 0);
2878  *
2879  * @endcode
2880  *
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.
2888  *
2889 
2890  *
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:
2899  *
2900 
2901  *
2902  * Fixed left hand side of the beam
2903  *
2904  * @code
2905  * {
2906  * const int boundary_id = 1;
2907  *
2908  * if (apply_dirichlet_bc == true)
2909  * VectorTools::interpolate_boundary_values(dof_handler_ref,
2910  * boundary_id,
2912  * constraints,
2913  * fe.component_mask(u_fe));
2914  * else
2915  * VectorTools::interpolate_boundary_values(dof_handler_ref,
2916  * boundary_id,
2918  * constraints,
2919  * fe.component_mask(u_fe));
2920  * }
2921  *
2922  * @endcode
2923  *
2924  * Zero Z-displacement through thickness direction
2925  * This corresponds to a plane strain condition being imposed on the beam
2926  *
2927  * @code
2928  * if (dim == 3)
2929  * {
2930  * const int boundary_id = 2;
2931  * const FEValuesExtractors::Scalar z_displacement(2);
2932  *
2933  * if (apply_dirichlet_bc == true)
2934  * VectorTools::interpolate_boundary_values(dof_handler_ref,
2935  * boundary_id,
2937  * constraints,
2938  * fe.component_mask(z_displacement));
2939  * else
2940  * VectorTools::interpolate_boundary_values(dof_handler_ref,
2941  * boundary_id,
2943  * constraints,
2944  * fe.component_mask(z_displacement));
2945  * }
2946  *
2947  * constraints.close();
2948  * }
2949  *
2950  * @endcode
2951  *
2952  *
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.
2957  *
2958  * @code
2959  * template <int dim,typename NumberType>
2960  * std::pair<unsigned int, double>
2961  * Solid<dim,NumberType>::solve_linear_system(BlockVector<double> &newton_update)
2962  * {
2963  * BlockVector<double> A(dofs_per_block);
2964  * BlockVector<double> B(dofs_per_block);
2965  *
2966  * unsigned int lin_it = 0;
2967  * double lin_res = 0.0;
2968  *
2969  * @endcode
2970  *
2971  * We solve for the incremental displacement @f$d\mathbf{u}@f$.
2972  *
2973  * @code
2974  * {
2975  * timer.enter_subsection("Linear solver");
2976  * std::cout << " SLV " << std::flush;
2977  * if (parameters.type_lin == "CG")
2978  * {
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
2982  * * system_rhs.block(u_dof).l2_norm();
2983  *
2984  * SolverControl solver_control(solver_its, tol_sol);
2985  *
2987  * SolverCG<Vector<double> > solver_CG(solver_control, GVM);
2988  *
2989  * @endcode
2990  *
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.
2996  *
2997  * @code
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));
3002  *
3003  * solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
3004  * newton_update.block(u_dof),
3005  * system_rhs.block(u_dof),
3006  * preconditioner);
3007  *
3008  * lin_it = solver_control.last_step();
3009  * lin_res = solver_control.last_value();
3010  * }
3011  * else if (parameters.type_lin == "Direct")
3012  * {
3013  * @endcode
3014  *
3015  * Otherwise if the problem is small
3016  * enough, a direct solver can be
3017  * utilised.
3018  *
3019  * @code
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));
3023  *
3024  * lin_it = 1;
3025  * lin_res = 0.0;
3026  * }
3027  * else
3028  * Assert (false, ExcMessage("Linear solver type not implemented"));
3029  *
3030  * timer.leave_subsection();
3031  * }
3032  *
3033  * @endcode
3034  *
3035  * Now that we have the displacement update, distribute the constraints
3036  * back to the Newton update:
3037  *
3038  * @code
3039  * constraints.distribute(newton_update);
3040  *
3041  * return std::make_pair(lin_it, lin_res);
3042  * }
3043  *
3044  * @endcode
3045  *
3046  *
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.
3052  *
3053  * @code
3054  * template <int dim,typename NumberType>
3055  * void Solid<dim,NumberType>::output_results() const
3056  * {
3057  * DataOut<dim> data_out;
3058  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3059  * data_component_interpretation(dim,
3060  * DataComponentInterpretation::component_is_part_of_vector);
3061  *
3062  * std::vector<std::string> solution_name(dim, "displacement");
3063  *
3064  * data_out.attach_dof_handler(dof_handler_ref);
3065  * data_out.add_data_vector(solution_n,
3066  * solution_name,
3067  * DataOut<dim>::type_dof_data,
3068  * data_component_interpretation);
3069  *
3070  * @endcode
3071  *
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.
3080  *
3081  * @code
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);
3087  *
3088  * std::ostringstream filename;
3089  * filename << "solution-" << time.get_timestep() << ".vtk";
3090  *
3091  * std::ofstream output(filename.str().c_str());
3092  * data_out.write_vtk(output);
3093  * }
3094  *
3095  * }
3096  *
3097  *
3098  * @endcode
3099  *
3100  *
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.
3105  *
3106  * @code
3107  * int main (int argc, char *argv[])
3108  * {
3109  * using namespace dealii;
3110  * using namespace Cook_Membrane;
3111  *
3112  * const unsigned int dim = 2;
3113  * try
3114  * {
3115  * deallog.depth_console(0);
3116  * Parameters::AllParameters parameters("parameters.prm");
3117  * if (parameters.automatic_differentiation_order == 0)
3118  * {
3119  * std::cout << "Assembly method: Residual and linearisation are computed manually." << std::endl;
3120  *
3121  * @endcode
3122  *
3123  * Allow multi-threading
3124  *
3125  * @code
3126  * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3127  * ::numbers::invalid_unsigned_int);
3128  *
3129  * typedef double NumberType;
3130  * Solid<dim,NumberType> solid_3d(parameters);
3131  * solid_3d.run();
3132  * }
3133  * #ifdef ENABLE_SACADO_FORMULATION
3134  * else if (parameters.automatic_differentiation_order == 1)
3135  * {
3136  * std::cout << "Assembly method: Residual computed manually; linearisation performed using AD." << std::endl;
3137  *
3138  * @endcode
3139  *
3140  * Allow multi-threading
3141  *
3142  * @code
3143  * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3144  * ::numbers::invalid_unsigned_int);
3145  *
3146  * typedef Sacado::Fad::DFad<double> NumberType;
3147  * Solid<dim,NumberType> solid_3d(parameters);
3148  * solid_3d.run();
3149  * }
3150  * else if (parameters.automatic_differentiation_order == 2)
3151  * {
3152  * std::cout << "Assembly method: Residual and linearisation computed using AD." << std::endl;
3153  *
3154  * @endcode
3155  *
3156  * Sacado Rad-Fad is not thread-safe, so disable threading.
3157  * Parallisation using MPI would be possible though.
3158  *
3159  * @code
3160  * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3161  * 1);
3162  *
3163  * typedef Sacado::Rad::ADvar< Sacado::Fad::DFad<double> > NumberType;
3164  * Solid<dim,NumberType> solid_3d(parameters);
3165  * solid_3d.run();
3166  * }
3167  * #endif
3168  * else
3169  * {
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."));
3174  * }
3175  * }
3176  * catch (std::exception &exc)
3177  * {
3178  * std::cerr << std::endl << std::endl
3179  * << "----------------------------------------------------"
3180  * << std::endl;
3181  * std::cerr << "Exception on processing: " << std::endl << exc.what()
3182  * << std::endl << "Aborting!" << std::endl
3183  * << "----------------------------------------------------"
3184  * << std::endl;
3185  *
3186  * return 1;
3187  * }
3188  * catch (...)
3189  * {
3190  * std::cerr << std::endl << std::endl
3191  * << "----------------------------------------------------"
3192  * << std::endl;
3193  * std::cerr << "Unknown exception!" << std::endl << "Aborting!"
3194  * << std::endl
3195  * << "----------------------------------------------------"
3196  * << std::endl;
3197  * return 1;
3198  * }
3199  *
3200  * return 0;
3201  * }
3202  * @endcode
3203 
3204 
3205 */
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
SolverCG
Definition: solver_cg.h:98
DEAL_II_VERSION_MAJOR
#define DEAL_II_VERSION_MAJOR
Definition: config.h:28
FE_Q
Definition: fe_q.h:554
std_cxx11
Definition: array.h:27
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:133
SymmetricTensor
Definition: symmetric_tensor.h:611
dealii
Definition: namespace_dealii.h:25
BlockVector< double >
DoFTools::count_dofs_per_block
void count_dofs_per_block(const DoFHandlerType &dof, std::vector< types::global_dof_index > &dofs_per_block, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition: dof_tools.cc:2001
mu
Definition: function_parser.h:32
Triangulation
Definition: tria.h:1109
FEValuesExtractors::Scalar
Definition: fe_values_extractors.h:95
TimerOutput::summary
@ summary
Definition: timer.h:605
SymmetricTensor::invert
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3467
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
DoFHandler::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
ComponentMask
Definition: component_mask.h:83
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
BlockDynamicSparsityPattern
Definition: block_sparsity_pattern.h:521
TimerOutput::wall_times
@ wall_times
Definition: timer.h:649
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
DoFTools::always
@ always
Definition: dof_tools.h:236
second
Point< 2 > second
Definition: grid_out.cc:4353
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler< dim >
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
Table
Definition: table.h:699
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
DoFHandler::distribute_dofs
virtual void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:1247
GrowingVectorMemory
Definition: vector_memory.h:320
BlockSparsityPattern::copy_from
void copy_from(const BlockDynamicSparsityPattern &dsp)
Definition: block_sparsity_pattern.cc:407
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
level
unsigned int level
Definition: grid_out.cc:4355
BlockVectorBase::collect_sizes
void collect_sizes()
TensorAccessors::extract
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
Definition: tensor_accessors.h:226
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
BlockMatrixBase::block
BlockType & block(const unsigned int row, const unsigned int column)
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
FEValuesExtractors::Vector
Definition: fe_values_extractors.h:150
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
Tensor
Definition: tensor.h:450
GridGenerator::subdivided_hyper_rectangle
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
BlockSparseMatrix::clear
void clear()
QGauss
Definition: quadrature_lib.h:40
DoFHandler::end
cell_iterator end() const
Definition: dof_handler.cc:951
BlockSparseMatrix::reinit
virtual void reinit(const BlockSparsityPattern &sparsity)
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
BlockVector::reinit
void reinit(const unsigned int n_blocks, const size_type block_size=0, const bool omit_zeroing_entries=false)
unsigned int
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AffineConstraints< double >
BlockSparsityPattern
Definition: block_sparsity_pattern.h:401
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
SymmetricTensor::symmetrize
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3547
SymmetricTensor::determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2645
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
GridTools::transform
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
DoFHandler::clear
virtual void clear()
Definition: dof_handler.cc:1352
Physics::Elasticity::Kinematics::F_iso
Tensor< 2, dim, Number > F_iso(const Tensor< 2, dim, Number > &F)
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
CellDataStorage
Definition: quadrature_point_data.h:64
DerivativeApproximation::internal::approximate
void approximate(SynchronousIterators< std::tuple< TriaActiveIterator< ::DoFCellAccessor< DoFHandlerType< dim, spacedim >, false >>, Vector< float >::iterator >> const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Definition: derivative_approximation.cc:924
Point< dim >
DoFRenumbering::Cuthill_McKee
void Cuthill_McKee(DoFHandlerType &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
Definition: dof_renumbering.cc:369
Functions::ZeroFunction
Definition: function.h:513
ParameterHandler
Definition: parameter_handler.h:845
Vector::l2_norm
real_type l2_norm() const
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
SolverControl
Definition: solver_control.h:67
Patterns::Selection
Definition: patterns.h:383
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
BlockVectorBase::block
BlockType & block(const unsigned int i)
first
Point< 2 > first
Definition: grid_out.cc:4352
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:935
DoFRenumbering::component_wise
void component_wise(DoFHandlerType &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
Definition: dof_renumbering.cc:633
FESystem
Definition: fe.h:44
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
ParameterHandler::parse_input
virtual void parse_input(std::istream &input, const std::string &filename="input file", const std::string &last_line="", const bool skip_undefined=false)
Definition: parameter_handler.cc:399
Patterns::Double
Definition: patterns.h:293
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
SparseMatrix::m
size_type m() const
Patterns::Integer
Definition: patterns.h:190
BlockSparseMatrix< double >