Reference documentation for deal.II version GIT relicensing-426-g7976cfd195 2024-04-18 21:10:01+00:00
\(\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
192 *  
193 *  
194 * @endcode
195 *
196 * We start by including all the necessary deal.II header files and some C++
197 * related ones. They have been discussed in detail in previous tutorial
198 * programs, so you need only refer to past tutorials for details.
199 *
200 * @code
201 *   #include <deal.II/base/function.h>
202 *   #include <deal.II/base/parameter_handler.h>
203 *   #include <deal.II/base/point.h>
204 *   #include <deal.II/base/quadrature_lib.h>
205 *   #include <deal.II/base/symmetric_tensor.h>
206 *   #include <deal.II/base/tensor.h>
207 *   #include <deal.II/base/timer.h>
208 *   #include <deal.II/base/work_stream.h>
209 *   #include <deal.II/base/quadrature_point_data.h>
210 *  
211 *   #include <deal.II/dofs/dof_renumbering.h>
212 *   #include <deal.II/dofs/dof_tools.h>
213 *  
214 *   #include <deal.II/grid/grid_generator.h>
215 *   #include <deal.II/grid/grid_tools.h>
216 *   #include <deal.II/grid/grid_in.h>
217 *   #include <deal.II/grid/tria.h>
218 *  
219 *   #include <deal.II/fe/fe_dgp_monomial.h>
220 *   #include <deal.II/fe/fe_q.h>
221 *   #include <deal.II/fe/fe_system.h>
222 *   #include <deal.II/fe/fe_tools.h>
223 *   #include <deal.II/fe/fe_values.h>
224 *   #include <deal.II/fe/mapping_q_eulerian.h>
225 *   #include <deal.II/fe/mapping_q.h>
226 *  
227 *   #include <deal.II/lac/affine_constraints.h>
228 *   #include <deal.II/lac/block_sparse_matrix.h>
229 *   #include <deal.II/lac/block_vector.h>
230 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
231 *   #include <deal.II/lac/full_matrix.h>
232 *   #include <deal.II/lac/precondition_selector.h>
233 *   #include <deal.II/lac/solver_cg.h>
234 *   #include <deal.II/lac/sparse_direct.h>
235 *  
236 *   #include <deal.II/numerics/data_out.h>
237 *   #include <deal.II/numerics/vector_tools.h>
238 *  
239 *   #include <deal.II/base/config.h>
240 *   #if DEAL_II_VERSION_MAJOR >= 9 && defined(DEAL_II_WITH_TRILINOS)
241 *   #include <deal.II/differentiation/ad.h>
242 *   #define ENABLE_SACADO_FORMULATION
243 *   #endif
244 *  
245 * @endcode
246 *
247 * These must be included below the AD headers so that
248 * their math functions are available for use in the
249 * definition of tensors and kinematic quantities
250 *
251 * @code
252 *   #include <deal.II/physics/elasticity/kinematics.h>
253 *   #include <deal.II/physics/elasticity/standard_tensors.h>
254 *  
255 *   #include <iostream>
256 *   #include <fstream>
257 *   #include <memory>
258 *  
259 *  
260 * @endcode
261 *
262 * We then stick everything that relates to this tutorial program into a
263 * namespace of its own, and import all the deal.II function and class names
264 * into it:
265 *
266 * @code
267 *   namespace Cook_Membrane
268 *   {
269 *   using namespace dealii;
270 *  
271 * @endcode
272 *
273 *
274 * <a name="cook_membrane.cc-Runtimeparameters"></a>
275 * <h3>Run-time parameters</h3>
276 *
277
278 *
279 * There are several parameters that can be set in the code so we set up a
280 * ParameterHandler object to read in the choices at run-time.
281 *
282 * @code
283 *   namespace Parameters
284 *   {
285 * @endcode
286 *
287 *
288 * <a name="cook_membrane.cc-Assemblymethod"></a>
289 * <h4>Assembly method</h4>
290 *
291
292 *
293 * Here we specify whether automatic differentiation is to be used to assemble
294 * the linear system, and if so then what order of differentiation is to be
295 * employed.
296 *
297 * @code
298 *   struct AssemblyMethod
299 *   {
300 *   unsigned int automatic_differentiation_order;
301 *  
302 *   static void
303 *   declare_parameters(ParameterHandler &prm);
304 *  
305 *   void
306 *   parse_parameters(ParameterHandler &prm);
307 *   };
308 *  
309 *  
310 *   void AssemblyMethod::declare_parameters(ParameterHandler &prm)
311 *   {
312 *   prm.enter_subsection("Assembly method");
313 *   {
314 *   prm.declare_entry("Automatic differentiation order", "0",
315 *   Patterns::Integer(0,2),
316 *   "The automatic differentiation order to be used in the assembly of the linear system.\n"
317 *   "# Order = 0: Both the residual and linearisation are computed manually.\n"
318 *   "# Order = 1: The residual is computed manually but the linearisation is performed using AD.\n"
319 *   "# Order = 2: Both the residual and linearisation are computed using AD.");
320 *   }
321 *   prm.leave_subsection();
322 *   }
323 *  
324 *   void AssemblyMethod::parse_parameters(ParameterHandler &prm)
325 *   {
326 *   prm.enter_subsection("Assembly method");
327 *   {
328 *   automatic_differentiation_order = prm.get_integer("Automatic differentiation order");
329 *   }
330 *   prm.leave_subsection();
331 *   }
332 *  
333 * @endcode
334 *
335 *
336 * <a name="cook_membrane.cc-FiniteElementsystem"></a>
337 * <h4>Finite Element system</h4>
338 *
339
340 *
341 * Here we specify the polynomial order used to approximate the solution.
342 * The quadrature order should be adjusted accordingly.
343 *
344 * @code
345 *   struct FESystem
346 *   {
347 *   unsigned int poly_degree;
348 *   unsigned int quad_order;
349 *  
350 *   static void
351 *   declare_parameters(ParameterHandler &prm);
352 *  
353 *   void
354 *   parse_parameters(ParameterHandler &prm);
355 *   };
356 *  
357 *  
358 *   void FESystem::declare_parameters(ParameterHandler &prm)
359 *   {
360 *   prm.enter_subsection("Finite element system");
361 *   {
362 *   prm.declare_entry("Polynomial degree", "2",
363 *   Patterns::Integer(0),
364 *   "Displacement system polynomial order");
365 *  
366 *   prm.declare_entry("Quadrature order", "3",
367 *   Patterns::Integer(0),
368 *   "Gauss quadrature order");
369 *   }
370 *   prm.leave_subsection();
371 *   }
372 *  
373 *   void FESystem::parse_parameters(ParameterHandler &prm)
374 *   {
375 *   prm.enter_subsection("Finite element system");
376 *   {
377 *   poly_degree = prm.get_integer("Polynomial degree");
378 *   quad_order = prm.get_integer("Quadrature order");
379 *   }
380 *   prm.leave_subsection();
381 *   }
382 *  
383 * @endcode
384 *
385 *
386 * <a name="cook_membrane.cc-Geometry"></a>
387 * <h4>Geometry</h4>
388 *
389
390 *
391 * Make adjustments to the problem geometry and its discretisation.
392 *
393 * @code
394 *   struct Geometry
395 *   {
396 *   unsigned int elements_per_edge;
397 *   double scale;
398 *  
399 *   static void
400 *   declare_parameters(ParameterHandler &prm);
401 *  
402 *   void
403 *   parse_parameters(ParameterHandler &prm);
404 *   };
405 *  
406 *   void Geometry::declare_parameters(ParameterHandler &prm)
407 *   {
408 *   prm.enter_subsection("Geometry");
409 *   {
410 *   prm.declare_entry("Elements per edge", "32",
411 *   Patterns::Integer(0),
412 *   "Number of elements per long edge of the beam");
413 *  
414 *   prm.declare_entry("Grid scale", "1e-3",
415 *   Patterns::Double(0.0),
416 *   "Global grid scaling factor");
417 *   }
418 *   prm.leave_subsection();
419 *   }
420 *  
421 *   void Geometry::parse_parameters(ParameterHandler &prm)
422 *   {
423 *   prm.enter_subsection("Geometry");
424 *   {
425 *   elements_per_edge = prm.get_integer("Elements per edge");
426 *   scale = prm.get_double("Grid scale");
427 *   }
428 *   prm.leave_subsection();
429 *   }
430 *  
431 * @endcode
432 *
433 *
434 * <a name="cook_membrane.cc-Materials"></a>
435 * <h4>Materials</h4>
436 *
437
438 *
439 * We also need the shear modulus @f$ \mu @f$ and Poisson ration @f$ \nu @f$ for the
440 * neo-Hookean material.
441 *
442 * @code
443 *   struct Materials
444 *   {
445 *   double nu;
446 *   double mu;
447 *  
448 *   static void
449 *   declare_parameters(ParameterHandler &prm);
450 *  
451 *   void
452 *   parse_parameters(ParameterHandler &prm);
453 *   };
454 *  
455 *   void Materials::declare_parameters(ParameterHandler &prm)
456 *   {
457 *   prm.enter_subsection("Material properties");
458 *   {
459 *   prm.declare_entry("Poisson's ratio", "0.3",
460 *   Patterns::Double(-1.0,0.5),
461 *   "Poisson's ratio");
462 *  
463 *   prm.declare_entry("Shear modulus", "0.4225e6",
464 *   Patterns::Double(),
465 *   "Shear modulus");
466 *   }
467 *   prm.leave_subsection();
468 *   }
469 *  
470 *   void Materials::parse_parameters(ParameterHandler &prm)
471 *   {
472 *   prm.enter_subsection("Material properties");
473 *   {
474 *   nu = prm.get_double("Poisson's ratio");
475 *   mu = prm.get_double("Shear modulus");
476 *   }
477 *   prm.leave_subsection();
478 *   }
479 *  
480 * @endcode
481 *
482 *
483 * <a name="cook_membrane.cc-Linearsolver"></a>
484 * <h4>Linear solver</h4>
485 *
486
487 *
488 * Next, we choose both solver and preconditioner settings. The use of an
489 * effective preconditioner is critical to ensure convergence when a large
490 * nonlinear motion occurs within a Newton increment.
491 *
492 * @code
493 *   struct LinearSolver
494 *   {
495 *   std::string type_lin;
496 *   double tol_lin;
497 *   double max_iterations_lin;
498 *   std::string preconditioner_type;
499 *   double preconditioner_relaxation;
500 *  
501 *   static void
502 *   declare_parameters(ParameterHandler &prm);
503 *  
504 *   void
505 *   parse_parameters(ParameterHandler &prm);
506 *   };
507 *  
508 *   void LinearSolver::declare_parameters(ParameterHandler &prm)
509 *   {
510 *   prm.enter_subsection("Linear solver");
511 *   {
512 *   prm.declare_entry("Solver type", "CG",
513 *   Patterns::Selection("CG|Direct"),
514 *   "Type of solver used to solve the linear system");
515 *  
516 *   prm.declare_entry("Residual", "1e-6",
517 *   Patterns::Double(0.0),
518 *   "Linear solver residual (scaled by residual norm)");
519 *  
520 *   prm.declare_entry("Max iteration multiplier", "1",
521 *   Patterns::Double(0.0),
522 *   "Linear solver iterations (multiples of the system matrix size)");
523 *  
524 *   prm.declare_entry("Preconditioner type", "ssor",
525 *   Patterns::Selection("jacobi|ssor"),
526 *   "Type of preconditioner");
527 *  
528 *   prm.declare_entry("Preconditioner relaxation", "0.65",
529 *   Patterns::Double(0.0),
530 *   "Preconditioner relaxation value");
531 *   }
532 *   prm.leave_subsection();
533 *   }
534 *  
535 *   void LinearSolver::parse_parameters(ParameterHandler &prm)
536 *   {
537 *   prm.enter_subsection("Linear solver");
538 *   {
539 *   type_lin = prm.get("Solver type");
540 *   tol_lin = prm.get_double("Residual");
541 *   max_iterations_lin = prm.get_double("Max iteration multiplier");
542 *   preconditioner_type = prm.get("Preconditioner type");
543 *   preconditioner_relaxation = prm.get_double("Preconditioner relaxation");
544 *   }
545 *   prm.leave_subsection();
546 *   }
547 *  
548 * @endcode
549 *
550 *
551 * <a name="cook_membrane.cc-Nonlinearsolver"></a>
552 * <h4>Nonlinear solver</h4>
553 *
554
555 *
556 * A Newton-Raphson scheme is used to solve the nonlinear system of governing
557 * equations. We now define the tolerances and the maximum number of
558 * iterations for the Newton-Raphson nonlinear solver.
559 *
560 * @code
561 *   struct NonlinearSolver
562 *   {
563 *   unsigned int max_iterations_NR;
564 *   double tol_f;
565 *   double tol_u;
566 *  
567 *   static void
568 *   declare_parameters(ParameterHandler &prm);
569 *  
570 *   void
571 *   parse_parameters(ParameterHandler &prm);
572 *   };
573 *  
574 *   void NonlinearSolver::declare_parameters(ParameterHandler &prm)
575 *   {
576 *   prm.enter_subsection("Nonlinear solver");
577 *   {
578 *   prm.declare_entry("Max iterations Newton-Raphson", "10",
579 *   Patterns::Integer(0),
580 *   "Number of Newton-Raphson iterations allowed");
581 *  
582 *   prm.declare_entry("Tolerance force", "1.0e-9",
583 *   Patterns::Double(0.0),
584 *   "Force residual tolerance");
585 *  
586 *   prm.declare_entry("Tolerance displacement", "1.0e-6",
587 *   Patterns::Double(0.0),
588 *   "Displacement error tolerance");
589 *   }
590 *   prm.leave_subsection();
591 *   }
592 *  
593 *   void NonlinearSolver::parse_parameters(ParameterHandler &prm)
594 *   {
595 *   prm.enter_subsection("Nonlinear solver");
596 *   {
597 *   max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
598 *   tol_f = prm.get_double("Tolerance force");
599 *   tol_u = prm.get_double("Tolerance displacement");
600 *   }
601 *   prm.leave_subsection();
602 *   }
603 *  
604 * @endcode
605 *
606 *
607 * <a name="cook_membrane.cc-Time"></a>
608 * <h4>Time</h4>
609 *
610
611 *
612 * Set the timestep size @f$ \varDelta t @f$ and the simulation end-time.
613 *
614 * @code
615 *   struct Time
616 *   {
617 *   double delta_t;
618 *   double end_time;
619 *  
620 *   static void
621 *   declare_parameters(ParameterHandler &prm);
622 *  
623 *   void
624 *   parse_parameters(ParameterHandler &prm);
625 *   };
626 *  
627 *   void Time::declare_parameters(ParameterHandler &prm)
628 *   {
629 *   prm.enter_subsection("Time");
630 *   {
631 *   prm.declare_entry("End time", "1",
632 *   Patterns::Double(),
633 *   "End time");
634 *  
635 *   prm.declare_entry("Time step size", "0.1",
636 *   Patterns::Double(),
637 *   "Time step size");
638 *   }
639 *   prm.leave_subsection();
640 *   }
641 *  
642 *   void Time::parse_parameters(ParameterHandler &prm)
643 *   {
644 *   prm.enter_subsection("Time");
645 *   {
646 *   end_time = prm.get_double("End time");
647 *   delta_t = prm.get_double("Time step size");
648 *   }
649 *   prm.leave_subsection();
650 *   }
651 *  
652 * @endcode
653 *
654 *
655 * <a name="cook_membrane.cc-Allparameters"></a>
656 * <h4>All parameters</h4>
657 *
658
659 *
660 * Finally we consolidate all of the above structures into a single container
661 * that holds all of our run-time selections.
662 *
663 * @code
664 *   struct AllParameters :
665 *   public AssemblyMethod,
666 *   public FESystem,
667 *   public Geometry,
668 *   public Materials,
669 *   public LinearSolver,
670 *   public NonlinearSolver,
671 *   public Time
672 *  
673 *   {
674 *   AllParameters(const std::string &input_file);
675 *  
676 *   static void
677 *   declare_parameters(ParameterHandler &prm);
678 *  
679 *   void
680 *   parse_parameters(ParameterHandler &prm);
681 *   };
682 *  
683 *   AllParameters::AllParameters(const std::string &input_file)
684 *   {
685 *   ParameterHandler prm;
686 *   declare_parameters(prm);
687 *   prm.parse_input(input_file);
688 *   parse_parameters(prm);
689 *   }
690 *  
691 *   void AllParameters::declare_parameters(ParameterHandler &prm)
692 *   {
693 *   AssemblyMethod::declare_parameters(prm);
694 *   FESystem::declare_parameters(prm);
695 *   Geometry::declare_parameters(prm);
696 *   Materials::declare_parameters(prm);
697 *   LinearSolver::declare_parameters(prm);
698 *   NonlinearSolver::declare_parameters(prm);
699 *   Time::declare_parameters(prm);
700 *   }
701 *  
702 *   void AllParameters::parse_parameters(ParameterHandler &prm)
703 *   {
704 *   AssemblyMethod::parse_parameters(prm);
705 *   FESystem::parse_parameters(prm);
706 *   Geometry::parse_parameters(prm);
707 *   Materials::parse_parameters(prm);
708 *   LinearSolver::parse_parameters(prm);
709 *   NonlinearSolver::parse_parameters(prm);
710 *   Time::parse_parameters(prm);
711 *   }
712 *   }
713 *  
714 *  
715 * @endcode
716 *
717 *
718 * <a name="cook_membrane.cc-Timeclass"></a>
719 * <h3>Time class</h3>
720 *
721
722 *
723 * A simple class to store time data. Its functioning is transparent so no
724 * discussion is necessary. For simplicity we assume a constant time step
725 * size.
726 *
727 * @code
728 *   class Time
729 *   {
730 *   public:
731 *   Time (const double time_end,
732 *   const double delta_t)
733 *   :
734 *   timestep(0),
735 *   time_current(0.0),
736 *   time_end(time_end),
737 *   delta_t(delta_t)
738 *   {}
739 *  
740 *   virtual ~Time()
741 *   {}
742 *  
743 *   double current() const
744 *   {
745 *   return time_current;
746 *   }
747 *   double end() const
748 *   {
749 *   return time_end;
750 *   }
751 *   double get_delta_t() const
752 *   {
753 *   return delta_t;
754 *   }
755 *   unsigned int get_timestep() const
756 *   {
757 *   return timestep;
758 *   }
759 *   void increment()
760 *   {
761 *   time_current += delta_t;
762 *   ++timestep;
763 *   }
764 *  
765 *   private:
766 *   unsigned int timestep;
767 *   double time_current;
768 *   const double time_end;
769 *   const double delta_t;
770 *   };
771 *  
772 * @endcode
773 *
774 *
775 * <a name="cook_membrane.cc-CompressibleneoHookeanmaterialwithinaonefieldformulation"></a>
776 * <h3>Compressible neo-Hookean material within a one-field formulation</h3>
777 *
778
779 *
780 * As discussed in the literature and @ref step_44 "step-44", Neo-Hookean materials are a type
781 * of hyperelastic materials. The entire domain is assumed to be composed of a
782 * compressible neo-Hookean material. This class defines the behaviour of
783 * this material within a one-field formulation. Compressible neo-Hookean
784 * materials can be described by a strain-energy function (SEF) @f$ \Psi =
785 * \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(J)
786 * @f$.
787 *
788
789 *
790 * The isochoric response is given by @f$
791 * \Psi_{\text{iso}}(\overline{\mathbf{b}}) = c_{1} [\overline{I}_{1} - 3] @f$
792 * where @f$ c_{1} = \frac{\mu}{2} @f$ and @f$\overline{I}_{1}@f$ is the first
793 * invariant of the left- or right-isochoric Cauchy-Green deformation tensors.
794 * That is @f$\overline{I}_1 :=\textrm{tr}(\overline{\mathbf{b}})@f$. In this
795 * example the SEF that governs the volumetric response is defined as @f$
796 * \Psi_{\text{vol}}(J) = \kappa \frac{1}{4} [ J^2 - 1
797 * - 2\textrm{ln}\; J ]@f$, where @f$\kappa:= \lambda + 2/3 \mu@f$ is
798 * the <a href="http://en.wikipedia.org/wiki/Bulk_modulus">bulk modulus</a>
799 * and @f$\lambda@f$ is <a
800 * href="http://en.wikipedia.org/wiki/Lam%C3%A9_parameters">Lame's first
801 * parameter</a>.
802 *
803
804 *
805 * The following class will be used to characterize the material we work with,
806 * and provides a central point that one would need to modify if one were to
807 * implement a different material model. For it to work, we will store one
808 * object of this type per quadrature point, and in each of these objects
809 * store the current state (characterized by the values or measures of the
810 * displacement field) so that we can compute the elastic coefficients
811 * linearized around the current state.
812 *
813 * @code
814 *   template <int dim,typename NumberType>
815 *   class Material_Compressible_Neo_Hook_One_Field
816 *   {
817 *   public:
818 *   Material_Compressible_Neo_Hook_One_Field(const double mu,
819 *   const double nu)
820 *   :
821 *   kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
822 *   c_1(mu / 2.0)
823 *   {
824 *   Assert(kappa > 0, ExcInternalError());
825 *   }
826 *  
827 *   ~Material_Compressible_Neo_Hook_One_Field()
828 *   {}
829 *  
830 * @endcode
831 *
832 * The first function is the total energy
833 * @f$\Psi = \Psi_{\textrm{iso}} + \Psi_{\textrm{vol}}@f$.
834 *
835 * @code
836 *   NumberType
837 *   get_Psi(const NumberType &det_F,
838 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
839 *   {
840 *   return get_Psi_vol(det_F) + get_Psi_iso(b_bar);
841 *   }
842 *  
843 * @endcode
844 *
845 * The second function determines the Kirchhoff stress @f$\boldsymbol{\tau}
846 * = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}@f$
847 *
848 * @code
849 *   SymmetricTensor<2,dim,NumberType>
850 *   get_tau(const NumberType &det_F,
851 *   const SymmetricTensor<2,dim,NumberType> &b_bar)
852 *   {
853 * @endcode
854 *
855 * See Holzapfel p231 eq6.98 onwards
856 *
857 * @code
858 *   return get_tau_vol(det_F) + get_tau_iso(b_bar);
859 *   }
860 *  
861 * @endcode
862 *
863 * The fourth-order elasticity tensor in the spatial setting
864 * @f$\mathfrak{c}@f$ is calculated from the SEF @f$\Psi@f$ as @f$ J
865 * \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}@f$
866 * where @f$ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial
867 * \mathbf{C} \partial \mathbf{C}}@f$
868 *
869 * @code
870 *   SymmetricTensor<4,dim,NumberType>
871 *   get_Jc(const NumberType &det_F,
872 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
873 *   {
874 *   return get_Jc_vol(det_F) + get_Jc_iso(b_bar);
875 *   }
876 *  
877 *   private:
878 * @endcode
879 *
880 * Define constitutive model parameters @f$\kappa@f$ (bulk modulus) and the
881 * neo-Hookean model parameter @f$c_1@f$:
882 *
883 * @code
884 *   const double kappa;
885 *   const double c_1;
886 *  
887 * @endcode
888 *
889 * Value of the volumetric free energy
890 *
891 * @code
892 *   NumberType
893 *   get_Psi_vol(const NumberType &det_F) const
894 *   {
895 *   return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
896 *   }
897 *  
898 * @endcode
899 *
900 * Value of the isochoric free energy
901 *
902 * @code
903 *   NumberType
904 *   get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
905 *   {
906 *   return c_1 * (trace(b_bar) - dim);
907 *   }
908 *  
909 * @endcode
910 *
911 * Derivative of the volumetric free energy with respect to
912 * @f$J@f$ return @f$\frac{\partial
913 * \Psi_{\text{vol}}(J)}{\partial J}@f$
914 *
915 * @code
916 *   NumberType
917 *   get_dPsi_vol_dJ(const NumberType &det_F) const
918 *   {
919 *   return (kappa / 2.0) * (det_F - 1.0 / det_F);
920 *   }
921 *  
922 * @endcode
923 *
924 * The following functions are used internally in determining the result
925 * of some of the public functions above. The first one determines the
926 * volumetric Kirchhoff stress @f$\boldsymbol{\tau}_{\textrm{vol}}@f$.
927 * Note the difference in its definition when compared to @ref step_44 "step-44".
928 *
929 * @code
930 *   SymmetricTensor<2,dim,NumberType>
931 *   get_tau_vol(const NumberType &det_F) const
932 *   {
933 *   return NumberType(get_dPsi_vol_dJ(det_F) * det_F) * Physics::Elasticity::StandardTensors<dim>::I;
934 *   }
935 *  
936 * @endcode
937 *
938 * Next, determine the isochoric Kirchhoff stress
939 * @f$\boldsymbol{\tau}_{\textrm{iso}} =
940 * \mathcal{P}:\overline{\boldsymbol{\tau}}@f$:
941 *
942 * @code
943 *   SymmetricTensor<2,dim,NumberType>
944 *   get_tau_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
945 *   {
946 *   return Physics::Elasticity::StandardTensors<dim>::dev_P * get_tau_bar(b_bar);
947 *   }
948 *  
949 * @endcode
950 *
951 * Then, determine the fictitious Kirchhoff stress
952 * @f$\overline{\boldsymbol{\tau}}@f$:
953 *
954 * @code
955 *   SymmetricTensor<2,dim,NumberType>
956 *   get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const
957 *   {
958 *   return 2.0 * c_1 * b_bar;
959 *   }
960 *  
961 * @endcode
962 *
963 * Second derivative of the volumetric free energy wrt @f$J@f$. We
964 * need the following computation explicitly in the tangent so we make it
965 * public. We calculate @f$\frac{\partial^2
966 * \Psi_{\textrm{vol}}(J)}{\partial J \partial
967 * J}@f$
968 *
969 * @code
970 *   NumberType
971 *   get_d2Psi_vol_dJ2(const NumberType &det_F) const
972 *   {
973 *   return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
974 *   }
975 *  
976 * @endcode
977 *
978 * Calculate the volumetric part of the tangent @f$J
979 * \mathfrak{c}_\textrm{vol}@f$. Again, note the difference in its
980 * definition when compared to @ref step_44 "step-44". The extra terms result from two
981 * quantities in @f$\boldsymbol{\tau}_{\textrm{vol}}@f$ being dependent on
982 * @f$\boldsymbol{F}@f$.
983 *
984 * @code
985 *   SymmetricTensor<4,dim,NumberType>
986 *   get_Jc_vol(const NumberType &det_F) const
987 *   {
988 * @endcode
989 *
990 * See Holzapfel p265
991 *
992 * @code
993 *   return det_F
994 *   * ( (get_dPsi_vol_dJ(det_F) + det_F * get_d2Psi_vol_dJ2(det_F))*Physics::Elasticity::StandardTensors<dim>::IxI
995 *   - (2.0 * get_dPsi_vol_dJ(det_F))*Physics::Elasticity::StandardTensors<dim>::S );
996 *   }
997 *  
998 * @endcode
999 *
1000 * Calculate the isochoric part of the tangent @f$J
1001 * \mathfrak{c}_\textrm{iso}@f$:
1002 *
1003 * @code
1004 *   SymmetricTensor<4,dim,NumberType>
1005 *   get_Jc_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
1006 *   {
1007 *   const SymmetricTensor<2, dim> tau_bar = get_tau_bar(b_bar);
1008 *   const SymmetricTensor<2, dim> tau_iso = get_tau_iso(b_bar);
1009 *   const SymmetricTensor<4, dim> tau_iso_x_I
1010 *   = outer_product(tau_iso,
1011 *   Physics::Elasticity::StandardTensors<dim>::I);
1012 *   const SymmetricTensor<4, dim> I_x_tau_iso
1013 *   = outer_product(Physics::Elasticity::StandardTensors<dim>::I,
1014 *   tau_iso);
1015 *   const SymmetricTensor<4, dim> c_bar = get_c_bar();
1016 *  
1017 *   return (2.0 / dim) * trace(tau_bar)
1018 *   * Physics::Elasticity::StandardTensors<dim>::dev_P
1019 *   - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso)
1020 *   + Physics::Elasticity::StandardTensors<dim>::dev_P * c_bar
1021 *   * Physics::Elasticity::StandardTensors<dim>::dev_P;
1022 *   }
1023 *  
1024 * @endcode
1025 *
1026 * Calculate the fictitious elasticity tensor @f$\overline{\mathfrak{c}}@f$.
1027 * For the material model chosen this is simply zero:
1028 *
1029 * @code
1030 *   SymmetricTensor<4,dim,double>
1031 *   get_c_bar() const
1032 *   {
1033 *   return SymmetricTensor<4, dim>();
1034 *   }
1035 *   };
1036 *  
1037 * @endcode
1038 *
1039 *
1040 * <a name="cook_membrane.cc-Quadraturepointhistory"></a>
1041 * <h3>Quadrature point history</h3>
1042 *
1043
1044 *
1045 * As seen in @ref step_18 "step-18", the <code> PointHistory </code> class offers a method
1046 * for storing data at the quadrature points. Here each quadrature point
1047 * holds a pointer to a material description. Thus, different material models
1048 * can be used in different regions of the domain. Among other data, we
1049 * choose to store the Kirchhoff stress @f$\boldsymbol{\tau}@f$ and the tangent
1050 * @f$J\mathfrak{c}@f$ for the quadrature points.
1051 *
1052 * @code
1053 *   template <int dim,typename NumberType>
1054 *   class PointHistory
1055 *   {
1056 *   public:
1057 *   PointHistory()
1058 *   {}
1059 *  
1060 *   virtual ~PointHistory()
1061 *   {}
1062 *  
1063 * @endcode
1064 *
1065 * The first function is used to create a material object and to
1066 * initialize all tensors correctly: The second one updates the stored
1067 * values and stresses based on the current deformation measure
1068 * @f$\textrm{Grad}\mathbf{u}_{\textrm{n}}@f$.
1069 *
1070 * @code
1071 *   void setup_lqp (const Parameters::AllParameters &parameters)
1072 *   {
1073 *   material.reset(new Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu,
1074 *   parameters.nu));
1075 *   }
1076 *  
1077 * @endcode
1078 *
1079 * We offer an interface to retrieve certain data.
1080 * This is the strain energy:
1081 *
1082 * @code
1083 *   NumberType
1084 *   get_Psi(const NumberType &det_F,
1085 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
1086 *   {
1087 *   return material->get_Psi(det_F,b_bar);
1088 *   }
1089 *  
1090 * @endcode
1091 *
1092 * Here are the kinetic variables. These are used in the material and
1093 * global tangent matrix and residual assembly operations:
1094 * First is the Kirchhoff stress:
1095 *
1096 * @code
1097 *   SymmetricTensor<2,dim,NumberType>
1098 *   get_tau(const NumberType &det_F,
1099 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
1100 *   {
1101 *   return material->get_tau(det_F,b_bar);
1102 *   }
1103 *  
1104 * @endcode
1105 *
1106 * And the tangent:
1107 *
1108 * @code
1109 *   SymmetricTensor<4,dim,NumberType>
1110 *   get_Jc(const NumberType &det_F,
1111 *   const SymmetricTensor<2,dim,NumberType> &b_bar) const
1112 *   {
1113 *   return material->get_Jc(det_F,b_bar);
1114 *   }
1115 *  
1116 * @endcode
1117 *
1118 * In terms of member functions, this class stores for the quadrature
1119 * point it represents a copy of a material type in case different
1120 * materials are used in different regions of the domain, as well as the
1121 * inverse of the deformation gradient...
1122 *
1123 * @code
1124 *   private:
1125 *   std::shared_ptr< Material_Compressible_Neo_Hook_One_Field<dim,NumberType> > material;
1126 *   };
1127 *  
1128 *  
1129 * @endcode
1130 *
1131 *
1132 * <a name="cook_membrane.cc-Quasistaticcompressiblefinitestrainsolid"></a>
1133 * <h3>Quasi-static compressible finite-strain solid</h3>
1134 *
1135
1136 *
1137 * Forward declarations for classes that will
1138 * perform assembly of the linear system.
1139 *
1140 * @code
1141 *   template <int dim,typename NumberType>
1142 *   struct Assembler_Base;
1143 *   template <int dim,typename NumberType>
1144 *   struct Assembler;
1145 *  
1146 * @endcode
1147 *
1148 * The Solid class is the central class in that it represents the problem at
1149 * hand. It follows the usual scheme in that all it really has is a
1150 * constructor, destructor and a <code>run()</code> function that dispatches
1151 * all the work to private functions of this class:
1152 *
1153 * @code
1154 *   template <int dim,typename NumberType>
1155 *   class Solid
1156 *   {
1157 *   public:
1158 *   Solid(const Parameters::AllParameters &parameters);
1159 *  
1160 *   virtual
1161 *   ~Solid();
1162 *  
1163 *   void
1164 *   run();
1165 *  
1166 *   private:
1167 *  
1168 * @endcode
1169 *
1170 * We start the collection of member functions with one that builds the
1171 * grid:
1172 *
1173 * @code
1174 *   void
1175 *   make_grid();
1176 *  
1177 * @endcode
1178 *
1179 * Set up the finite element system to be solved:
1180 *
1181 * @code
1182 *   void
1183 *   system_setup();
1184 *  
1185 * @endcode
1186 *
1187 * Several functions to assemble the system and right hand side matrices
1188 * using multithreading. Each of them comes as a wrapper function, one
1189 * that is executed to do the work in the WorkStream model on one cell,
1190 * and one that copies the work done on this one cell into the global
1191 * object that represents it:
1192 *
1193 * @code
1194 *   void
1195 *   assemble_system(const BlockVector<double> &solution_delta);
1196 *  
1197 * @endcode
1198 *
1199 * We use a separate data structure to perform the assembly. It needs access
1200 * to some low-level data, so we simply befriend the class instead of
1201 * creating a complex interface to provide access as necessary.
1202 *
1203 * @code
1204 *   friend struct Assembler_Base<dim,NumberType>;
1205 *   friend struct Assembler<dim,NumberType>;
1206 *  
1207 * @endcode
1208 *
1209 * Apply Dirichlet boundary conditions on the displacement field
1210 *
1211 * @code
1212 *   void
1213 *   make_constraints(const int &it_nr);
1214 *  
1215 * @endcode
1216 *
1217 * Create and update the quadrature points. Here, no data needs to be
1218 * copied into a global object, so the copy_local_to_global function is
1219 * empty:
1220 *
1221 * @code
1222 *   void
1223 *   setup_qph();
1224 *  
1225 * @endcode
1226 *
1227 * Solve for the displacement using a Newton-Raphson method. We break this
1228 * function into the nonlinear loop and the function that solves the
1229 * linearized Newton-Raphson step:
1230 *
1231 * @code
1232 *   void
1233 *   solve_nonlinear_timestep(BlockVector<double> &solution_delta);
1234 *  
1235 *   std::pair<unsigned int, double>
1236 *   solve_linear_system(BlockVector<double> &newton_update);
1237 *  
1238 * @endcode
1239 *
1240 * Solution retrieval as well as post-processing and writing data to file :
1241 *
1242 * @code
1243 *   BlockVector<double>
1244 *   get_total_solution(const BlockVector<double> &solution_delta) const;
1245 *  
1246 *   void
1247 *   output_results() const;
1248 *  
1249 * @endcode
1250 *
1251 * Finally, some member variables that describe the current state: A
1252 * collection of the parameters used to describe the problem setup...
1253 *
1254 * @code
1255 *   const Parameters::AllParameters &parameters;
1256 *  
1257 * @endcode
1258 *
1259 * ...the volume of the reference and current configurations...
1260 *
1261 * @code
1262 *   double vol_reference;
1263 *   double vol_current;
1264 *  
1265 * @endcode
1266 *
1267 * ...and description of the geometry on which the problem is solved:
1268 *
1269 * @code
1270 *   Triangulation<dim> triangulation;
1271 *  
1272 * @endcode
1273 *
1274 * Also, keep track of the current time and the time spent evaluating
1275 * certain functions
1276 *
1277 * @code
1278 *   Time time;
1279 *   TimerOutput timer;
1280 *  
1281 * @endcode
1282 *
1283 * A storage object for quadrature point information. As opposed to
1284 * @ref step_18 "step-18", deal.II's native quadrature point data manager is employed here.
1285 *
1286 * @code
1288 *   PointHistory<dim,NumberType> > quadrature_point_history;
1289 *  
1290 * @endcode
1291 *
1292 * A description of the finite-element system including the displacement
1293 * polynomial degree, the degree-of-freedom handler, number of DoFs per
1294 * cell and the extractor objects used to retrieve information from the
1295 * solution vectors:
1296 *
1297 * @code
1298 *   const unsigned int degree;
1299 *   const FESystem<dim> fe;
1300 *   DoFHandler<dim> dof_handler_ref;
1301 *   const unsigned int dofs_per_cell;
1302 *   const FEValuesExtractors::Vector u_fe;
1303 *  
1304 * @endcode
1305 *
1306 * Description of how the block-system is arranged. There is just 1 block,
1307 * that contains a vector DOF @f$\mathbf{u}@f$.
1308 * There are two reasons that we retain the block system in this problem.
1309 * The first is pure laziness to perform further modifications to the
1310 * code from which this work originated. The second is that a block system
1311 * would typically necessary when extending this code to multiphysics
1312 * problems.
1313 *
1314 * @code
1315 *   static const unsigned int n_blocks = 1;
1316 *   static const unsigned int n_components = dim;
1317 *   static const unsigned int first_u_component = 0;
1318 *  
1319 *   enum
1320 *   {
1321 *   u_dof = 0
1322 *   };
1323 *  
1324 *   std::vector<types::global_dof_index> dofs_per_block;
1325 *  
1326 * @endcode
1327 *
1328 * Rules for Gauss-quadrature on both the cell and faces. The number of
1329 * quadrature points on both cells and faces is recorded.
1330 *
1331 * @code
1332 *   const QGauss<dim> qf_cell;
1333 *   const QGauss<dim - 1> qf_face;
1334 *   const unsigned int n_q_points;
1335 *   const unsigned int n_q_points_f;
1336 *  
1337 * @endcode
1338 *
1339 * Objects that store the converged solution and right-hand side vectors,
1340 * as well as the tangent matrix. There is a AffineConstraints object used
1341 * to keep track of constraints. We make use of a sparsity pattern
1342 * designed for a block system.
1343 *
1344 * @code
1345 *   AffineConstraints<double> constraints;
1346 *   BlockSparsityPattern sparsity_pattern;
1347 *   BlockSparseMatrix<double> tangent_matrix;
1348 *   BlockVector<double> system_rhs;
1349 *   BlockVector<double> solution_n;
1350 *  
1351 * @endcode
1352 *
1353 * Then define a number of variables to store norms and update norms and
1354 * normalisation factors.
1355 *
1356 * @code
1357 *   struct Errors
1358 *   {
1359 *   Errors()
1360 *   :
1361 *   norm(1.0), u(1.0)
1362 *   {}
1363 *  
1364 *   void reset()
1365 *   {
1366 *   norm = 1.0;
1367 *   u = 1.0;
1368 *   }
1369 *   void normalise(const Errors &rhs)
1370 *   {
1371 *   if (rhs.norm != 0.0)
1372 *   norm /= rhs.norm;
1373 *   if (rhs.u != 0.0)
1374 *   u /= rhs.u;
1375 *   }
1376 *  
1377 *   double norm, u;
1378 *   };
1379 *  
1380 *   Errors error_residual, error_residual_0, error_residual_norm, error_update,
1381 *   error_update_0, error_update_norm;
1382 *  
1383 * @endcode
1384 *
1385 * Methods to calculate error measures
1386 *
1387 * @code
1388 *   void
1389 *   get_error_residual(Errors &error_residual);
1390 *  
1391 *   void
1392 *   get_error_update(const BlockVector<double> &newton_update,
1393 *   Errors &error_update);
1394 *  
1395 * @endcode
1396 *
1397 * Print information to screen in a pleasing way...
1398 *
1399 * @code
1400 *   static
1401 *   void
1402 *   print_conv_header();
1403 *  
1404 *   void
1405 *   print_conv_footer();
1406 *  
1407 *   void
1408 *   print_vertical_tip_displacement();
1409 *   };
1410 *  
1411 * @endcode
1412 *
1413 *
1414 * <a name="cook_membrane.cc-ImplementationofthecodeSolidcodeclass"></a>
1415 * <h3>Implementation of the <code>Solid</code> class</h3>
1416 *
1417
1418 *
1419 *
1420 * <a name="cook_membrane.cc-Publicinterface"></a>
1421 * <h4>Public interface</h4>
1422 *
1423
1424 *
1425 * We initialise the Solid class using data extracted from the parameter file.
1426 *
1427 * @code
1428 *   template <int dim,typename NumberType>
1429 *   Solid<dim,NumberType>::Solid(const Parameters::AllParameters &parameters)
1430 *   :
1431 *   parameters(parameters),
1432 *   vol_reference (0.0),
1433 *   vol_current (0.0),
1435 *   time(parameters.end_time, parameters.delta_t),
1436 *   timer(std::cout,
1439 *   degree(parameters.poly_degree),
1440 * @endcode
1441 *
1442 * The Finite Element System is composed of dim continuous displacement
1443 * DOFs.
1444 *
1445 * @code
1446 *   fe(FE_Q<dim>(parameters.poly_degree), dim), // displacement
1447 *   dof_handler_ref(triangulation),
1448 *   dofs_per_cell (fe.dofs_per_cell),
1449 *   u_fe(first_u_component),
1450 *   dofs_per_block(n_blocks),
1451 *   qf_cell(parameters.quad_order),
1452 *   qf_face(parameters.quad_order),
1453 *   n_q_points (qf_cell.size()),
1454 *   n_q_points_f (qf_face.size())
1455 *   {
1456 *  
1457 *   }
1458 *  
1459 * @endcode
1460 *
1461 * The class destructor simply clears the data held by the DOFHandler
1462 *
1463 * @code
1464 *   template <int dim,typename NumberType>
1465 *   Solid<dim,NumberType>::~Solid()
1466 *   {
1467 *   dof_handler_ref.clear();
1468 *   }
1469 *  
1470 *  
1471 * @endcode
1472 *
1473 * In solving the quasi-static problem, the time becomes a loading parameter,
1474 * i.e. we increasing the loading linearly with time, making the two concepts
1475 * interchangeable. We choose to increment time linearly using a constant time
1476 * step size.
1477 *
1478
1479 *
1480 * We start the function with preprocessing, and then output the initial grid
1481 * before starting the simulation proper with the first time (and loading)
1482 * increment.
1483 *
1484
1485 *
1486 *
1487 * @code
1488 *   template <int dim,typename NumberType>
1489 *   void Solid<dim,NumberType>::run()
1490 *   {
1491 *   make_grid();
1492 *   system_setup();
1493 *   output_results();
1494 *   time.increment();
1495 *  
1496 * @endcode
1497 *
1498 * We then declare the incremental solution update @f$\varDelta
1499 * \mathbf{\Xi}:= \{\varDelta \mathbf{u}\}@f$ and start the loop over the
1500 * time domain.
1501 *
1502
1503 *
1504 * At the beginning, we reset the solution update for this time step...
1505 *
1506 * @code
1507 *   BlockVector<double> solution_delta(dofs_per_block);
1508 *   while (time.current() <= time.end())
1509 *   {
1510 *   solution_delta = 0.0;
1511 *  
1512 * @endcode
1513 *
1514 * ...solve the current time step and update total solution vector
1515 * @f$\mathbf{\Xi}_{\textrm{n}} = \mathbf{\Xi}_{\textrm{n-1}} +
1516 * \varDelta \mathbf{\Xi}@f$...
1517 *
1518 * @code
1519 *   solve_nonlinear_timestep(solution_delta);
1520 *   solution_n += solution_delta;
1521 *  
1522 * @endcode
1523 *
1524 * ...and plot the results before moving on happily to the next time
1525 * step:
1526 *
1527 * @code
1528 *   output_results();
1529 *   time.increment();
1530 *   }
1531 *  
1532 * @endcode
1533 *
1534 * Lastly, we print the vertical tip displacement of the Cook cantilever
1535 * after the full load is applied
1536 *
1537 * @code
1538 *   print_vertical_tip_displacement();
1539 *   }
1540 *  
1541 *  
1542 * @endcode
1543 *
1544 *
1545 * <a name="cook_membrane.cc-Privateinterface"></a>
1546 * <h3>Private interface</h3>
1547 *
1548
1549 *
1550 *
1551 * <a name="cook_membrane.cc-Solidmake_grid"></a>
1552 * <h4>Solid::make_grid</h4>
1553 *
1554
1555 *
1556 * On to the first of the private member functions. Here we create the
1557 * triangulation of the domain, for which we choose a scaled an anisotripically
1558 * discretised rectangle which is subsequently transformed into the correct
1559 * of the Cook cantilever. Each relevant boundary face is then given a boundary
1560 * ID number.
1561 *
1562
1563 *
1564 * We then determine the volume of the reference configuration and print it
1565 * for comparison.
1566 *
1567
1568 *
1569 *
1570 * @code
1571 *   template <int dim>
1572 *   Point<dim> grid_y_transform (const Point<dim> &pt_in)
1573 *   {
1574 *   const double &x = pt_in[0];
1575 *   const double &y = pt_in[1];
1576 *  
1577 *   const double y_upper = 44.0 + (16.0/48.0)*x; // Line defining upper edge of beam
1578 *   const double y_lower = 0.0 + (44.0/48.0)*x; // Line defining lower edge of beam
1579 *   const double theta = y/44.0; // Fraction of height along left side of beam
1580 *   const double y_transform = (1-theta)*y_lower + theta*y_upper; // Final transformation
1581 *  
1582 *   Point<dim> pt_out = pt_in;
1583 *   pt_out[1] = y_transform;
1584 *  
1585 *   return pt_out;
1586 *   }
1587 *  
1588 *   template <int dim,typename NumberType>
1589 *   void Solid<dim,NumberType>::make_grid()
1590 *   {
1591 * @endcode
1592 *
1593 * Divide the beam, but only along the x- and y-coordinate directions
1594 *
1595 * @code
1596 *   std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge);
1597 * @endcode
1598 *
1599 * Only allow one element through the thickness
1600 * (modelling a plane strain condition)
1601 *
1602 * @code
1603 *   if (dim == 3)
1604 *   repetitions[dim-1] = 1;
1605 *  
1606 *   const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, -0.5) : Point<dim>(0.0, 0.0));
1607 *   const Point<dim> top_right = (dim == 3 ? Point<dim>(48.0, 44.0, 0.5) : Point<dim>(48.0, 44.0));
1608 *  
1610 *   repetitions,
1611 *   bottom_left,
1612 *   top_right);
1613 *  
1614 * @endcode
1615 *
1616 * Since we wish to apply a Neumann BC to the right-hand surface, we
1617 * must find the cell faces in this part of the domain and mark them with
1618 * a distinct boundary ID number. The faces we are looking for are on the
1619 * +x surface and will get boundary ID 11.
1620 * Dirichlet boundaries exist on the left-hand face of the beam (this fixed
1621 * boundary will get ID 1) and on the +Z and -Z faces (which correspond to
1622 * ID 2 and we will use to impose the plane strain condition)
1623 *
1624 * @code
1625 *   const double tol_boundary = 1e-6;
1626 *   typename Triangulation<dim>::active_cell_iterator cell =
1627 *   triangulation.begin_active(), endc = triangulation.end();
1628 *   for (; cell != endc; ++cell)
1629 *   for (unsigned int face = 0;
1630 *   face < GeometryInfo<dim>::faces_per_cell; ++face)
1631 *   if (cell->face(face)->at_boundary() == true)
1632 *   {
1633 *   if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary)
1634 *   cell->face(face)->set_boundary_id(1); // -X faces
1635 *   else if (std::abs(cell->face(face)->center()[0] - 48.0) < tol_boundary)
1636 *   cell->face(face)->set_boundary_id(11); // +X faces
1637 *   else if (dim == 3 && std::abs(std::abs(cell->face(face)->center()[2]) - 0.5) < tol_boundary)
1638 *   cell->face(face)->set_boundary_id(2); // +Z and -Z faces
1639 *   }
1640 *  
1641 * @endcode
1642 *
1643 * Transform the hyper-rectangle into the beam shape
1644 *
1645 * @code
1646 *   GridTools::transform(&grid_y_transform<dim>, triangulation);
1647 *  
1648 *   GridTools::scale(parameters.scale, triangulation);
1649 *  
1650 *   vol_reference = GridTools::volume(triangulation);
1651 *   vol_current = vol_reference;
1652 *   std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
1653 *   }
1654 *  
1655 *  
1656 * @endcode
1657 *
1658 *
1659 * <a name="cook_membrane.cc-Solidsystem_setup"></a>
1660 * <h4>Solid::system_setup</h4>
1661 *
1662
1663 *
1664 * Next we describe how the FE system is setup. We first determine the number
1665 * of components per block. Since the displacement is a vector component, the
1666 * first dim components belong to it.
1667 *
1668 * @code
1669 *   template <int dim,typename NumberType>
1670 *   void Solid<dim,NumberType>::system_setup()
1671 *   {
1672 *   timer.enter_subsection("Setup system");
1673 *  
1674 *   std::vector<unsigned int> block_component(n_components, u_dof); // Displacement
1675 *  
1676 * @endcode
1677 *
1678 * The DOF handler is then initialised and we renumber the grid in an
1679 * efficient manner. We also record the number of DOFs per block.
1680 *
1681 * @code
1682 *   dof_handler_ref.distribute_dofs(fe);
1683 *   DoFRenumbering::Cuthill_McKee(dof_handler_ref);
1684 *   DoFRenumbering::component_wise(dof_handler_ref, block_component);
1685 *   dofs_per_block = DoFTools::count_dofs_per_fe_block(dof_handler_ref, block_component);
1686 *  
1687 *   std::cout << "Triangulation:"
1688 *   << "\n\t Number of active cells: " << triangulation.n_active_cells()
1689 *   << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs()
1690 *   << std::endl;
1691 *  
1692 * @endcode
1693 *
1694 * Setup the sparsity pattern and tangent matrix
1695 *
1696 * @code
1697 *   tangent_matrix.clear();
1698 *   {
1699 *   const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
1700 *  
1701 *   BlockDynamicSparsityPattern csp(n_blocks, n_blocks);
1702 *  
1703 *   csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
1704 *   csp.collect_sizes();
1705 *  
1706 * @endcode
1707 *
1708 * Naturally, for a one-field vector-valued problem, all of the
1709 * components of the system are coupled.
1710 *
1711 * @code
1712 *   Table<2, DoFTools::Coupling> coupling(n_components, n_components);
1713 *   for (unsigned int ii = 0; ii < n_components; ++ii)
1714 *   for (unsigned int jj = 0; jj < n_components; ++jj)
1715 *   coupling[ii][jj] = DoFTools::always;
1716 *   DoFTools::make_sparsity_pattern(dof_handler_ref,
1717 *   coupling,
1718 *   csp,
1719 *   constraints,
1720 *   false);
1721 *   sparsity_pattern.copy_from(csp);
1722 *   }
1723 *  
1724 *   tangent_matrix.reinit(sparsity_pattern);
1725 *  
1726 * @endcode
1727 *
1728 * We then set up storage vectors
1729 *
1730 * @code
1731 *   system_rhs.reinit(dofs_per_block);
1732 *   system_rhs.collect_sizes();
1733 *  
1734 *   solution_n.reinit(dofs_per_block);
1735 *   solution_n.collect_sizes();
1736 *  
1737 * @endcode
1738 *
1739 * ...and finally set up the quadrature
1740 * point history:
1741 *
1742 * @code
1743 *   setup_qph();
1744 *  
1745 *   timer.leave_subsection();
1746 *   }
1747 *  
1748 *  
1749 * @endcode
1750 *
1751 *
1752 * <a name="cook_membrane.cc-Solidsetup_qph"></a>
1753 * <h4>Solid::setup_qph</h4>
1754 * The method used to store quadrature information is already described in
1755 * @ref step_18 "step-18" and @ref step_44 "step-44". Here we implement a similar setup for a SMP machine.
1756 *
1757
1758 *
1759 * Firstly the actual QPH data objects are created. This must be done only
1760 * once the grid is refined to its finest level.
1761 *
1762 * @code
1763 *   template <int dim,typename NumberType>
1764 *   void Solid<dim,NumberType>::setup_qph()
1765 *   {
1766 *   std::cout << " Setting up quadrature point data..." << std::endl;
1767 *  
1768 *   quadrature_point_history.initialize(triangulation.begin_active(),
1769 *   triangulation.end(),
1770 *   n_q_points);
1771 *  
1772 * @endcode
1773 *
1774 * Next we setup the initial quadrature point data. Note that when
1775 * the quadrature point data is retrieved, it is returned as a vector
1776 * of smart pointers.
1777 *
1778 * @code
1779 *   for (typename Triangulation<dim>::active_cell_iterator cell =
1780 *   triangulation.begin_active(); cell != triangulation.end(); ++cell)
1781 *   {
1782 *   const std::vector<std::shared_ptr<PointHistory<dim,NumberType> > > lqph =
1783 *   quadrature_point_history.get_data(cell);
1784 *   Assert(lqph.size() == n_q_points, ExcInternalError());
1785 *  
1786 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1787 *   lqph[q_point]->setup_lqp(parameters);
1788 *   }
1789 *   }
1790 *  
1791 *  
1792 * @endcode
1793 *
1794 *
1795 * <a name="cook_membrane.cc-Solidsolve_nonlinear_timestep"></a>
1796 * <h4>Solid::solve_nonlinear_timestep</h4>
1797 *
1798
1799 *
1800 * The next function is the driver method for the Newton-Raphson scheme. At
1801 * its top we create a new vector to store the current Newton update step,
1802 * reset the error storage objects and print solver header.
1803 *
1804 * @code
1805 *   template <int dim,typename NumberType>
1806 *   void
1807 *   Solid<dim,NumberType>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
1808 *   {
1809 *   std::cout << std::endl << "Timestep " << time.get_timestep() << " @ "
1810 *   << time.current() << "s" << std::endl;
1811 *  
1812 *   BlockVector<double> newton_update(dofs_per_block);
1813 *  
1814 *   error_residual.reset();
1815 *   error_residual_0.reset();
1816 *   error_residual_norm.reset();
1817 *   error_update.reset();
1818 *   error_update_0.reset();
1819 *   error_update_norm.reset();
1820 *  
1821 *   print_conv_header();
1822 *  
1823 * @endcode
1824 *
1825 * We now perform a number of Newton iterations to iteratively solve the
1826 * nonlinear problem. Since the problem is fully nonlinear and we are
1827 * using a full Newton method, the data stored in the tangent matrix and
1828 * right-hand side vector is not reusable and must be cleared at each
1829 * Newton step. We then initially build the right-hand side vector to
1830 * check for convergence (and store this value in the first iteration).
1831 * The unconstrained DOFs of the rhs vector hold the out-of-balance
1832 * forces. The building is done before assembling the system matrix as the
1833 * latter is an expensive operation and we can potentially avoid an extra
1834 * assembly process by not assembling the tangent matrix when convergence
1835 * is attained.
1836 *
1837 * @code
1838 *   unsigned int newton_iteration = 0;
1839 *   for (; newton_iteration < parameters.max_iterations_NR;
1840 *   ++newton_iteration)
1841 *   {
1842 *   std::cout << " " << std::setw(2) << newton_iteration << " " << std::flush;
1843 *  
1844 * @endcode
1845 *
1846 * If we have decided that we want to continue with the iteration, we
1847 * assemble the tangent, make and impose the Dirichlet constraints,
1848 * and do the solve of the linearized system:
1849 *
1850 * @code
1851 *   make_constraints(newton_iteration);
1852 *   assemble_system(solution_delta);
1853 *  
1854 *   get_error_residual(error_residual);
1855 *  
1856 *   if (newton_iteration == 0)
1857 *   error_residual_0 = error_residual;
1858 *  
1859 * @endcode
1860 *
1861 * We can now determine the normalised residual error and check for
1862 * solution convergence:
1863 *
1864 * @code
1865 *   error_residual_norm = error_residual;
1866 *   error_residual_norm.normalise(error_residual_0);
1867 *  
1868 *   if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u
1869 *   && error_residual_norm.u <= parameters.tol_f)
1870 *   {
1871 *   std::cout << " CONVERGED! " << std::endl;
1872 *   print_conv_footer();
1873 *  
1874 *   break;
1875 *   }
1876 *  
1877 *   const std::pair<unsigned int, double>
1878 *   lin_solver_output = solve_linear_system(newton_update);
1879 *  
1880 *   get_error_update(newton_update, error_update);
1881 *   if (newton_iteration == 0)
1882 *   error_update_0 = error_update;
1883 *  
1884 * @endcode
1885 *
1886 * We can now determine the normalised Newton update error, and
1887 * perform the actual update of the solution increment for the current
1888 * time step, update all quadrature point information pertaining to
1889 * this new displacement and stress state and continue iterating:
1890 *
1891 * @code
1892 *   error_update_norm = error_update;
1893 *   error_update_norm.normalise(error_update_0);
1894 *  
1895 *   solution_delta += newton_update;
1896 *  
1897 *   std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
1898 *   << std::scientific << lin_solver_output.first << " "
1899 *   << lin_solver_output.second << " " << error_residual_norm.norm
1900 *   << " " << error_residual_norm.u << " "
1901 *   << " " << error_update_norm.norm << " " << error_update_norm.u
1902 *   << " " << std::endl;
1903 *   }
1904 *  
1905 * @endcode
1906 *
1907 * At the end, if it turns out that we have in fact done more iterations
1908 * than the parameter file allowed, we raise an exception that can be
1909 * caught in the main() function. The call <code>AssertThrow(condition,
1910 * exc_object)</code> is in essence equivalent to <code>if (!cond) throw
1911 * exc_object;</code> but the former form fills certain fields in the
1912 * exception object that identify the location (filename and line number)
1913 * where the exception was raised to make it simpler to identify where the
1914 * problem happened.
1915 *
1916 * @code
1917 *   AssertThrow (newton_iteration <= parameters.max_iterations_NR,
1918 *   ExcMessage("No convergence in nonlinear solver!"));
1919 *   }
1920 *  
1921 *  
1922 * @endcode
1923 *
1924 *
1925 * <a name="cook_membrane.cc-Solidprint_conv_headerSolidprint_conv_footerandSolidprint_vertical_tip_displacement"></a>
1926 * <h4>Solid::print_conv_header, Solid::print_conv_footer and Solid::print_vertical_tip_displacement</h4>
1927 *
1928
1929 *
1930 * This program prints out data in a nice table that is updated
1931 * on a per-iteration basis. The next two functions set up the table
1932 * header and footer:
1933 *
1934 * @code
1935 *   template <int dim,typename NumberType>
1936 *   void Solid<dim,NumberType>::print_conv_header()
1937 *   {
1938 *   static const unsigned int l_width = 87;
1939 *  
1940 *   for (unsigned int i = 0; i < l_width; ++i)
1941 *   std::cout << "_";
1942 *   std::cout << std::endl;
1943 *  
1944 *   std::cout << " SOLVER STEP "
1945 *   << " | LIN_IT LIN_RES RES_NORM "
1946 *   << " RES_U NU_NORM "
1947 *   << " NU_U " << std::endl;
1948 *  
1949 *   for (unsigned int i = 0; i < l_width; ++i)
1950 *   std::cout << "_";
1951 *   std::cout << std::endl;
1952 *   }
1953 *  
1954 *  
1955 *  
1956 *   template <int dim,typename NumberType>
1957 *   void Solid<dim,NumberType>::print_conv_footer()
1958 *   {
1959 *   static const unsigned int l_width = 87;
1960 *  
1961 *   for (unsigned int i = 0; i < l_width; ++i)
1962 *   std::cout << "_";
1963 *   std::cout << std::endl;
1964 *  
1965 *   std::cout << "Relative errors:" << std::endl
1966 *   << "Displacement:\t" << error_update.u / error_update_0.u << std::endl
1967 *   << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
1968 *   << "v / V_0:\t" << vol_current << " / " << vol_reference
1969 *   << std::endl;
1970 *   }
1971 *  
1972 * @endcode
1973 *
1974 * At the end we also output the result that can be compared to that found in
1975 * the literature, namely the displacement at the upper right corner of the
1976 * beam.
1977 *
1978 * @code
1979 *   template <int dim,typename NumberType>
1980 *   void Solid<dim,NumberType>::print_vertical_tip_displacement()
1981 *   {
1982 *   static const unsigned int l_width = 87;
1983 *  
1984 *   for (unsigned int i = 0; i < l_width; ++i)
1985 *   std::cout << "_";
1986 *   std::cout << std::endl;
1987 *  
1988 * @endcode
1989 *
1990 * The measurement point, as stated in the reference paper, is at the midway
1991 * point of the surface on which the traction is applied.
1992 *
1993 * @code
1994 *   const Point<dim> soln_pt = (dim == 3 ?
1995 *   Point<dim>(48.0*parameters.scale, 52.0*parameters.scale, 0.5*parameters.scale) :
1996 *   Point<dim>(48.0*parameters.scale, 52.0*parameters.scale));
1997 *   double vertical_tip_displacement = 0.0;
1998 *   double vertical_tip_displacement_check = 0.0;
1999 *  
2000 *   typename DoFHandler<dim>::active_cell_iterator cell =
2001 *   dof_handler_ref.begin_active(), endc = dof_handler_ref.end();
2002 *   for (; cell != endc; ++cell)
2003 *   {
2004 * @endcode
2005 *
2006 * if (cell->point_inside(soln_pt) == true)
2007 *
2008 * @code
2009 *   for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
2010 *   if (cell->vertex(v).distance(soln_pt) < 1e-6)
2011 *   {
2012 * @endcode
2013 *
2014 * Extract y-component of solution at the given point
2015 * This point is coindicent with a vertex, so we can
2016 * extract it directly as we're using FE_Q finite elements
2017 * that have support at the vertices
2018 *
2019 * @code
2020 *   vertical_tip_displacement = solution_n(cell->vertex_dof_index(v,u_dof+1));
2021 *  
2022 * @endcode
2023 *
2024 * Sanity check using alternate method to extract the solution
2025 * at the given point. To do this, we must create an FEValues instance
2026 * to help us extract the solution value at the desired point
2027 *
2028 * @code
2029 *   const MappingQ<dim> mapping (parameters.poly_degree);
2030 *   const Point<dim> qp_unit = mapping.transform_real_to_unit_cell(cell,soln_pt);
2031 *   const Quadrature<dim> soln_qrule (qp_unit);
2032 *   AssertThrow(soln_qrule.size() == 1, ExcInternalError());
2033 *   FEValues<dim> fe_values_soln (fe, soln_qrule, update_values);
2034 *   fe_values_soln.reinit(cell);
2035 *  
2036 * @endcode
2037 *
2038 * Extract y-component of solution at given point
2039 *
2040 * @code
2041 *   std::vector< Tensor<1,dim> > soln_values (soln_qrule.size());
2042 *   fe_values_soln[u_fe].get_function_values(solution_n,
2043 *   soln_values);
2044 *   vertical_tip_displacement_check = soln_values[0][u_dof+1];
2045 *  
2046 *   break;
2047 *   }
2048 *   }
2049 *   AssertThrow(vertical_tip_displacement > 0.0, ExcMessage("Found no cell with point inside!"));
2050 *  
2051 *   std::cout << "Vertical tip displacement: " << vertical_tip_displacement
2052 *   << "\t Check: " << vertical_tip_displacement_check
2053 *   << std::endl;
2054 *   }
2055 *  
2056 *  
2057 * @endcode
2058 *
2059 *
2060 * <a name="cook_membrane.cc-Solidget_error_residual"></a>
2061 * <h4>Solid::get_error_residual</h4>
2062 *
2063
2064 *
2065 * Determine the true residual error for the problem. That is, determine the
2066 * error in the residual for the unconstrained degrees of freedom. Note that to
2067 * do so, we need to ignore constrained DOFs by setting the residual in these
2068 * vector components to zero.
2069 *
2070 * @code
2071 *   template <int dim,typename NumberType>
2072 *   void Solid<dim,NumberType>::get_error_residual(Errors &error_residual)
2073 *   {
2074 *   BlockVector<double> error_res(dofs_per_block);
2075 *  
2076 *   for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
2077 *   if (!constraints.is_constrained(i))
2078 *   error_res(i) = system_rhs(i);
2079 *  
2080 *   error_residual.norm = error_res.l2_norm();
2081 *   error_residual.u = error_res.block(u_dof).l2_norm();
2082 *   }
2083 *  
2084 *  
2085 * @endcode
2086 *
2087 *
2088 * <a name="cook_membrane.cc-Solidget_error_udpate"></a>
2089 * <h4>Solid::get_error_udpate</h4>
2090 *
2091
2092 *
2093 * Determine the true Newton update error for the problem
2094 *
2095 * @code
2096 *   template <int dim,typename NumberType>
2097 *   void Solid<dim,NumberType>::get_error_update(const BlockVector<double> &newton_update,
2098 *   Errors &error_update)
2099 *   {
2100 *   BlockVector<double> error_ud(dofs_per_block);
2101 *   for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
2102 *   if (!constraints.is_constrained(i))
2103 *   error_ud(i) = newton_update(i);
2104 *  
2105 *   error_update.norm = error_ud.l2_norm();
2106 *   error_update.u = error_ud.block(u_dof).l2_norm();
2107 *   }
2108 *  
2109 *  
2110 *  
2111 * @endcode
2112 *
2113 *
2114 * <a name="cook_membrane.cc-Solidget_total_solution"></a>
2115 * <h4>Solid::get_total_solution</h4>
2116 *
2117
2118 *
2119 * This function provides the total solution, which is valid at any Newton step.
2120 * This is required as, to reduce computational error, the total solution is
2121 * only updated at the end of the timestep.
2122 *
2123 * @code
2124 *   template <int dim,typename NumberType>
2125 *   BlockVector<double>
2126 *   Solid<dim,NumberType>::get_total_solution(const BlockVector<double> &solution_delta) const
2127 *   {
2128 *   BlockVector<double> solution_total(solution_n);
2129 *   solution_total += solution_delta;
2130 *   return solution_total;
2131 *   }
2132 *  
2133 *  
2134 * @endcode
2135 *
2136 *
2137 * <a name="cook_membrane.cc-Solidassemble_system"></a>
2138 * <h4>Solid::assemble_system</h4>
2139 *
2140
2141 *
2142 *
2143 * @code
2144 *   template <int dim,typename NumberType>
2145 *   struct Assembler_Base
2146 *   {
2147 *   virtual ~Assembler_Base() {}
2148 *  
2149 * @endcode
2150 *
2151 * Here we deal with the tangent matrix assembly structures. The
2152 * PerTaskData object stores local contributions.
2153 *
2154 * @code
2155 *   struct PerTaskData_ASM
2156 *   {
2157 *   const Solid<dim,NumberType> *solid;
2158 *   FullMatrix<double> cell_matrix;
2159 *   Vector<double> cell_rhs;
2160 *   std::vector<types::global_dof_index> local_dof_indices;
2161 *  
2162 *   PerTaskData_ASM(const Solid<dim,NumberType> *solid)
2163 *   :
2164 *   solid (solid),
2165 *   cell_matrix(solid->dofs_per_cell, solid->dofs_per_cell),
2166 *   cell_rhs(solid->dofs_per_cell),
2167 *   local_dof_indices(solid->dofs_per_cell)
2168 *   {}
2169 *  
2170 *   void reset()
2171 *   {
2172 *   cell_matrix = 0.0;
2173 *   cell_rhs = 0.0;
2174 *   }
2175 *   };
2176 *  
2177 * @endcode
2178 *
2179 * On the other hand, the ScratchData object stores the larger objects such as
2180 * the shape-function values array (<code>Nx</code>) and a shape function
2181 * gradient and symmetric gradient vector which we will use during the
2182 * assembly.
2183 *
2184 * @code
2185 *   struct ScratchData_ASM
2186 *   {
2187 *   const BlockVector<double> &solution_total;
2188 *   std::vector<Tensor<2, dim,NumberType> > solution_grads_u_total;
2189 *  
2190 *   FEValues<dim> fe_values_ref;
2191 *   FEFaceValues<dim> fe_face_values_ref;
2192 *  
2193 *   std::vector<std::vector<Tensor<2, dim,NumberType> > > grad_Nx;
2194 *   std::vector<std::vector<SymmetricTensor<2,dim,NumberType> > > symm_grad_Nx;
2195 *  
2196 *   ScratchData_ASM(const FiniteElement<dim> &fe_cell,
2197 *   const QGauss<dim> &qf_cell,
2198 *   const UpdateFlags uf_cell,
2199 *   const QGauss<dim-1> & qf_face,
2200 *   const UpdateFlags uf_face,
2201 *   const BlockVector<double> &solution_total)
2202 *   :
2203 *   solution_total(solution_total),
2204 *   solution_grads_u_total(qf_cell.size()),
2205 *   fe_values_ref(fe_cell, qf_cell, uf_cell),
2206 *   fe_face_values_ref(fe_cell, qf_face, uf_face),
2207 *   grad_Nx(qf_cell.size(),
2208 *   std::vector<Tensor<2,dim,NumberType> >(fe_cell.dofs_per_cell)),
2209 *   symm_grad_Nx(qf_cell.size(),
2210 *   std::vector<SymmetricTensor<2,dim,NumberType> >
2211 *   (fe_cell.dofs_per_cell))
2212 *   {}
2213 *  
2214 *   ScratchData_ASM(const ScratchData_ASM &rhs)
2215 *   :
2216 *   solution_total (rhs.solution_total),
2217 *   solution_grads_u_total(rhs.solution_grads_u_total),
2218 *   fe_values_ref(rhs.fe_values_ref.get_fe(),
2219 *   rhs.fe_values_ref.get_quadrature(),
2220 *   rhs.fe_values_ref.get_update_flags()),
2221 *   fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
2222 *   rhs.fe_face_values_ref.get_quadrature(),
2223 *   rhs.fe_face_values_ref.get_update_flags()),
2224 *   grad_Nx(rhs.grad_Nx),
2225 *   symm_grad_Nx(rhs.symm_grad_Nx)
2226 *   {}
2227 *  
2228 *   void reset()
2229 *   {
2230 *   const unsigned int n_q_points = fe_values_ref.get_quadrature().size();
2231 *   const unsigned int n_dofs_per_cell = fe_values_ref.dofs_per_cell;
2232 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2233 *   {
2234 *   Assert( grad_Nx[q_point].size() == n_dofs_per_cell,
2235 *   ExcInternalError());
2236 *   Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell,
2237 *   ExcInternalError());
2238 *  
2239 *   solution_grads_u_total[q_point] = Tensor<2,dim,NumberType>();
2240 *   for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
2241 *   {
2242 *   grad_Nx[q_point][k] = Tensor<2,dim,NumberType>();
2243 *   symm_grad_Nx[q_point][k] = SymmetricTensor<2,dim,NumberType>();
2244 *   }
2245 *   }
2246 *   }
2247 *  
2248 *   };
2249 *  
2250 * @endcode
2251 *
2252 * Of course, we still have to define how we assemble the tangent matrix
2253 * contribution for a single cell.
2254 *
2255 * @code
2256 *   void
2257 *   assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2258 *   ScratchData_ASM &scratch,
2259 *   PerTaskData_ASM &data)
2260 *   {
2261 * @endcode
2262 *
2263 * Due to the C++ specialization rules, we need one more
2264 * level of indirection in order to define the assembly
2265 * routine for all different number. The next function call
2266 * is specialized for each NumberType, but to prevent having
2267 * to specialize the whole class along with it we have inlined
2268 * the definition of the other functions that are common to
2269 * all implementations.
2270 *
2271 * @code
2272 *   assemble_system_tangent_residual_one_cell(cell, scratch, data);
2273 *   assemble_neumann_contribution_one_cell(cell, scratch, data);
2274 *   }
2275 *  
2276 * @endcode
2277 *
2278 * This function adds the local contribution to the system matrix.
2279 *
2280 * @code
2281 *   void
2282 *   copy_local_to_global_ASM(const PerTaskData_ASM &data)
2283 *   {
2284 *   const AffineConstraints<double> &constraints = data.solid->constraints;
2285 *   BlockSparseMatrix<double> &tangent_matrix = const_cast<Solid<dim,NumberType> *>(data.solid)->tangent_matrix;
2286 *   BlockVector<double> &system_rhs = const_cast<Solid<dim,NumberType> *>(data.solid)->system_rhs;
2287 *  
2288 *   constraints.distribute_local_to_global(
2289 *   data.cell_matrix, data.cell_rhs,
2290 *   data.local_dof_indices,
2291 *   tangent_matrix, system_rhs);
2292 *   }
2293 *  
2294 *   protected:
2295 *  
2296 * @endcode
2297 *
2298 * This function needs to exist in the base class for
2299 * Workstream to work with a reference to the base class.
2300 *
2301 * @code
2302 *   virtual void
2303 *   assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &/*cell*/,
2304 *   ScratchData_ASM &/*scratch*/,
2305 *   PerTaskData_ASM &/*data*/)
2306 *   {
2307 *   AssertThrow(false, ExcPureFunctionCalled());
2308 *   }
2309 *  
2310 *   void
2311 *   assemble_neumann_contribution_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2312 *   ScratchData_ASM &scratch,
2313 *   PerTaskData_ASM &data)
2314 *   {
2315 * @endcode
2316 *
2317 * Aliases for data referenced from the Solid class
2318 *
2319 * @code
2320 *   const unsigned int &n_q_points_f = data.solid->n_q_points_f;
2321 *   const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2322 *   const Parameters::AllParameters &parameters = data.solid->parameters;
2323 *   const Time &time = data.solid->time;
2324 *   const FESystem<dim> &fe = data.solid->fe;
2325 *   const unsigned int &u_dof = data.solid->u_dof;
2326 *  
2327 * @endcode
2328 *
2329 * Next we assemble the Neumann contribution. We first check to see it the
2330 * cell face exists on a boundary on which a traction is applied and add
2331 * the contribution if this is the case.
2332 *
2333 * @code
2334 *   for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2335 *   ++face)
2336 *   if (cell->face(face)->at_boundary() == true
2337 *   && cell->face(face)->boundary_id() == 11)
2338 *   {
2339 *   scratch.fe_face_values_ref.reinit(cell, face);
2340 *  
2341 *   for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
2342 *   ++f_q_point)
2343 *   {
2344 * @endcode
2345 *
2346 * We specify the traction in reference configuration.
2347 * For this problem, a defined total vertical force is applied
2348 * in the reference configuration.
2349 * The direction of the applied traction is assumed not to
2350 * evolve with the deformation of the domain.
2351 *
2352
2353 *
2354 * Note that the contributions to the right hand side vector we
2355 * compute here only exist in the displacement components of the
2356 * vector.
2357 *
2358 * @code
2359 *   const double time_ramp = (time.current() / time.end());
2360 *   const double magnitude = (1.0/(16.0*parameters.scale*1.0*parameters.scale))*time_ramp; // (Total force) / (RHS surface area)
2361 *   Tensor<1,dim> dir;
2362 *   dir[1] = 1.0;
2363 *   const Tensor<1, dim> traction = magnitude*dir;
2364 *  
2365 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2366 *   {
2367 *   const unsigned int i_group =
2368 *   fe.system_to_base_index(i).first.first;
2369 *  
2370 *   if (i_group == u_dof)
2371 *   {
2372 *   const unsigned int component_i =
2373 *   fe.system_to_component_index(i).first;
2374 *   const double Ni =
2375 *   scratch.fe_face_values_ref.shape_value(i,
2376 *   f_q_point);
2377 *   const double JxW = scratch.fe_face_values_ref.JxW(
2378 *   f_q_point);
2379 *  
2380 *   data.cell_rhs(i) += (Ni * traction[component_i])
2381 *   * JxW;
2382 *   }
2383 *   }
2384 *   }
2385 *   }
2386 *   }
2387 *  
2388 *   };
2389 *  
2390 *   template <int dim>
2391 *   struct Assembler<dim,double> : Assembler_Base<dim,double>
2392 *   {
2393 *   typedef double NumberType;
2394 *   using typename Assembler_Base<dim,NumberType>::ScratchData_ASM;
2395 *   using typename Assembler_Base<dim,NumberType>::PerTaskData_ASM;
2396 *  
2397 *   virtual ~Assembler() {}
2398 *  
2399 *   virtual void
2400 *   assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2401 *   ScratchData_ASM &scratch,
2402 *   PerTaskData_ASM &data) override
2403 *   {
2404 * @endcode
2405 *
2406 * Aliases for data referenced from the Solid class
2407 *
2408 * @code
2409 *   const unsigned int &n_q_points = data.solid->n_q_points;
2410 *   const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2411 *   const FESystem<dim> &fe = data.solid->fe;
2412 *   const unsigned int &u_dof = data.solid->u_dof;
2413 *   const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
2414 *  
2415 *   data.reset();
2416 *   scratch.reset();
2417 *   scratch.fe_values_ref.reinit(cell);
2418 *   cell->get_dof_indices(data.local_dof_indices);
2419 *  
2420 *   const std::vector<std::shared_ptr<const PointHistory<dim,NumberType> > > lqph =
2421 *   const_cast<const Solid<dim,NumberType> *>(data.solid)->quadrature_point_history.get_data(cell);
2422 *   Assert(lqph.size() == n_q_points, ExcInternalError());
2423 *  
2424 * @endcode
2425 *
2426 * We first need to find the solution gradients at quadrature points
2427 * inside the current cell and then we update each local QP using the
2428 * displacement gradient:
2429 *
2430 * @code
2431 *   scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
2432 *   scratch.solution_grads_u_total);
2433 *  
2434 * @endcode
2435 *
2436 * Now we build the local cell stiffness matrix. Since the global and
2437 * local system matrices are symmetric, we can exploit this property by
2438 * building only the lower half of the local matrix and copying the values
2439 * to the upper half.
2440 *
2441
2442 *
2443 * In doing so, we first extract some configuration dependent variables
2444 * from our QPH history objects for the current quadrature point.
2445 *
2446 * @code
2447 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2448 *   {
2449 *   const Tensor<2,dim,NumberType> &grad_u = scratch.solution_grads_u_total[q_point];
2450 *   const Tensor<2,dim,NumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
2451 *   const NumberType det_F = determinant(F);
2452 *   const Tensor<2,dim,NumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
2453 *   const SymmetricTensor<2,dim,NumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
2454 *   const Tensor<2,dim,NumberType> F_inv = invert(F);
2455 *   Assert(det_F > NumberType(0.0), ExcInternalError());
2456 *  
2457 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
2458 *   {
2459 *   const unsigned int k_group = fe.system_to_base_index(k).first.first;
2460 *  
2461 *   if (k_group == u_dof)
2462 *   {
2463 *   scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
2464 *   scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
2465 *   }
2466 *   else
2467 *   Assert(k_group <= u_dof, ExcInternalError());
2468 *   }
2469 *  
2470 *   const SymmetricTensor<2,dim,NumberType> tau = lqph[q_point]->get_tau(det_F,b_bar);
2471 *   const SymmetricTensor<4,dim,NumberType> Jc = lqph[q_point]->get_Jc(det_F,b_bar);
2472 *   const Tensor<2,dim,NumberType> tau_ns (tau);
2473 *  
2474 * @endcode
2475 *
2476 * Next we define some aliases to make the assembly process easier to
2477 * follow
2478 *
2479 * @code
2480 *   const std::vector<SymmetricTensor<2, dim> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
2481 *   const std::vector<Tensor<2, dim> > &grad_Nx = scratch.grad_Nx[q_point];
2482 *   const double JxW = scratch.fe_values_ref.JxW(q_point);
2483 *  
2484 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2485 *   {
2486 *   const unsigned int component_i = fe.system_to_component_index(i).first;
2487 *   const unsigned int i_group = fe.system_to_base_index(i).first.first;
2488 *  
2489 *   if (i_group == u_dof)
2490 *   data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
2491 *   else
2492 *   Assert(i_group <= u_dof, ExcInternalError());
2493 *  
2494 *   for (unsigned int j = 0; j <= i; ++j)
2495 *   {
2496 *   const unsigned int component_j = fe.system_to_component_index(j).first;
2497 *   const unsigned int j_group = fe.system_to_base_index(j).first.first;
2498 *  
2499 * @endcode
2500 *
2501 * This is the @f$\mathsf{\mathbf{k}}_{\mathbf{u} \mathbf{u}}@f$
2502 * contribution. It comprises a material contribution, and a
2503 * geometrical stress contribution which is only added along
2504 * the local matrix diagonals:
2505 *
2506 * @code
2507 *   if ((i_group == j_group) && (i_group == u_dof))
2508 *   {
2509 *   data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
2510 *   * symm_grad_Nx[j] * JxW;
2511 *   if (component_i == component_j) // geometrical stress contribution
2512 *   data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau_ns
2513 *   * grad_Nx[j][component_j] * JxW;
2514 *   }
2515 *   else
2516 *   Assert((i_group <= u_dof) && (j_group <= u_dof),
2517 *   ExcInternalError());
2518 *   }
2519 *   }
2520 *   }
2521 *  
2522 *  
2523 * @endcode
2524 *
2525 * Finally, we need to copy the lower half of the local matrix into the
2526 * upper half:
2527 *
2528 * @code
2529 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2530 *   for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
2531 *   data.cell_matrix(i, j) = data.cell_matrix(j, i);
2532 *   }
2533 *  
2534 *   };
2535 *  
2536 *   #ifdef ENABLE_SACADO_FORMULATION
2537 *  
2538 *  
2539 *   template <int dim>
2540 *   struct Assembler<dim,Sacado::Fad::DFad<double> > : Assembler_Base<dim,Sacado::Fad::DFad<double> >
2541 *   {
2542 *   typedef Sacado::Fad::DFad<double> ADNumberType;
2543 *   using typename Assembler_Base<dim,ADNumberType>::ScratchData_ASM;
2544 *   using typename Assembler_Base<dim,ADNumberType>::PerTaskData_ASM;
2545 *  
2546 *   virtual ~Assembler() {}
2547 *  
2548 *   virtual void
2549 *   assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2550 *   ScratchData_ASM &scratch,
2551 *   PerTaskData_ASM &data) override
2552 *   {
2553 * @endcode
2554 *
2555 * Aliases for data referenced from the Solid class
2556 *
2557 * @code
2558 *   const unsigned int &n_q_points = data.solid->n_q_points;
2559 *   const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2560 *   const FESystem<dim> &fe = data.solid->fe;
2561 *   const unsigned int &u_dof = data.solid->u_dof;
2562 *   const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
2563 *  
2564 *   data.reset();
2565 *   scratch.reset();
2566 *   scratch.fe_values_ref.reinit(cell);
2567 *   cell->get_dof_indices(data.local_dof_indices);
2568 *  
2569 *   const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
2570 *   const_cast<const Solid<dim,ADNumberType> *>(data.solid)->quadrature_point_history.get_data(cell);
2571 *   Assert(lqph.size() == n_q_points, ExcInternalError());
2572 *  
2573 *   const unsigned int n_independent_variables = data.local_dof_indices.size();
2574 *   std::vector<double> local_dof_values(n_independent_variables);
2575 *   cell->get_dof_values(scratch.solution_total,
2576 *   local_dof_values.begin(),
2577 *   local_dof_values.end());
2578 *  
2579 * @endcode
2580 *
2581 * We now retrieve a set of degree-of-freedom values that
2582 * have the operations that are performed with them tracked.
2583 *
2584 * @code
2585 *   std::vector<ADNumberType> local_dof_values_ad (n_independent_variables);
2586 *   for (unsigned int i=0; i<n_independent_variables; ++i)
2587 *   local_dof_values_ad[i] = ADNumberType(n_independent_variables, i, local_dof_values[i]);
2588 *  
2589 * @endcode
2590 *
2591 * Compute all values, gradients etc. based on sensitive
2592 * AD degree-of-freedom values.
2593 *
2594 * @code
2595 *   scratch.fe_values_ref[u_fe].get_function_gradients_from_local_dof_values(
2596 *   local_dof_values_ad,
2597 *   scratch.solution_grads_u_total);
2598 *  
2599 * @endcode
2600 *
2601 * Accumulate the residual value for each degree of freedom.
2602 * Note: Its important that the vectors is initialised (zero'd) correctly.
2603 *
2604 * @code
2605 *   std::vector<ADNumberType> residual_ad (dofs_per_cell, ADNumberType(0.0));
2606 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2607 *   {
2608 *   const Tensor<2,dim,ADNumberType> &grad_u = scratch.solution_grads_u_total[q_point];
2610 *   const ADNumberType det_F = determinant(F);
2613 *   const Tensor<2,dim,ADNumberType> F_inv = invert(F);
2614 *   Assert(det_F > ADNumberType(0.0), ExcInternalError());
2615 *  
2616 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
2617 *   {
2618 *   const unsigned int k_group = fe.system_to_base_index(k).first.first;
2619 *  
2620 *   if (k_group == u_dof)
2621 *   {
2622 *   scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
2623 *   scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
2624 *   }
2625 *   else
2626 *   Assert(k_group <= u_dof, ExcInternalError());
2627 *   }
2628 *  
2629 *   const SymmetricTensor<2,dim,ADNumberType> tau = lqph[q_point]->get_tau(det_F,b_bar);
2630 *  
2631 * @endcode
2632 *
2633 * Next we define some position-dependent aliases, again to
2634 * make the assembly process easier to follow.
2635 *
2636 * @code
2637 *   const std::vector<SymmetricTensor<2, dim,ADNumberType> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
2638 *   const double JxW = scratch.fe_values_ref.JxW(q_point);
2639 *  
2640 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2641 *   {
2642 *   const unsigned int i_group = fe.system_to_base_index(i).first.first;
2643 *  
2644 *   if (i_group == u_dof)
2645 *   residual_ad[i] += (symm_grad_Nx[i] * tau) * JxW;
2646 *   else
2647 *   Assert(i_group <= u_dof, ExcInternalError());
2648 *   }
2649 *   }
2650 *  
2651 *   for (unsigned int I=0; I<n_independent_variables; ++I)
2652 *   {
2653 *   const ADNumberType &residual_I = residual_ad[I];
2654 *   data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
2655 *   for (unsigned int J=0; J<n_independent_variables; ++J)
2656 *   {
2657 * @endcode
2658 *
2659 * Compute the gradients of the residual entry [forward-mode]
2660 *
2661 * @code
2662 *   data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
2663 *   }
2664 *   }
2665 *   }
2666 *  
2667 *   };
2668 *  
2669 *  
2670 *   template <int dim>
2671 *   struct Assembler<dim,Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > > : Assembler_Base<dim,Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > >
2672 *   {
2673 *   typedef Sacado::Fad::DFad<double> ADDerivType;
2674 *   typedef Sacado::Rad::ADvar<ADDerivType> ADNumberType;
2675 *   using typename Assembler_Base<dim,ADNumberType>::ScratchData_ASM;
2676 *   using typename Assembler_Base<dim,ADNumberType>::PerTaskData_ASM;
2677 *  
2678 *   virtual ~Assembler() {}
2679 *  
2680 *   virtual void
2681 *   assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2682 *   ScratchData_ASM &scratch,
2683 *   PerTaskData_ASM &data) override
2684 *   {
2685 * @endcode
2686 *
2687 * Aliases for data referenced from the Solid class
2688 *
2689 * @code
2690 *   const unsigned int &n_q_points = data.solid->n_q_points;
2691 *   const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
2692 *  
2693 *   data.reset();
2694 *   scratch.reset();
2695 *   scratch.fe_values_ref.reinit(cell);
2696 *   cell->get_dof_indices(data.local_dof_indices);
2697 *  
2698 *   const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
2699 *   data.solid->quadrature_point_history.get_data(cell);
2700 *   Assert(lqph.size() == n_q_points, ExcInternalError());
2701 *  
2702 *   const unsigned int n_independent_variables = data.local_dof_indices.size();
2703 *   std::vector<double> local_dof_values(n_independent_variables);
2704 *   cell->get_dof_values(scratch.solution_total,
2705 *   local_dof_values.begin(),
2706 *   local_dof_values.end());
2707 *  
2708 * @endcode
2709 *
2710 * We now retrieve a set of degree-of-freedom values that
2711 * have the operations that are performed with them tracked.
2712 *
2713 * @code
2714 *   std::vector<ADNumberType> local_dof_values_ad (n_independent_variables);
2715 *   for (unsigned int i=0; i<n_independent_variables; ++i)
2716 *   local_dof_values_ad[i] = ADDerivType(n_independent_variables, i, local_dof_values[i]);
2717 *  
2718 * @endcode
2719 *
2720 * Compute all values, gradients etc. based on sensitive
2721 * AD degree-of-freedom values.
2722 *
2723 * @code
2724 *   scratch.fe_values_ref[u_fe].get_function_gradients_from_local_dof_values(
2725 *   local_dof_values_ad,
2726 *   scratch.solution_grads_u_total);
2727 *  
2728 * @endcode
2729 *
2730 * Next we compute the total potential energy of the element.
2731 * This is defined as follows:
2732 * Total energy = (internal - external) energies
2733 * Note: Its important that this value is initialised (zero'd) correctly.
2734 *
2735 * @code
2736 *   ADNumberType cell_energy_ad = ADNumberType(0.0);
2737 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2738 *   {
2739 *   const Tensor<2,dim,ADNumberType> &grad_u = scratch.solution_grads_u_total[q_point];
2740 *   const Tensor<2,dim,ADNumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
2741 *   const ADNumberType det_F = determinant(F);
2742 *   const Tensor<2,dim,ADNumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
2743 *   const SymmetricTensor<2,dim,ADNumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
2744 *   Assert(det_F > ADNumberType(0.0), ExcInternalError());
2745 *  
2746 * @endcode
2747 *
2748 * Next we define some position-dependent aliases, again to
2749 * make the assembly process easier to follow.
2750 *
2751 * @code
2752 *   const double JxW = scratch.fe_values_ref.JxW(q_point);
2753 *  
2754 *   const ADNumberType Psi = lqph[q_point]->get_Psi(det_F,b_bar);
2755 *  
2756 * @endcode
2757 *
2758 * We extract the configuration-dependent material energy
2759 * from our QPH history objects for the current quadrature point
2760 * and integrate its contribution to increment the total
2761 * cell energy.
2762 *
2763 * @code
2764 *   cell_energy_ad += Psi * JxW;
2765 *   }
2766 *  
2767 * @endcode
2768 *
2769 * Compute derivatives of reverse-mode AD variables
2770 *
2771 * @code
2772 *   ADNumberType::Gradcomp();
2773 *  
2774 *   for (unsigned int I=0; I<n_independent_variables; ++I)
2775 *   {
2776 * @endcode
2777 *
2778 * This computes the adjoint df/dX_{i} [reverse-mode]
2779 *
2780 * @code
2781 *   const ADDerivType residual_I = local_dof_values_ad[I].adj();
2782 *   data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
2783 *   for (unsigned int J=0; J<n_independent_variables; ++J)
2784 *   {
2785 * @endcode
2786 *
2787 * Compute the gradients of the residual entry [forward-mode]
2788 *
2789 * @code
2790 *   data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
2791 *   }
2792 *   }
2793 *   }
2794 *  
2795 *   };
2796 *  
2797 *  
2798 *   #endif
2799 *  
2800 *  
2801 * @endcode
2802 *
2803 * Since we use TBB for assembly, we simply setup a copy of the
2804 * data structures required for the process and pass them, along
2805 * with the memory addresses of the assembly functions to the
2806 * WorkStream object for processing. Note that we must ensure that
2807 * the matrix is reset before any assembly operations can occur.
2808 *
2809 * @code
2810 *   template <int dim,typename NumberType>
2811 *   void Solid<dim,NumberType>::assemble_system(const BlockVector<double> &solution_delta)
2812 *   {
2813 *   timer.enter_subsection("Assemble linear system");
2814 *   std::cout << " ASM " << std::flush;
2815 *  
2816 *   tangent_matrix = 0.0;
2817 *   system_rhs = 0.0;
2818 *  
2819 *   const UpdateFlags uf_cell(update_gradients |
2820 *   update_JxW_values);
2821 *   const UpdateFlags uf_face(update_values |
2822 *   update_JxW_values);
2823 *  
2824 *   const BlockVector<double> solution_total(get_total_solution(solution_delta));
2825 *   typename Assembler_Base<dim,NumberType>::PerTaskData_ASM per_task_data(this);
2826 *   typename Assembler_Base<dim,NumberType>::ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total);
2827 *   Assembler<dim,NumberType> assembler;
2828 *  
2829 *   WorkStream::run(dof_handler_ref.begin_active(),
2830 *   dof_handler_ref.end(),
2831 *   static_cast<Assembler_Base<dim,NumberType>&>(assembler),
2832 *   &Assembler_Base<dim,NumberType>::assemble_system_one_cell,
2833 *   &Assembler_Base<dim,NumberType>::copy_local_to_global_ASM,
2834 *   scratch_data,
2835 *   per_task_data);
2836 *  
2837 *   timer.leave_subsection();
2838 *   }
2839 *  
2840 *  
2841 * @endcode
2842 *
2843 *
2844 * <a name="cook_membrane.cc-Solidmake_constraints"></a>
2845 * <h4>Solid::make_constraints</h4>
2846 * The constraints for this problem are simple to describe.
2847 * However, since we are dealing with an iterative Newton method,
2848 * it should be noted that any displacement constraints should only
2849 * be specified at the zeroth iteration and subsequently no
2850 * additional contributions are to be made since the constraints
2851 * are already exactly satisfied.
2852 *
2853 * @code
2854 *   template <int dim,typename NumberType>
2855 *   void Solid<dim,NumberType>::make_constraints(const int &it_nr)
2856 *   {
2857 *   std::cout << " CST " << std::flush;
2858 *  
2859 * @endcode
2860 *
2861 * Since the constraints are different at different Newton iterations, we
2862 * need to clear the constraints matrix and completely rebuild
2863 * it. However, after the first iteration, the constraints remain the same
2864 * and we can simply skip the rebuilding step if we do not clear it.
2865 *
2866 * @code
2867 *   if (it_nr > 1)
2868 *   return;
2869 *   const bool apply_dirichlet_bc = (it_nr == 0);
2870 *  
2871 * @endcode
2872 *
2873 * The boundary conditions for the indentation problem are as follows: On
2874 * the -x face (ID = 1), we set up a zero-displacement condition, -y and +y traction
2875 * free boundary condition (don't need to take care); -z and +z faces (ID = 2) are
2876 * not allowed to move along z axis so that it is a plane strain problem.
2877 * Finally, as described earlier, +x face (ID = 11) has an the applied
2878 * distributed shear force (converted by total force per unit area) which
2879 * needs to be taken care as an inhomogeneous Newmann boundary condition.
2880 *
2881
2882 *
2883 * In the following, we will have to tell the function interpolation
2884 * boundary values which components of the solution vector should be
2885 * constrained (i.e., whether it's the x-, y-, z-displacements or
2886 * combinations thereof). This is done using ComponentMask objects (see
2887 * @ref GlossComponentMask) which we can get from the finite element if we
2888 * provide it with an extractor object for the component we wish to
2889 * select. To this end we first set up such extractor objects and later
2890 * use it when generating the relevant component masks:
2891 *
2892
2893 *
2894 *
2895 * @code
2896 *   if (apply_dirichlet_bc)
2897 *   {
2898 *   constraints.clear();
2899 *  
2900 * @endcode
2901 *
2902 * Fixed left hand side of the beam
2903 *
2904 * @code
2905 *   {
2906 *   const int boundary_id = 1;
2907 *   VectorTools::interpolate_boundary_values(dof_handler_ref,
2908 *   boundary_id,
2909 *   Functions::ZeroFunction<dim>(n_components),
2910 *   constraints,
2911 *   fe.component_mask(u_fe));
2912 *   }
2913 *  
2914 * @endcode
2915 *
2916 * Zero Z-displacement through thickness direction
2917 * This corresponds to a plane strain condition being imposed on the beam
2918 *
2919 * @code
2920 *   if (dim == 3)
2921 *   {
2922 *   const int boundary_id = 2;
2923 *   const FEValuesExtractors::Scalar z_displacement(2);
2924 *   VectorTools::interpolate_boundary_values(dof_handler_ref,
2925 *   boundary_id,
2926 *   Functions::ZeroFunction<dim>(n_components),
2927 *   constraints,
2928 *   fe.component_mask(z_displacement));
2929 *   }
2930 *   }
2931 *   else
2932 *   {
2933 *   if (constraints.has_inhomogeneities())
2934 *   {
2935 *   AffineConstraints<double> homogeneous_constraints(constraints);
2936 *   for (unsigned int dof = 0; dof != dof_handler_ref.n_dofs(); ++dof)
2937 *   if (homogeneous_constraints.is_inhomogeneously_constrained(dof))
2938 *   homogeneous_constraints.set_inhomogeneity(dof, 0.0);
2939 *   constraints.clear();
2940 *   constraints.copy_from(homogeneous_constraints);
2941 *   }
2942 *   }
2943 *  
2944 *   constraints.close();
2945 *   }
2946 *  
2947 * @endcode
2948 *
2949 *
2950 * <a name="cook_membrane.cc-Solidsolve_linear_system"></a>
2951 * <h4>Solid::solve_linear_system</h4>
2952 * As the system is composed of a single block, defining a solution scheme
2953 * for the linear problem is straight-forward.
2954 *
2955 * @code
2956 *   template <int dim,typename NumberType>
2957 *   std::pair<unsigned int, double>
2958 *   Solid<dim,NumberType>::solve_linear_system(BlockVector<double> &newton_update)
2959 *   {
2960 *   BlockVector<double> A(dofs_per_block);
2961 *   BlockVector<double> B(dofs_per_block);
2962 *  
2963 *   unsigned int lin_it = 0;
2964 *   double lin_res = 0.0;
2965 *  
2966 * @endcode
2967 *
2968 * We solve for the incremental displacement @f$d\mathbf{u}@f$.
2969 *
2970 * @code
2971 *   {
2972 *   timer.enter_subsection("Linear solver");
2973 *   std::cout << " SLV " << std::flush;
2974 *   if (parameters.type_lin == "CG")
2975 *   {
2976 *   const int solver_its = static_cast<unsigned int>(
2977 *   tangent_matrix.block(u_dof, u_dof).m()
2978 *   * parameters.max_iterations_lin);
2979 *   const double tol_sol = parameters.tol_lin
2980 *   * system_rhs.block(u_dof).l2_norm();
2981 *  
2982 *   SolverControl solver_control(solver_its, tol_sol);
2983 *  
2985 *   SolverCG<Vector<double> > solver_CG(solver_control, GVM);
2986 *  
2987 * @endcode
2988 *
2989 * We've chosen by default a SSOR preconditioner as it appears to
2990 * provide the fastest solver convergence characteristics for this
2991 * problem on a single-thread machine. However, for multicore
2992 * computing, the Jacobi preconditioner which is multithreaded may
2993 * converge quicker for larger linear systems.
2994 *
2995 * @code
2996 *   PreconditionSelector<SparseMatrix<double>, Vector<double> >
2997 *   preconditioner (parameters.preconditioner_type,
2998 *   parameters.preconditioner_relaxation);
2999 *   preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
3000 *  
3001 *   solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
3002 *   newton_update.block(u_dof),
3003 *   system_rhs.block(u_dof),
3004 *   preconditioner);
3005 *  
3006 *   lin_it = solver_control.last_step();
3007 *   lin_res = solver_control.last_value();
3008 *   }
3009 *   else if (parameters.type_lin == "Direct")
3010 *   {
3011 * @endcode
3012 *
3013 * Otherwise if the problem is small
3014 * enough, a direct solver can be
3015 * utilised.
3016 *
3017 * @code
3018 *   SparseDirectUMFPACK A_direct;
3019 *   A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
3020 *   A_direct.vmult(newton_update.block(u_dof), system_rhs.block(u_dof));
3021 *  
3022 *   lin_it = 1;
3023 *   lin_res = 0.0;
3024 *   }
3025 *   else
3026 *   Assert (false, ExcMessage("Linear solver type not implemented"));
3027 *  
3028 *   timer.leave_subsection();
3029 *   }
3030 *  
3031 * @endcode
3032 *
3033 * Now that we have the displacement update, distribute the constraints
3034 * back to the Newton update:
3035 *
3036 * @code
3037 *   constraints.distribute(newton_update);
3038 *  
3039 *   return std::make_pair(lin_it, lin_res);
3040 *   }
3041 *  
3042 * @endcode
3043 *
3044 *
3045 * <a name="cook_membrane.cc-Solidoutput_results"></a>
3046 * <h4>Solid::output_results</h4>
3047 * Here we present how the results are written to file to be viewed
3048 * using ParaView or Visit. The method is similar to that shown in the
3049 * tutorials so will not be discussed in detail.
3050 *
3051 * @code
3052 *   template <int dim,typename NumberType>
3053 *   void Solid<dim,NumberType>::output_results() const
3054 *   {
3055 *   DataOut<dim> data_out;
3056 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
3057 *   data_component_interpretation(dim,
3058 *   DataComponentInterpretation::component_is_part_of_vector);
3059 *  
3060 *   std::vector<std::string> solution_name(dim, "displacement");
3061 *  
3062 *   data_out.attach_dof_handler(dof_handler_ref);
3063 *   data_out.add_data_vector(solution_n,
3064 *   solution_name,
3065 *   DataOut<dim>::type_dof_data,
3066 *   data_component_interpretation);
3067 *  
3068 * @endcode
3069 *
3070 * Since we are dealing with a large deformation problem, it would be nice
3071 * to display the result on a displaced grid! The MappingQEulerian class
3072 * linked with the DataOut class provides an interface through which this
3073 * can be achieved without physically moving the grid points in the
3074 * Triangulation object ourselves. We first need to copy the solution to
3075 * a temporary vector and then create the Eulerian mapping. We also
3076 * specify the polynomial degree to the DataOut object in order to produce
3077 * a more refined output data set when higher order polynomials are used.
3078 *
3079 * @code
3080 *   Vector<double> soln(solution_n.size());
3081 *   for (unsigned int i = 0; i < soln.size(); ++i)
3082 *   soln(i) = solution_n(i);
3083 *   MappingQEulerian<dim> q_mapping(degree, dof_handler_ref, soln);
3084 *   data_out.build_patches(q_mapping, degree);
3085 *  
3086 *   std::ostringstream filename;
3087 *   filename << "solution-" << time.get_timestep() << ".vtk";
3088 *  
3089 *   std::ofstream output(filename.str().c_str());
3090 *   data_out.write_vtk(output);
3091 *   }
3092 *  
3093 *   }
3094 *  
3095 *  
3096 * @endcode
3097 *
3098 *
3099 * <a name="cook_membrane.cc-Mainfunction"></a>
3100 * <h3>Main function</h3>
3101 * Lastly we provide the main driver function which appears
3102 * no different to the other tutorials.
3103 *
3104 * @code
3105 *   int main (int argc, char *argv[])
3106 *   {
3107 *   using namespace dealii;
3108 *   using namespace Cook_Membrane;
3109 *  
3110 *   const unsigned int dim = 3;
3111 *   try
3112 *   {
3113 *   deallog.depth_console(0);
3114 *   Parameters::AllParameters parameters("parameters.prm");
3115 *   if (parameters.automatic_differentiation_order == 0)
3116 *   {
3117 *   std::cout << "Assembly method: Residual and linearisation are computed manually." << std::endl;
3118 *  
3119 * @endcode
3120 *
3121 * Allow multi-threading
3122 *
3123 * @code
3124 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3125 *   ::numbers::invalid_unsigned_int);
3126 *  
3127 *   typedef double NumberType;
3128 *   Solid<dim,NumberType> solid_3d(parameters);
3129 *   solid_3d.run();
3130 *   }
3131 *   #ifdef ENABLE_SACADO_FORMULATION
3132 *   else if (parameters.automatic_differentiation_order == 1)
3133 *   {
3134 *   std::cout << "Assembly method: Residual computed manually; linearisation performed using AD." << std::endl;
3135 *  
3136 * @endcode
3137 *
3138 * Allow multi-threading
3139 *
3140 * @code
3141 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3142 *   ::numbers::invalid_unsigned_int);
3143 *  
3144 *   typedef Sacado::Fad::DFad<double> NumberType;
3145 *   Solid<dim,NumberType> solid_3d(parameters);
3146 *   solid_3d.run();
3147 *   }
3148 *   else if (parameters.automatic_differentiation_order == 2)
3149 *   {
3150 *   std::cout << "Assembly method: Residual and linearisation computed using AD." << std::endl;
3151 *  
3152 * @endcode
3153 *
3154 * Sacado Rad-Fad is not thread-safe, so disable threading.
3155 * Parallisation using MPI would be possible though.
3156 *
3157 * @code
3158 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3159 *   1);
3160 *  
3161 *   typedef Sacado::Rad::ADvar< Sacado::Fad::DFad<double> > NumberType;
3162 *   Solid<dim,NumberType> solid_3d(parameters);
3163 *   solid_3d.run();
3164 *   }
3165 *   #endif
3166 *   else
3167 *   {
3168 *   AssertThrow(false,
3169 *   ExcMessage("The selected assembly method is not supported. "
3170 *   "You need deal.II 9.0 and Trilinos with the Sacado package "
3171 *   "to enable assembly using automatic differentiation."));
3172 *   }
3173 *   }
3174 *   catch (std::exception &exc)
3175 *   {
3176 *   std::cerr << std::endl << std::endl
3177 *   << "----------------------------------------------------"
3178 *   << std::endl;
3179 *   std::cerr << "Exception on processing: " << std::endl << exc.what()
3180 *   << std::endl << "Aborting!" << std::endl
3181 *   << "----------------------------------------------------"
3182 *   << std::endl;
3183 *  
3184 *   return 1;
3185 *   }
3186 *   catch (...)
3187 *   {
3188 *   std::cerr << std::endl << std::endl
3189 *   << "----------------------------------------------------"
3190 *   << std::endl;
3191 *   std::cerr << "Unknown exception!" << std::endl << "Aborting!"
3192 *   << std::endl
3193 *   << "----------------------------------------------------"
3194 *   << std::endl;
3195 *   return 1;
3196 *   }
3197 *  
3198 *   return 0;
3199 *   }
3200 * @endcode
3201
3202
3203*/
Definition fe_q.h:550
Definition point.h:111
numbers::NumberTraits< Number >::real_type norm() const
constexpr void clear()
@ wall_times
Definition timer.h:651
#define DEAL_II_VERSION_MAJOR
Definition config.h:27
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
__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(IteratorType begin, std_cxx20::type_identity_t< IteratorType > 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, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:442
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
const Event initial
Definition event.cc:64
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &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)
const std::vector< bool > & used
const Triangulation< dim, spacedim > & tria
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:471
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:49
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
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)
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)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
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={})
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:70
const InputIterator OutputIterator out
Definition parallel.h:167
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int boundary_id
Definition types.h:144
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 > &)