Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
Quasi_static_Finite_strain_Compressible_Elasticity.h
Go to the documentation of this file.
1
194 *  
195 *   /*
196 *   * Authors: Jean-Paul Pelteret, University of Erlangen-Nuremberg,
197 *   * Andrew McBride, University of Cape Town, 2015, 2017
198 *   */
199 *  
200 *  
201 * @endcode
202 *
203 * We start by including all the necessary deal.II header files and some C++
204 * related ones. They have been discussed in detail in previous tutorial
205 * programs, so you need only refer to past tutorials for details.
206 *
207 * @code
208 *   #include <deal.II/base/function.h>
209 *   #include <deal.II/base/parameter_handler.h>
210 *   #include <deal.II/base/point.h>
211 *   #include <deal.II/base/quadrature_lib.h>
212 *   #include <deal.II/base/symmetric_tensor.h>
213 *   #include <deal.II/base/tensor.h>
214 *   #include <deal.II/base/timer.h>
215 *   #include <deal.II/base/work_stream.h>
216 *   #include <deal.II/base/quadrature_point_data.h>
217 *  
218 *   #include <deal.II/dofs/dof_renumbering.h>
219 *   #include <deal.II/dofs/dof_tools.h>
220 *  
221 *   #include <deal.II/grid/grid_generator.h>
222 *   #include <deal.II/grid/grid_tools.h>
223 *   #include <deal.II/grid/grid_in.h>
224 *   #include <deal.II/grid/tria.h>
225 *  
226 *   #include <deal.II/fe/fe_dgp_monomial.h>
227 *   #include <deal.II/fe/fe_q.h>
228 *   #include <deal.II/fe/fe_system.h>
229 *   #include <deal.II/fe/fe_tools.h>
230 *   #include <deal.II/fe/fe_values.h>
231 *   #include <deal.II/fe/mapping_q_eulerian.h>
232 *   #include <deal.II/fe/mapping_q.h>
233 *  
234 *   #include <deal.II/lac/affine_constraints.h>
235 *   #include <deal.II/lac/block_sparse_matrix.h>
236 *   #include <deal.II/lac/block_vector.h>
237 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
238 *   #include <deal.II/lac/full_matrix.h>
239 *   #include <deal.II/lac/precondition_selector.h>
240 *   #include <deal.II/lac/solver_cg.h>
241 *   #include <deal.II/lac/sparse_direct.h>
242 *  
243 *   #include <deal.II/numerics/data_out.h>
244 *   #include <deal.II/numerics/vector_tools.h>
245 *  
246 *   #include <deal.II/base/config.h>
247 *   #if DEAL_II_VERSION_MAJOR >= 9 && defined(DEAL_II_WITH_TRILINOS)
248 *   #include <deal.II/differentiation/ad.h>
249 *   #define ENABLE_SACADO_FORMULATION
250 *   #endif
251 *  
252 * @endcode
253 *
254 * These must be included below the AD headers so that
255 * their math functions are available for use in the
256 * definition of tensors and kinematic quantities
257 *
258 * @code
259 *   #include <deal.II/physics/elasticity/kinematics.h>
260 *   #include <deal.II/physics/elasticity/standard_tensors.h>
261 *  
262 *   #include <iostream>
263 *   #include <fstream>
264 *   #include <memory>
265 *  
266 *  
267 * @endcode
268 *
269 * We then stick everything that relates to this tutorial program into a
270 * namespace of its own, and import all the deal.II function and class names
271 * into it:
272 *
273 * @code
274 *   namespace Cook_Membrane
275 *   {
276 *   using namespace dealii;
277 *  
278 * @endcode
279 *
280 *
281 * <a name="Runtimeparameters"></a>
282 * <h3>Run-time parameters</h3>
283 *
284
285 *
286 * There are several parameters that can be set in the code so we set up a
287 * ParameterHandler object to read in the choices at run-time.
288 *
289 * @code
290 *   namespace Parameters
291 *   {
292 * @endcode
293 *
294 *
295 * <a name="Assemblymethod"></a>
296 * <h4>Assembly method</h4>
297 *
298
299 *
300 * Here we specify whether automatic differentiation is to be used to assemble
301 * the linear system, and if so then what order of differentiation is to be
302 * employed.
303 *
304 * @code
305 *   struct AssemblyMethod
306 *   {
307 *   unsigned int automatic_differentiation_order;
308 *  
309 *   static void
310 *   declare_parameters(ParameterHandler &prm);
311 *  
312 *   void
313 *   parse_parameters(ParameterHandler &prm);
314 *   };
315 *  
316 *  
317 *   void AssemblyMethod::declare_parameters(ParameterHandler &prm)
318 *   {
319 *   prm.enter_subsection("Assembly method");
320 *   {
321 *   prm.declare_entry("Automatic differentiation order", "0",
322 *   Patterns::Integer(0,2),
323 *   "The automatic differentiation order to be used in the assembly of the linear system.\n"
324 *   "# Order = 0: Both the residual and linearisation are computed manually.\n"
325 *   "# Order = 1: The residual is computed manually but the linearisation is performed using AD.\n"
326 *   "# Order = 2: Both the residual and linearisation are computed using AD.");
327 *   }
328 *   prm.leave_subsection();
329 *   }
330 *  
331 *   void AssemblyMethod::parse_parameters(ParameterHandler &prm)
332 *   {
333 *   prm.enter_subsection("Assembly method");
334 *   {
335 *   automatic_differentiation_order = prm.get_integer("Automatic differentiation order");
336 *   }
337 *   prm.leave_subsection();
338 *   }
339 *  
340 * @endcode
341 *
342 *
343 * <a name="FiniteElementsystem"></a>
344 * <h4>Finite Element system</h4>
345 *
346
347 *
348 * Here we specify the polynomial order used to approximate the solution.
349 * The quadrature order should be adjusted accordingly.
350 *
351 * @code
352 *   struct FESystem
353 *   {
354 *   unsigned int poly_degree;
355 *   unsigned int quad_order;
356 *  
357 *   static void
358 *   declare_parameters(ParameterHandler &prm);
359 *  
360 *   void
361 *   parse_parameters(ParameterHandler &prm);
362 *   };
363 *  
364 *  
365 *   void FESystem::declare_parameters(ParameterHandler &prm)
366 *   {
367 *   prm.enter_subsection("Finite element system");
368 *   {
369 *   prm.declare_entry("Polynomial degree", "2",
370 *   Patterns::Integer(0),
371 *   "Displacement system polynomial order");
372 *  
373 *   prm.declare_entry("Quadrature order", "3",
374 *   Patterns::Integer(0),
375 *   "Gauss quadrature order");
376 *   }
377 *   prm.leave_subsection();
378 *   }
379 *  
380 *   void FESystem::parse_parameters(ParameterHandler &prm)
381 *   {
382 *   prm.enter_subsection("Finite element system");
383 *   {
384 *   poly_degree = prm.get_integer("Polynomial degree");
385 *   quad_order = prm.get_integer("Quadrature order");
386 *   }
387 *   prm.leave_subsection();
388 *   }
389 *  
390 * @endcode
391 *
392 *
393 * <a name="Geometry"></a>
394 * <h4>Geometry</h4>
395 *
396
397 *
398 * Make adjustments to the problem geometry and its discretisation.
399 *
400 * @code
401 *   struct Geometry
402 *   {
403 *   unsigned int elements_per_edge;
404 *   double scale;
405 *  
406 *   static void
407 *   declare_parameters(ParameterHandler &prm);
408 *  
409 *   void
410 *   parse_parameters(ParameterHandler &prm);
411 *   };
412 *  
413 *   void Geometry::declare_parameters(ParameterHandler &prm)
414 *   {
415 *   prm.enter_subsection("Geometry");
416 *   {
417 *   prm.declare_entry("Elements per edge", "32",
418 *   Patterns::Integer(0),
419 *   "Number of elements per long edge of the beam");
420 *  
421 *   prm.declare_entry("Grid scale", "1e-3",
422 *   Patterns::Double(0.0),
423 *   "Global grid scaling factor");
424 *   }
425 *   prm.leave_subsection();
426 *   }
427 *  
428 *   void Geometry::parse_parameters(ParameterHandler &prm)
429 *   {
430 *   prm.enter_subsection("Geometry");
431 *   {
432 *   elements_per_edge = prm.get_integer("Elements per edge");
433 *   scale = prm.get_double("Grid scale");
434 *   }
435 *   prm.leave_subsection();
436 *   }
437 *  
438 * @endcode
439 *
440 *
441 * <a name="Materials"></a>
442 * <h4>Materials</h4>
443 *
444
445 *
446 * We also need the shear modulus @f$ \mu @f$ and Poisson ration @f$ \nu @f$ for the
447 * neo-Hookean material.
448 *
449 * @code
450 *   struct Materials
451 *   {
452 *   double nu;
453 *   double mu;
454 *  
455 *   static void
456 *   declare_parameters(ParameterHandler &prm);
457 *  
458 *   void
459 *   parse_parameters(ParameterHandler &prm);
460 *   };
461 *  
462 *   void Materials::declare_parameters(ParameterHandler &prm)
463 *   {
464 *   prm.enter_subsection("Material properties");
465 *   {
466 *   prm.declare_entry("Poisson's ratio", "0.3",
467 *   Patterns::Double(-1.0,0.5),
468 *   "Poisson's ratio");
469 *  
470 *   prm.declare_entry("Shear modulus", "0.4225e6",
471 *   Patterns::Double(),
472 *   "Shear modulus");
473 *   }
474 *   prm.leave_subsection();
475 *   }
476 *  
477 *   void Materials::parse_parameters(ParameterHandler &prm)
478 *   {
479 *   prm.enter_subsection("Material properties");
480 *   {
481 *   nu = prm.get_double("Poisson's ratio");
482 *   mu = prm.get_double("Shear modulus");
483 *   }
484 *   prm.leave_subsection();
485 *   }
486 *  
487 * @endcode
488 *
489 *
490 * <a name="Linearsolver"></a>
491 * <h4>Linear solver</h4>
492 *
493
494 *
495 * Next, we choose both solver and preconditioner settings. The use of an
496 * effective preconditioner is critical to ensure convergence when a large
497 * nonlinear motion occurs within a Newton increment.
498 *
499 * @code
500 *   struct LinearSolver
501 *   {
502 *   std::string type_lin;
503 *   double tol_lin;
504 *   double max_iterations_lin;
505 *   std::string preconditioner_type;
506 *   double preconditioner_relaxation;
507 *  
508 *   static void
509 *   declare_parameters(ParameterHandler &prm);
510 *  
511 *   void
512 *   parse_parameters(ParameterHandler &prm);
513 *   };
514 *  
515 *   void LinearSolver::declare_parameters(ParameterHandler &prm)
516 *   {
517 *   prm.enter_subsection("Linear solver");
518 *   {
519 *   prm.declare_entry("Solver type", "CG",
520 *   Patterns::Selection("CG|Direct"),
521 *   "Type of solver used to solve the linear system");
522 *  
523 *   prm.declare_entry("Residual", "1e-6",
524 *   Patterns::Double(0.0),
525 *   "Linear solver residual (scaled by residual norm)");
526 *  
527 *   prm.declare_entry("Max iteration multiplier", "1",
528 *   Patterns::Double(0.0),
529 *   "Linear solver iterations (multiples of the system matrix size)");
530 *  
531 *   prm.declare_entry("Preconditioner type", "ssor",
532 *   Patterns::Selection("jacobi|ssor"),
533 *   "Type of preconditioner");
534 *  
535 *   prm.declare_entry("Preconditioner relaxation", "0.65",
536 *   Patterns::Double(0.0),
537 *   "Preconditioner relaxation value");
538 *   }
539 *   prm.leave_subsection();
540 *   }
541 *  
542 *   void LinearSolver::parse_parameters(ParameterHandler &prm)
543 *   {
544 *   prm.enter_subsection("Linear solver");
545 *   {
546 *   type_lin = prm.get("Solver type");
547 *   tol_lin = prm.get_double("Residual");
548 *   max_iterations_lin = prm.get_double("Max iteration multiplier");
549 *   preconditioner_type = prm.get("Preconditioner type");
550 *   preconditioner_relaxation = prm.get_double("Preconditioner relaxation");
551 *   }
552 *   prm.leave_subsection();
553 *   }
554 *  
555 * @endcode
556 *
557 *
558 * <a name="Nonlinearsolver"></a>
559 * <h4>Nonlinear solver</h4>
560 *
561
562 *
563 * A Newton-Raphson scheme is used to solve the nonlinear system of governing
564 * equations. We now define the tolerances and the maximum number of
565 * iterations for the Newton-Raphson nonlinear solver.
566 *
567 * @code
568 *   struct NonlinearSolver
569 *   {
570 *   unsigned int max_iterations_NR;
571 *   double tol_f;
572 *   double tol_u;
573 *  
574 *   static void
575 *   declare_parameters(ParameterHandler &prm);
576 *  
577 *   void
578 *   parse_parameters(ParameterHandler &prm);
579 *   };
580 *  
581 *   void NonlinearSolver::declare_parameters(ParameterHandler &prm)
582 *   {
583 *   prm.enter_subsection("Nonlinear solver");
584 *   {
585 *   prm.declare_entry("Max iterations Newton-Raphson", "10",
586 *   Patterns::Integer(0),
587 *   "Number of Newton-Raphson iterations allowed");
588 *  
589 *   prm.declare_entry("Tolerance force", "1.0e-9",
590 *   Patterns::Double(0.0),
591 *   "Force residual tolerance");
592 *  
593 *   prm.declare_entry("Tolerance displacement", "1.0e-6",
594 *   Patterns::Double(0.0),
595 *   "Displacement error tolerance");
596 *   }
597 *   prm.leave_subsection();
598 *   }
599 *  
600 *   void NonlinearSolver::parse_parameters(ParameterHandler &prm)
601 *   {
602 *   prm.enter_subsection("Nonlinear solver");
603 *   {
604 *   max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
605 *   tol_f = prm.get_double("Tolerance force");
606 *   tol_u = prm.get_double("Tolerance displacement");
607 *   }
608 *   prm.leave_subsection();
609 *   }
610 *  
611 * @endcode
612 *
613 *
614 * <a name="Time"></a>
615 * <h4>Time</h4>
616 *
617
618 *
619 * Set the timestep size @f$ \varDelta t @f$ and the simulation end-time.
620 *
621 * @code
622 *   struct Time
623 *   {
624 *   double delta_t;
625 *   double end_time;
626 *  
627 *   static void
628 *   declare_parameters(ParameterHandler &prm);
629 *  
630 *   void
631 *   parse_parameters(ParameterHandler &prm);
632 *   };
633 *  
634 *   void Time::declare_parameters(ParameterHandler &prm)
635 *   {
636 *   prm.enter_subsection("Time");
637 *   {
638 *   prm.declare_entry("End time", "1",
639 *   Patterns::Double(),
640 *   "End time");
641 *  
642 *   prm.declare_entry("Time step size", "0.1",
643 *   Patterns::Double(),
644 *   "Time step size");
645 *   }
646 *   prm.leave_subsection();
647 *   }
648 *  
649 *   void Time::parse_parameters(ParameterHandler &prm)
650 *   {
651 *   prm.enter_subsection("Time");
652 *   {
653 *   end_time = prm.get_double("End time");
654 *   delta_t = prm.get_double("Time step size");
655 *   }
656 *   prm.leave_subsection();
657 *   }
658 *  
659 * @endcode
660 *
661 *
662 * <a name="Allparameters"></a>
663 * <h4>All parameters</h4>
664 *
665
666 *
667 * Finally we consolidate all of the above structures into a single container
668 * that holds all of our run-time selections.
669 *
670 * @code
671 *   struct AllParameters :
672 *   public AssemblyMethod,
673 *   public FESystem,
674 *   public Geometry,
675 *   public Materials,
676 *   public LinearSolver,
677 *   public NonlinearSolver,
678 *   public Time
679 *  
680 *   {
681 *   AllParameters(const std::string &input_file);
682 *  
683 *   static void
684 *   declare_parameters(ParameterHandler &prm);
685 *  
686 *   void
687 *   parse_parameters(ParameterHandler &prm);
688 *   };
689 *  
690 *   AllParameters::AllParameters(const std::string &input_file)
691 *   {
692 *   ParameterHandler prm;
693 *   declare_parameters(prm);
694 *   prm.parse_input(input_file);
695 *   parse_parameters(prm);
696 *   }
697 *  
698 *   void AllParameters::declare_parameters(ParameterHandler &prm)
699 *   {
700 *   AssemblyMethod::declare_parameters(prm);
701 *   FESystem::declare_parameters(prm);
702 *   Geometry::declare_parameters(prm);
703 *   Materials::declare_parameters(prm);
704 *   LinearSolver::declare_parameters(prm);
705 *   NonlinearSolver::declare_parameters(prm);
706 *   Time::declare_parameters(prm);
707 *   }
708 *  
709 *   void AllParameters::parse_parameters(ParameterHandler &prm)
710 *   {
711 *   AssemblyMethod::parse_parameters(prm);
712 *   FESystem::parse_parameters(prm);
713 *   Geometry::parse_parameters(prm);
714 *   Materials::parse_parameters(prm);
715 *   LinearSolver::parse_parameters(prm);
716 *   NonlinearSolver::parse_parameters(prm);
717 *   Time::parse_parameters(prm);
718 *   }
719 *   }
720 *  
721 *  
722 * @endcode
723 *
724 *
725 * <a name="Timeclass"></a>
726 * <h3>Time class</h3>
727 *
728
729 *
730 * A simple class to store time data. Its functioning is transparent so no
731 * discussion is necessary. For simplicity we assume a constant time step
732 * size.
733 *
734 * @code
735 *   class Time
736 *   {
737 *   public:
738 *   Time (const double time_end,
739 *   const double delta_t)
740 *   :
741 *   timestep(0),
742 *   time_current(0.0),
743 *   time_end(time_end),
744 *   delta_t(delta_t)
745 *   {}
746 *  
747 *   virtual ~Time()
748 *   {}
749 *  
750 *   double current() const
751 *   {
752 *   return time_current;
753 *   }
754 *   double end() const
755 *   {
756 *   return time_end;
757 *   }
758 *   double get_delta_t() const
759 *   {
760 *   return delta_t;
761 *   }
762 *   unsigned int get_timestep() const
763 *   {
764 *   return timestep;
765 *   }
766 *   void increment()
767 *   {
768 *   time_current += delta_t;
769 *   ++timestep;
770 *   }
771 *  
772 *   private:
773 *   unsigned int timestep;
774 *   double time_current;
775 *   const double time_end;
776 *   const double delta_t;
777 *   };
778 *  
779 * @endcode
780 *
781 *
782 * <a name="CompressibleneoHookeanmaterialwithinaonefieldformulation"></a>
783 * <h3>Compressible neo-Hookean material within a one-field formulation</h3>
784 *
785
786 *
787 * As discussed in the literature and @ref step_44 "step-44", Neo-Hookean materials are a type
788 * of hyperelastic materials. The entire domain is assumed to be composed of a
789 * compressible neo-Hookean material. This class defines the behaviour of
790 * this material within a one-field formulation. Compressible neo-Hookean
791 * materials can be described by a strain-energy function (SEF) @f$ \Psi =
792 * \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(J)
793 * @f$.
794 *
795
796 *
797 * The isochoric response is given by @f$
798 * \Psi_{\text{iso}}(\overline{\mathbf{b}}) = c_{1} [\overline{I}_{1} - 3] @f$
799 * where @f$ c_{1} = \frac{\mu}{2} @f$ and @f$\overline{I}_{1}@f$ is the first
800 * invariant of the left- or right-isochoric Cauchy-Green deformation tensors.
801 * That is @f$\overline{I}_1 :=\textrm{tr}(\overline{\mathbf{b}})@f$. In this
802 * example the SEF that governs the volumetric response is defined as @f$
803 * \Psi_{\text{vol}}(J) = \kappa \frac{1}{4} [ J^2 - 1
804 * - 2\textrm{ln}\; J ]@f$, where @f$\kappa:= \lambda + 2/3 \mu@f$ is
805 * the <a href="http://en.wikipedia.org/wiki/Bulk_modulus">bulk modulus</a>
806 * and @f$\lambda@f$ is <a
807 * href="http://en.wikipedia.org/wiki/Lam%C3%A9_parameters">Lame's first
808 * parameter</a>.
809 *
810
811 *
812 * The following class will be used to characterize the material we work with,
813 * and provides a central point that one would need to modify if one were to
814 * implement a different material model. For it to work, we will store one
815 * object of this type per quadrature point, and in each of these objects
816 * store the current state (characterized by the values or measures of the
817 * displacement field) so that we can compute the elastic coefficients
818 * linearized around the current state.
819 *
820 * @code
821 *   template <int dim,typename NumberType>
822 *   class Material_Compressible_Neo_Hook_One_Field
823 *   {
824 *   public:
825 *   Material_Compressible_Neo_Hook_One_Field(const double mu,
826 *   const double nu)
827 *   :
828 *   kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
829 *   c_1(mu / 2.0)
830 *   {
831 *   Assert(kappa > 0, ExcInternalError());
832 *   }
833 *  
834 *   ~Material_Compressible_Neo_Hook_One_Field()
835 *   {}
836 *  
837 * @endcode
838 *
839 * The first function is the total energy
840 * @f$\Psi = \Psi_{\textrm{iso}} + \Psi_{\textrm{vol}}@f$.
841 *
842 * @code
843 *   NumberType
844 *   get_Psi(const NumberType &det_F,
845 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
846 *   {
847 *   return get_Psi_vol(det_F) + get_Psi_iso(b_bar);
848 *   }
849 *  
850 * @endcode
851 *
852 * The second function determines the Kirchhoff stress @f$\boldsymbol{\tau}
853 * = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}@f$
854 *
855 * @code
856 *   SymmetricTensor<2,dim,NumberType>
857 *   get_tau(const NumberType &det_F,
858 *   const SymmetricTensor<2,dim,NumberType> &b_bar)
859 *   {
860 * @endcode
861 *
862 * See Holzapfel p231 eq6.98 onwards
863 *
864 * @code
865 *   return get_tau_vol(det_F) + get_tau_iso(b_bar);
866 *   }
867 *  
868 * @endcode
869 *
870 * The fourth-order elasticity tensor in the spatial setting
871 * @f$\mathfrak{c}@f$ is calculated from the SEF @f$\Psi@f$ as @f$ J
872 * \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}@f$
873 * where @f$ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial
874 * \mathbf{C} \partial \mathbf{C}}@f$
875 *
876 * @code
877 *   SymmetricTensor<4,dim,NumberType>
878 *   get_Jc(const NumberType &det_F,
879 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
880 *   {
881 *   return get_Jc_vol(det_F) + get_Jc_iso(b_bar);
882 *   }
883 *  
884 *   private:
885 * @endcode
886 *
887 * Define constitutive model parameters @f$\kappa@f$ (bulk modulus) and the
888 * neo-Hookean model parameter @f$c_1@f$:
889 *
890 * @code
891 *   const double kappa;
892 *   const double c_1;
893 *  
894 * @endcode
895 *
896 * Value of the volumetric free energy
897 *
898 * @code
899 *   NumberType
900 *   get_Psi_vol(const NumberType &det_F) const
901 *   {
902 *   return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
903 *   }
904 *  
905 * @endcode
906 *
907 * Value of the isochoric free energy
908 *
909 * @code
910 *   NumberType
911 *   get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
912 *   {
913 *   return c_1 * (trace(b_bar) - dim);
914 *   }
915 *  
916 * @endcode
917 *
918 * Derivative of the volumetric free energy with respect to
919 * @f$J@f$ return @f$\frac{\partial
920 * \Psi_{\text{vol}}(J)}{\partial J}@f$
921 *
922 * @code
923 *   NumberType
924 *   get_dPsi_vol_dJ(const NumberType &det_F) const
925 *   {
926 *   return (kappa / 2.0) * (det_F - 1.0 / det_F);
927 *   }
928 *  
929 * @endcode
930 *
931 * The following functions are used internally in determining the result
932 * of some of the public functions above. The first one determines the
933 * volumetric Kirchhoff stress @f$\boldsymbol{\tau}_{\textrm{vol}}@f$.
934 * Note the difference in its definition when compared to @ref step_44 "step-44".
935 *
936 * @code
937 *   SymmetricTensor<2,dim,NumberType>
938 *   get_tau_vol(const NumberType &det_F) const
939 *   {
940 *   return NumberType(get_dPsi_vol_dJ(det_F) * det_F) * Physics::Elasticity::StandardTensors<dim>::I;
941 *   }
942 *  
943 * @endcode
944 *
945 * Next, determine the isochoric Kirchhoff stress
946 * @f$\boldsymbol{\tau}_{\textrm{iso}} =
947 * \mathcal{P}:\overline{\boldsymbol{\tau}}@f$:
948 *
949 * @code
950 *   SymmetricTensor<2,dim,NumberType>
951 *   get_tau_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
952 *   {
953 *   return Physics::Elasticity::StandardTensors<dim>::dev_P * get_tau_bar(b_bar);
954 *   }
955 *  
956 * @endcode
957 *
958 * Then, determine the fictitious Kirchhoff stress
959 * @f$\overline{\boldsymbol{\tau}}@f$:
960 *
961 * @code
962 *   SymmetricTensor<2,dim,NumberType>
963 *   get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const
964 *   {
965 *   return 2.0 * c_1 * b_bar;
966 *   }
967 *  
968 * @endcode
969 *
970 * Second derivative of the volumetric free energy wrt @f$J@f$. We
971 * need the following computation explicitly in the tangent so we make it
972 * public. We calculate @f$\frac{\partial^2
973 * \Psi_{\textrm{vol}}(J)}{\partial J \partial
974 * J}@f$
975 *
976 * @code
977 *   NumberType
978 *   get_d2Psi_vol_dJ2(const NumberType &det_F) const
979 *   {
980 *   return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
981 *   }
982 *  
983 * @endcode
984 *
985 * Calculate the volumetric part of the tangent @f$J
986 * \mathfrak{c}_\textrm{vol}@f$. Again, note the difference in its
987 * definition when compared to @ref step_44 "step-44". The extra terms result from two
988 * quantities in @f$\boldsymbol{\tau}_{\textrm{vol}}@f$ being dependent on
989 * @f$\boldsymbol{F}@f$.
990 *
991 * @code
992 *   SymmetricTensor<4,dim,NumberType>
993 *   get_Jc_vol(const NumberType &det_F) const
994 *   {
995 * @endcode
996 *
997 * See Holzapfel p265
998 *
999 * @code
1000 *   return det_F
1001 *   * ( (get_dPsi_vol_dJ(det_F) + det_F * get_d2Psi_vol_dJ2(det_F))*Physics::Elasticity::StandardTensors<dim>::IxI
1002 *   - (2.0 * get_dPsi_vol_dJ(det_F))*Physics::Elasticity::StandardTensors<dim>::S );
1003 *   }
1004 *  
1005 * @endcode
1006 *
1007 * Calculate the isochoric part of the tangent @f$J
1008 * \mathfrak{c}_\textrm{iso}@f$:
1009 *
1010 * @code
1011 *   SymmetricTensor<4,dim,NumberType>
1012 *   get_Jc_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
1013 *   {
1014 *   const SymmetricTensor<2, dim> tau_bar = get_tau_bar(b_bar);
1015 *   const SymmetricTensor<2, dim> tau_iso = get_tau_iso(b_bar);
1016 *   const SymmetricTensor<4, dim> tau_iso_x_I
1017 *   = outer_product(tau_iso,
1018 *   Physics::Elasticity::StandardTensors<dim>::I);
1019 *   const SymmetricTensor<4, dim> I_x_tau_iso
1020 *   = outer_product(Physics::Elasticity::StandardTensors<dim>::I,
1021 *   tau_iso);
1022 *   const SymmetricTensor<4, dim> c_bar = get_c_bar();
1023 *  
1024 *   return (2.0 / dim) * trace(tau_bar)
1025 *   * Physics::Elasticity::StandardTensors<dim>::dev_P
1026 *   - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso)
1027 *   + Physics::Elasticity::StandardTensors<dim>::dev_P * c_bar
1028 *   * Physics::Elasticity::StandardTensors<dim>::dev_P;
1029 *   }
1030 *  
1031 * @endcode
1032 *
1033 * Calculate the fictitious elasticity tensor @f$\overline{\mathfrak{c}}@f$.
1034 * For the material model chosen this is simply zero:
1035 *
1036 * @code
1037 *   SymmetricTensor<4,dim,double>
1038 *   get_c_bar() const
1039 *   {
1040 *   return SymmetricTensor<4, dim>();
1041 *   }
1042 *   };
1043 *  
1044 * @endcode
1045 *
1046 *
1047 * <a name="Quadraturepointhistory"></a>
1048 * <h3>Quadrature point history</h3>
1049 *
1050
1051 *
1052 * As seen in @ref step_18 "step-18", the <code> PointHistory </code> class offers a method
1053 * for storing data at the quadrature points. Here each quadrature point
1054 * holds a pointer to a material description. Thus, different material models
1055 * can be used in different regions of the domain. Among other data, we
1056 * choose to store the Kirchhoff stress @f$\boldsymbol{\tau}@f$ and the tangent
1057 * @f$J\mathfrak{c}@f$ for the quadrature points.
1058 *
1059 * @code
1060 *   template <int dim,typename NumberType>
1061 *   class PointHistory
1062 *   {
1063 *   public:
1064 *   PointHistory()
1065 *   {}
1066 *  
1067 *   virtual ~PointHistory()
1068 *   {}
1069 *  
1070 * @endcode
1071 *
1072 * The first function is used to create a material object and to
1073 * initialize all tensors correctly: The second one updates the stored
1074 * values and stresses based on the current deformation measure
1075 * @f$\textrm{Grad}\mathbf{u}_{\textrm{n}}@f$.
1076 *
1077 * @code
1078 *   void setup_lqp (const Parameters::AllParameters &parameters)
1079 *   {
1080 *   material.reset(new Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu,
1081 *   parameters.nu));
1082 *   }
1083 *  
1084 * @endcode
1085 *
1086 * We offer an interface to retrieve certain data.
1087 * This is the strain energy:
1088 *
1089 * @code
1090 *   NumberType
1091 *   get_Psi(const NumberType &det_F,
1092 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
1093 *   {
1094 *   return material->get_Psi(det_F,b_bar);
1095 *   }
1096 *  
1097 * @endcode
1098 *
1099 * Here are the kinetic variables. These are used in the material and
1100 * global tangent matrix and residual assembly operations:
1101 * First is the Kirchhoff stress:
1102 *
1103 * @code
1104 *   SymmetricTensor<2,dim,NumberType>
1105 *   get_tau(const NumberType &det_F,
1106 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
1107 *   {
1108 *   return material->get_tau(det_F,b_bar);
1109 *   }
1110 *  
1111 * @endcode
1112 *
1113 * And the tangent:
1114 *
1115 * @code
1116 *   SymmetricTensor<4,dim,NumberType>
1117 *   get_Jc(const NumberType &det_F,
1118 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
1119 *   {
1120 *   return material->get_Jc(det_F,b_bar);
1121 *   }
1122 *  
1123 * @endcode
1124 *
1125 * In terms of member functions, this class stores for the quadrature
1126 * point it represents a copy of a material type in case different
1127 * materials are used in different regions of the domain, as well as the
1128 * inverse of the deformation gradient...
1129 *
1130 * @code
1131 *   private:
1132 *   std::shared_ptr< Material_Compressible_Neo_Hook_One_Field<dim,NumberType> > material;
1133 *   };
1134 *  
1135 *  
1136 * @endcode
1137 *
1138 *
1139 * <a name="Quasistaticcompressiblefinitestrainsolid"></a>
1140 * <h3>Quasi-static compressible finite-strain solid</h3>
1141 *
1142
1143 *
1144 * Forward declarations for classes that will
1145 * perform assembly of the linear system.
1146 *
1147 * @code
1148 *   template <int dim,typename NumberType>
1149 *   struct Assembler_Base;
1150 *   template <int dim,typename NumberType>
1151 *   struct Assembler;
1152 *  
1153 * @endcode
1154 *
1155 * The Solid class is the central class in that it represents the problem at
1156 * hand. It follows the usual scheme in that all it really has is a
1157 * constructor, destructor and a <code>run()</code> function that dispatches
1158 * all the work to private functions of this class:
1159 *
1160 * @code
1161 *   template <int dim,typename NumberType>
1162 *   class Solid
1163 *   {
1164 *   public:
1165 *   Solid(const Parameters::AllParameters &parameters);
1166 *  
1167 *   virtual
1168 *   ~Solid();
1169 *  
1170 *   void
1171 *   run();
1172 *  
1173 *   private:
1174 *  
1175 * @endcode
1176 *
1177 * We start the collection of member functions with one that builds the
1178 * grid:
1179 *
1180 * @code
1181 *   void
1182 *   make_grid();
1183 *  
1184 * @endcode
1185 *
1186 * Set up the finite element system to be solved:
1187 *
1188 * @code
1189 *   void
1190 *   system_setup();
1191 *  
1192 * @endcode
1193 *
1194 * Several functions to assemble the system and right hand side matrices
1195 * using multithreading. Each of them comes as a wrapper function, one
1196 * that is executed to do the work in the WorkStream model on one cell,
1197 * and one that copies the work done on this one cell into the global
1198 * object that represents it:
1199 *
1200 * @code
1201 *   void
1202 *   assemble_system(const BlockVector<double> &solution_delta);
1203 *  
1204 * @endcode
1205 *
1206 * We use a separate data structure to perform the assembly. It needs access
1207 * to some low-level data, so we simply befriend the class instead of
1208 * creating a complex interface to provide access as necessary.
1209 *
1210 * @code
1211 *   friend struct Assembler_Base<dim,NumberType>;
1212 *   friend struct Assembler<dim,NumberType>;
1213 *  
1214 * @endcode
1215 *
1216 * Apply Dirichlet boundary conditions on the displacement field
1217 *
1218 * @code
1219 *   void
1220 *   make_constraints(const int &it_nr);
1221 *  
1222 * @endcode
1223 *
1224 * Create and update the quadrature points. Here, no data needs to be
1225 * copied into a global object, so the copy_local_to_global function is
1226 * empty:
1227 *
1228 * @code
1229 *   void
1230 *   setup_qph();
1231 *  
1232 * @endcode
1233 *
1234 * Solve for the displacement using a Newton-Raphson method. We break this
1235 * function into the nonlinear loop and the function that solves the
1236 * linearized Newton-Raphson step:
1237 *
1238 * @code
1239 *   void
1240 *   solve_nonlinear_timestep(BlockVector<double> &solution_delta);
1241 *  
1242 *   std::pair<unsigned int, double>
1243 *   solve_linear_system(BlockVector<double> &newton_update);
1244 *  
1245 * @endcode
1246 *
1247 * Solution retrieval as well as post-processing and writing data to file :
1248 *
1249 * @code
1250 *   BlockVector<double>
1251 *   get_total_solution(const BlockVector<double> &solution_delta) const;
1252 *  
1253 *   void
1254 *   output_results() const;
1255 *  
1256 * @endcode
1257 *
1258 * Finally, some member variables that describe the current state: A
1259 * collection of the parameters used to describe the problem setup...
1260 *
1261 * @code
1262 *   const Parameters::AllParameters &parameters;
1263 *  
1264 * @endcode
1265 *
1266 * ...the volume of the reference and current configurations...
1267 *
1268 * @code
1269 *   double vol_reference;
1270 *   double vol_current;
1271 *  
1272 * @endcode
1273 *
1274 * ...and description of the geometry on which the problem is solved:
1275 *
1276 * @code
1277 *   Triangulation<dim> triangulation;
1278 *  
1279 * @endcode
1280 *
1281 * Also, keep track of the current time and the time spent evaluating
1282 * certain functions
1283 *
1284 * @code
1285 *   Time time;
1286 *   TimerOutput timer;
1287 *  
1288 * @endcode
1289 *
1290 * A storage object for quadrature point information. As opposed to
1291 * @ref step_18 "step-18", deal.II's native quadrature point data manager is employed here.
1292 *
1293 * @code
1295 *   PointHistory<dim,NumberType> > quadrature_point_history;
1296 *  
1297 * @endcode
1298 *
1299 * A description of the finite-element system including the displacement
1300 * polynomial degree, the degree-of-freedom handler, number of DoFs per
1301 * cell and the extractor objects used to retrieve information from the
1302 * solution vectors:
1303 *
1304 * @code
1305 *   const unsigned int degree;
1306 *   const FESystem<dim> fe;
1307 *   DoFHandler<dim> dof_handler_ref;
1308 *   const unsigned int dofs_per_cell;
1309 *   const FEValuesExtractors::Vector u_fe;
1310 *  
1311 * @endcode
1312 *
1313 * Description of how the block-system is arranged. There is just 1 block,
1314 * that contains a vector DOF @f$\mathbf{u}@f$.
1315 * There are two reasons that we retain the block system in this problem.
1316 * The first is pure laziness to perform further modifications to the
1317 * code from which this work originated. The second is that a block system
1318 * would typically necessary when extending this code to multiphysics
1319 * problems.
1320 *
1321 * @code
1322 *   static const unsigned int n_blocks = 1;
1323 *   static const unsigned int n_components = dim;
1324 *   static const unsigned int first_u_component = 0;
1325 *  
1326 *   enum
1327 *   {
1328 *   u_dof = 0
1329 *   };
1330 *  
1331 *   std::vector<types::global_dof_index> dofs_per_block;
1332 *  
1333 * @endcode
1334 *
1335 * Rules for Gauss-quadrature on both the cell and faces. The number of
1336 * quadrature points on both cells and faces is recorded.
1337 *
1338 * @code
1339 *   const QGauss<dim> qf_cell;
1340 *   const QGauss<dim - 1> qf_face;
1341 *   const unsigned int n_q_points;
1342 *   const unsigned int n_q_points_f;
1343 *  
1344 * @endcode
1345 *
1346 * Objects that store the converged solution and right-hand side vectors,
1347 * as well as the tangent matrix. There is a AffineConstraints object used
1348 * to keep track of constraints. We make use of a sparsity pattern
1349 * designed for a block system.
1350 *
1351 * @code
1352 *   AffineConstraints<double> constraints;
1353 *   BlockSparsityPattern sparsity_pattern;
1354 *   BlockSparseMatrix<double> tangent_matrix;
1355 *   BlockVector<double> system_rhs;
1356 *   BlockVector<double> solution_n;
1357 *  
1358 * @endcode
1359 *
1360 * Then define a number of variables to store norms and update norms and
1361 * normalisation factors.
1362 *
1363 * @code
1364 *   struct Errors
1365 *   {
1366 *   Errors()
1367 *   :
1368 *   norm(1.0), u(1.0)
1369 *   {}
1370 *  
1371 *   void reset()
1372 *   {
1373 *   norm = 1.0;
1374 *   u = 1.0;
1375 *   }
1376 *   void normalise(const Errors &rhs)
1377 *   {
1378 *   if (rhs.norm != 0.0)
1379 *   norm /= rhs.norm;
1380 *   if (rhs.u != 0.0)
1381 *   u /= rhs.u;
1382 *   }
1383 *  
1384 *   double norm, u;
1385 *   };
1386 *  
1387 *   Errors error_residual, error_residual_0, error_residual_norm, error_update,
1388 *   error_update_0, error_update_norm;
1389 *  
1390 * @endcode
1391 *
1392 * Methods to calculate error measures
1393 *
1394 * @code
1395 *   void
1396 *   get_error_residual(Errors &error_residual);
1397 *  
1398 *   void
1399 *   get_error_update(const BlockVector<double> &newton_update,
1400 *   Errors &error_update);
1401 *  
1402 * @endcode
1403 *
1404 * Print information to screen in a pleasing way...
1405 *
1406 * @code
1407 *   static
1408 *   void
1409 *   print_conv_header();
1410 *  
1411 *   void
1412 *   print_conv_footer();
1413 *  
1414 *   void
1415 *   print_vertical_tip_displacement();
1416 *   };
1417 *  
1418 * @endcode
1419 *
1420 *
1421 * <a name="ImplementationofthecodeSolidcodeclass"></a>
1422 * <h3>Implementation of the <code>Solid</code> class</h3>
1423 *
1424
1425 *
1426 *
1427 * <a name="Publicinterface"></a>
1428 * <h4>Public interface</h4>
1429 *
1430
1431 *
1432 * We initialise the Solid class using data extracted from the parameter file.
1433 *
1434 * @code
1435 *   template <int dim,typename NumberType>
1436 *   Solid<dim,NumberType>::Solid(const Parameters::AllParameters &parameters)
1437 *   :
1438 *   parameters(parameters),
1439 *   vol_reference (0.0),
1440 *   vol_current (0.0),
1442 *   time(parameters.end_time, parameters.delta_t),
1443 *   timer(std::cout,
1446 *   degree(parameters.poly_degree),
1447 * @endcode
1448 *
1449 * The Finite Element System is composed of dim continuous displacement
1450 * DOFs.
1451 *
1452 * @code
1453 *   fe(FE_Q<dim>(parameters.poly_degree), dim), // displacement
1454 *   dof_handler_ref(triangulation),
1455 *   dofs_per_cell (fe.dofs_per_cell),
1456 *   u_fe(first_u_component),
1457 *   dofs_per_block(n_blocks),
1458 *   qf_cell(parameters.quad_order),
1459 *   qf_face(parameters.quad_order),
1460 *   n_q_points (qf_cell.size()),
1461 *   n_q_points_f (qf_face.size())
1462 *   {
1463 *  
1464 *   }
1465 *  
1466 * @endcode
1467 *
1468 * The class destructor simply clears the data held by the DOFHandler
1469 *
1470 * @code
1471 *   template <int dim,typename NumberType>
1472 *   Solid<dim,NumberType>::~Solid()
1473 *   {
1474 *   dof_handler_ref.clear();
1475 *   }
1476 *  
1477 *  
1478 * @endcode
1479 *
1480 * In solving the quasi-static problem, the time becomes a loading parameter,
1481 * i.e. we increasing the loading linearly with time, making the two concepts
1482 * interchangeable. We choose to increment time linearly using a constant time
1483 * step size.
1484 *
1485
1486 *
1487 * We start the function with preprocessing, and then output the initial grid
1488 * before starting the simulation proper with the first time (and loading)
1489 * increment.
1490 *
1491
1492 *
1493 *
1494 * @code
1495 *   template <int dim,typename NumberType>
1496 *   void Solid<dim,NumberType>::run()
1497 *   {
1498 *   make_grid();
1499 *   system_setup();
1500 *   output_results();
1501 *   time.increment();
1502 *  
1503 * @endcode
1504 *
1505 * We then declare the incremental solution update @f$\varDelta
1506 * \mathbf{\Xi}:= \{\varDelta \mathbf{u}\}@f$ and start the loop over the
1507 * time domain.
1508 *
1509
1510 *
1511 * At the beginning, we reset the solution update for this time step...
1512 *
1513 * @code
1514 *   BlockVector<double> solution_delta(dofs_per_block);
1515 *   while (time.current() <= time.end())
1516 *   {
1517 *   solution_delta = 0.0;
1518 *  
1519 * @endcode
1520 *
1521 * ...solve the current time step and update total solution vector
1522 * @f$\mathbf{\Xi}_{\textrm{n}} = \mathbf{\Xi}_{\textrm{n-1}} +
1523 * \varDelta \mathbf{\Xi}@f$...
1524 *
1525 * @code
1526 *   solve_nonlinear_timestep(solution_delta);
1527 *   solution_n += solution_delta;
1528 *  
1529 * @endcode
1530 *
1531 * ...and plot the results before moving on happily to the next time
1532 * step:
1533 *
1534 * @code
1535 *   output_results();
1536 *   time.increment();
1537 *   }
1538 *  
1539 * @endcode
1540 *
1541 * Lastly, we print the vertical tip displacement of the Cook cantilever
1542 * after the full load is applied
1543 *
1544 * @code
1545 *   print_vertical_tip_displacement();
1546 *   }
1547 *  
1548 *  
1549 * @endcode
1550 *
1551 *
1552 * <a name="Privateinterface"></a>
1553 * <h3>Private interface</h3>
1554 *
1555
1556 *
1557 *
1558 * <a name="Solidmake_grid"></a>
1559 * <h4>Solid::make_grid</h4>
1560 *
1561
1562 *
1563 * On to the first of the private member functions. Here we create the
1564 * triangulation of the domain, for which we choose a scaled an anisotripically
1565 * discretised rectangle which is subsequently transformed into the correct
1566 * of the Cook cantilever. Each relevant boundary face is then given a boundary
1567 * ID number.
1568 *
1569
1570 *
1571 * We then determine the volume of the reference configuration and print it
1572 * for comparison.
1573 *
1574
1575 *
1576 *
1577 * @code
1578 *   template <int dim>
1579 *   Point<dim> grid_y_transform (const Point<dim> &pt_in)
1580 *   {
1581 *   const double &x = pt_in[0];
1582 *   const double &y = pt_in[1];
1583 *  
1584 *   const double y_upper = 44.0 + (16.0/48.0)*x; // Line defining upper edge of beam
1585 *   const double y_lower = 0.0 + (44.0/48.0)*x; // Line defining lower edge of beam
1586 *   const double theta = y/44.0; // Fraction of height along left side of beam
1587 *   const double y_transform = (1-theta)*y_lower + theta*y_upper; // Final transformation
1588 *  
1589 *   Point<dim> pt_out = pt_in;
1590 *   pt_out[1] = y_transform;
1591 *  
1592 *   return pt_out;
1593 *   }
1594 *  
1595 *   template <int dim,typename NumberType>
1596 *   void Solid<dim,NumberType>::make_grid()
1597 *   {
1598 * @endcode
1599 *
1600 * Divide the beam, but only along the x- and y-coordinate directions
1601 *
1602 * @code
1603 *   std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge);
1604 * @endcode
1605 *
1606 * Only allow one element through the thickness
1607 * (modelling a plane strain condition)
1608 *
1609 * @code
1610 *   if (dim == 3)
1611 *   repetitions[dim-1] = 1;
1612 *  
1613 *   const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, -0.5) : Point<dim>(0.0, 0.0));
1614 *   const Point<dim> top_right = (dim == 3 ? Point<dim>(48.0, 44.0, 0.5) : Point<dim>(48.0, 44.0));
1615 *  
1617 *   repetitions,
1618 *   bottom_left,
1619 *   top_right);
1620 *  
1621 * @endcode
1622 *
1623 * Since we wish to apply a Neumann BC to the right-hand surface, we
1624 * must find the cell faces in this part of the domain and mark them with
1625 * a distinct boundary ID number. The faces we are looking for are on the
1626 * +x surface and will get boundary ID 11.
1627 * Dirichlet boundaries exist on the left-hand face of the beam (this fixed
1628 * boundary will get ID 1) and on the +Z and -Z faces (which correspond to
1629 * ID 2 and we will use to impose the plane strain condition)
1630 *
1631 * @code
1632 *   const double tol_boundary = 1e-6;
1633 *   typename Triangulation<dim>::active_cell_iterator cell =
1634 *   triangulation.begin_active(), endc = triangulation.end();
1635 *   for (; cell != endc; ++cell)
1636 *   for (unsigned int face = 0;
1637 *   face < GeometryInfo<dim>::faces_per_cell; ++face)
1638 *   if (cell->face(face)->at_boundary() == true)
1639 *   {
1640 *   if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary)
1641 *   cell->face(face)->set_boundary_id(1); // -X faces
1642 *   else if (std::abs(cell->face(face)->center()[0] - 48.0) < tol_boundary)
1643 *   cell->face(face)->set_boundary_id(11); // +X faces
1644 *   else if (dim == 3 && std::abs(std::abs(cell->face(face)->center()[2]) - 0.5) < tol_boundary)
1645 *   cell->face(face)->set_boundary_id(2); // +Z and -Z faces
1646 *   }
1647 *  
1648 * @endcode
1649 *
1650 * Transform the hyper-rectangle into the beam shape
1651 *
1652 * @code
1653 *   GridTools::transform(&grid_y_transform<dim>, triangulation);
1654 *  
1655 *   GridTools::scale(parameters.scale, triangulation);
1656 *  
1657 *   vol_reference = GridTools::volume(triangulation);
1658 *   vol_current = vol_reference;
1659 *   std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
1660 *   }
1661 *  
1662 *  
1663 * @endcode
1664 *
1665 *
1666 * <a name="Solidsystem_setup"></a>
1667 * <h4>Solid::system_setup</h4>
1668 *
1669
1670 *
1671 * Next we describe how the FE system is setup. We first determine the number
1672 * of components per block. Since the displacement is a vector component, the
1673 * first dim components belong to it.
1674 *
1675 * @code
1676 *   template <int dim,typename NumberType>
1677 *   void Solid<dim,NumberType>::system_setup()
1678 *   {
1679 *   timer.enter_subsection("Setup system");
1680 *  
1681 *   std::vector<unsigned int> block_component(n_components, u_dof); // Displacement
1682 *  
1683 * @endcode
1684 *
1685 * The DOF handler is then initialised and we renumber the grid in an
1686 * efficient manner. We also record the number of DOFs per block.
1687 *
1688 * @code
1689 *   dof_handler_ref.distribute_dofs(fe);
1690 *   DoFRenumbering::Cuthill_McKee(dof_handler_ref);
1691 *   DoFRenumbering::component_wise(dof_handler_ref, block_component);
1692 *   dofs_per_block = DoFTools::count_dofs_per_fe_block(dof_handler_ref, block_component);
1693 *  
1694 *   std::cout << "Triangulation:"
1695 *   << "\n\t Number of active cells: " << triangulation.n_active_cells()
1696 *   << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs()
1697 *   << std::endl;
1698 *  
1699 * @endcode
1700 *
1701 * Setup the sparsity pattern and tangent matrix
1702 *
1703 * @code
1704 *   tangent_matrix.clear();
1705 *   {
1706 *   const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
1707 *  
1708 *   BlockDynamicSparsityPattern csp(n_blocks, n_blocks);
1709 *  
1710 *   csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
1711 *   csp.collect_sizes();
1712 *  
1713 * @endcode
1714 *
1715 * Naturally, for a one-field vector-valued problem, all of the
1716 * components of the system are coupled.
1717 *
1718 * @code
1719 *   Table<2, DoFTools::Coupling> coupling(n_components, n_components);
1720 *   for (unsigned int ii = 0; ii < n_components; ++ii)
1721 *   for (unsigned int jj = 0; jj < n_components; ++jj)
1722 *   coupling[ii][jj] = DoFTools::always;
1723 *   DoFTools::make_sparsity_pattern(dof_handler_ref,
1724 *   coupling,
1725 *   csp,
1726 *   constraints,
1727 *   false);
1728 *   sparsity_pattern.copy_from(csp);
1729 *   }
1730 *  
1731 *   tangent_matrix.reinit(sparsity_pattern);
1732 *  
1733 * @endcode
1734 *
1735 * We then set up storage vectors
1736 *
1737 * @code
1738 *   system_rhs.reinit(dofs_per_block);
1739 *   system_rhs.collect_sizes();
1740 *  
1741 *   solution_n.reinit(dofs_per_block);
1742 *   solution_n.collect_sizes();
1743 *  
1744 * @endcode
1745 *
1746 * ...and finally set up the quadrature
1747 * point history:
1748 *
1749 * @code
1750 *   setup_qph();
1751 *  
1752 *   timer.leave_subsection();
1753 *   }
1754 *  
1755 *  
1756 * @endcode
1757 *
1758 *
1759 * <a name="Solidsetup_qph"></a>
1760 * <h4>Solid::setup_qph</h4>
1761 * The method used to store quadrature information is already described in
1762 * @ref step_18 "step-18" and @ref step_44 "step-44". Here we implement a similar setup for a SMP machine.
1763 *
1764
1765 *
1766 * Firstly the actual QPH data objects are created. This must be done only
1767 * once the grid is refined to its finest level.
1768 *
1769 * @code
1770 *   template <int dim,typename NumberType>
1771 *   void Solid<dim,NumberType>::setup_qph()
1772 *   {
1773 *   std::cout << " Setting up quadrature point data..." << std::endl;
1774 *  
1775 *   quadrature_point_history.initialize(triangulation.begin_active(),
1776 *   triangulation.end(),
1777 *   n_q_points);
1778 *  
1779 * @endcode
1780 *
1781 * Next we setup the initial quadrature point data. Note that when
1782 * the quadrature point data is retrieved, it is returned as a vector
1783 * of smart pointers.
1784 *
1785 * @code
1786 *   for (typename Triangulation<dim>::active_cell_iterator cell =
1787 *   triangulation.begin_active(); cell != triangulation.end(); ++cell)
1788 *   {
1789 *   const std::vector<std::shared_ptr<PointHistory<dim,NumberType> > > lqph =
1790 *   quadrature_point_history.get_data(cell);
1791 *   Assert(lqph.size() == n_q_points, ExcInternalError());
1792 *  
1793 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1794 *   lqph[q_point]->setup_lqp(parameters);
1795 *   }
1796 *   }
1797 *  
1798 *  
1799 * @endcode
1800 *
1801 *
1802 * <a name="Solidsolve_nonlinear_timestep"></a>
1803 * <h4>Solid::solve_nonlinear_timestep</h4>
1804 *
1805
1806 *
1807 * The next function is the driver method for the Newton-Raphson scheme. At
1808 * its top we create a new vector to store the current Newton update step,
1809 * reset the error storage objects and print solver header.
1810 *
1811 * @code
1812 *   template <int dim,typename NumberType>
1813 *   void
1814 *   Solid<dim,NumberType>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
1815 *   {
1816 *   std::cout << std::endl << "Timestep " << time.get_timestep() << " @ "
1817 *   << time.current() << "s" << std::endl;
1818 *  
1819 *   BlockVector<double> newton_update(dofs_per_block);
1820 *  
1821 *   error_residual.reset();
1822 *   error_residual_0.reset();
1823 *   error_residual_norm.reset();
1824 *   error_update.reset();
1825 *   error_update_0.reset();
1826 *   error_update_norm.reset();
1827 *  
1828 *   print_conv_header();
1829 *  
1830 * @endcode
1831 *
1832 * We now perform a number of Newton iterations to iteratively solve the
1833 * nonlinear problem. Since the problem is fully nonlinear and we are
1834 * using a full Newton method, the data stored in the tangent matrix and
1835 * right-hand side vector is not reusable and must be cleared at each
1836 * Newton step. We then initially build the right-hand side vector to
1837 * check for convergence (and store this value in the first iteration).
1838 * The unconstrained DOFs of the rhs vector hold the out-of-balance
1839 * forces. The building is done before assembling the system matrix as the
1840 * latter is an expensive operation and we can potentially avoid an extra
1841 * assembly process by not assembling the tangent matrix when convergence
1842 * is attained.
1843 *
1844 * @code
1845 *   unsigned int newton_iteration = 0;
1846 *   for (; newton_iteration < parameters.max_iterations_NR;
1847 *   ++newton_iteration)
1848 *   {
1849 *   std::cout << " " << std::setw(2) << newton_iteration << " " << std::flush;
1850 *  
1851 * @endcode
1852 *
1853 * If we have decided that we want to continue with the iteration, we
1854 * assemble the tangent, make and impose the Dirichlet constraints,
1855 * and do the solve of the linearized system:
1856 *
1857 * @code
1858 *   make_constraints(newton_iteration);
1859 *   assemble_system(solution_delta);
1860 *  
1861 *   get_error_residual(error_residual);
1862 *  
1863 *   if (newton_iteration == 0)
1864 *   error_residual_0 = error_residual;
1865 *  
1866 * @endcode
1867 *
1868 * We can now determine the normalised residual error and check for
1869 * solution convergence:
1870 *
1871 * @code
1872 *   error_residual_norm = error_residual;
1873 *   error_residual_norm.normalise(error_residual_0);
1874 *  
1875 *   if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u
1876 *   && error_residual_norm.u <= parameters.tol_f)
1877 *   {
1878 *   std::cout << " CONVERGED! " << std::endl;
1879 *   print_conv_footer();
1880 *  
1881 *   break;
1882 *   }
1883 *  
1884 *   const std::pair<unsigned int, double>
1885 *   lin_solver_output = solve_linear_system(newton_update);
1886 *  
1887 *   get_error_update(newton_update, error_update);
1888 *   if (newton_iteration == 0)
1889 *   error_update_0 = error_update;
1890 *  
1891 * @endcode
1892 *
1893 * We can now determine the normalised Newton update error, and
1894 * perform the actual update of the solution increment for the current
1895 * time step, update all quadrature point information pertaining to
1896 * this new displacement and stress state and continue iterating:
1897 *
1898 * @code
1899 *   error_update_norm = error_update;
1900 *   error_update_norm.normalise(error_update_0);
1901 *  
1902 *   solution_delta += newton_update;
1903 *  
1904 *   std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
1905 *   << std::scientific << lin_solver_output.first << " "
1906 *   << lin_solver_output.second << " " << error_residual_norm.norm
1907 *   << " " << error_residual_norm.u << " "
1908 *   << " " << error_update_norm.norm << " " << error_update_norm.u
1909 *   << " " << std::endl;
1910 *   }
1911 *  
1912 * @endcode
1913 *
1914 * At the end, if it turns out that we have in fact done more iterations
1915 * than the parameter file allowed, we raise an exception that can be
1916 * caught in the main() function. The call <code>AssertThrow(condition,
1917 * exc_object)</code> is in essence equivalent to <code>if (!cond) throw
1918 * exc_object;</code> but the former form fills certain fields in the
1919 * exception object that identify the location (filename and line number)
1920 * where the exception was raised to make it simpler to identify where the
1921 * problem happened.
1922 *
1923 * @code
1924 *   AssertThrow (newton_iteration <= parameters.max_iterations_NR,
1925 *   ExcMessage("No convergence in nonlinear solver!"));
1926 *   }
1927 *  
1928 *  
1929 * @endcode
1930 *
1931 *
1932 * <a name="Solidprint_conv_headerSolidprint_conv_footerandSolidprint_vertical_tip_displacement"></a>
1933 * <h4>Solid::print_conv_header, Solid::print_conv_footer and Solid::print_vertical_tip_displacement</h4>
1934 *
1935
1936 *
1937 * This program prints out data in a nice table that is updated
1938 * on a per-iteration basis. The next two functions set up the table
1939 * header and footer:
1940 *
1941 * @code
1942 *   template <int dim,typename NumberType>
1943 *   void Solid<dim,NumberType>::print_conv_header()
1944 *   {
1945 *   static const unsigned int l_width = 87;
1946 *  
1947 *   for (unsigned int i = 0; i < l_width; ++i)
1948 *   std::cout << "_";
1949 *   std::cout << std::endl;
1950 *  
1951 *   std::cout << " SOLVER STEP "
1952 *   << " | LIN_IT LIN_RES RES_NORM "
1953 *   << " RES_U NU_NORM "
1954 *   << " NU_U " << std::endl;
1955 *  
1956 *   for (unsigned int i = 0; i < l_width; ++i)
1957 *   std::cout << "_";
1958 *   std::cout << std::endl;
1959 *   }
1960 *  
1961 *  
1962 *  
1963 *   template <int dim,typename NumberType>
1964 *   void Solid<dim,NumberType>::print_conv_footer()
1965 *   {
1966 *   static const unsigned int l_width = 87;
1967 *  
1968 *   for (unsigned int i = 0; i < l_width; ++i)
1969 *   std::cout << "_";
1970 *   std::cout << std::endl;
1971 *  
1972 *   std::cout << "Relative errors:" << std::endl
1973 *   << "Displacement:\t" << error_update.u / error_update_0.u << std::endl
1974 *   << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
1975 *   << "v / V_0:\t" << vol_current << " / " << vol_reference
1976 *   << std::endl;
1977 *   }
1978 *  
1979 * @endcode
1980 *
1981 * At the end we also output the result that can be compared to that found in
1982 * the literature, namely the displacement at the upper right corner of the
1983 * beam.
1984 *
1985 * @code
1986 *   template <int dim,typename NumberType>
1987 *   void Solid<dim,NumberType>::print_vertical_tip_displacement()
1988 *   {
1989 *   static const unsigned int l_width = 87;
1990 *  
1991 *   for (unsigned int i = 0; i < l_width; ++i)
1992 *   std::cout << "_";
1993 *   std::cout << std::endl;
1994 *  
1995 * @endcode
1996 *
1997 * The measurement point, as stated in the reference paper, is at the midway
1998 * point of the surface on which the traction is applied.
1999 *
2000 * @code
2001 *   const Point<dim> soln_pt = (dim == 3 ?
2002 *   Point<dim>(48.0*parameters.scale, 52.0*parameters.scale, 0.5*parameters.scale) :
2003 *   Point<dim>(48.0*parameters.scale, 52.0*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 AffineConstraints<double> &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) override
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) override
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) override
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 *   const bool apply_dirichlet_bc = (it_nr == 0);
2877 *  
2878 * @endcode
2879 *
2880 * The boundary conditions for the indentation problem are as follows: On
2881 * the -x face (ID = 1), we set up a zero-displacement condition, -y and +y traction
2882 * free boundary condition (don't need to take care); -z and +z faces (ID = 2) are
2883 * not allowed to move along z axis so that it is a plane strain problem.
2884 * Finally, as described earlier, +x face (ID = 11) has an the applied
2885 * distributed shear force (converted by total force per unit area) which
2886 * needs to be taken care as an inhomogeneous Newmann boundary condition.
2887 *
2888
2889 *
2890 * In the following, we will have to tell the function interpolation
2891 * boundary values which components of the solution vector should be
2892 * constrained (i.e., whether it's the x-, y-, z-displacements or
2893 * combinations thereof). This is done using ComponentMask objects (see
2894 * @ref GlossComponentMask) which we can get from the finite element if we
2895 * provide it with an extractor object for the component we wish to
2896 * select. To this end we first set up such extractor objects and later
2897 * use it when generating the relevant component masks:
2898 *
2899
2900 *
2901 *
2902 * @code
2903 *   if (apply_dirichlet_bc)
2904 *   {
2905 *   constraints.clear();
2906 *  
2907 * @endcode
2908 *
2909 * Fixed left hand side of the beam
2910 *
2911 * @code
2912 *   {
2913 *   const int boundary_id = 1;
2914 *   VectorTools::interpolate_boundary_values(dof_handler_ref,
2915 *   boundary_id,
2916 *   Functions::ZeroFunction<dim>(n_components),
2917 *   constraints,
2918 *   fe.component_mask(u_fe));
2919 *   }
2920 *  
2921 * @endcode
2922 *
2923 * Zero Z-displacement through thickness direction
2924 * This corresponds to a plane strain condition being imposed on the beam
2925 *
2926 * @code
2927 *   if (dim == 3)
2928 *   {
2929 *   const int boundary_id = 2;
2930 *   const FEValuesExtractors::Scalar z_displacement(2);
2931 *   VectorTools::interpolate_boundary_values(dof_handler_ref,
2932 *   boundary_id,
2933 *   Functions::ZeroFunction<dim>(n_components),
2934 *   constraints,
2935 *   fe.component_mask(z_displacement));
2936 *   }
2937 *   }
2938 *   else
2939 *   {
2940 *   if (constraints.has_inhomogeneities())
2941 *   {
2942 *   AffineConstraints<double> homogeneous_constraints(constraints);
2943 *   for (unsigned int dof = 0; dof != dof_handler_ref.n_dofs(); ++dof)
2944 *   if (homogeneous_constraints.is_inhomogeneously_constrained(dof))
2945 *   homogeneous_constraints.set_inhomogeneity(dof, 0.0);
2946 *   constraints.clear();
2947 *   constraints.copy_from(homogeneous_constraints);
2948 *   }
2949 *   }
2950 *  
2951 *   constraints.close();
2952 *   }
2953 *  
2954 * @endcode
2955 *
2956 *
2957 * <a name="Solidsolve_linear_system"></a>
2958 * <h4>Solid::solve_linear_system</h4>
2959 * As the system is composed of a single block, defining a solution scheme
2960 * for the linear problem is straight-forward.
2961 *
2962 * @code
2963 *   template <int dim,typename NumberType>
2964 *   std::pair<unsigned int, double>
2965 *   Solid<dim,NumberType>::solve_linear_system(BlockVector<double> &newton_update)
2966 *   {
2967 *   BlockVector<double> A(dofs_per_block);
2968 *   BlockVector<double> B(dofs_per_block);
2969 *  
2970 *   unsigned int lin_it = 0;
2971 *   double lin_res = 0.0;
2972 *  
2973 * @endcode
2974 *
2975 * We solve for the incremental displacement @f$d\mathbf{u}@f$.
2976 *
2977 * @code
2978 *   {
2979 *   timer.enter_subsection("Linear solver");
2980 *   std::cout << " SLV " << std::flush;
2981 *   if (parameters.type_lin == "CG")
2982 *   {
2983 *   const int solver_its = static_cast<unsigned int>(
2984 *   tangent_matrix.block(u_dof, u_dof).m()
2985 *   * parameters.max_iterations_lin);
2986 *   const double tol_sol = parameters.tol_lin
2987 *   * system_rhs.block(u_dof).l2_norm();
2988 *  
2989 *   SolverControl solver_control(solver_its, tol_sol);
2990 *  
2992 *   SolverCG<Vector<double> > solver_CG(solver_control, GVM);
2993 *  
2994 * @endcode
2995 *
2996 * We've chosen by default a SSOR preconditioner as it appears to
2997 * provide the fastest solver convergence characteristics for this
2998 * problem on a single-thread machine. However, for multicore
2999 * computing, the Jacobi preconditioner which is multithreaded may
3000 * converge quicker for larger linear systems.
3001 *
3002 * @code
3003 *   PreconditionSelector<SparseMatrix<double>, Vector<double> >
3004 *   preconditioner (parameters.preconditioner_type,
3005 *   parameters.preconditioner_relaxation);
3006 *   preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
3007 *  
3008 *   solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
3009 *   newton_update.block(u_dof),
3010 *   system_rhs.block(u_dof),
3011 *   preconditioner);
3012 *  
3013 *   lin_it = solver_control.last_step();
3014 *   lin_res = solver_control.last_value();
3015 *   }
3016 *   else if (parameters.type_lin == "Direct")
3017 *   {
3018 * @endcode
3019 *
3020 * Otherwise if the problem is small
3021 * enough, a direct solver can be
3022 * utilised.
3023 *
3024 * @code
3025 *   SparseDirectUMFPACK A_direct;
3026 *   A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
3027 *   A_direct.vmult(newton_update.block(u_dof), system_rhs.block(u_dof));
3028 *  
3029 *   lin_it = 1;
3030 *   lin_res = 0.0;
3031 *   }
3032 *   else
3033 *   Assert (false, ExcMessage("Linear solver type not implemented"));
3034 *  
3035 *   timer.leave_subsection();
3036 *   }
3037 *  
3038 * @endcode
3039 *
3040 * Now that we have the displacement update, distribute the constraints
3041 * back to the Newton update:
3042 *
3043 * @code
3044 *   constraints.distribute(newton_update);
3045 *  
3046 *   return std::make_pair(lin_it, lin_res);
3047 *   }
3048 *  
3049 * @endcode
3050 *
3051 *
3052 * <a name="Solidoutput_results"></a>
3053 * <h4>Solid::output_results</h4>
3054 * Here we present how the results are written to file to be viewed
3055 * using ParaView or Visit. The method is similar to that shown in the
3056 * tutorials so will not be discussed in detail.
3057 *
3058 * @code
3059 *   template <int dim,typename NumberType>
3060 *   void Solid<dim,NumberType>::output_results() const
3061 *   {
3062 *   DataOut<dim> data_out;
3063 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
3064 *   data_component_interpretation(dim,
3065 *   DataComponentInterpretation::component_is_part_of_vector);
3066 *  
3067 *   std::vector<std::string> solution_name(dim, "displacement");
3068 *  
3069 *   data_out.attach_dof_handler(dof_handler_ref);
3070 *   data_out.add_data_vector(solution_n,
3071 *   solution_name,
3072 *   DataOut<dim>::type_dof_data,
3073 *   data_component_interpretation);
3074 *  
3075 * @endcode
3076 *
3077 * Since we are dealing with a large deformation problem, it would be nice
3078 * to display the result on a displaced grid! The MappingQEulerian class
3079 * linked with the DataOut class provides an interface through which this
3080 * can be achieved without physically moving the grid points in the
3081 * Triangulation object ourselves. We first need to copy the solution to
3082 * a temporary vector and then create the Eulerian mapping. We also
3083 * specify the polynomial degree to the DataOut object in order to produce
3084 * a more refined output data set when higher order polynomials are used.
3085 *
3086 * @code
3087 *   Vector<double> soln(solution_n.size());
3088 *   for (unsigned int i = 0; i < soln.size(); ++i)
3089 *   soln(i) = solution_n(i);
3090 *   MappingQEulerian<dim> q_mapping(degree, dof_handler_ref, soln);
3091 *   data_out.build_patches(q_mapping, degree);
3092 *  
3093 *   std::ostringstream filename;
3094 *   filename << "solution-" << time.get_timestep() << ".vtk";
3095 *  
3096 *   std::ofstream output(filename.str().c_str());
3097 *   data_out.write_vtk(output);
3098 *   }
3099 *  
3100 *   }
3101 *  
3102 *  
3103 * @endcode
3104 *
3105 *
3106 * <a name="Mainfunction"></a>
3107 * <h3>Main function</h3>
3108 * Lastly we provide the main driver function which appears
3109 * no different to the other tutorials.
3110 *
3111 * @code
3112 *   int main (int argc, char *argv[])
3113 *   {
3114 *   using namespace dealii;
3115 *   using namespace Cook_Membrane;
3116 *  
3117 *   const unsigned int dim = 3;
3118 *   try
3119 *   {
3120 *   deallog.depth_console(0);
3121 *   Parameters::AllParameters parameters("parameters.prm");
3122 *   if (parameters.automatic_differentiation_order == 0)
3123 *   {
3124 *   std::cout << "Assembly method: Residual and linearisation are computed manually." << std::endl;
3125 *  
3126 * @endcode
3127 *
3128 * Allow multi-threading
3129 *
3130 * @code
3131 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3132 *   ::numbers::invalid_unsigned_int);
3133 *  
3134 *   typedef double NumberType;
3135 *   Solid<dim,NumberType> solid_3d(parameters);
3136 *   solid_3d.run();
3137 *   }
3138 *   #ifdef ENABLE_SACADO_FORMULATION
3139 *   else if (parameters.automatic_differentiation_order == 1)
3140 *   {
3141 *   std::cout << "Assembly method: Residual computed manually; linearisation performed using AD." << std::endl;
3142 *  
3143 * @endcode
3144 *
3145 * Allow multi-threading
3146 *
3147 * @code
3148 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3149 *   ::numbers::invalid_unsigned_int);
3150 *  
3151 *   typedef Sacado::Fad::DFad<double> NumberType;
3152 *   Solid<dim,NumberType> solid_3d(parameters);
3153 *   solid_3d.run();
3154 *   }
3155 *   else if (parameters.automatic_differentiation_order == 2)
3156 *   {
3157 *   std::cout << "Assembly method: Residual and linearisation computed using AD." << std::endl;
3158 *  
3159 * @endcode
3160 *
3161 * Sacado Rad-Fad is not thread-safe, so disable threading.
3162 * Parallisation using MPI would be possible though.
3163 *
3164 * @code
3165 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3166 *   1);
3167 *  
3168 *   typedef Sacado::Rad::ADvar< Sacado::Fad::DFad<double> > NumberType;
3169 *   Solid<dim,NumberType> solid_3d(parameters);
3170 *   solid_3d.run();
3171 *   }
3172 *   #endif
3173 *   else
3174 *   {
3175 *   AssertThrow(false,
3176 *   ExcMessage("The selected assembly method is not supported. "
3177 *   "You need deal.II 9.0 and Trilinos with the Sacado package "
3178 *   "to enable assembly using automatic differentiation."));
3179 *   }
3180 *   }
3181 *   catch (std::exception &exc)
3182 *   {
3183 *   std::cerr << std::endl << std::endl
3184 *   << "----------------------------------------------------"
3185 *   << std::endl;
3186 *   std::cerr << "Exception on processing: " << std::endl << exc.what()
3187 *   << std::endl << "Aborting!" << std::endl
3188 *   << "----------------------------------------------------"
3189 *   << std::endl;
3190 *  
3191 *   return 1;
3192 *   }
3193 *   catch (...)
3194 *   {
3195 *   std::cerr << std::endl << std::endl
3196 *   << "----------------------------------------------------"
3197 *   << std::endl;
3198 *   std::cerr << "Unknown exception!" << std::endl << "Aborting!"
3199 *   << std::endl
3200 *   << "----------------------------------------------------"
3201 *   << std::endl;
3202 *   return 1;
3203 *   }
3204 *  
3205 *   return 0;
3206 *   }
3207 * @endcode
3208
3209
3210*/
Definition fe_q.h:551
Definition point.h:112
numbers::NumberTraits< Number >::real_type norm() const
constexpr void clear()
@ wall_times
Definition timer.h:652
#define DEAL_II_VERSION_MAJOR
Definition config.h:28
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
unsigned int level
Definition grid_out.cc:4618
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > 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:439
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
const Event initial
Definition event.cc:65
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
void approximate(SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &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 >())
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
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)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
double volume(const Triangulation< dim, spacedim > &tria)
@ matrix
Contents is actually a matrix.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:472
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:50
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F_iso(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< 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())
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:71
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int boundary_id
Definition types.h:141
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
const TriangulationDescription::Settings settings
const ::Triangulation< dim, spacedim > & tria