Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
Nonlinear_PoroViscoelasticity.h
Go to the documentation of this file.
1
598 *  
599 *   /* Authors: Ester Comellas and Jean-Paul Pelteret,
600 *   * University of Erlangen-Nuremberg, 2018
601 *   */
602 *  
603 * @endcode
604 *
605 * We start by including all the necessary deal.II header files and some C++
606 * related ones. They have been discussed in detail in previous tutorial
607 * programs, so you need only refer to past tutorials for details.
608 *
609
610 *
611 *
612 * @code
613 *   #include <deal.II/base/function.h>
614 *   #include <deal.II/base/parameter_handler.h>
615 *   #include <deal.II/base/point.h>
616 *   #include <deal.II/base/quadrature_lib.h>
617 *   #include <deal.II/base/symmetric_tensor.h>
618 *   #include <deal.II/base/tensor.h>
619 *   #include <deal.II/base/timer.h>
620 *   #include <deal.II/base/work_stream.h>
621 *   #include <deal.II/base/mpi.h>
622 *   #include <deal.II/base/quadrature_point_data.h>
623 *  
624 *   #include <deal.II/differentiation/ad.h>
625 *  
626 *   #include <deal.II/distributed/shared_tria.h>
627 *  
628 *   #include <deal.II/dofs/dof_renumbering.h>
629 *   #include <deal.II/dofs/dof_tools.h>
630 *   #include <deal.II/dofs/dof_accessor.h>
631 *  
632 *   #include <deal.II/grid/filtered_iterator.h>
633 *   #include <deal.II/grid/grid_generator.h>
634 *   #include <deal.II/grid/grid_tools.h>
635 *   #include <deal.II/grid/grid_in.h>
636 *   #include <deal.II/grid/grid_out.h>
637 *   #include <deal.II/grid/manifold_lib.h>
638 *   #include <deal.II/grid/tria_accessor.h>
639 *   #include <deal.II/grid/tria_iterator.h>
640 *  
641 *   #include <deal.II/fe/fe_dgp_monomial.h>
642 *   #include <deal.II/fe/fe_q.h>
643 *   #include <deal.II/fe/fe_system.h>
644 *   #include <deal.II/fe/fe_tools.h>
645 *   #include <deal.II/fe/fe_values.h>
646 *  
647 *   #include <deal.II/lac/block_sparsity_pattern.h>
648 *   #include <deal.II/lac/affine_constraints.h>
649 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
650 *   #include <deal.II/lac/full_matrix.h>
651 *   #include <deal.II/lac/linear_operator.h>
652 *   #include <deal.II/lac/packaged_operation.h>
653 *  
654 *   #include <deal.II/lac/trilinos_block_sparse_matrix.h>
655 *   #include <deal.II/lac/trilinos_linear_operator.h>
656 *   #include <deal.II/lac/trilinos_parallel_block_vector.h>
657 *   #include <deal.II/lac/trilinos_precondition.h>
658 *   #include <deal.II/lac/trilinos_sparse_matrix.h>
659 *   #include <deal.II/lac/trilinos_sparsity_pattern.h>
660 *   #include <deal.II/lac/trilinos_solver.h>
661 *   #include <deal.II/lac/trilinos_vector.h>
662 *  
663 *   #include <deal.II/lac/block_vector.h>
664 *   #include <deal.II/lac/vector.h>
665 *  
666 *   #include <deal.II/numerics/data_postprocessor.h>
667 *   #include <deal.II/numerics/data_out.h>
668 *   #include <deal.II/numerics/data_out_faces.h>
669 *   #include <deal.II/numerics/fe_field_function.h>
670 *   #include <deal.II/numerics/vector_tools.h>
671 *  
672 *   #include <deal.II/physics/transformations.h>
673 *   #include <deal.II/physics/elasticity/kinematics.h>
674 *   #include <deal.II/physics/elasticity/standard_tensors.h>
675 *  
676 *   #include <iostream>
677 *   #include <fstream>
678 *   #include <numeric>
679 *   #include <iomanip>
680 *  
681 *  
682 * @endcode
683 *
684 * We create a namespace for everything that relates to
685 * the nonlinear poro-viscoelastic formulation,
686 * and import all the deal.II function and class names into it:
687 *
688 * @code
689 *   namespace NonLinearPoroViscoElasticity
690 *   {
691 *   using namespace dealii;
692 *  
693 * @endcode
694 *
695 *
696 * <a name="Runtimeparameters"></a>
697 * <h3>Run-time parameters</h3>
698 *
699
700 *
701 * Set up a ParameterHandler object to read in the parameter choices at run-time
702 * introduced by the user through the file "parameters.prm"
703 *
704 * @code
705 *   namespace Parameters
706 *   {
707 * @endcode
708 *
709 *
710 * <a name="FiniteElementsystem"></a>
711 * <h4>Finite Element system</h4>
712 * Here we specify the polynomial order used to approximate the solution,
713 * both for the displacements and pressure unknowns.
714 * The quadrature order should be adjusted accordingly.
715 *
716 * @code
717 *   struct FESystem
718 *   {
719 *   unsigned int poly_degree_displ;
720 *   unsigned int poly_degree_pore;
721 *   unsigned int quad_order;
722 *  
723 *   static void
724 *   declare_parameters(ParameterHandler &prm);
725 *  
726 *   void
727 *   parse_parameters(ParameterHandler &prm);
728 *   };
729 *  
730 *   void FESystem::declare_parameters(ParameterHandler &prm)
731 *   {
732 *   prm.enter_subsection("Finite element system");
733 *   {
734 *   prm.declare_entry("Polynomial degree displ", "2",
735 *   Patterns::Integer(0),
736 *   "Displacement system polynomial order");
737 *  
738 *   prm.declare_entry("Polynomial degree pore", "1",
739 *   Patterns::Integer(0),
740 *   "Pore pressure system polynomial order");
741 *  
742 *   prm.declare_entry("Quadrature order", "3",
743 *   Patterns::Integer(0),
744 *   "Gauss quadrature order");
745 *   }
746 *   prm.leave_subsection();
747 *   }
748 *  
749 *   void FESystem::parse_parameters(ParameterHandler &prm)
750 *   {
751 *   prm.enter_subsection("Finite element system");
752 *   {
753 *   poly_degree_displ = prm.get_integer("Polynomial degree displ");
754 *   poly_degree_pore = prm.get_integer("Polynomial degree pore");
755 *   quad_order = prm.get_integer("Quadrature order");
756 *   }
757 *   prm.leave_subsection();
758 *   }
759 *  
760 * @endcode
761 *
762 *
763 * <a name="Geometry"></a>
764 * <h4>Geometry</h4>
765 * These parameters are related to the geometry definition and mesh generation.
766 * We select the type of problem to solve and introduce the desired load values.
767 *
768 * @code
769 *   struct Geometry
770 *   {
771 *   std::string geom_type;
772 *   unsigned int global_refinement;
773 *   double scale;
774 *   std::string load_type;
775 *   double load;
776 *   unsigned int num_cycle_sets;
777 *   double fluid_flow;
778 *   double drained_pressure;
779 *  
780 *   static void
781 *   declare_parameters(ParameterHandler &prm);
782 *  
783 *   void
784 *   parse_parameters(ParameterHandler &prm);
785 *   };
786 *  
787 *   void Geometry::declare_parameters(ParameterHandler &prm)
788 *   {
789 *   prm.enter_subsection("Geometry");
790 *   {
791 *   prm.declare_entry("Geometry type", "Ehlers_tube_step_load",
792 *   Patterns::Selection("Ehlers_tube_step_load"
793 *   "|Ehlers_tube_increase_load"
794 *   "|Ehlers_cube_consolidation"
795 *   "|Franceschini_consolidation"
796 *   "|Budday_cube_tension_compression"
797 *   "|Budday_cube_tension_compression_fully_fixed"
798 *   "|Budday_cube_shear_fully_fixed"),
799 *   "Type of geometry used. "
800 *   "For Ehlers verification examples see Ehlers and Eipper (1999). "
801 *   "For Franceschini brain consolidation see Franceschini et al. (2006)"
802 *   "For Budday brain examples see Budday et al. (2017)");
803 *  
804 *   prm.declare_entry("Global refinement", "1",
805 *   Patterns::Integer(0),
806 *   "Global refinement level");
807 *  
808 *   prm.declare_entry("Grid scale", "1.0",
809 *   Patterns::Double(0.0),
810 *   "Global grid scaling factor");
811 *  
812 *   prm.declare_entry("Load type", "pressure",
813 *   Patterns::Selection("pressure|displacement|none"),
814 *   "Type of loading");
815 *  
816 *   prm.declare_entry("Load value", "-7.5e+6",
817 *   Patterns::Double(),
818 *   "Loading value");
819 *  
820 *   prm.declare_entry("Number of cycle sets", "1",
821 *   Patterns::Integer(1,2),
822 *   "Number of times each set of 3 cycles is repeated, only for "
823 *   "Budday_cube_tension_compression and Budday_cube_tension_compression_fully_fixed. "
824 *   "Load value is doubled in second set, load rate is kept constant."
825 *   "Final time indicates end of second cycle set.");
826 *  
827 *   prm.declare_entry("Fluid flow value", "0.0",
828 *   Patterns::Double(),
829 *   "Prescribed fluid flow. Not implemented in any example yet.");
830 *  
831 *   prm.declare_entry("Drained pressure", "0.0",
832 *   Patterns::Double(),
833 *   "Increase of pressure value at drained boundary w.r.t the atmospheric pressure.");
834 *   }
835 *   prm.leave_subsection();
836 *   }
837 *  
838 *   void Geometry::parse_parameters(ParameterHandler &prm)
839 *   {
840 *   prm.enter_subsection("Geometry");
841 *   {
842 *   geom_type = prm.get("Geometry type");
843 *   global_refinement = prm.get_integer("Global refinement");
844 *   scale = prm.get_double("Grid scale");
845 *   load_type = prm.get("Load type");
846 *   load = prm.get_double("Load value");
847 *   num_cycle_sets = prm.get_integer("Number of cycle sets");
848 *   fluid_flow = prm.get_double("Fluid flow value");
849 *   drained_pressure = prm.get_double("Drained pressure");
850 *   }
851 *   prm.leave_subsection();
852 *   }
853 *  
854 * @endcode
855 *
856 *
857 * <a name="Materials"></a>
858 * <h4>Materials</h4>
859 *
860
861 *
862 * Here we select the type of material for the solid component
863 * and define the corresponding material parameters.
864 * Then we define he fluid data, including the type of
865 * seepage velocity definition to use.
866 *
867 * @code
868 *   struct Materials
869 *   {
870 *   std::string mat_type;
871 *   double lambda;
872 *   double mu;
873 *   double mu1_infty;
874 *   double mu2_infty;
875 *   double mu3_infty;
876 *   double alpha1_infty;
877 *   double alpha2_infty;
878 *   double alpha3_infty;
879 *   double mu1_mode_1;
880 *   double mu2_mode_1;
881 *   double mu3_mode_1;
882 *   double alpha1_mode_1;
883 *   double alpha2_mode_1;
884 *   double alpha3_mode_1;
885 *   double viscosity_mode_1;
886 *   std::string fluid_type;
887 *   double solid_vol_frac;
888 *   double kappa_darcy;
889 *   double init_intrinsic_perm;
890 *   double viscosity_FR;
891 *   double init_darcy_coef;
892 *   double weight_FR;
893 *   bool gravity_term;
894 *   int gravity_direction;
895 *   double gravity_value;
896 *   double density_FR;
897 *   double density_SR;
898 *   enum SymmetricTensorEigenvectorMethod eigen_solver;
899 *  
900 *   static void
901 *   declare_parameters(ParameterHandler &prm);
902 *  
903 *   void
904 *   parse_parameters(ParameterHandler &prm);
905 *   };
906 *  
907 *   void Materials::declare_parameters(ParameterHandler &prm)
908 *   {
909 *   prm.enter_subsection("Material properties");
910 *   {
911 *   prm.declare_entry("material", "Neo-Hooke",
912 *   Patterns::Selection("Neo-Hooke|Ogden|visco-Ogden"),
913 *   "Type of material used in the problem");
914 *  
915 *   prm.declare_entry("lambda", "8.375e6",
916 *   Patterns::Double(0,1e100),
917 *   "First Lamé parameter for extension function related to compactation point in solid material [Pa].");
918 *  
919 *   prm.declare_entry("shear modulus", "5.583e6",
920 *   Patterns::Double(0,1e100),
921 *   "shear modulus for Neo-Hooke materials [Pa].");
922 *  
923 *   prm.declare_entry("eigen solver", "QL Implicit Shifts",
924 *   Patterns::Selection("QL Implicit Shifts|Jacobi"),
925 *   "The type of eigen solver to be used for Ogden and visco-Ogden models.");
926 *  
927 *   prm.declare_entry("mu1", "0.0",
928 *   Patterns::Double(),
929 *   "Shear material parameter 'mu1' for Ogden material [Pa].");
930 *  
931 *   prm.declare_entry("mu2", "0.0",
932 *   Patterns::Double(),
933 *   "Shear material parameter 'mu2' for Ogden material [Pa].");
934 *  
935 *   prm.declare_entry("mu3", "0.0",
936 *   Patterns::Double(),
937 *   "Shear material parameter 'mu1' for Ogden material [Pa].");
938 *  
939 *   prm.declare_entry("alpha1", "1.0",
940 *   Patterns::Double(),
941 *   "Stiffness material parameter 'alpha1' for Ogden material [-].");
942 *  
943 *   prm.declare_entry("alpha2", "1.0",
944 *   Patterns::Double(),
945 *   "Stiffness material parameter 'alpha2' for Ogden material [-].");
946 *  
947 *   prm.declare_entry("alpha3", "1.0",
948 *   Patterns::Double(),
949 *   "Stiffness material parameter 'alpha3' for Ogden material [-].");
950 *  
951 *   prm.declare_entry("mu1_1", "0.0",
952 *   Patterns::Double(),
953 *   "Shear material parameter 'mu1' for first viscous mode in Ogden material [Pa].");
954 *  
955 *   prm.declare_entry("mu2_1", "0.0",
956 *   Patterns::Double(),
957 *   "Shear material parameter 'mu2' for first viscous mode in Ogden material [Pa].");
958 *  
959 *   prm.declare_entry("mu3_1", "0.0",
960 *   Patterns::Double(),
961 *   "Shear material parameter 'mu1' for first viscous mode in Ogden material [Pa].");
962 *  
963 *   prm.declare_entry("alpha1_1", "1.0",
964 *   Patterns::Double(),
965 *   "Stiffness material parameter 'alpha1' for first viscous mode in Ogden material [-].");
966 *  
967 *   prm.declare_entry("alpha2_1", "1.0",
968 *   Patterns::Double(),
969 *   "Stiffness material parameter 'alpha2' for first viscous mode in Ogden material [-].");
970 *  
971 *   prm.declare_entry("alpha3_1", "1.0",
972 *   Patterns::Double(),
973 *   "Stiffness material parameter 'alpha3' for first viscous mode in Ogden material [-].");
974 *  
975 *   prm.declare_entry("viscosity_1", "1e-10",
976 *   Patterns::Double(1e-10,1e100),
977 *   "Deformation-independent viscosity parameter 'eta_1' for first viscous mode in Ogden material [-].");
978 *  
979 *   prm.declare_entry("seepage definition", "Ehlers",
980 *   Patterns::Selection("Markert|Ehlers"),
981 *   "Type of formulation used to define the seepage velocity in the problem. "
982 *   "Choose between Markert formulation of deformation-dependent intrinsic permeability "
983 *   "and Ehlers formulation of deformation-dependent Darcy flow coefficient.");
984 *  
985 *   prm.declare_entry("initial solid volume fraction", "0.67",
986 *   Patterns::Double(0.001,0.999),
987 *   "Initial porosity (solid volume fraction, 0 < n_0s < 1)");
988 *  
989 *   prm.declare_entry("kappa", "0.0",
990 *   Patterns::Double(0,100),
991 *   "Deformation-dependency control parameter for specific permeability (kappa >= 0)");
992 *  
993 *   prm.declare_entry("initial intrinsic permeability", "0.0",
994 *   Patterns::Double(0,1e100),
995 *   "Initial intrinsic permeability parameter [m^2] (isotropic permeability). To be used with Markert formulation.");
996 *  
997 *   prm.declare_entry("fluid viscosity", "0.0",
998 *   Patterns::Double(0, 1e100),
999 *   "Effective shear viscosity parameter of the fluid [Pa·s, (N·s)/m^2]. To be used with Markert formulation.");
1000 *  
1001 *   prm.declare_entry("initial Darcy coefficient", "1.0e-4",
1002 *   Patterns::Double(0,1e100),
1003 *   "Initial Darcy flow coefficient [m/s] (isotropic permeability). To be used with Ehlers formulation.");
1004 *  
1005 *   prm.declare_entry("fluid weight", "1.0e4",
1006 *   Patterns::Double(0, 1e100),
1007 *   "Effective weight of the fluid [N/m^3]. To be used with Ehlers formulation.");
1008 *  
1009 *   prm.declare_entry("gravity term", "false",
1010 *   Patterns::Bool(),
1011 *   "Gravity term considered (true) or neglected (false)");
1012 *  
1013 *   prm.declare_entry("fluid density", "1.0",
1014 *   Patterns::Double(0,1e100),
1015 *   "Real (or effective) density of the fluid");
1016 *  
1017 *   prm.declare_entry("solid density", "1.0",
1018 *   Patterns::Double(0,1e100),
1019 *   "Real (or effective) density of the solid");
1020 *  
1021 *   prm.declare_entry("gravity direction", "2",
1022 *   Patterns::Integer(0,2),
1023 *   "Direction of gravity (unit vector 0 for x, 1 for y, 2 for z)");
1024 *  
1025 *   prm.declare_entry("gravity value", "-9.81",
1026 *   Patterns::Double(),
1027 *   "Value of gravity (be careful to have consistent units!)");
1028 *   }
1029 *   prm.leave_subsection();
1030 *   }
1031 *  
1032 *   void Materials::parse_parameters(ParameterHandler &prm)
1033 *   {
1034 *   prm.enter_subsection("Material properties");
1035 *   {
1036 * @endcode
1037 *
1038 * Solid
1039 *
1040 * @code
1041 *   mat_type = prm.get("material");
1042 *   lambda = prm.get_double("lambda");
1043 *   mu = prm.get_double("shear modulus");
1044 *   mu1_infty = prm.get_double("mu1");
1045 *   mu2_infty = prm.get_double("mu2");
1046 *   mu3_infty = prm.get_double("mu3");
1047 *   alpha1_infty = prm.get_double("alpha1");
1048 *   alpha2_infty = prm.get_double("alpha2");
1049 *   alpha3_infty = prm.get_double("alpha3");
1050 *   mu1_mode_1 = prm.get_double("mu1_1");
1051 *   mu2_mode_1 = prm.get_double("mu2_1");
1052 *   mu3_mode_1 = prm.get_double("mu3_1");
1053 *   alpha1_mode_1 = prm.get_double("alpha1_1");
1054 *   alpha2_mode_1 = prm.get_double("alpha2_1");
1055 *   alpha3_mode_1 = prm.get_double("alpha3_1");
1056 *   viscosity_mode_1 = prm.get_double("viscosity_1");
1057 * @endcode
1058 *
1059 * Fluid
1060 *
1061 * @code
1062 *   fluid_type = prm.get("seepage definition");
1063 *   solid_vol_frac = prm.get_double("initial solid volume fraction");
1064 *   kappa_darcy = prm.get_double("kappa");
1065 *   init_intrinsic_perm = prm.get_double("initial intrinsic permeability");
1066 *   viscosity_FR = prm.get_double("fluid viscosity");
1067 *   init_darcy_coef = prm.get_double("initial Darcy coefficient");
1068 *   weight_FR = prm.get_double("fluid weight");
1069 * @endcode
1070 *
1071 * Gravity effects
1072 *
1073 * @code
1074 *   gravity_term = prm.get_bool("gravity term");
1075 *   density_FR = prm.get_double("fluid density");
1076 *   density_SR = prm.get_double("solid density");
1077 *   gravity_direction = prm.get_integer("gravity direction");
1078 *   gravity_value = prm.get_double("gravity value");
1079 *  
1080 *   if ( (fluid_type == "Markert") && ((init_intrinsic_perm == 0.0) || (viscosity_FR == 0.0)) )
1081 *   AssertThrow(false, ExcMessage("Markert seepage velocity formulation requires the definition of "
1082 *   "'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0."));
1083 *  
1084 *   if ( (fluid_type == "Ehlers") && ((init_darcy_coef == 0.0) || (weight_FR == 0.0)) )
1085 *   AssertThrow(false, ExcMessage("Ehler seepage velocity formulation requires the definition of "
1086 *   "'initial Darcy coefficient' and 'fluid weight' greater than 0.0."));
1087 *  
1088 *   const std::string eigen_solver_type = prm.get("eigen solver");
1089 *   if (eigen_solver_type == "QL Implicit Shifts")
1090 *   eigen_solver = SymmetricTensorEigenvectorMethod::ql_implicit_shifts;
1091 *   else if (eigen_solver_type == "Jacobi")
1092 *   eigen_solver = SymmetricTensorEigenvectorMethod::jacobi;
1093 *   else
1094 *   {
1095 *   AssertThrow(false, ExcMessage("Unknown eigen solver selected."));
1096 *   }
1097 *   }
1098 *   prm.leave_subsection();
1099 *   }
1100 *  
1101 * @endcode
1102 *
1103 *
1104 * <a name="Nonlinearsolver"></a>
1105 * <h4>Nonlinear solver</h4>
1106 *
1107
1108 *
1109 * We now define the tolerances and the maximum number of iterations for the
1110 * Newton-Raphson scheme used to solve the nonlinear system of governing equations.
1111 *
1112 * @code
1113 *   struct NonlinearSolver
1114 *   {
1115 *   unsigned int max_iterations_NR;
1116 *   double tol_f;
1117 *   double tol_u;
1118 *   double tol_p_fluid;
1119 *  
1120 *   static void
1121 *   declare_parameters(ParameterHandler &prm);
1122 *  
1123 *   void
1124 *   parse_parameters(ParameterHandler &prm);
1125 *   };
1126 *  
1127 *   void NonlinearSolver::declare_parameters(ParameterHandler &prm)
1128 *   {
1129 *   prm.enter_subsection("Nonlinear solver");
1130 *   {
1131 *   prm.declare_entry("Max iterations Newton-Raphson", "15",
1132 *   Patterns::Integer(0),
1133 *   "Number of Newton-Raphson iterations allowed");
1134 *  
1135 *   prm.declare_entry("Tolerance force", "1.0e-8",
1136 *   Patterns::Double(0.0),
1137 *   "Force residual tolerance");
1138 *  
1139 *   prm.declare_entry("Tolerance displacement", "1.0e-6",
1140 *   Patterns::Double(0.0),
1141 *   "Displacement error tolerance");
1142 *  
1143 *   prm.declare_entry("Tolerance pore pressure", "1.0e-6",
1144 *   Patterns::Double(0.0),
1145 *   "Pore pressure error tolerance");
1146 *   }
1147 *   prm.leave_subsection();
1148 *   }
1149 *  
1150 *   void NonlinearSolver::parse_parameters(ParameterHandler &prm)
1151 *   {
1152 *   prm.enter_subsection("Nonlinear solver");
1153 *   {
1154 *   max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
1155 *   tol_f = prm.get_double("Tolerance force");
1156 *   tol_u = prm.get_double("Tolerance displacement");
1157 *   tol_p_fluid = prm.get_double("Tolerance pore pressure");
1158 *   }
1159 *   prm.leave_subsection();
1160 *   }
1161 *  
1162 * @endcode
1163 *
1164 *
1165 * <a name="Time"></a>
1166 * <h4>Time</h4>
1167 * Here we set the timestep size @f$ \varDelta t @f$ and the simulation end-time.
1168 *
1169 * @code
1170 *   struct Time
1171 *   {
1172 *   double end_time;
1173 *   double delta_t;
1174 *   static void
1175 *   declare_parameters(ParameterHandler &prm);
1176 *  
1177 *   void
1178 *   parse_parameters(ParameterHandler &prm);
1179 *   };
1180 *  
1181 *   void Time::declare_parameters(ParameterHandler &prm)
1182 *   {
1183 *   prm.enter_subsection("Time");
1184 *   {
1185 *   prm.declare_entry("End time", "10.0",
1186 *   Patterns::Double(),
1187 *   "End time");
1188 *  
1189 *   prm.declare_entry("Time step size", "0.002",
1190 *   Patterns::Double(1.0e-6),
1191 *   "Time step size. The value must be larger than the displacement error tolerance defined.");
1192 *   }
1193 *   prm.leave_subsection();
1194 *   }
1195 *  
1196 *   void Time::parse_parameters(ParameterHandler &prm)
1197 *   {
1198 *   prm.enter_subsection("Time");
1199 *   {
1200 *   end_time = prm.get_double("End time");
1201 *   delta_t = prm.get_double("Time step size");
1202 *   }
1203 *   prm.leave_subsection();
1204 *   }
1205 *  
1206 *  
1207 * @endcode
1208 *
1209 *
1210 * <a name="Output"></a>
1211 * <h4>Output</h4>
1212 * We can choose the frequency of the data for the output files.
1213 *
1214 * @code
1215 *   struct OutputParam
1216 *   {
1217 *  
1218 *   std::string outfiles_requested;
1219 *   unsigned int timestep_output;
1220 *   std::string outtype;
1221 *  
1222 *   static void
1223 *   declare_parameters(ParameterHandler &prm);
1224 *  
1225 *   void
1226 *   parse_parameters(ParameterHandler &prm);
1227 *   };
1228 *  
1229 *   void OutputParam::declare_parameters(ParameterHandler &prm)
1230 *   {
1231 *   prm.enter_subsection("Output parameters");
1232 *   {
1233 *   prm.declare_entry("Output files", "true",
1234 *   Patterns::Selection("true|false"),
1235 *   "Paraview output files to generate.");
1236 *   prm.declare_entry("Time step number output", "1",
1237 *   Patterns::Integer(0),
1238 *   "Output data for time steps multiple of the given "
1239 *   "integer value.");
1240 *   prm.declare_entry("Averaged results", "nodes",
1241 *   Patterns::Selection("elements|nodes"),
1242 *   "Output data associated with integration point values"
1243 *   " averaged on elements or on nodes.");
1244 *   }
1245 *   prm.leave_subsection();
1246 *   }
1247 *  
1248 *   void OutputParam::parse_parameters(ParameterHandler &prm)
1249 *   {
1250 *   prm.enter_subsection("Output parameters");
1251 *   {
1252 *   outfiles_requested = prm.get("Output files");
1253 *   timestep_output = prm.get_integer("Time step number output");
1254 *   outtype = prm.get("Averaged results");
1255 *   }
1256 *   prm.leave_subsection();
1257 *   }
1258 *  
1259 * @endcode
1260 *
1261 *
1262 * <a name="Allparameters"></a>
1263 * <h4>All parameters</h4>
1264 * We finally consolidate all of the above structures into a single container that holds all the run-time selections.
1265 *
1266 * @code
1267 *   struct AllParameters : public FESystem,
1268 *   public Geometry,
1269 *   public Materials,
1270 *   public NonlinearSolver,
1271 *   public Time,
1272 *   public OutputParam
1273 *   {
1274 *   AllParameters(const std::string &input_file);
1275 *  
1276 *   static void
1277 *   declare_parameters(ParameterHandler &prm);
1278 *  
1279 *   void
1280 *   parse_parameters(ParameterHandler &prm);
1281 *   };
1282 *  
1283 *   AllParameters::AllParameters(const std::string &input_file)
1284 *   {
1285 *   ParameterHandler prm;
1286 *   declare_parameters(prm);
1287 *   prm.parse_input(input_file);
1288 *   parse_parameters(prm);
1289 *   }
1290 *  
1291 *   void AllParameters::declare_parameters(ParameterHandler &prm)
1292 *   {
1293 *   FESystem::declare_parameters(prm);
1294 *   Geometry::declare_parameters(prm);
1295 *   Materials::declare_parameters(prm);
1296 *   NonlinearSolver::declare_parameters(prm);
1297 *   Time::declare_parameters(prm);
1298 *   OutputParam::declare_parameters(prm);
1299 *   }
1300 *  
1301 *   void AllParameters::parse_parameters(ParameterHandler &prm)
1302 *   {
1303 *   FESystem::parse_parameters(prm);
1304 *   Geometry::parse_parameters(prm);
1305 *   Materials::parse_parameters(prm);
1306 *   NonlinearSolver::parse_parameters(prm);
1307 *   Time::parse_parameters(prm);
1308 *   OutputParam::parse_parameters(prm);
1309 *   }
1310 *   }
1311 *  
1312 * @endcode
1313 *
1314 *
1315 * <a name="Timeclass"></a>
1316 * <h3>Time class</h3>
1317 * A simple class to store time data.
1318 * For simplicity we assume a constant time step size.
1319 *
1320 * @code
1321 *   class Time
1322 *   {
1323 *   public:
1324 *   Time (const double time_end,
1325 *   const double delta_t)
1326 *   :
1327 *   timestep(0),
1328 *   time_current(0.0),
1329 *   time_end(time_end),
1330 *   delta_t(delta_t)
1331 *   {}
1332 *  
1333 *   virtual ~Time()
1334 *   {}
1335 *  
1336 *   double get_current() const
1337 *   {
1338 *   return time_current;
1339 *   }
1340 *   double get_end() const
1341 *   {
1342 *   return time_end;
1343 *   }
1344 *   double get_delta_t() const
1345 *   {
1346 *   return delta_t;
1347 *   }
1348 *   unsigned int get_timestep() const
1349 *   {
1350 *   return timestep;
1351 *   }
1352 *   void increment_time ()
1353 *   {
1354 *   time_current += delta_t;
1355 *   ++timestep;
1356 *   }
1357 *  
1358 *   private:
1359 *   unsigned int timestep;
1360 *   double time_current;
1361 *   double time_end;
1362 *   const double delta_t;
1363 *   };
1364 *  
1365 * @endcode
1366 *
1367 *
1368 * <a name="Constitutiveequationforthesolidcomponentofthebiphasicmaterial"></a>
1369 * <h3>Constitutive equation for the solid component of the biphasic material</h3>
1370 *
1371
1372 *
1373 *
1374 * <a name="Baseclassgenerichyperelasticmaterial"></a>
1375 * <h4>Base class: generic hyperelastic material</h4>
1376 * The ``extra" Kirchhoff stress in the solid component is the sum of isochoric
1377 * and a volumetric part.
1378 * @f$\mathbf{\tau} = \mathbf{\tau}_E^{(\bullet)} + \mathbf{\tau}^{\textrm{vol}}@f$
1379 * The deviatoric part changes depending on the type of material model selected:
1380 * Neo-Hooken hyperelasticity, Ogden hyperelasticiy,
1381 * or a single-mode finite viscoelasticity based on the Ogden hyperelastic model.
1382 * In this base class we declare it as a virtual function,
1383 * and it will be defined for each model type in the corresponding derived class.
1384 * We define here the volumetric component, which depends on the
1385 * extension function @f$U(J_S)@f$ selected, and in this case is the same for all models.
1386 * We use the function proposed by
1387 * Ehlers & Eipper 1999 doi:10.1023/A:1006565509095
1388 * We also define some public functions to access and update the internal variables.
1389 *
1390 * @code
1391 *   template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1392 *   class Material_Hyperelastic
1393 *   {
1394 *   public:
1395 *   Material_Hyperelastic(const Parameters::AllParameters &parameters,
1396 *   const Time &time)
1397 *   :
1398 *   n_OS (parameters.solid_vol_frac),
1399 *   lambda (parameters.lambda),
1400 *   time(time),
1401 *   det_F (1.0),
1402 *   det_F_converged (1.0),
1403 *   eigen_solver (parameters.eigen_solver)
1404 *   {}
1405 *   ~Material_Hyperelastic()
1406 *   {}
1407 *  
1408 *   SymmetricTensor<2, dim, NumberType>
1409 *   get_tau_E(const Tensor<2,dim, NumberType> &F) const
1410 *   {
1411 *   return ( get_tau_E_base(F) + get_tau_E_ext_func(F) );
1412 *   }
1413 *  
1414 *   SymmetricTensor<2, dim, NumberType>
1415 *   get_Cauchy_E(const Tensor<2, dim, NumberType> &F) const
1416 *   {
1417 *   const NumberType det_F = determinant(F);
1418 *   Assert(det_F > 0, ExcInternalError());
1419 *   return get_tau_E(F)*NumberType(1/det_F);
1420 *   }
1421 *  
1422 *   double
1423 *   get_converged_det_F() const
1424 *   {
1425 *   return det_F_converged;
1426 *   }
1427 *  
1428 *   virtual void
1429 *   update_end_timestep()
1430 *   {
1431 *   det_F_converged = det_F;
1432 *   }
1433 *  
1434 *   virtual void
1435 *   update_internal_equilibrium( const Tensor<2, dim, NumberType> &F )
1436 *   {
1437 *   det_F = Tensor<0,dim,double>(determinant(F));
1438 *   }
1439 *  
1440 *   virtual double
1441 *   get_viscous_dissipation( ) const = 0;
1442 *  
1443 *   const double n_OS;
1444 *   const double lambda;
1445 *   const Time &time;
1446 *   double det_F;
1447 *   double det_F_converged;
1448 *   const enum SymmetricTensorEigenvectorMethod eigen_solver;
1449 *  
1450 *   protected:
1451 *   SymmetricTensor<2, dim, NumberType>
1452 *   get_tau_E_ext_func(const Tensor<2,dim, NumberType> &F) const
1453 *   {
1454 *   const NumberType det_F = determinant(F);
1455 *   Assert(det_F > 0, ExcInternalError());
1456 *  
1457 *   static const SymmetricTensor< 2, dim, double>
1458 *   I (Physics::Elasticity::StandardTensors<dim>::I);
1459 *   return ( NumberType(lambda * (1.0-n_OS)*(1.0-n_OS)
1460 *   * (det_F/(1.0-n_OS) - det_F/(det_F-n_OS))) * I );
1461 *   }
1462 *  
1463 *   virtual SymmetricTensor<2, dim, NumberType>
1464 *   get_tau_E_base(const Tensor<2,dim, NumberType> &F) const = 0;
1465 *   };
1466 *  
1467 * @endcode
1468 *
1469 *
1470 * <a name="DerivedclassNeoHookeanhyperelasticmaterial"></a>
1471 * <h4>Derived class: Neo-Hookean hyperelastic material</h4>
1472 *
1473 * @code
1474 *   template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1475 *   class NeoHooke : public Material_Hyperelastic < dim, NumberType >
1476 *   {
1477 *   public:
1478 *   NeoHooke(const Parameters::AllParameters &parameters,
1479 *   const Time &time)
1480 *   :
1481 *   Material_Hyperelastic< dim, NumberType > (parameters,time),
1482 *   mu(parameters.mu)
1483 *   {}
1484 *   virtual ~NeoHooke()
1485 *   {}
1486 *  
1487 *   double
1488 *   get_viscous_dissipation() const override
1489 *   {
1490 *   return 0.0;
1491 *   }
1492 *  
1493 *   protected:
1494 *   const double mu;
1495 *  
1496 *   SymmetricTensor<2, dim, NumberType>
1497 *   get_tau_E_base(const Tensor<2,dim, NumberType> &F) const override
1498 *   {
1499 *   static const SymmetricTensor< 2, dim, double>
1500 *   I (Physics::Elasticity::StandardTensors<dim>::I);
1501 *  
1502 *   const bool use_standard_model = true;
1503 *  
1504 *   if (use_standard_model)
1505 *   {
1506 * @endcode
1507 *
1508 * Standard Neo-Hooke
1509 *
1510 * @code
1511 *   return ( mu * ( symmetrize(F * transpose(F)) - I ) );
1512 *   }
1513 *   else
1514 *   {
1515 * @endcode
1516 *
1517 * Neo-Hooke in terms of principal stretches
1518 *
1519 * @code
1520 *   const SymmetricTensor<2, dim, NumberType>
1521 *   B = symmetrize(F * transpose(F));
1522 *   const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1523 *   eigen_B = eigenvectors(B, this->eigen_solver);
1524 *  
1525 *   SymmetricTensor<2, dim, NumberType> B_ev;
1526 *   for (unsigned int d=0; d<dim; ++d)
1527 *   B_ev += eigen_B[d].first*symmetrize(outer_product(eigen_B[d].second,eigen_B[d].second));
1528 *  
1529 *   return ( mu*(B_ev-I) );
1530 *   }
1531 *   }
1532 *   };
1533 *  
1534 * @endcode
1535 *
1536 *
1537 * <a name="DerivedclassOgdenhyperelasticmaterial"></a>
1538 * <h4>Derived class: Ogden hyperelastic material</h4>
1539 *
1540 * @code
1541 *   template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1542 *   class Ogden : public Material_Hyperelastic < dim, NumberType >
1543 *   {
1544 *   public:
1545 *   Ogden(const Parameters::AllParameters &parameters,
1546 *   const Time &time)
1547 *   :
1548 *   Material_Hyperelastic< dim, NumberType > (parameters,time),
1549 *   mu({parameters.mu1_infty,
1550 *   parameters.mu2_infty,
1551 *   parameters.mu3_infty}),
1552 *   alpha({parameters.alpha1_infty,
1553 *   parameters.alpha2_infty,
1554 *   parameters.alpha3_infty})
1555 *   {}
1556 *   virtual ~Ogden()
1557 *   {}
1558 *  
1559 *   double
1560 *   get_viscous_dissipation() const override
1561 *   {
1562 *   return 0.0;
1563 *   }
1564 *  
1565 *   protected:
1566 *   std::vector<double> mu;
1567 *   std::vector<double> alpha;
1568 *  
1569 *   SymmetricTensor<2, dim, NumberType>
1570 *   get_tau_E_base(const Tensor<2,dim, NumberType> &F) const override
1571 *   {
1572 *   const SymmetricTensor<2, dim, NumberType>
1573 *   B = symmetrize(F * transpose(F));
1574 *  
1575 *   const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1576 *   eigen_B = eigenvectors(B, this->eigen_solver);
1577 *  
1578 *   SymmetricTensor<2, dim, NumberType> tau;
1579 *   static const SymmetricTensor< 2, dim, double>
1580 *   I (Physics::Elasticity::StandardTensors<dim>::I);
1581 *  
1582 *   for (unsigned int i = 0; i < 3; ++i)
1583 *   {
1584 *   for (unsigned int A = 0; A < dim; ++A)
1585 *   {
1586 *   SymmetricTensor<2, dim, NumberType> tau_aux1 = symmetrize(
1587 *   outer_product(eigen_B[A].second,eigen_B[A].second));
1588 *   tau_aux1 *= mu[i]*std::pow(eigen_B[A].first, (alpha[i]/2.) );
1589 *   tau += tau_aux1;
1590 *   }
1591 *   SymmetricTensor<2, dim, NumberType> tau_aux2 (I);
1592 *   tau_aux2 *= mu[i];
1593 *   tau -= tau_aux2;
1594 *   }
1595 *   return tau;
1596 *   }
1597 *   };
1598 *  
1599 * @endcode
1600 *
1601 *
1602 * <a name="DerivedclassSinglemodeOgdenviscoelasticmaterial"></a>
1603 * <h4>Derived class: Single-mode Ogden viscoelastic material</h4>
1604 * We use the finite viscoelastic model described in
1605 * Reese & Govindjee (1998) doi:10.1016/S0020-7683(97)00217-5
1606 * The algorithm for the implicit exponential time integration is given in
1607 * Budday et al. (2017) doi: 10.1016/j.actbio.2017.06.024
1608 *
1609 * @code
1610 *   template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1611 *   class visco_Ogden : public Material_Hyperelastic < dim, NumberType >
1612 *   {
1613 *   public:
1614 *   visco_Ogden(const Parameters::AllParameters &parameters,
1615 *   const Time &time)
1616 *   :
1617 *   Material_Hyperelastic< dim, NumberType > (parameters,time),
1618 *   mu_infty({parameters.mu1_infty,
1619 *   parameters.mu2_infty,
1620 *   parameters.mu3_infty}),
1621 *   alpha_infty({parameters.alpha1_infty,
1622 *   parameters.alpha2_infty,
1623 *   parameters.alpha3_infty}),
1624 *   mu_mode_1({parameters.mu1_mode_1,
1625 *   parameters.mu2_mode_1,
1626 *   parameters.mu3_mode_1}),
1627 *   alpha_mode_1({parameters.alpha1_mode_1,
1628 *   parameters.alpha2_mode_1,
1629 *   parameters.alpha3_mode_1}),
1630 *   viscosity_mode_1(parameters.viscosity_mode_1),
1631 *   Cinv_v_1(Physics::Elasticity::StandardTensors<dim>::I),
1632 *   Cinv_v_1_converged(Physics::Elasticity::StandardTensors<dim>::I)
1633 *   {}
1634 *   virtual ~visco_Ogden()
1635 *   {}
1636 *  
1637 *   void
1638 *   update_internal_equilibrium( const Tensor<2, dim, NumberType> &F ) override
1639 *   {
1640 *   Material_Hyperelastic < dim, NumberType >::update_internal_equilibrium(F);
1641 *  
1642 *   this->Cinv_v_1 = this->Cinv_v_1_converged;
1643 *   SymmetricTensor<2, dim, NumberType> B_e_1_tr = symmetrize(F * this->Cinv_v_1 * transpose(F));
1644 *  
1645 *   const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1646 *   eigen_B_e_1_tr = eigenvectors(B_e_1_tr, this->eigen_solver);
1647 *  
1648 *   Tensor< 1, dim, NumberType > lambdas_e_1_tr;
1649 *   Tensor< 1, dim, NumberType > epsilon_e_1_tr;
1650 *   for (int a = 0; a < dim; ++a)
1651 *   {
1652 *   lambdas_e_1_tr[a] = std::sqrt(eigen_B_e_1_tr[a].first);
1653 *   epsilon_e_1_tr[a] = std::log(lambdas_e_1_tr[a]);
1654 *   }
1655 *  
1656 *   const double tolerance = 1e-8;
1657 *   double residual_check = tolerance*10.0;
1658 *   Tensor< 1, dim, NumberType > residual;
1659 *   Tensor< 2, dim, NumberType > tangent;
1660 *   static const SymmetricTensor< 2, dim, double> I(Physics::Elasticity::StandardTensors<dim>::I);
1661 *   NumberType J_e_1 = std::sqrt(determinant(B_e_1_tr));
1662 *  
1663 *   std::vector<NumberType> lambdas_e_1_iso(dim);
1664 *   SymmetricTensor<2, dim, NumberType> B_e_1;
1665 *   int iteration = 0;
1666 *  
1667 *   Tensor< 1, dim, NumberType > lambdas_e_1;
1668 *   Tensor< 1, dim, NumberType > epsilon_e_1;
1669 *   epsilon_e_1 = epsilon_e_1_tr;
1670 *  
1671 *   while(residual_check > tolerance)
1672 *   {
1673 *   NumberType aux_J_e_1 = 1.0;
1674 *   for (unsigned int a = 0; a < dim; ++a)
1675 *   {
1676 *   lambdas_e_1[a] = std::exp(epsilon_e_1[a]);
1677 *   aux_J_e_1 *= lambdas_e_1[a];
1678 *   }
1679 *  
1680 *   J_e_1 = aux_J_e_1;
1681 *  
1682 *   for (unsigned int a = 0; a < dim; ++a)
1683 *   lambdas_e_1_iso[a] = lambdas_e_1[a]*std::pow(J_e_1,-1.0/dim);
1684 *  
1685 *   for (unsigned int a = 0; a < dim; ++a)
1686 *   {
1687 *   residual[a] = get_beta_mode_1(lambdas_e_1_iso, a);
1688 *   residual[a] *= this->time.get_delta_t()/(2.0*viscosity_mode_1);
1689 *   residual[a] += epsilon_e_1[a];
1690 *   residual[a] -= epsilon_e_1_tr[a];
1691 *  
1692 *   for (unsigned int b = 0; b < dim; ++b)
1693 *   {
1694 *   tangent[a][b] = get_gamma_mode_1(lambdas_e_1_iso, a, b);
1695 *   tangent[a][b] *= this->time.get_delta_t()/(2.0*viscosity_mode_1);
1696 *   tangent[a][b] += I[a][b];
1697 *   }
1698 *  
1699 *   }
1700 *   epsilon_e_1 -= invert(tangent)*residual;
1701 *  
1702 *   residual_check = 0.0;
1703 *   for (unsigned int a = 0; a < dim; ++a)
1704 *   {
1705 *   if ( std::abs(residual[a]) > residual_check)
1706 *   residual_check = std::abs(Tensor<0,dim,double>(residual[a]));
1707 *   }
1708 *   iteration += 1;
1709 *   if (iteration > 15 )
1710 *   AssertThrow(false, ExcMessage("No convergence in local Newton iteration for the "
1711 *   "viscoelastic exponential time integration algorithm."));
1712 *   }
1713 *  
1714 *   NumberType aux_J_e_1 = 1.0;
1715 *   for (unsigned int a = 0; a < dim; ++a)
1716 *   {
1717 *   lambdas_e_1[a] = std::exp(epsilon_e_1[a]);
1718 *   aux_J_e_1 *= lambdas_e_1[a];
1719 *   }
1720 *   J_e_1 = aux_J_e_1;
1721 *  
1722 *   for (unsigned int a = 0; a < dim; ++a)
1723 *   lambdas_e_1_iso[a] = lambdas_e_1[a]*std::pow(J_e_1,-1.0/dim);
1724 *  
1725 *   for (unsigned int a = 0; a < dim; ++a)
1726 *   {
1727 *   SymmetricTensor<2, dim, NumberType>
1728 *   B_e_1_aux = symmetrize(outer_product(eigen_B_e_1_tr[a].second,eigen_B_e_1_tr[a].second));
1729 *   B_e_1_aux *= lambdas_e_1[a] * lambdas_e_1[a];
1730 *   B_e_1 += B_e_1_aux;
1731 *   }
1732 *  
1733 *   Tensor<2, dim, NumberType>Cinv_v_1_AD = symmetrize(invert(F) * B_e_1 * invert(transpose(F)));
1734 *  
1735 *   this->tau_neq_1 = 0;
1736 *   for (unsigned int a = 0; a < dim; ++a)
1737 *   {
1738 *   SymmetricTensor<2, dim, NumberType>
1739 *   tau_neq_1_aux = symmetrize(outer_product(eigen_B_e_1_tr[a].second,eigen_B_e_1_tr[a].second));
1740 *   tau_neq_1_aux *= get_beta_mode_1(lambdas_e_1_iso, a);
1741 *   this->tau_neq_1 += tau_neq_1_aux;
1742 *   }
1743 *  
1744 * @endcode
1745 *
1746 * Store history
1747 *
1748 * @code
1749 *   for (unsigned int a = 0; a < dim; ++a)
1750 *   for (unsigned int b = 0; b < dim; ++b)
1751 *   this->Cinv_v_1[a][b]= Tensor<0,dim,double>(Cinv_v_1_AD[a][b]);
1752 *   }
1753 *  
1754 *   void update_end_timestep() override
1755 *   {
1756 *   Material_Hyperelastic < dim, NumberType >::update_end_timestep();
1757 *   this->Cinv_v_1_converged = this->Cinv_v_1;
1758 *   }
1759 *  
1760 *   double get_viscous_dissipation() const override
1761 *   {
1762 *   NumberType dissipation_term = get_tau_E_neq() * get_tau_E_neq(); //Double contract the two SymmetricTensor
1763 *   dissipation_term /= (2*viscosity_mode_1);
1764 *  
1765 *   return dissipation_term.val();
1766 *   }
1767 *  
1768 *   protected:
1769 *   std::vector<double> mu_infty;
1770 *   std::vector<double> alpha_infty;
1771 *   std::vector<double> mu_mode_1;
1772 *   std::vector<double> alpha_mode_1;
1773 *   double viscosity_mode_1;
1774 *   SymmetricTensor<2, dim, double> Cinv_v_1;
1775 *   SymmetricTensor<2, dim, double> Cinv_v_1_converged;
1776 *   SymmetricTensor<2, dim, NumberType> tau_neq_1;
1777 *  
1778 *   SymmetricTensor<2, dim, NumberType>
1779 *   get_tau_E_base(const Tensor<2,dim, NumberType> &F) const override
1780 *   {
1781 *   return ( get_tau_E_neq() + get_tau_E_eq(F) );
1782 *   }
1783 *  
1784 *   SymmetricTensor<2, dim, NumberType>
1785 *   get_tau_E_eq(const Tensor<2,dim, NumberType> &F) const
1786 *   {
1787 *   const SymmetricTensor<2, dim, NumberType> B = symmetrize(F * transpose(F));
1788 *  
1789 *   std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim > eigen_B;
1790 *   eigen_B = eigenvectors(B, this->eigen_solver);
1791 *  
1792 *   SymmetricTensor<2, dim, NumberType> tau;
1793 *   static const SymmetricTensor< 2, dim, double>
1794 *   I (Physics::Elasticity::StandardTensors<dim>::I);
1795 *  
1796 *   for (unsigned int i = 0; i < 3; ++i)
1797 *   {
1798 *   for (unsigned int A = 0; A < dim; ++A)
1799 *   {
1800 *   SymmetricTensor<2, dim, NumberType> tau_aux1 = symmetrize(
1801 *   outer_product(eigen_B[A].second,eigen_B[A].second));
1802 *   tau_aux1 *= mu_infty[i]*std::pow(eigen_B[A].first, (alpha_infty[i]/2.) );
1803 *   tau += tau_aux1;
1804 *   }
1805 *   SymmetricTensor<2, dim, NumberType> tau_aux2 (I);
1806 *   tau_aux2 *= mu_infty[i];
1807 *   tau -= tau_aux2;
1808 *   }
1809 *   return tau;
1810 *   }
1811 *  
1812 *   SymmetricTensor<2, dim, NumberType>
1813 *   get_tau_E_neq() const
1814 *   {
1815 *   return tau_neq_1;
1816 *   }
1817 *  
1818 *   NumberType
1819 *   get_beta_mode_1(std::vector< NumberType > &lambda, const int &A) const
1820 *   {
1821 *   NumberType beta = 0.0;
1822 *  
1823 *   for (unsigned int i = 0; i < 3; ++i) //3rd-order Ogden model
1824 *   {
1825 *  
1826 *   NumberType aux = 0.0;
1827 *   for (int p = 0; p < dim; ++p)
1828 *   aux += std::pow(lambda[p],alpha_mode_1[i]);
1829 *  
1830 *   aux *= -1.0/dim;
1831 *   aux += std::pow(lambda[A], alpha_mode_1[i]);
1832 *   aux *= mu_mode_1[i];
1833 *  
1834 *   beta += aux;
1835 *   }
1836 *   return beta;
1837 *   }
1838 *  
1839 *   NumberType
1840 *   get_gamma_mode_1(std::vector< NumberType > &lambda,
1841 *   const int &A,
1842 *   const int &B ) const
1843 *   {
1844 *   NumberType gamma = 0.0;
1845 *  
1846 *   if (A==B)
1847 *   {
1848 *   for (unsigned int i = 0; i < 3; ++i)
1849 *   {
1850 *   NumberType aux = 0.0;
1851 *   for (int p = 0; p < dim; ++p)
1852 *   aux += std::pow(lambda[p],alpha_mode_1[i]);
1853 *  
1854 *   aux *= 1.0/(dim*dim);
1855 *   aux += 1.0/dim * std::pow(lambda[A], alpha_mode_1[i]);
1856 *   aux *= mu_mode_1[i]*alpha_mode_1[i];
1857 *  
1858 *   gamma += aux;
1859 *   }
1860 *   }
1861 *   else
1862 *   {
1863 *   for (unsigned int i = 0; i < 3; ++i)
1864 *   {
1865 *   NumberType aux = 0.0;
1866 *   for (int p = 0; p < dim; ++p)
1867 *   aux += std::pow(lambda[p],alpha_mode_1[i]);
1868 *  
1869 *   aux *= 1.0/(dim*dim);
1870 *   aux -= 1.0/dim * std::pow(lambda[A], alpha_mode_1[i]);
1871 *   aux -= 1.0/dim * std::pow(lambda[B], alpha_mode_1[i]);
1872 *   aux *= mu_mode_1[i]*alpha_mode_1[i];
1873 *  
1874 *   gamma += aux;
1875 *   }
1876 *   }
1877 *  
1878 *   return gamma;
1879 *   }
1880 *   };
1881 *  
1882 *  
1883 * @endcode
1884 *
1885 *
1886 * <a name="Constitutiveequationforthefluidcomponentofthebiphasicmaterial"></a>
1887 * <h3>Constitutive equation for the fluid component of the biphasic material</h3>
1888 * We consider two slightly different definitions to define the seepage velocity with a Darcy-like law.
1889 * Ehlers & Eipper 1999, doi:10.1023/A:1006565509095
1890 * Markert 2007, doi:10.1007/s11242-007-9107-6
1891 * The selection of one or another is made by the user via the parameters file.
1892 *
1893 * @code
1894 *   template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1895 *   class Material_Darcy_Fluid
1896 *   {
1897 *   public:
1898 *   Material_Darcy_Fluid(const Parameters::AllParameters &parameters)
1899 *   :
1900 *   fluid_type(parameters.fluid_type),
1901 *   n_OS(parameters.solid_vol_frac),
1902 *   initial_intrinsic_permeability(parameters.init_intrinsic_perm),
1903 *   viscosity_FR(parameters.viscosity_FR),
1904 *   initial_darcy_coefficient(parameters.init_darcy_coef),
1905 *   weight_FR(parameters.weight_FR),
1906 *   kappa_darcy(parameters.kappa_darcy),
1907 *   gravity_term(parameters.gravity_term),
1908 *   density_FR(parameters.density_FR),
1909 *   gravity_direction(parameters.gravity_direction),
1910 *   gravity_value(parameters.gravity_value)
1911 *   {
1912 *   Assert(kappa_darcy >= 0, ExcInternalError());
1913 *   }
1914 *   ~Material_Darcy_Fluid()
1915 *   {}
1916 *  
1917 *   Tensor<1, dim, NumberType> get_seepage_velocity_current
1918 *   (const Tensor<2,dim, NumberType> &F,
1919 *   const Tensor<1,dim, NumberType> &grad_p_fluid) const
1920 *   {
1921 *   const NumberType det_F = determinant(F);
1922 *   Assert(det_F > 0.0, ExcInternalError());
1923 *  
1924 *   Tensor<2, dim, NumberType> permeability_term;
1925 *  
1926 *   if (fluid_type == "Markert")
1927 *   permeability_term = get_instrinsic_permeability_current(F) / viscosity_FR;
1928 *  
1929 *   else if (fluid_type == "Ehlers")
1930 *   permeability_term = get_darcy_flow_current(F) / weight_FR;
1931 *  
1932 *   else
1933 *   AssertThrow(false, ExcMessage(
1934 *   "Material_Darcy_Fluid --> Only Markert "
1935 *   "and Ehlers formulations have been implemented."));
1936 *  
1937 *   return ( -1.0 * permeability_term * det_F
1938 *   * (grad_p_fluid - get_body_force_FR_current()) );
1939 *   }
1940 *  
1941 *   double get_porous_dissipation(const Tensor<2,dim, NumberType> &F,
1942 *   const Tensor<1,dim, NumberType> &grad_p_fluid) const
1943 *   {
1944 *   NumberType dissipation_term;
1945 *   Tensor<1, dim, NumberType> seepage_velocity;
1946 *   Tensor<2, dim, NumberType> permeability_term;
1947 *  
1948 *   const NumberType det_F = determinant(F);
1949 *   Assert(det_F > 0.0, ExcInternalError());
1950 *  
1951 *   if (fluid_type == "Markert")
1952 *   {
1953 *   permeability_term = get_instrinsic_permeability_current(F) / viscosity_FR;
1954 *   seepage_velocity = get_seepage_velocity_current(F,grad_p_fluid);
1955 *   }
1956 *   else if (fluid_type == "Ehlers")
1957 *   {
1958 *   permeability_term = get_darcy_flow_current(F) / weight_FR;
1959 *   seepage_velocity = get_seepage_velocity_current(F,grad_p_fluid);
1960 *   }
1961 *   else
1962 *   AssertThrow(false, ExcMessage(
1963 *   "Material_Darcy_Fluid --> Only Markert and Ehlers "
1964 *   "formulations have been implemented."));
1965 *  
1966 *   dissipation_term = ( invert(permeability_term) * seepage_velocity ) * seepage_velocity;
1967 *   dissipation_term *= 1.0/(det_F*det_F);
1968 *   return Tensor<0,dim,double>(dissipation_term);
1969 *   }
1970 *  
1971 *   protected:
1972 *   const std::string fluid_type;
1973 *   const double n_OS;
1974 *   const double initial_intrinsic_permeability;
1975 *   const double viscosity_FR;
1976 *   const double initial_darcy_coefficient;
1977 *   const double weight_FR;
1978 *   const double kappa_darcy;
1979 *   const bool gravity_term;
1980 *   const double density_FR;
1981 *   const int gravity_direction;
1982 *   const double gravity_value;
1983 *  
1984 *   Tensor<2, dim, NumberType>
1985 *   get_instrinsic_permeability_current(const Tensor<2,dim, NumberType> &F) const
1986 *   {
1987 *   static const SymmetricTensor< 2, dim, double>
1988 *   I (Physics::Elasticity::StandardTensors<dim>::I);
1989 *   const Tensor<2, dim, NumberType> initial_instrinsic_permeability_tensor
1990 *   = Tensor<2, dim, double>(initial_intrinsic_permeability * I);
1991 *  
1992 *   const NumberType det_F = determinant(F);
1993 *   Assert(det_F > 0.0, ExcInternalError());
1994 *  
1995 *   const NumberType fraction = (det_F - n_OS)/(1 - n_OS);
1996 *   return ( NumberType (std::pow(fraction, kappa_darcy))
1997 *   * initial_instrinsic_permeability_tensor );
1998 *   }
1999 *  
2000 *   Tensor<2, dim, NumberType>
2001 *   get_darcy_flow_current(const Tensor<2,dim, NumberType> &F) const
2002 *   {
2003 *   static const SymmetricTensor< 2, dim, double>
2004 *   I (Physics::Elasticity::StandardTensors<dim>::I);
2005 *   const Tensor<2, dim, NumberType> initial_darcy_flow_tensor
2006 *   = Tensor<2, dim, double>(initial_darcy_coefficient * I);
2007 *  
2008 *   const NumberType det_F = determinant(F);
2009 *   Assert(det_F > 0.0, ExcInternalError());
2010 *  
2011 *   const NumberType fraction = (1.0 - (n_OS / det_F) )/(1.0 - n_OS);
2012 *   return ( NumberType (std::pow(fraction, kappa_darcy))
2013 *   * initial_darcy_flow_tensor);
2014 *   }
2015 *  
2016 *   Tensor<1, dim, NumberType>
2017 *   get_body_force_FR_current() const
2018 *   {
2019 *   Tensor<1, dim, NumberType> body_force_FR_current;
2020 *  
2021 *   if (gravity_term == true)
2022 *   {
2023 *   Tensor<1, dim, NumberType> gravity_vector;
2024 *   gravity_vector[gravity_direction] = gravity_value;
2025 *   body_force_FR_current = density_FR * gravity_vector;
2026 *   }
2027 *   return body_force_FR_current;
2028 *   }
2029 *   };
2030 *  
2031 * @endcode
2032 *
2033 *
2034 * <a name="Quadraturepointhistory"></a>
2035 * <h3>Quadrature point history</h3>
2036 * As seen in @ref step_18 "step-18", the <code> PointHistory </code> class offers a method
2037 * for storing data at the quadrature points. Here each quadrature point
2038 * holds a pointer to a material description. Thus, different material models
2039 * can be used in different regions of the domain. Among other data, we
2040 * choose to store the ``extra" Kirchhoff stress @f$\boldsymbol{\tau}_E@f$ and
2041 * the dissipation values @f$\mathcal{D}_p@f$ and @f$\mathcal{D}_v@f$.
2042 *
2043 * @code
2044 *   template <int dim, typename NumberType = Sacado::Fad::DFad<double> > //double>
2045 *   class PointHistory
2046 *   {
2047 *   public:
2048 *   PointHistory()
2049 *   {}
2050 *  
2051 *   virtual ~PointHistory()
2052 *   {}
2053 *  
2054 *   void setup_lqp (const Parameters::AllParameters &parameters,
2055 *   const Time &time)
2056 *   {
2057 *   if (parameters.mat_type == "Neo-Hooke")
2058 *   solid_material.reset(new NeoHooke<dim,NumberType>(parameters,time));
2059 *   else if (parameters.mat_type == "Ogden")
2060 *   solid_material.reset(new Ogden<dim,NumberType>(parameters,time));
2061 *   else if (parameters.mat_type == "visco-Ogden")
2062 *   solid_material.reset(new visco_Ogden<dim,NumberType>(parameters,time));
2063 *   else
2064 *   Assert (false, ExcMessage("Material type not implemented"));
2065 *  
2066 *   fluid_material.reset(new Material_Darcy_Fluid<dim,NumberType>(parameters));
2067 *   }
2068 *  
2069 *   SymmetricTensor<2, dim, NumberType>
2070 *   get_tau_E(const Tensor<2, dim, NumberType> &F) const
2071 *   {
2072 *   return solid_material->get_tau_E(F);
2073 *   }
2074 *  
2075 *   SymmetricTensor<2, dim, NumberType>
2076 *   get_Cauchy_E(const Tensor<2, dim, NumberType> &F) const
2077 *   {
2078 *   return solid_material->get_Cauchy_E(F);
2079 *   }
2080 *  
2081 *   double
2082 *   get_converged_det_F() const
2083 *   {
2084 *   return solid_material->get_converged_det_F();
2085 *   }
2086 *  
2087 *   void
2088 *   update_end_timestep()
2089 *   {
2090 *   solid_material->update_end_timestep();
2091 *   }
2092 *  
2093 *   void
2094 *   update_internal_equilibrium(const Tensor<2, dim, NumberType> &F )
2095 *   {
2096 *   solid_material->update_internal_equilibrium(F);
2097 *   }
2098 *  
2099 *   double
2100 *   get_viscous_dissipation() const
2101 *   {
2102 *   return solid_material->get_viscous_dissipation();
2103 *   }
2104 *  
2105 *   Tensor<1,dim, NumberType>
2106 *   get_seepage_velocity_current (const Tensor<2,dim, NumberType> &F,
2107 *   const Tensor<1,dim, NumberType> &grad_p_fluid) const
2108 *   {
2109 *   return fluid_material->get_seepage_velocity_current(F, grad_p_fluid);
2110 *   }
2111 *  
2112 *   double
2113 *   get_porous_dissipation(const Tensor<2,dim, NumberType> &F,
2114 *   const Tensor<1,dim, NumberType> &grad_p_fluid) const
2115 *   {
2116 *   return fluid_material->get_porous_dissipation(F, grad_p_fluid);
2117 *   }
2118 *  
2119 *   Tensor<1, dim, NumberType>
2120 *   get_overall_body_force (const Tensor<2,dim, NumberType> &F,
2121 *   const Parameters::AllParameters &parameters) const
2122 *   {
2123 *   Tensor<1, dim, NumberType> body_force;
2124 *  
2125 *   if (parameters.gravity_term == true)
2126 *   {
2127 *   const NumberType det_F_AD = determinant(F);
2128 *   Assert(det_F_AD > 0.0, ExcInternalError());
2129 *  
2130 *   const NumberType overall_density_ref
2131 *   = parameters.density_SR * parameters.solid_vol_frac
2132 *   + parameters.density_FR
2133 *   * (det_F_AD - parameters.solid_vol_frac);
2134 *  
2135 *   Tensor<1, dim, NumberType> gravity_vector;
2136 *   gravity_vector[parameters.gravity_direction] = parameters.gravity_value;
2137 *   body_force = overall_density_ref * gravity_vector;
2138 *   }
2139 *  
2140 *   return body_force;
2141 *   }
2142 *   private:
2143 *   std::shared_ptr< Material_Hyperelastic<dim, NumberType> > solid_material;
2144 *   std::shared_ptr< Material_Darcy_Fluid<dim, NumberType> > fluid_material;
2145 *   };
2146 *  
2147 * @endcode
2148 *
2149 *
2150 * <a name="Nonlinearporoviscoelasticsolid"></a>
2151 * <h3>Nonlinear poro-viscoelastic solid</h3>
2152 * The Solid class is the central class as it represents the problem at hand:
2153 * the nonlinear poro-viscoelastic solid
2154 *
2155 * @code
2156 *   template <int dim>
2157 *   class Solid
2158 *   {
2159 *   public:
2160 *   Solid(const Parameters::AllParameters &parameters);
2161 *   virtual ~Solid();
2162 *   void run();
2163 *  
2164 *   protected:
2165 *   using ADNumberType = Sacado::Fad::DFad<double>;
2166 *  
2167 *   std::ofstream outfile;
2168 *   std::ofstream pointfile;
2169 *  
2170 *   struct PerTaskData_ASM;
2171 *   template<typename NumberType = double> struct ScratchData_ASM;
2172 *  
2173 * @endcode
2174 *
2175 * Generate mesh
2176 *
2177 * @code
2178 *   virtual void make_grid() = 0;
2179 *  
2180 * @endcode
2181 *
2182 * Define points for post-processing
2183 *
2184 * @code
2185 *   virtual void define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) = 0;
2186 *  
2187 * @endcode
2188 *
2189 * Set up the finite element system to be solved:
2190 *
2191 * @code
2192 *   void system_setup(TrilinosWrappers::MPI::BlockVector &solution_delta_OUT);
2193 *  
2194 * @endcode
2195 *
2196 * Extract sub-blocks from the global matrix
2197 *
2198 * @code
2199 *   void determine_component_extractors();
2200 *  
2201 * @endcode
2202 *
2203 * Several functions to assemble the system and right hand side matrices using multithreading.
2204 *
2205 * @code
2206 *   void assemble_system
2207 *   (const TrilinosWrappers::MPI::BlockVector &solution_delta_OUT );
2208 *   void assemble_system_one_cell
2209 *   (const typename DoFHandler<dim>::active_cell_iterator &cell,
2210 *   ScratchData_ASM<ADNumberType> &scratch,
2211 *   PerTaskData_ASM &data) const;
2212 *   void copy_local_to_global_system(const PerTaskData_ASM &data);
2213 *  
2214 * @endcode
2215 *
2216 * Define boundary conditions
2217 *
2218 * @code
2219 *   virtual void make_constraints(const int &it_nr);
2220 *   virtual void make_dirichlet_constraints(AffineConstraints<double> &constraints) = 0;
2221 *   virtual Tensor<1,dim> get_neumann_traction
2222 *   (const types::boundary_id &boundary_id,
2223 *   const Point<dim> &pt,
2224 *   const Tensor<1,dim> &N) const = 0;
2225 *   virtual double get_prescribed_fluid_flow
2226 *   (const types::boundary_id &boundary_id,
2227 *   const Point<dim> &pt) const = 0;
2228 *   virtual types::boundary_id
2229 *   get_reaction_boundary_id_for_output () const = 0;
2230 *   virtual std::pair<types::boundary_id,types::boundary_id>
2231 *   get_drained_boundary_id_for_output () const = 0;
2232 *   virtual std::vector<double> get_dirichlet_load
2233 *   (const types::boundary_id &boundary_id,
2234 *   const int &direction) const = 0;
2235 *  
2236 * @endcode
2237 *
2238 * Create and update the quadrature points.
2239 *
2240 * @code
2241 *   void setup_qph();
2242 *  
2243 * @endcode
2244 *
2245 * Solve non-linear system using a Newton-Raphson scheme
2246 *
2247 * @code
2248 *   void solve_nonlinear_timestep(TrilinosWrappers::MPI::BlockVector &solution_delta_OUT);
2249 *  
2250 * @endcode
2251 *
2252 * Solve the linearized equations using a direct solver
2253 *
2254 * @code
2255 *   void solve_linear_system ( TrilinosWrappers::MPI::BlockVector &newton_update_OUT);
2256 *  
2257 * @endcode
2258 *
2259 * Retrieve the solution
2260 *
2261 * @code
2262 *   TrilinosWrappers::MPI::BlockVector
2263 *   get_total_solution(const TrilinosWrappers::MPI::BlockVector &solution_delta_IN) const;
2264 *  
2265 * @endcode
2266 *
2267 * Store the converged values of the internal variables at the end of each timestep
2268 *
2269 * @code
2270 *   void update_end_timestep();
2271 *  
2272 * @endcode
2273 *
2274 * Post-processing and writing data to files
2275 *
2276 * @code
2277 *   void output_results_to_vtu(const unsigned int timestep,
2278 *   const double current_time,
2279 *   TrilinosWrappers::MPI::BlockVector solution) const;
2280 *   void output_results_to_plot(const unsigned int timestep,
2281 *   const double current_time,
2282 *   TrilinosWrappers::MPI::BlockVector solution,
2283 *   std::vector<Point<dim> > &tracked_vertices,
2284 *   std::ofstream &pointfile) const;
2285 *  
2286 * @endcode
2287 *
2288 * Headers and footer for the output files
2289 *
2290 * @code
2291 *   void print_console_file_header( std::ofstream &outfile) const;
2292 *   void print_plot_file_header(std::vector<Point<dim> > &tracked_vertices,
2293 *   std::ofstream &pointfile) const;
2294 *   void print_console_file_footer(std::ofstream &outfile) const;
2295 *   void print_plot_file_footer( std::ofstream &pointfile) const;
2296 *  
2297 * @endcode
2298 *
2299 * For parallel communication
2300 *
2301 * @code
2302 *   MPI_Comm mpi_communicator;
2303 *   const unsigned int n_mpi_processes;
2304 *   const unsigned int this_mpi_process;
2305 *   mutable ConditionalOStream pcout;
2306 *  
2307 * @endcode
2308 *
2309 * A collection of the parameters used to describe the problem setup
2310 *
2311 * @code
2312 *   const Parameters::AllParameters &parameters;
2313 *  
2314 * @endcode
2315 *
2316 * Declare an instance of dealii Triangulation class (mesh)
2317 *
2318 * @code
2319 *   parallel::shared::Triangulation<dim> triangulation;
2320 *  
2321 * @endcode
2322 *
2323 * Keep track of the current time and the time spent evaluating certain functions
2324 *
2325 * @code
2326 *   Time time;
2327 *   TimerOutput timerconsole;
2328 *   TimerOutput timerfile;
2329 *  
2330 * @endcode
2331 *
2332 * A storage object for quadrature point information.
2333 *
2334 * @code
2335 *   CellDataStorage<typename Triangulation<dim>::cell_iterator, PointHistory<dim,ADNumberType> > quadrature_point_history;
2336 *  
2337 * @endcode
2338 *
2339 * Integers to store polynomial degree (needed for output)
2340 *
2341 * @code
2342 *   const unsigned int degree_displ;
2343 *   const unsigned int degree_pore;
2344 *  
2345 * @endcode
2346 *
2347 * Declare an instance of dealii FESystem class (finite element definition)
2348 *
2349 * @code
2350 *   const FESystem<dim> fe;
2351 *  
2352 * @endcode
2353 *
2354 * Declare an instance of dealii DoFHandler class (assign DoFs to mesh)
2355 *
2356 * @code
2357 *   DoFHandler<dim> dof_handler_ref;
2358 *  
2359 * @endcode
2360 *
2361 * Integer to store DoFs per element (this value will be used often)
2362 *
2363 * @code
2364 *   const unsigned int dofs_per_cell;
2365 *  
2366 * @endcode
2367 *
2368 * Declare an instance of dealii Extractor objects used to retrieve information from the solution vectors
2369 * We will use "u_fe" and "p_fluid_fe"as subscript in operator [] expressions on FEValues and FEFaceValues
2370 * objects to extract the components of the displacement vector and fluid pressure, respectively.
2371 *
2372 * @code
2373 *   const FEValuesExtractors::Vector u_fe;
2374 *   const FEValuesExtractors::Scalar p_fluid_fe;
2375 *  
2376 * @endcode
2377 *
2378 * Description of how the block-system is arranged. There are 3 blocks:
2379 * 0 - vector DOF displacements u
2380 * 1 - scalar DOF fluid pressure p_fluid
2381 *
2382 * @code
2383 *   static const unsigned int n_blocks = 2;
2384 *   static const unsigned int n_components = dim+1;
2385 *   static const unsigned int first_u_component = 0;
2386 *   static const unsigned int p_fluid_component = dim;
2387 *  
2388 *   enum
2389 *   {
2390 *   u_block = 0,
2391 *   p_fluid_block = 1
2392 *   };
2393 *  
2394 * @endcode
2395 *
2396 * Extractors
2397 *
2398 * @code
2399 *   const FEValuesExtractors::Scalar x_displacement;
2400 *   const FEValuesExtractors::Scalar y_displacement;
2401 *   const FEValuesExtractors::Scalar z_displacement;
2402 *   const FEValuesExtractors::Scalar pressure;
2403 *  
2404 * @endcode
2405 *
2406 * Block data
2407 *
2408 * @code
2409 *   std::vector<unsigned int> block_component;
2410 *  
2411 * @endcode
2412 *
2413 * DoF index data
2414 *
2415 * @code
2416 *   std::vector<IndexSet> all_locally_owned_dofs;
2417 *   IndexSet locally_owned_dofs;
2418 *   IndexSet locally_relevant_dofs;
2419 *   std::vector<IndexSet> locally_owned_partitioning;
2420 *   std::vector<IndexSet> locally_relevant_partitioning;
2421 *  
2422 *   std::vector<types::global_dof_index> dofs_per_block;
2423 *   std::vector<types::global_dof_index> element_indices_u;
2424 *   std::vector<types::global_dof_index> element_indices_p_fluid;
2425 *  
2426 * @endcode
2427 *
2428 * Declare an instance of dealii QGauss class (The Gauss-Legendre family of quadrature rules for numerical integration)
2429 * Gauss Points in element, with n quadrature points (in each space direction <dim> )
2430 *
2431 * @code
2432 *   const QGauss<dim> qf_cell;
2433 * @endcode
2434 *
2435 * Gauss Points on element faces (used for definition of BCs)
2436 *
2437 * @code
2438 *   const QGauss<dim - 1> qf_face;
2439 * @endcode
2440 *
2441 * Integer to store num GPs per element (this value will be used often)
2442 *
2443 * @code
2444 *   const unsigned int n_q_points;
2445 * @endcode
2446 *
2447 * Integer to store num GPs per face (this value will be used often)
2448 *
2449 * @code
2450 *   const unsigned int n_q_points_f;
2451 *  
2452 * @endcode
2453 *
2454 * Declare an instance of dealii AffineConstraints class (linear constraints on DoFs due to hanging nodes or BCs)
2455 *
2456 * @code
2457 *   AffineConstraints<double> constraints;
2458 *  
2459 * @endcode
2460 *
2461 * Declare an instance of dealii classes necessary for FE system set-up and assembly
2462 * Store elements of tangent matrix (indicated by SparsityPattern class) as sparse matrix (more efficient)
2463 *
2464 * @code
2465 *   TrilinosWrappers::BlockSparseMatrix tangent_matrix;
2466 *   TrilinosWrappers::BlockSparseMatrix tangent_matrix_preconditioner;
2467 * @endcode
2468 *
2469 * Right hand side vector of forces
2470 *
2471 * @code
2472 *   TrilinosWrappers::MPI::BlockVector system_rhs;
2473 * @endcode
2474 *
2475 * Total displacement values + pressure (accumulated solution to FE system)
2476 *
2477 * @code
2478 *   TrilinosWrappers::MPI::BlockVector solution_n;
2479 *  
2480 * @endcode
2481 *
2482 * Non-block system for the direct solver. We will copy the block system into these to solve the linearized system of equations.
2483 *
2484 * @code
2485 *   TrilinosWrappers::SparseMatrix tangent_matrix_nb;
2486 *   TrilinosWrappers::MPI::Vector system_rhs_nb;
2487 *  
2488 * @endcode
2489 *
2490 * We define variables to store norms and update norms and normalisation factors.
2491 *
2492 * @code
2493 *   struct Errors
2494 *   {
2495 *   Errors()
2496 *   :
2497 *   norm(1.0), u(1.0), p_fluid(1.0)
2498 *   {}
2499 *  
2500 *   void reset()
2501 *   {
2502 *   norm = 1.0;
2503 *   u = 1.0;
2504 *   p_fluid = 1.0;
2505 *   }
2506 *   void normalise(const Errors &rhs)
2507 *   {
2508 *   if (rhs.norm != 0.0)
2509 *   norm /= rhs.norm;
2510 *   if (rhs.u != 0.0)
2511 *   u /= rhs.u;
2512 *   if (rhs.p_fluid != 0.0)
2513 *   p_fluid /= rhs.p_fluid;
2514 *   }
2515 *  
2516 *   double norm, u, p_fluid;
2517 *   };
2518 *  
2519 * @endcode
2520 *
2521 * Declare several instances of the "Error" structure
2522 *
2523 * @code
2524 *   Errors error_residual, error_residual_0, error_residual_norm, error_update,
2525 *   error_update_0, error_update_norm;
2526 *  
2527 * @endcode
2528 *
2529 * Methods to calculate error measures
2530 *
2531 * @code
2532 *   void get_error_residual(Errors &error_residual_OUT);
2533 *   void get_error_update
2534 *   (const TrilinosWrappers::MPI::BlockVector &newton_update_IN,
2535 *   Errors &error_update_OUT);
2536 *  
2537 * @endcode
2538 *
2539 * Print information to screen
2540 *
2541 * @code
2542 *   void print_conv_header();
2543 *   void print_conv_footer();
2544 *  
2545 * @endcode
2546 *
2547 * NOTE: In all functions, we pass by reference (&), so these functions work on the original copy (not a clone copy),
2548 * modifying the input variables inside the functions will change them outside the function.
2549 *
2550 * @code
2551 *   };
2552 *  
2553 * @endcode
2554 *
2555 *
2556 * <a name="ImplementationofthecodeSolidcodeclass"></a>
2557 * <h3>Implementation of the <code>Solid</code> class</h3>
2558 *
2559 * <a name="Publicinterface"></a>
2560 * <h4>Public interface</h4>
2561 * We initialise the Solid class using data extracted from the parameter file.
2562 *
2563 * @code
2564 *   template <int dim>
2565 *   Solid<dim>::Solid(const Parameters::AllParameters &parameters)
2566 *   :
2567 *   mpi_communicator(MPI_COMM_WORLD),
2568 *   n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
2569 *   this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
2570 *   pcout(std::cout, this_mpi_process == 0),
2571 *   parameters(parameters),
2572 *   triangulation(mpi_communicator,Triangulation<dim>::maximum_smoothing),
2573 *   time(parameters.end_time, parameters.delta_t),
2574 *   timerconsole( mpi_communicator,
2575 *   pcout,
2576 *   TimerOutput::summary,
2577 *   TimerOutput::wall_times),
2578 *   timerfile( mpi_communicator,
2579 *   outfile,
2580 *   TimerOutput::summary,
2581 *   TimerOutput::wall_times),
2582 *   degree_displ(parameters.poly_degree_displ),
2583 *   degree_pore(parameters.poly_degree_pore),
2584 *   fe( FE_Q<dim>(parameters.poly_degree_displ), dim,
2585 *   FE_Q<dim>(parameters.poly_degree_pore), 1 ),
2586 *   dof_handler_ref(triangulation),
2587 *   dofs_per_cell (fe.dofs_per_cell),
2588 *   u_fe(first_u_component),
2589 *   p_fluid_fe(p_fluid_component),
2590 *   x_displacement(first_u_component),
2591 *   y_displacement(first_u_component+1),
2592 *   z_displacement(first_u_component+2),
2593 *   pressure(p_fluid_component),
2594 *   dofs_per_block(n_blocks),
2595 *   qf_cell(parameters.quad_order),
2596 *   qf_face(parameters.quad_order),
2597 *   n_q_points (qf_cell.size()),
2598 *   n_q_points_f (qf_face.size())
2599 *   {
2600 *   Assert(dim==3, ExcMessage("This problem only works in 3 space dimensions."));
2601 *   determine_component_extractors();
2602 *   }
2603 *  
2604 * @endcode
2605 *
2606 * The class destructor simply clears the data held by the DOFHandler
2607 *
2608 * @code
2609 *   template <int dim>
2610 *   Solid<dim>::~Solid()
2611 *   {
2612 *   dof_handler_ref.clear();
2613 *   }
2614 *  
2615 * @endcode
2616 *
2617 * Runs the 3D solid problem
2618 *
2619 * @code
2620 *   template <int dim>
2621 *   void Solid<dim>::run()
2622 *   {
2623 * @endcode
2624 *
2625 * The current solution increment is defined as a block vector to reflect the structure
2626 * of the PDE system, with multiple solution components
2627 *
2628 * @code
2629 *   TrilinosWrappers::MPI::BlockVector solution_delta;
2630 *  
2631 * @endcode
2632 *
2633 * Open file
2634 *
2635 * @code
2636 *   if (this_mpi_process == 0)
2637 *   {
2638 *   outfile.open("console-output.sol");
2639 *   print_console_file_header(outfile);
2640 *   }
2641 *  
2642 * @endcode
2643 *
2644 * Generate mesh
2645 *
2646 * @code
2647 *   make_grid();
2648 *  
2649 * @endcode
2650 *
2651 * Assign DOFs and create the stiffness and right-hand-side force vector
2652 *
2653 * @code
2654 *   system_setup(solution_delta);
2655 *  
2656 * @endcode
2657 *
2658 * Define points for post-processing
2659 *
2660 * @code
2661 *   std::vector<Point<dim> > tracked_vertices (2);
2662 *   define_tracked_vertices(tracked_vertices);
2663 *   std::vector<Point<dim>> reaction_force;
2664 *  
2665 *   if (this_mpi_process == 0)
2666 *   {
2667 *   pointfile.open("data-for-gnuplot.sol");
2668 *   print_plot_file_header(tracked_vertices, pointfile);
2669 *   }
2670 *  
2671 * @endcode
2672 *
2673 * Print results to output file
2674 *
2675 * @code
2676 *   if (parameters.outfiles_requested == "true")
2677 *   {
2678 *   output_results_to_vtu(time.get_timestep(),
2679 *   time.get_current(),
2680 *   solution_n );
2681 *   }
2682 *  
2683 *   output_results_to_plot(time.get_timestep(),
2684 *   time.get_current(),
2685 *   solution_n,
2686 *   tracked_vertices,
2687 *   pointfile);
2688 *  
2689 * @endcode
2690 *
2691 * Increment time step (=load step)
2692 * NOTE: In solving the quasi-static problem, the time becomes a loading parameter,
2693 * i.e. we increase the loading linearly with time, making the two concepts interchangeable.
2694 *
2695 * @code
2696 *   time.increment_time();
2697 *  
2698 * @endcode
2699 *
2700 * Print information on screen
2701 *
2702 * @code
2703 *   pcout << "\nSolver:";
2704 *   pcout << "\n CST = make constraints";
2705 *   pcout << "\n ASM_SYS = assemble system";
2706 *   pcout << "\n SLV = linear solver \n";
2707 *  
2708 * @endcode
2709 *
2710 * Print information on file
2711 *
2712 * @code
2713 *   outfile << "\nSolver:";
2714 *   outfile << "\n CST = make constraints";
2715 *   outfile << "\n ASM_SYS = assemble system";
2716 *   outfile << "\n SLV = linear solver \n";
2717 *  
2718 *   while ( (time.get_end() - time.get_current()) > -1.0*parameters.tol_u )
2719 *   {
2720 * @endcode
2721 *
2722 * Initialize the current solution increment to zero
2723 *
2724 * @code
2725 *   solution_delta = 0.0;
2726 *  
2727 * @endcode
2728 *
2729 * Solve the non-linear system using a Newton-Rapshon scheme
2730 *
2731 * @code
2732 *   solve_nonlinear_timestep(solution_delta);
2733 *  
2734 * @endcode
2735 *
2736 * Add the computed solution increment to total solution
2737 *
2738 * @code
2739 *   solution_n += solution_delta;
2740 *  
2741 * @endcode
2742 *
2743 * Store the converged values of the internal variables
2744 *
2745 * @code
2746 *   update_end_timestep();
2747 *  
2748 * @endcode
2749 *
2750 * Output results
2751 *
2752 * @code
2753 *   if (( (time.get_timestep()%parameters.timestep_output) == 0 )
2754 *   && (parameters.outfiles_requested == "true") )
2755 *   {
2756 *   output_results_to_vtu(time.get_timestep(),
2757 *   time.get_current(),
2758 *   solution_n );
2759 *   }
2760 *  
2761 *   output_results_to_plot(time.get_timestep(),
2762 *   time.get_current(),
2763 *   solution_n,
2764 *   tracked_vertices,
2765 *   pointfile);
2766 *  
2767 * @endcode
2768 *
2769 * Increment the time step (=load step)
2770 *
2771 * @code
2772 *   time.increment_time();
2773 *   }
2774 *  
2775 * @endcode
2776 *
2777 * Print the footers and close files
2778 *
2779 * @code
2780 *   if (this_mpi_process == 0)
2781 *   {
2782 *   print_plot_file_footer(pointfile);
2783 *   pointfile.close ();
2784 *   print_console_file_footer(outfile);
2785 *  
2786 * @endcode
2787 *
2788 * NOTE: ideally, we should close the outfile here [ >> outfile.close (); ]
2789 * But if we do, then the timer output will not be printed. That is why we leave it open.
2790 *
2791 * @code
2792 *   }
2793 *   }
2794 *  
2795 * @endcode
2796 *
2797 *
2798 * <a name="Privateinterface"></a>
2799 * <h4>Private interface</h4>
2800 * We define the structures needed for parallelization with Threading Building Blocks (TBB)
2801 * Tangent matrix and right-hand side force vector assembly structures.
2802 * PerTaskData_ASM stores local contributions
2803 *
2804 * @code
2805 *   template <int dim>
2806 *   struct Solid<dim>::PerTaskData_ASM
2807 *   {
2808 *   FullMatrix<double> cell_matrix;
2809 *   Vector<double> cell_rhs;
2810 *   std::vector<types::global_dof_index> local_dof_indices;
2811 *  
2812 *   PerTaskData_ASM(const unsigned int dofs_per_cell)
2813 *   :
2814 *   cell_matrix(dofs_per_cell, dofs_per_cell),
2815 *   cell_rhs(dofs_per_cell),
2816 *   local_dof_indices(dofs_per_cell)
2817 *   {}
2818 *  
2819 *   void reset()
2820 *   {
2821 *   cell_matrix = 0.0;
2822 *   cell_rhs = 0.0;
2823 *   }
2824 *   };
2825 *  
2826 * @endcode
2827 *
2828 * ScratchData_ASM stores larger objects used during the assembly
2829 *
2830 * @code
2831 *   template <int dim>
2832 *   template <typename NumberType>
2833 *   struct Solid<dim>::ScratchData_ASM
2834 *   {
2835 *   const TrilinosWrappers::MPI::BlockVector &solution_total;
2836 *  
2837 * @endcode
2838 *
2839 * Integration helper
2840 *
2841 * @code
2842 *   FEValues<dim> fe_values_ref;
2843 *   FEFaceValues<dim> fe_face_values_ref;
2844 *  
2845 * @endcode
2846 *
2847 * Quadrature point solution
2848 *
2849 * @code
2850 *   std::vector<NumberType> local_dof_values;
2851 *   std::vector<Tensor<2, dim, NumberType> > solution_grads_u_total;
2852 *   std::vector<NumberType> solution_values_p_fluid_total;
2853 *   std::vector<Tensor<1, dim, NumberType> > solution_grads_p_fluid_total;
2854 *   std::vector<Tensor<1, dim, NumberType> > solution_grads_face_p_fluid_total;
2855 *  
2856 * @endcode
2857 *
2858 * shape function values
2859 *
2860 * @code
2861 *   std::vector<std::vector<Tensor<1,dim>>> Nx;
2862 *   std::vector<std::vector<double>> Nx_p_fluid;
2863 * @endcode
2864 *
2865 * shape function gradients
2866 *
2867 * @code
2868 *   std::vector<std::vector<Tensor<2,dim, NumberType>>> grad_Nx;
2869 *   std::vector<std::vector<SymmetricTensor<2,dim, NumberType>>> symm_grad_Nx;
2870 *   std::vector<std::vector<Tensor<1,dim, NumberType>>> grad_Nx_p_fluid;
2871 *  
2872 *   ScratchData_ASM(const FiniteElement<dim> &fe_cell,
2873 *   const QGauss<dim> &qf_cell, const UpdateFlags uf_cell,
2874 *   const QGauss<dim - 1> & qf_face, const UpdateFlags uf_face,
2875 *   const TrilinosWrappers::MPI::BlockVector &solution_total )
2876 *   :
2877 *   solution_total (solution_total),
2878 *   fe_values_ref(fe_cell, qf_cell, uf_cell),
2879 *   fe_face_values_ref(fe_cell, qf_face, uf_face),
2880 *   local_dof_values(fe_cell.dofs_per_cell),
2881 *   solution_grads_u_total(qf_cell.size()),
2882 *   solution_values_p_fluid_total(qf_cell.size()),
2883 *   solution_grads_p_fluid_total(qf_cell.size()),
2884 *   solution_grads_face_p_fluid_total(qf_face.size()),
2885 *   Nx(qf_cell.size(), std::vector<Tensor<1,dim>>(fe_cell.dofs_per_cell)),
2886 *   Nx_p_fluid(qf_cell.size(), std::vector<double>(fe_cell.dofs_per_cell)),
2887 *   grad_Nx(qf_cell.size(), std::vector<Tensor<2, dim, NumberType>>(fe_cell.dofs_per_cell)),
2888 *   symm_grad_Nx(qf_cell.size(), std::vector<SymmetricTensor<2, dim, NumberType>> (fe_cell.dofs_per_cell)),
2889 *   grad_Nx_p_fluid(qf_cell.size(), std::vector<Tensor<1, dim, NumberType>>(fe_cell.dofs_per_cell))
2890 *   {}
2891 *  
2892 *   ScratchData_ASM(const ScratchData_ASM &rhs)
2893 *   :
2894 *   solution_total (rhs.solution_total),
2895 *   fe_values_ref(rhs.fe_values_ref.get_fe(),
2896 *   rhs.fe_values_ref.get_quadrature(),
2897 *   rhs.fe_values_ref.get_update_flags()),
2898 *   fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
2899 *   rhs.fe_face_values_ref.get_quadrature(),
2900 *   rhs.fe_face_values_ref.get_update_flags()),
2901 *   local_dof_values(rhs.local_dof_values),
2902 *   solution_grads_u_total(rhs.solution_grads_u_total),
2903 *   solution_values_p_fluid_total(rhs.solution_values_p_fluid_total),
2904 *   solution_grads_p_fluid_total(rhs.solution_grads_p_fluid_total),
2905 *   solution_grads_face_p_fluid_total(rhs.solution_grads_face_p_fluid_total),
2906 *   Nx(rhs.Nx),
2907 *   Nx_p_fluid(rhs.Nx_p_fluid),
2908 *   grad_Nx(rhs.grad_Nx),
2909 *   symm_grad_Nx(rhs.symm_grad_Nx),
2910 *   grad_Nx_p_fluid(rhs.grad_Nx_p_fluid)
2911 *   {}
2912 *  
2913 *   void reset()
2914 *   {
2915 *   const unsigned int n_q_points = Nx_p_fluid.size();
2916 *   const unsigned int n_dofs_per_cell = Nx_p_fluid[0].size();
2917 *  
2918 *   Assert(local_dof_values.size() == n_dofs_per_cell, ExcInternalError());
2919 *  
2920 *   for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
2921 *   {
2922 *   local_dof_values[k] = 0.0;
2923 *   }
2924 *  
2925 *   Assert(solution_grads_u_total.size() == n_q_points, ExcInternalError());
2926 *   Assert(solution_values_p_fluid_total.size() == n_q_points, ExcInternalError());
2927 *   Assert(solution_grads_p_fluid_total.size() == n_q_points, ExcInternalError());
2928 *  
2929 *   Assert(Nx.size() == n_q_points, ExcInternalError());
2930 *   Assert(grad_Nx.size() == n_q_points, ExcInternalError());
2931 *   Assert(symm_grad_Nx.size() == n_q_points, ExcInternalError());
2932 *  
2933 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2934 *   {
2935 *   Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
2936 *   Assert( grad_Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
2937 *   Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
2938 *  
2939 *   solution_grads_u_total[q_point] = 0.0;
2940 *   solution_values_p_fluid_total[q_point] = 0.0;
2941 *   solution_grads_p_fluid_total[q_point] = 0.0;
2942 *  
2943 *   for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
2944 *   {
2945 *   Nx[q_point][k] = 0.0;
2946 *   Nx_p_fluid[q_point][k] = 0.0;
2947 *   grad_Nx[q_point][k] = 0.0;
2948 *   symm_grad_Nx[q_point][k] = 0.0;
2949 *   grad_Nx_p_fluid[q_point][k] = 0.0;
2950 *   }
2951 *   }
2952 *  
2953 *   const unsigned int n_f_q_points = solution_grads_face_p_fluid_total.size();
2954 *   Assert(solution_grads_face_p_fluid_total.size() == n_f_q_points, ExcInternalError());
2955 *  
2956 *   for (unsigned int f_q_point = 0; f_q_point < n_f_q_points; ++f_q_point)
2957 *   solution_grads_face_p_fluid_total[f_q_point] = 0.0;
2958 *   }
2959 *   };
2960 *  
2961 * @endcode
2962 *
2963 * Define the boundary conditions on the mesh
2964 *
2965 * @code
2966 *   template <int dim>
2967 *   void Solid<dim>::make_constraints(const int &it_nr_IN)
2968 *   {
2969 *   pcout << " CST " << std::flush;
2970 *   outfile << " CST " << std::flush;
2971 *  
2972 *   if (it_nr_IN > 1) return;
2973 *  
2974 *   const bool apply_dirichlet_bc = (it_nr_IN == 0);
2975 *  
2976 *   if (apply_dirichlet_bc)
2977 *   {
2978 *   constraints.clear();
2979 *   make_dirichlet_constraints(constraints);
2980 *   }
2981 *   else
2982 *   {
2983 *   for (unsigned int i=0; i<dof_handler_ref.n_dofs(); ++i)
2984 *   if (constraints.is_inhomogeneously_constrained(i) == true)
2985 *   constraints.set_inhomogeneity(i,0.0);
2986 *   }
2987 *   constraints.close();
2988 *   }
2989 *  
2990 * @endcode
2991 *
2992 * Set-up the FE system
2993 *
2994 * @code
2995 *   template <int dim>
2996 *   void Solid<dim>::system_setup(TrilinosWrappers::MPI::BlockVector &solution_delta_OUT)
2997 *   {
2998 *   timerconsole.enter_subsection("Setup system");
2999 *   timerfile.enter_subsection("Setup system");
3000 *  
3001 * @endcode
3002 *
3003 * Determine number of components per block
3004 *
3005 * @code
3006 *   std::vector<unsigned int> block_component(n_components, u_block);
3007 *   block_component[p_fluid_component] = p_fluid_block;
3008 *  
3009 * @endcode
3010 *
3011 * The DOF handler is initialised and we renumber the grid in an efficient manner.
3012 *
3013 * @code
3014 *   dof_handler_ref.distribute_dofs(fe);
3015 *   DoFRenumbering::Cuthill_McKee(dof_handler_ref);
3016 *   DoFRenumbering::component_wise(dof_handler_ref, block_component);
3017 *  
3018 * @endcode
3019 *
3020 * Count the number of DoFs in each block
3021 *
3022 * @code
3023 *   dofs_per_block = DoFTools::count_dofs_per_fe_block(dof_handler_ref, block_component);
3024 *  
3025 * @endcode
3026 *
3027 * Setup the sparsity pattern and tangent matrix
3028 *
3029 * @code
3030 *   all_locally_owned_dofs = DoFTools::locally_owned_dofs_per_subdomain (dof_handler_ref);
3031 *   std::vector<IndexSet> all_locally_relevant_dofs
3032 *   = DoFTools::locally_relevant_dofs_per_subdomain (dof_handler_ref);
3033 *  
3034 *   locally_owned_dofs.clear();
3035 *   locally_owned_partitioning.clear();
3036 *   Assert(all_locally_owned_dofs.size() > this_mpi_process, ExcInternalError());
3037 *   locally_owned_dofs = all_locally_owned_dofs[this_mpi_process];
3038 *  
3039 *   locally_relevant_dofs.clear();
3040 *   locally_relevant_partitioning.clear();
3041 *   Assert(all_locally_relevant_dofs.size() > this_mpi_process, ExcInternalError());
3042 *   locally_relevant_dofs = all_locally_relevant_dofs[this_mpi_process];
3043 *  
3044 *   locally_owned_partitioning.reserve(n_blocks);
3045 *   locally_relevant_partitioning.reserve(n_blocks);
3046 *  
3047 *   for (unsigned int b=0; b<n_blocks; ++b)
3048 *   {
3049 *   const types::global_dof_index idx_begin
3050 *   = std::accumulate(dofs_per_block.begin(),
3051 *   std::next(dofs_per_block.begin(),b), 0);
3052 *   const types::global_dof_index idx_end
3053 *   = std::accumulate(dofs_per_block.begin(),
3054 *   std::next(dofs_per_block.begin(),b+1), 0);
3055 *   locally_owned_partitioning.push_back(locally_owned_dofs.get_view(idx_begin, idx_end));
3056 *   locally_relevant_partitioning.push_back(locally_relevant_dofs.get_view(idx_begin, idx_end));
3057 *   }
3058 *  
3059 * @endcode
3060 *
3061 * Print information on screen
3062 *
3063 * @code
3064 *   pcout << "\nTriangulation:\n"
3065 *   << " Number of active cells: "
3066 *   << triangulation.n_active_cells()
3067 *   << " (by partition:";
3068 *   for (unsigned int p=0; p<n_mpi_processes; ++p)
3069 *   pcout << (p==0 ? ' ' : '+')
3070 *   << (GridTools::count_cells_with_subdomain_association (triangulation,p));
3071 *   pcout << ")"
3072 *   << std::endl;
3073 *   pcout << " Number of degrees of freedom: "
3074 *   << dof_handler_ref.n_dofs()
3075 *   << " (by partition:";
3076 *   for (unsigned int p=0; p<n_mpi_processes; ++p)
3077 *   pcout << (p==0 ? ' ' : '+')
3078 *   << (DoFTools::count_dofs_with_subdomain_association (dof_handler_ref,p));
3079 *   pcout << ")"
3080 *   << std::endl;
3081 *   pcout << " Number of degrees of freedom per block: "
3082 *   << "[n_u, n_p_fluid] = ["
3083 *   << dofs_per_block[u_block]
3084 *   << ", "
3085 *   << dofs_per_block[p_fluid_block]
3086 *   << "]"
3087 *   << std::endl;
3088 *  
3089 * @endcode
3090 *
3091 * Print information to file
3092 *
3093 * @code
3094 *   outfile << "\nTriangulation:\n"
3095 *   << " Number of active cells: "
3096 *   << triangulation.n_active_cells()
3097 *   << " (by partition:";
3098 *   for (unsigned int p=0; p<n_mpi_processes; ++p)
3099 *   outfile << (p==0 ? ' ' : '+')
3100 *   << (GridTools::count_cells_with_subdomain_association (triangulation,p));
3101 *   outfile << ")"
3102 *   << std::endl;
3103 *   outfile << " Number of degrees of freedom: "
3104 *   << dof_handler_ref.n_dofs()
3105 *   << " (by partition:";
3106 *   for (unsigned int p=0; p<n_mpi_processes; ++p)
3107 *   outfile << (p==0 ? ' ' : '+')
3108 *   << (DoFTools::count_dofs_with_subdomain_association (dof_handler_ref,p));
3109 *   outfile << ")"
3110 *   << std::endl;
3111 *   outfile << " Number of degrees of freedom per block: "
3112 *   << "[n_u, n_p_fluid] = ["
3113 *   << dofs_per_block[u_block]
3114 *   << ", "
3115 *   << dofs_per_block[p_fluid_block]
3116 *   << "]"
3117 *   << std::endl;
3118 *  
3119 * @endcode
3120 *
3121 * We optimise the sparsity pattern to reflect this structure and prevent
3122 * unnecessary data creation for the right-diagonal block components.
3123 *
3124 * @code
3125 *   Table<2, DoFTools::Coupling> coupling(n_components, n_components);
3126 *   for (unsigned int ii = 0; ii < n_components; ++ii)
3127 *   for (unsigned int jj = 0; jj < n_components; ++jj)
3128 *  
3129 * @endcode
3130 *
3131 * Identify "zero" matrix components of FE-system (The two components do not couple)
3132 *
3133 * @code
3134 *   if (((ii == p_fluid_component) && (jj < p_fluid_component))
3135 *   || ((ii < p_fluid_component) && (jj == p_fluid_component)) )
3136 *   coupling[ii][jj] = DoFTools::none;
3137 *  
3138 * @endcode
3139 *
3140 * The rest of components always couple
3141 *
3142 * @code
3143 *   else
3144 *   coupling[ii][jj] = DoFTools::always;
3145 *  
3146 *   TrilinosWrappers::BlockSparsityPattern bsp (locally_owned_partitioning,
3147 *   mpi_communicator);
3148 *  
3149 *   DoFTools::make_sparsity_pattern (dof_handler_ref, bsp, constraints,
3150 *   false, this_mpi_process);
3151 *   bsp.compress();
3152 *  
3153 * @endcode
3154 *
3155 * Reinitialize the (sparse) tangent matrix with the given sparsity pattern.
3156 *
3157 * @code
3158 *   tangent_matrix.reinit (bsp);
3159 *  
3160 * @endcode
3161 *
3162 * Initialize the right hand side and solution vectors with number of DoFs
3163 *
3164 * @code
3165 *   system_rhs.reinit(locally_owned_partitioning, mpi_communicator);
3166 *   solution_n.reinit(locally_owned_partitioning, mpi_communicator);
3167 *   solution_delta_OUT.reinit(locally_owned_partitioning, mpi_communicator);
3168 *  
3169 * @endcode
3170 *
3171 * Non-block system
3172 *
3173 * @code
3174 *   TrilinosWrappers::SparsityPattern sp (locally_owned_dofs,
3175 *   mpi_communicator);
3176 *   DoFTools::make_sparsity_pattern (dof_handler_ref, sp, constraints,
3177 *   false, this_mpi_process);
3178 *   sp.compress();
3179 *   tangent_matrix_nb.reinit (sp);
3180 *   system_rhs_nb.reinit(locally_owned_dofs, mpi_communicator);
3181 *  
3182 * @endcode
3183 *
3184 * Set up the quadrature point history
3185 *
3186 * @code
3187 *   setup_qph();
3188 *  
3189 *   timerconsole.leave_subsection();
3190 *   timerfile.leave_subsection();
3191 *   }
3192 *  
3193 * @endcode
3194 *
3195 * Component extractors: used to extract sub-blocks from the global matrix
3196 * Description of which local element DOFs are attached to which block component
3197 *
3198 * @code
3199 *   template <int dim>
3200 *   void Solid<dim>::determine_component_extractors()
3201 *   {
3202 *   element_indices_u.clear();
3203 *   element_indices_p_fluid.clear();
3204 *  
3205 *   for (unsigned int k = 0; k < fe.dofs_per_cell; ++k)
3206 *   {
3207 *   const unsigned int k_group = fe.system_to_base_index(k).first.first;
3208 *   if (k_group == u_block)
3209 *   element_indices_u.push_back(k);
3210 *   else if (k_group == p_fluid_block)
3211 *   element_indices_p_fluid.push_back(k);
3212 *   else
3213 *   {
3214 *   Assert(k_group <= p_fluid_block, ExcInternalError());
3215 *   }
3216 *   }
3217 *   }
3218 *  
3219 * @endcode
3220 *
3221 * Set-up quadrature point history (QPH) data objects
3222 *
3223 * @code
3224 *   template <int dim>
3225 *   void Solid<dim>::setup_qph()
3226 *   {
3227 *   pcout << "\nSetting up quadrature point data..." << std::endl;
3228 *   outfile << "\nSetting up quadrature point data..." << std::endl;
3229 *  
3230 * @endcode
3231 *
3232 * Create QPH data objects.
3233 *
3234 * @code
3235 *   quadrature_point_history.initialize(triangulation.begin_active(),
3236 *   triangulation.end(), n_q_points);
3237 *  
3238 * @endcode
3239 *
3240 * Setup the initial quadrature point data using the info stored in parameters
3241 *
3242 * @code
3243 *   FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
3244 *   cell (IteratorFilters::LocallyOwnedCell(),
3245 *   dof_handler_ref.begin_active()),
3246 *   endc (IteratorFilters::LocallyOwnedCell(),
3247 *   dof_handler_ref.end());
3248 *   for (; cell!=endc; ++cell)
3249 *   {
3250 *   Assert(cell->is_locally_owned(), ExcInternalError());
3251 *   Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
3252 *  
3253 *   const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >
3254 *   lqph = quadrature_point_history.get_data(cell);
3255 *   Assert(lqph.size() == n_q_points, ExcInternalError());
3256 *  
3257 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3258 *   lqph[q_point]->setup_lqp(parameters, time);
3259 *   }
3260 *   }
3261 *  
3262 * @endcode
3263 *
3264 * Solve the non-linear system using a Newton-Raphson scheme
3265 *
3266 * @code
3267 *   template <int dim>
3268 *   void Solid<dim>::solve_nonlinear_timestep(TrilinosWrappers::MPI::BlockVector &solution_delta_OUT)
3269 *   {
3270 * @endcode
3271 *
3272 * Print the load step
3273 *
3274 * @code
3275 *   pcout << std::endl
3276 *   << "\nTimestep "
3277 *   << time.get_timestep()
3278 *   << " @ "
3279 *   << time.get_current()
3280 *   << "s"
3281 *   << std::endl;
3282 *   outfile << std::endl
3283 *   << "\nTimestep "
3284 *   << time.get_timestep()
3285 *   << " @ "
3286 *   << time.get_current()
3287 *   << "s"
3288 *   << std::endl;
3289 *  
3290 * @endcode
3291 *
3292 * Declare newton_update vector (solution of a Newton iteration),
3293 * which must have as many positions as global DoFs.
3294 *
3295 * @code
3296 *   TrilinosWrappers::MPI::BlockVector newton_update
3297 *   (locally_owned_partitioning, mpi_communicator);
3298 *  
3299 * @endcode
3300 *
3301 * Reset the error storage objects
3302 *
3303 * @code
3304 *   error_residual.reset();
3305 *   error_residual_0.reset();
3306 *   error_residual_norm.reset();
3307 *   error_update.reset();
3308 *   error_update_0.reset();
3309 *   error_update_norm.reset();
3310 *  
3311 *   print_conv_header();
3312 *  
3313 * @endcode
3314 *
3315 * Declare and initialize iterator for the Newton-Raphson algorithm steps
3316 *
3317 * @code
3318 *   unsigned int newton_iteration = 0;
3319 *  
3320 * @endcode
3321 *
3322 * Iterate until error is below tolerance or max number iterations are reached
3323 *
3324 * @code
3325 *   while(newton_iteration < parameters.max_iterations_NR)
3326 *   {
3327 *   pcout << " " << std::setw(2) << newton_iteration << " " << std::flush;
3328 *   outfile << " " << std::setw(2) << newton_iteration << " " << std::flush;
3329 *  
3330 * @endcode
3331 *
3332 * Initialize global stiffness matrix and global force vector to zero
3333 *
3334 * @code
3335 *   tangent_matrix = 0.0;
3336 *   system_rhs = 0.0;
3337 *  
3338 *   tangent_matrix_nb = 0.0;
3339 *   system_rhs_nb = 0.0;
3340 *  
3341 * @endcode
3342 *
3343 * Apply boundary conditions
3344 *
3345 * @code
3346 *   make_constraints(newton_iteration);
3347 *   assemble_system(solution_delta_OUT);
3348 *  
3349 * @endcode
3350 *
3351 * Compute the rhs residual (error between external and internal forces in FE system)
3352 *
3353 * @code
3354 *   get_error_residual(error_residual);
3355 *  
3356 * @endcode
3357 *
3358 * error_residual in first iteration is stored to normalize posterior error measures
3359 *
3360 * @code
3361 *   if (newton_iteration == 0)
3362 *   error_residual_0 = error_residual;
3363 *  
3364 * @endcode
3365 *
3366 * Determine the normalised residual error
3367 *
3368 * @code
3369 *   error_residual_norm = error_residual;
3370 *   error_residual_norm.normalise(error_residual_0);
3371 *  
3372 * @endcode
3373 *
3374 * If both errors are below the tolerances, exit the loop.
3375 * We need to check the residual vector directly for convergence
3376 * in the load steps where no external forces or displacements are imposed.
3377 *
3378 * @code
3379 *   if ( ((newton_iteration > 0)
3380 *   && (error_update_norm.u <= parameters.tol_u)
3381 *   && (error_update_norm.p_fluid <= parameters.tol_p_fluid)
3382 *   && (error_residual_norm.u <= parameters.tol_f)
3383 *   && (error_residual_norm.p_fluid <= parameters.tol_f))
3384 *   || ( (newton_iteration > 0)
3385 *   && system_rhs.l2_norm() <= parameters.tol_f) )
3386 *   {
3387 *   pcout << "\n ***** CONVERGED! ***** "
3388 *   << system_rhs.l2_norm() << " "
3389 *   << " " << error_residual_norm.norm
3390 *   << " " << error_residual_norm.u
3391 *   << " " << error_residual_norm.p_fluid
3392 *   << " " << error_update_norm.norm
3393 *   << " " << error_update_norm.u
3394 *   << " " << error_update_norm.p_fluid
3395 *   << " " << std::endl;
3396 *   outfile << "\n ***** CONVERGED! ***** "
3397 *   << system_rhs.l2_norm() << " "
3398 *   << " " << error_residual_norm.norm
3399 *   << " " << error_residual_norm.u
3400 *   << " " << error_residual_norm.p_fluid
3401 *   << " " << error_update_norm.norm
3402 *   << " " << error_update_norm.u
3403 *   << " " << error_update_norm.p_fluid
3404 *   << " " << std::endl;
3405 *   print_conv_footer();
3406 *  
3407 *   break;
3408 *   }
3409 *  
3410 * @endcode
3411 *
3412 * Solve the linearized system
3413 *
3414 * @code
3415 *   solve_linear_system(newton_update);
3416 *   constraints.distribute(newton_update);
3417 *  
3418 * @endcode
3419 *
3420 * Compute the displacement error
3421 *
3422 * @code
3423 *   get_error_update(newton_update, error_update);
3424 *  
3425 * @endcode
3426 *
3427 * error_update in first iteration is stored to normalize posterior error measures
3428 *
3429 * @code
3430 *   if (newton_iteration == 0)
3431 *   error_update_0 = error_update;
3432 *  
3433 * @endcode
3434 *
3435 * Determine the normalised Newton update error
3436 *
3437 * @code
3438 *   error_update_norm = error_update;
3439 *   error_update_norm.normalise(error_update_0);
3440 *  
3441 * @endcode
3442 *
3443 * Determine the normalised residual error
3444 *
3445 * @code
3446 *   error_residual_norm = error_residual;
3447 *   error_residual_norm.normalise(error_residual_0);
3448 *  
3449 * @endcode
3450 *
3451 * Print error values
3452 *
3453 * @code
3454 *   pcout << " | " << std::fixed << std::setprecision(3)
3455 *   << std::setw(7) << std::scientific
3456 *   << system_rhs.l2_norm()
3457 *   << " " << error_residual_norm.norm
3458 *   << " " << error_residual_norm.u
3459 *   << " " << error_residual_norm.p_fluid
3460 *   << " " << error_update_norm.norm
3461 *   << " " << error_update_norm.u
3462 *   << " " << error_update_norm.p_fluid
3463 *   << " " << std::endl;
3464 *  
3465 *   outfile << " | " << std::fixed << std::setprecision(3)
3466 *   << std::setw(7) << std::scientific
3467 *   << system_rhs.l2_norm()
3468 *   << " " << error_residual_norm.norm
3469 *   << " " << error_residual_norm.u
3470 *   << " " << error_residual_norm.p_fluid
3471 *   << " " << error_update_norm.norm
3472 *   << " " << error_update_norm.u
3473 *   << " " << error_update_norm.p_fluid
3474 *   << " " << std::endl;
3475 *  
3476 * @endcode
3477 *
3478 * Update
3479 *
3480 * @code
3481 *   solution_delta_OUT += newton_update;
3482 *   newton_update = 0.0;
3483 *   newton_iteration++;
3484 *   }
3485 *  
3486 * @endcode
3487 *
3488 * If maximum allowed number of iterations for Newton algorithm are reached, print non-convergence message and abort program
3489 *
3490 * @code
3491 *   AssertThrow (newton_iteration < parameters.max_iterations_NR, ExcMessage("No convergence in nonlinear solver!"));
3492 *   }
3493 *  
3494 * @endcode
3495 *
3496 * Prints the header for convergence info on console
3497 *
3498 * @code
3499 *   template <int dim>
3500 *   void Solid<dim>::print_conv_header()
3501 *   {
3502 *   static const unsigned int l_width = 120;
3503 *  
3504 *   for (unsigned int i = 0; i < l_width; ++i)
3505 *   {
3506 *   pcout << "_";
3507 *   outfile << "_";
3508 *   }
3509 *  
3510 *   pcout << std::endl;
3511 *   outfile << std::endl;
3512 *  
3513 *   pcout << "\n SOLVER STEP | SYS_RES "
3514 *   << "RES_NORM RES_U RES_P "
3515 *   << "NU_NORM NU_U NU_P " << std::endl;
3516 *   outfile << "\n SOLVER STEP | SYS_RES "
3517 *   << "RES_NORM RES_U RES_P "
3518 *   << "NU_NORM NU_U NU_P " << std::endl;
3519 *  
3520 *   for (unsigned int i = 0; i < l_width; ++i)
3521 *   {
3522 *   pcout << "_";
3523 *   outfile << "_";
3524 *   }
3525 *   pcout << std::endl << std::endl;
3526 *   outfile << std::endl << std::endl;
3527 *   }
3528 *  
3529 * @endcode
3530 *
3531 * Prints the footer for convergence info on console
3532 *
3533 * @code
3534 *   template <int dim>
3535 *   void Solid<dim>::print_conv_footer()
3536 *   {
3537 *   static const unsigned int l_width = 120;
3538 *  
3539 *   for (unsigned int i = 0; i < l_width; ++i)
3540 *   {
3541 *   pcout << "_";
3542 *   outfile << "_";
3543 *   }
3544 *   pcout << std::endl << std::endl;
3545 *   outfile << std::endl << std::endl;
3546 *  
3547 *   pcout << "Relative errors:" << std::endl
3548 *   << "Displacement: "
3549 *   << error_update.u / error_update_0.u << std::endl
3550 *   << "Force (displ): "
3551 *   << error_residual.u / error_residual_0.u << std::endl
3552 *   << "Pore pressure: "
3553 *   << error_update.p_fluid / error_update_0.p_fluid << std::endl
3554 *   << "Force (pore): "
3555 *   << error_residual.p_fluid / error_residual_0.p_fluid << std::endl;
3556 *   outfile << "Relative errors:" << std::endl
3557 *   << "Displacement: "
3558 *   << error_update.u / error_update_0.u << std::endl
3559 *   << "Force (displ): "
3560 *   << error_residual.u / error_residual_0.u << std::endl
3561 *   << "Pore pressure: "
3562 *   << error_update.p_fluid / error_update_0.p_fluid << std::endl
3563 *   << "Force (pore): "
3564 *   << error_residual.p_fluid / error_residual_0.p_fluid << std::endl;
3565 *   }
3566 *  
3567 * @endcode
3568 *
3569 * Determine the true residual error for the problem
3570 *
3571 * @code
3572 *   template <int dim>
3573 *   void Solid<dim>::get_error_residual(Errors &error_residual_OUT)
3574 *   {
3575 *   TrilinosWrappers::MPI::BlockVector error_res(system_rhs);
3576 *   constraints.set_zero(error_res);
3577 *  
3578 *   error_residual_OUT.norm = error_res.l2_norm();
3579 *   error_residual_OUT.u = error_res.block(u_block).l2_norm();
3580 *   error_residual_OUT.p_fluid = error_res.block(p_fluid_block).l2_norm();
3581 *   }
3582 *  
3583 * @endcode
3584 *
3585 * Determine the true Newton update error for the problem
3586 *
3587 * @code
3588 *   template <int dim>
3589 *   void Solid<dim>::get_error_update
3590 *   (const TrilinosWrappers::MPI::BlockVector &newton_update_IN,
3591 *   Errors &error_update_OUT)
3592 *   {
3593 *   TrilinosWrappers::MPI::BlockVector error_ud(newton_update_IN);
3594 *   constraints.set_zero(error_ud);
3595 *  
3596 *   error_update_OUT.norm = error_ud.l2_norm();
3597 *   error_update_OUT.u = error_ud.block(u_block).l2_norm();
3598 *   error_update_OUT.p_fluid = error_ud.block(p_fluid_block).l2_norm();
3599 *   }
3600 *  
3601 * @endcode
3602 *
3603 * Compute the total solution, which is valid at any Newton step. This is required as, to reduce
3604 * computational error, the total solution is only updated at the end of the timestep.
3605 *
3606 * @code
3607 *   template <int dim>
3608 *   TrilinosWrappers::MPI::BlockVector
3609 *   Solid<dim>::get_total_solution(const TrilinosWrappers::MPI::BlockVector &solution_delta_IN) const
3610 *   {
3611 * @endcode
3612 *
3613 * Cell interpolation -> Ghosted vector
3614 *
3615 * @code
3616 *   TrilinosWrappers::MPI::BlockVector
3617 *   solution_total (locally_owned_partitioning,
3618 *   locally_relevant_partitioning,
3619 *   mpi_communicator,
3620 *   /*vector_writable = */ false);
3621 *   TrilinosWrappers::MPI::BlockVector tmp (solution_total);
3622 *   solution_total = solution_n;
3623 *   tmp = solution_delta_IN;
3624 *   solution_total += tmp;
3625 *   return solution_total;
3626 *   }
3627 *  
3628 * @endcode
3629 *
3630 * Compute elemental stiffness tensor and right-hand side force vector, and assemble into global ones
3631 *
3632 * @code
3633 *   template <int dim>
3634 *   void Solid<dim>::assemble_system( const TrilinosWrappers::MPI::BlockVector &solution_delta )
3635 *   {
3636 *   timerconsole.enter_subsection("Assemble system");
3637 *   timerfile.enter_subsection("Assemble system");
3638 *   pcout << " ASM_SYS " << std::flush;
3639 *   outfile << " ASM_SYS " << std::flush;
3640 *  
3641 *   const TrilinosWrappers::MPI::BlockVector solution_total(get_total_solution(solution_delta));
3642 *  
3643 * @endcode
3644 *
3645 * Info given to FEValues and FEFaceValues constructors, to indicate which data will be needed at each element.
3646 *
3647 * @code
3648 *   const UpdateFlags uf_cell(update_values |
3649 *   update_gradients |
3650 *   update_JxW_values);
3651 *   const UpdateFlags uf_face(update_values |
3652 *   update_gradients |
3653 *   update_normal_vectors |
3654 *   update_quadrature_points |
3655 *   update_JxW_values );
3656 *  
3657 * @endcode
3658 *
3659 * Setup a copy of the data structures required for the process and pass them, along with the
3660 * memory addresses of the assembly functions to the WorkStream object for processing
3661 *
3662 * @code
3663 *   PerTaskData_ASM per_task_data(dofs_per_cell);
3664 *   ScratchData_ASM<ADNumberType> scratch_data(fe, qf_cell, uf_cell,
3665 *   qf_face, uf_face,
3666 *   solution_total);
3667 *  
3668 *   FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
3669 *   cell (IteratorFilters::LocallyOwnedCell(),
3670 *   dof_handler_ref.begin_active()),
3671 *   endc (IteratorFilters::LocallyOwnedCell(),
3672 *   dof_handler_ref.end());
3673 *   for (; cell != endc; ++cell)
3674 *   {
3675 *   Assert(cell->is_locally_owned(), ExcInternalError());
3676 *   Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
3677 *  
3678 *   assemble_system_one_cell(cell, scratch_data, per_task_data);
3679 *   copy_local_to_global_system(per_task_data);
3680 *   }
3681 *   tangent_matrix.compress(VectorOperation::add);
3682 *   system_rhs.compress(VectorOperation::add);
3683 *  
3684 *   tangent_matrix_nb.compress(VectorOperation::add);
3685 *   system_rhs_nb.compress(VectorOperation::add);
3686 *  
3687 *   timerconsole.leave_subsection();
3688 *   timerfile.leave_subsection();
3689 *   }
3690 *  
3691 * @endcode
3692 *
3693 * Add the local elemental contribution to the global stiffness tensor
3694 * We do it twice, for the block and the non-block systems
3695 *
3696 * @code
3697 *   template <int dim>
3698 *   void Solid<dim>::copy_local_to_global_system (const PerTaskData_ASM &data)
3699 *   {
3700 *   constraints.distribute_local_to_global(data.cell_matrix,
3701 *   data.cell_rhs,
3702 *   data.local_dof_indices,
3703 *   tangent_matrix,
3704 *   system_rhs);
3705 *  
3706 *   constraints.distribute_local_to_global(data.cell_matrix,
3707 *   data.cell_rhs,
3708 *   data.local_dof_indices,
3709 *   tangent_matrix_nb,
3710 *   system_rhs_nb);
3711 *   }
3712 *  
3713 * @endcode
3714 *
3715 * Compute stiffness matrix and corresponding rhs for one element
3716 *
3717 * @code
3718 *   template <int dim>
3719 *   void Solid<dim>::assemble_system_one_cell
3720 *   (const typename DoFHandler<dim>::active_cell_iterator &cell,
3721 *   ScratchData_ASM<ADNumberType> &scratch,
3722 *   PerTaskData_ASM &data) const
3723 *   {
3724 *   Assert(cell->is_locally_owned(), ExcInternalError());
3725 *  
3726 *   data.reset();
3727 *   scratch.reset();
3728 *   scratch.fe_values_ref.reinit(cell);
3729 *   cell->get_dof_indices(data.local_dof_indices);
3730 *  
3731 * @endcode
3732 *
3733 * Setup automatic differentiation
3734 *
3735 * @code
3736 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
3737 *   {
3738 * @endcode
3739 *
3740 * Initialise the dofs for the cell using the current solution.
3741 *
3742 * @code
3743 *   scratch.local_dof_values[k] = scratch.solution_total[data.local_dof_indices[k]];
3744 * @endcode
3745 *
3746 * Mark this cell DoF as an independent variable
3747 *
3748 * @code
3749 *   scratch.local_dof_values[k].diff(k, dofs_per_cell);
3750 *   }
3751 *  
3752 * @endcode
3753 *
3754 * Update the quadrature point solution
3755 * Compute the values and gradients of the solution in terms of the AD variables
3756 *
3757 * @code
3758 *   for (unsigned int q = 0; q < n_q_points; ++q)
3759 *   {
3760 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
3761 *   {
3762 *   const unsigned int k_group = fe.system_to_base_index(k).first.first;
3763 *   if (k_group == u_block)
3764 *   {
3765 *   const Tensor<2, dim> Grad_Nx_u =
3766 *   scratch.fe_values_ref[u_fe].gradient(k, q);
3767 *   for (unsigned int dd = 0; dd < dim; ++dd)
3768 *   {
3769 *   for (unsigned int ee = 0; ee < dim; ++ee)
3770 *   {
3771 *   scratch.solution_grads_u_total[q][dd][ee]
3772 *   += scratch.local_dof_values[k] * Grad_Nx_u[dd][ee];
3773 *   }
3774 *   }
3775 *   }
3776 *   else if (k_group == p_fluid_block)
3777 *   {
3778 *   const double Nx_p = scratch.fe_values_ref[p_fluid_fe].value(k, q);
3779 *   const Tensor<1, dim> Grad_Nx_p =
3780 *   scratch.fe_values_ref[p_fluid_fe].gradient(k, q);
3781 *  
3782 *   scratch.solution_values_p_fluid_total[q]
3783 *   += scratch.local_dof_values[k] * Nx_p;
3784 *   for (unsigned int dd = 0; dd < dim; ++dd)
3785 *   {
3786 *   scratch.solution_grads_p_fluid_total[q][dd]
3787 *   += scratch.local_dof_values[k] * Grad_Nx_p[dd];
3788 *   }
3789 *   }
3790 *   else
3791 *   Assert(k_group <= p_fluid_block, ExcInternalError());
3792 *   }
3793 *   }
3794 *  
3795 * @endcode
3796 *
3797 * Set up pointer "lgph" to the PointHistory object of this element
3798 *
3799 * @code
3800 *   const std::vector<std::shared_ptr<const PointHistory<dim, ADNumberType> > >
3801 *   lqph = quadrature_point_history.get_data(cell);
3802 *   Assert(lqph.size() == n_q_points, ExcInternalError());
3803 *  
3804 *  
3805 * @endcode
3806 *
3807 * Precalculate the element shape function values and gradients
3808 *
3809 * @code
3810 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3811 *   {
3812 *   Tensor<2, dim, ADNumberType> F_AD = scratch.solution_grads_u_total[q_point];
3813 *   F_AD += Tensor<2, dim, double>(Physics::Elasticity::StandardTensors<dim>::I);
3814 *   Assert(determinant(F_AD) > 0, ExcMessage("Invalid deformation map"));
3815 *   const Tensor<2, dim, ADNumberType> F_inv_AD = invert(F_AD);
3816 *  
3817 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
3818 *   {
3819 *   const unsigned int i_group = fe.system_to_base_index(i).first.first;
3820 *  
3821 *   if (i_group == u_block)
3822 *   {
3823 *   scratch.Nx[q_point][i] =
3824 *   scratch.fe_values_ref[u_fe].value(i, q_point);
3825 *   scratch.grad_Nx[q_point][i] =
3826 *   scratch.fe_values_ref[u_fe].gradient(i, q_point)*F_inv_AD;
3827 *   scratch.symm_grad_Nx[q_point][i] =
3828 *   symmetrize(scratch.grad_Nx[q_point][i]);
3829 *   }
3830 *   else if (i_group == p_fluid_block)
3831 *   {
3832 *   scratch.Nx_p_fluid[q_point][i] =
3833 *   scratch.fe_values_ref[p_fluid_fe].value(i, q_point);
3834 *   scratch.grad_Nx_p_fluid[q_point][i] =
3835 *   scratch.fe_values_ref[p_fluid_fe].gradient(i, q_point)*F_inv_AD;
3836 *   }
3837 *   else
3838 *   Assert(i_group <= p_fluid_block, ExcInternalError());
3839 *   }
3840 *   }
3841 *  
3842 * @endcode
3843 *
3844 * Assemble the stiffness matrix and rhs vector
3845 *
3846 * @code
3847 *   std::vector<ADNumberType> residual_ad (dofs_per_cell, ADNumberType(0.0));
3848 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3849 *   {
3850 *   Tensor<2, dim, ADNumberType> F_AD = scratch.solution_grads_u_total[q_point];
3851 *   F_AD += Tensor<2, dim,double>(Physics::Elasticity::StandardTensors<dim>::I);
3852 *   const ADNumberType det_F_AD = determinant(F_AD);
3853 *  
3854 *   Assert(det_F_AD > 0, ExcInternalError());
3855 *   const Tensor<2, dim, ADNumberType> F_inv_AD = invert(F_AD); //inverse of def. gradient tensor
3856 *  
3857 *   const ADNumberType p_fluid = scratch.solution_values_p_fluid_total[q_point];
3858 *  
3859 *   {
3860 *   PointHistory<dim, ADNumberType> *lqph_q_point_nc =
3861 *   const_cast<PointHistory<dim, ADNumberType>*>(lqph[q_point].get());
3862 *   lqph_q_point_nc->update_internal_equilibrium(F_AD);
3863 *   }
3864 *  
3865 * @endcode
3866 *
3867 * Get some info from constitutive model of solid
3868 *
3869 * @code
3870 *   static const SymmetricTensor< 2, dim, double>
3871 *   I (Physics::Elasticity::StandardTensors<dim>::I);
3872 *   const SymmetricTensor<2, dim, ADNumberType>
3873 *   tau_E = lqph[q_point]->get_tau_E(F_AD);
3874 *   SymmetricTensor<2, dim, ADNumberType> tau_fluid_vol (I);
3875 *   tau_fluid_vol *= -1.0 * p_fluid * det_F_AD;
3876 *  
3877 * @endcode
3878 *
3879 * Get some info from constitutive model of fluid
3880 *
3881 * @code
3882 *   const ADNumberType det_F_aux = lqph[q_point]->get_converged_det_F();
3883 *   const double det_F_converged = Tensor<0,dim,double>(det_F_aux); //Needs to be double, not AD number
3884 *   const Tensor<1, dim, ADNumberType> overall_body_force
3885 *   = lqph[q_point]->get_overall_body_force(F_AD, parameters);
3886 *  
3887 * @endcode
3888 *
3889 * Define some aliases to make the assembly process easier to follow
3890 *
3891 * @code
3892 *   const std::vector<Tensor<1,dim>> &Nu = scratch.Nx[q_point];
3893 *   const std::vector<SymmetricTensor<2, dim, ADNumberType>>
3894 *   &symm_grad_Nu = scratch.symm_grad_Nx[q_point];
3895 *   const std::vector<double> &Np = scratch.Nx_p_fluid[q_point];
3896 *   const std::vector<Tensor<1, dim, ADNumberType> > &grad_Np
3897 *   = scratch.grad_Nx_p_fluid[q_point];
3898 *   const Tensor<1, dim, ADNumberType> grad_p
3899 *   = scratch.solution_grads_p_fluid_total[q_point]*F_inv_AD;
3900 *   const double JxW = scratch.fe_values_ref.JxW(q_point);
3901 *  
3902 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
3903 *   {
3904 *   const unsigned int i_group = fe.system_to_base_index(i).first.first;
3905 *  
3906 *   if (i_group == u_block)
3907 *   {
3908 *   residual_ad[i] += symm_grad_Nu[i] * ( tau_E + tau_fluid_vol ) * JxW;
3909 *   residual_ad[i] -= Nu[i] * overall_body_force * JxW;
3910 *   }
3911 *   else if (i_group == p_fluid_block)
3912 *   {
3913 *   const Tensor<1, dim, ADNumberType> seepage_vel_current
3914 *   = lqph[q_point]->get_seepage_velocity_current(F_AD, grad_p);
3915 *   residual_ad[i] += Np[i] * (det_F_AD - det_F_converged) * JxW;
3916 *   residual_ad[i] -= time.get_delta_t() * grad_Np[i]
3917 *   * seepage_vel_current * JxW;
3918 *   }
3919 *   else
3920 *   Assert(i_group <= p_fluid_block, ExcInternalError());
3921 *   }
3922 *   }
3923 *  
3924 * @endcode
3925 *
3926 * Assemble the Neumann contribution (external force contribution).
3927 *
3928 * @code
3929 *   for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face) //Loop over faces in element
3930 *   {
3931 *   if (cell->face(face)->at_boundary() == true)
3932 *   {
3933 *   scratch.fe_face_values_ref.reinit(cell, face);
3934 *  
3935 *   for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; ++f_q_point)
3936 *   {
3937 *   const Tensor<1, dim> &N
3938 *   = scratch.fe_face_values_ref.normal_vector(f_q_point);
3939 *   const Point<dim> &pt
3940 *   = scratch.fe_face_values_ref.quadrature_point(f_q_point);
3941 *   const Tensor<1, dim> traction
3942 *   = get_neumann_traction(cell->face(face)->boundary_id(), pt, N);
3943 *   const double flow
3944 *   = get_prescribed_fluid_flow(cell->face(face)->boundary_id(), pt);
3945 *  
3946 *   if ( (traction.norm() < 1e-12) && (std::abs(flow) < 1e-12) ) continue;
3947 *  
3948 *   const double JxW_f = scratch.fe_face_values_ref.JxW(f_q_point);
3949 *  
3950 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
3951 *   {
3952 *   const unsigned int i_group = fe.system_to_base_index(i).first.first;
3953 *  
3954 *   if ((i_group == u_block) && (traction.norm() > 1e-12))
3955 *   {
3956 *   const unsigned int component_i
3957 *   = fe.system_to_component_index(i).first;
3958 *   const double Nu_f
3959 *   = scratch.fe_face_values_ref.shape_value(i, f_q_point);
3960 *   residual_ad[i] -= (Nu_f * traction[component_i]) * JxW_f;
3961 *   }
3962 *   if ((i_group == p_fluid_block) && (std::abs(flow) > 1e-12))
3963 *   {
3964 *   const double Nu_p
3965 *   = scratch.fe_face_values_ref.shape_value(i, f_q_point);
3966 *   residual_ad[i] -= (Nu_p * flow) * JxW_f;
3967 *   }
3968 *   }
3969 *   }
3970 *   }
3971 *   }
3972 *  
3973 * @endcode
3974 *
3975 * Linearise the residual
3976 *
3977 * @code
3978 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
3979 *   {
3980 *   const ADNumberType &R_i = residual_ad[i];
3981 *  
3982 *   data.cell_rhs(i) -= R_i.val();
3983 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
3984 *   data.cell_matrix(i,j) += R_i.fastAccessDx(j);
3985 *   }
3986 *   }
3987 *  
3988 * @endcode
3989 *
3990 * Store the converged values of the internal variables
3991 *
3992 * @code
3993 *   template <int dim>
3994 *   void Solid<dim>::update_end_timestep()
3995 *   {
3996 *   FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
3997 *   cell (IteratorFilters::LocallyOwnedCell(),
3998 *   dof_handler_ref.begin_active()),
3999 *   endc (IteratorFilters::LocallyOwnedCell(),
4000 *   dof_handler_ref.end());
4001 *   for (; cell!=endc; ++cell)
4002 *   {
4003 *   Assert(cell->is_locally_owned(), ExcInternalError());
4004 *   Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
4005 *  
4006 *   const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >
4007 *   lqph = quadrature_point_history.get_data(cell);
4008 *   Assert(lqph.size() == n_q_points, ExcInternalError());
4009 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
4010 *   lqph[q_point]->update_end_timestep();
4011 *   }
4012 *   }
4013 *  
4014 *  
4015 * @endcode
4016 *
4017 * Solve the linearized equations
4018 *
4019 * @code
4020 *   template <int dim>
4021 *   void Solid<dim>::solve_linear_system( TrilinosWrappers::MPI::BlockVector &newton_update_OUT)
4022 *   {
4023 *  
4024 *   timerconsole.enter_subsection("Linear solver");
4025 *   timerfile.enter_subsection("Linear solver");
4026 *   pcout << " SLV " << std::flush;
4027 *   outfile << " SLV " << std::flush;
4028 *  
4029 *   TrilinosWrappers::MPI::Vector newton_update_nb;
4030 *   newton_update_nb.reinit(locally_owned_dofs, mpi_communicator);
4031 *  
4032 *   SolverControl solver_control (tangent_matrix_nb.m(),
4033 *   1.0e-6 * system_rhs_nb.l2_norm());
4034 *   TrilinosWrappers::SolverDirect solver (solver_control);
4035 *   solver.solve(tangent_matrix_nb, newton_update_nb, system_rhs_nb);
4036 *  
4037 * @endcode
4038 *
4039 * Copy the non-block solution back to block system
4040 *
4041 * @code
4042 *   for (unsigned int i=0; i<locally_owned_dofs.n_elements(); ++i)
4043 *   {
4044 *   const types::global_dof_index idx_i
4045 *   = locally_owned_dofs.nth_index_in_set(i);
4046 *   newton_update_OUT(idx_i) = newton_update_nb(idx_i);
4047 *   }
4048 *   newton_update_OUT.compress(VectorOperation::insert);
4049 *  
4050 *   timerconsole.leave_subsection();
4051 *   timerfile.leave_subsection();
4052 *   }
4053 *  
4054 * @endcode
4055 *
4056 * Class to compute gradient of the pressure
4057 *
4058 * @code
4059 *   template <int dim>
4060 *   class GradientPostprocessor : public DataPostprocessorVector<dim>
4061 *   {
4062 *   public:
4063 *   GradientPostprocessor (const unsigned int p_fluid_component)
4064 *   :
4065 *   DataPostprocessorVector<dim> ("grad_p",
4066 *   update_gradients),
4067 *   p_fluid_component (p_fluid_component)
4068 *   {}
4069 *  
4070 *   virtual ~GradientPostprocessor(){}
4071 *  
4072 *   virtual void
4073 *   evaluate_vector_field
4074 *   (const DataPostprocessorInputs::Vector<dim> &input_data,
4075 *   std::vector<Vector<double> > &computed_quantities) const override
4076 *   {
4077 *   AssertDimension (input_data.solution_gradients.size(),
4078 *   computed_quantities.size());
4079 *   for (unsigned int p=0; p<input_data.solution_gradients.size(); ++p)
4080 *   {
4081 *   AssertDimension (computed_quantities[p].size(), dim);
4082 *   for (unsigned int d=0; d<dim; ++d)
4083 *   computed_quantities[p][d]
4084 *   = input_data.solution_gradients[p][p_fluid_component][d];
4085 *   }
4086 *   }
4087 *  
4088 *   private:
4089 *   const unsigned int p_fluid_component;
4090 *   };
4091 *  
4092 *  
4093 * @endcode
4094 *
4095 * Print results to vtu file
4096 *
4097 * @code
4098 *   template <int dim> void Solid<dim>::output_results_to_vtu
4099 *   (const unsigned int timestep,
4100 *   const double current_time,
4101 *   TrilinosWrappers::MPI::BlockVector solution_IN) const
4102 *   {
4103 *   TrilinosWrappers::MPI::BlockVector solution_total(locally_owned_partitioning,
4104 *   locally_relevant_partitioning,
4105 *   mpi_communicator,
4106 *   false);
4107 *   solution_total = solution_IN;
4108 *   Vector<double> material_id;
4109 *   material_id.reinit(triangulation.n_active_cells());
4110 *   std::vector<types::subdomain_id> partition_int(triangulation.n_active_cells());
4111 *   GradientPostprocessor<dim> gradient_postprocessor(p_fluid_component);
4112 *  
4113 * @endcode
4114 *
4115 * Declare local variables with number of stress components
4116 * & assign value according to "dim" value
4117 *
4118 * @code
4119 *   unsigned int num_comp_symm_tensor = 6;
4120 *  
4121 * @endcode
4122 *
4123 * Declare local vectors to store values
4124 * OUTPUT AVERAGED ON ELEMENTS -------------------------------------------
4125 *
4126 * @code
4127 *   std::vector<Vector<double>>cauchy_stresses_total_elements
4128 *   (num_comp_symm_tensor,
4129 *   Vector<double> (triangulation.n_active_cells()));
4130 *   std::vector<Vector<double>>cauchy_stresses_E_elements
4131 *   (num_comp_symm_tensor,
4132 *   Vector<double> (triangulation.n_active_cells()));
4133 *   std::vector<Vector<double>>stretches_elements
4134 *   (dim,
4135 *   Vector<double> (triangulation.n_active_cells()));
4136 *   std::vector<Vector<double>>seepage_velocity_elements
4137 *   (dim,
4138 *   Vector<double> (triangulation.n_active_cells()));
4139 *   Vector<double> porous_dissipation_elements
4140 *   (triangulation.n_active_cells());
4141 *   Vector<double> viscous_dissipation_elements
4142 *   (triangulation.n_active_cells());
4143 *   Vector<double> solid_vol_fraction_elements
4144 *   (triangulation.n_active_cells());
4145 *  
4146 * @endcode
4147 *
4148 * OUTPUT AVERAGED ON NODES ----------------------------------------------
4149 * We need to create a new FE space with a single dof per node to avoid
4150 * duplication of the output on nodes for our problem with dim+1 dofs.
4151 *
4152 * @code
4153 *   FE_Q<dim> fe_vertex(1);
4154 *   DoFHandler<dim> vertex_handler_ref(triangulation);
4155 *   vertex_handler_ref.distribute_dofs(fe_vertex);
4156 *   AssertThrow(vertex_handler_ref.n_dofs() == triangulation.n_vertices(),
4157 *   ExcDimensionMismatch(vertex_handler_ref.n_dofs(),
4158 *   triangulation.n_vertices()));
4159 *  
4160 *   Vector<double> counter_on_vertices_mpi
4161 *   (vertex_handler_ref.n_dofs());
4162 *   Vector<double> sum_counter_on_vertices
4163 *   (vertex_handler_ref.n_dofs());
4164 *  
4165 *   std::vector<Vector<double>>cauchy_stresses_total_vertex_mpi
4166 *   (num_comp_symm_tensor,
4167 *   Vector<double>(vertex_handler_ref.n_dofs()));
4168 *   std::vector<Vector<double>>sum_cauchy_stresses_total_vertex
4169 *   (num_comp_symm_tensor,
4170 *   Vector<double>(vertex_handler_ref.n_dofs()));
4171 *   std::vector<Vector<double>>cauchy_stresses_E_vertex_mpi
4172 *   (num_comp_symm_tensor,
4173 *   Vector<double>(vertex_handler_ref.n_dofs()));
4174 *   std::vector<Vector<double>>sum_cauchy_stresses_E_vertex
4175 *   (num_comp_symm_tensor,
4176 *   Vector<double>(vertex_handler_ref.n_dofs()));
4177 *   std::vector<Vector<double>>stretches_vertex_mpi
4178 *   (dim,
4179 *   Vector<double>(vertex_handler_ref.n_dofs()));
4180 *   std::vector<Vector<double>>sum_stretches_vertex
4181 *   (dim,
4182 *   Vector<double>(vertex_handler_ref.n_dofs()));
4183 *   Vector<double> porous_dissipation_vertex_mpi(vertex_handler_ref.n_dofs());
4184 *   Vector<double> sum_porous_dissipation_vertex(vertex_handler_ref.n_dofs());
4185 *   Vector<double> viscous_dissipation_vertex_mpi(vertex_handler_ref.n_dofs());
4186 *   Vector<double> sum_viscous_dissipation_vertex(vertex_handler_ref.n_dofs());
4187 *   Vector<double> solid_vol_fraction_vertex_mpi(vertex_handler_ref.n_dofs());
4188 *   Vector<double> sum_solid_vol_fraction_vertex(vertex_handler_ref.n_dofs());
4189 *  
4190 * @endcode
4191 *
4192 * We need to create a new FE space with a dim dof per node to
4193 * be able to ouput data on nodes in vector form
4194 *
4195 * @code
4196 *   FESystem<dim> fe_vertex_vec(FE_Q<dim>(1),dim);
4197 *   DoFHandler<dim> vertex_vec_handler_ref(triangulation);
4198 *   vertex_vec_handler_ref.distribute_dofs(fe_vertex_vec);
4199 *   AssertThrow(vertex_vec_handler_ref.n_dofs() == (dim*triangulation.n_vertices()),
4200 *   ExcDimensionMismatch(vertex_vec_handler_ref.n_dofs(),
4201 *   (dim*triangulation.n_vertices())));
4202 *  
4203 *   Vector<double> seepage_velocity_vertex_vec_mpi(vertex_vec_handler_ref.n_dofs());
4204 *   Vector<double> sum_seepage_velocity_vertex_vec(vertex_vec_handler_ref.n_dofs());
4205 *   Vector<double> counter_on_vertices_vec_mpi(vertex_vec_handler_ref.n_dofs());
4206 *   Vector<double> sum_counter_on_vertices_vec(vertex_vec_handler_ref.n_dofs());
4207 * @endcode
4208 *
4209 * -----------------------------------------------------------------------
4210 *
4211
4212 *
4213 * Declare and initialize local unit vectors (to construct tensor basis)
4214 *
4215 * @code
4216 *   std::vector<Tensor<1,dim>> basis_vectors (dim, Tensor<1,dim>() );
4217 *   for (unsigned int i=0; i<dim; ++i)
4218 *   basis_vectors[i][i] = 1;
4219 *  
4220 * @endcode
4221 *
4222 * Declare an instance of the material class object
4223 *
4224 * @code
4225 *   if (parameters.mat_type == "Neo-Hooke")
4226 *   NeoHooke<dim,ADNumberType> material(parameters,time);
4227 *   else if (parameters.mat_type == "Ogden")
4228 *   Ogden<dim,ADNumberType> material(parameters,time);
4229 *   else if (parameters.mat_type == "visco-Ogden")
4230 *   visco_Ogden <dim,ADNumberType>material(parameters,time);
4231 *   else
4232 *   Assert (false, ExcMessage("Material type not implemented"));
4233 *  
4234 * @endcode
4235 *
4236 * Define a local instance of FEValues to compute updated values required
4237 * to calculate stresses
4238 *
4239 * @code
4240 *   const UpdateFlags uf_cell(update_values | update_gradients |
4241 *   update_JxW_values);
4242 *   FEValues<dim> fe_values_ref (fe, qf_cell, uf_cell);
4243 *  
4244 * @endcode
4245 *
4246 * Iterate through elements (cells) and Gauss Points
4247 *
4248 * @code
4249 *   FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
4250 *   cell(IteratorFilters::LocallyOwnedCell(),
4251 *   dof_handler_ref.begin_active()),
4252 *   endc(IteratorFilters::LocallyOwnedCell(),
4253 *   dof_handler_ref.end()),
4254 *   cell_v(IteratorFilters::LocallyOwnedCell(),
4255 *   vertex_handler_ref.begin_active()),
4256 *   cell_v_vec(IteratorFilters::LocallyOwnedCell(),
4257 *   vertex_vec_handler_ref.begin_active());
4258 * @endcode
4259 *
4260 * start cell loop
4261 *
4262 * @code
4263 *   for (; cell!=endc; ++cell, ++cell_v, ++cell_v_vec)
4264 *   {
4265 *   Assert(cell->is_locally_owned(), ExcInternalError());
4266 *   Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
4267 *  
4268 *   material_id(cell->active_cell_index())=
4269 *   static_cast<int>(cell->material_id());
4270 *  
4271 *   fe_values_ref.reinit(cell);
4272 *  
4273 *   std::vector<Tensor<2,dim>> solution_grads_u(n_q_points);
4274 *   fe_values_ref[u_fe].get_function_gradients(solution_total,
4275 *   solution_grads_u);
4276 *  
4277 *   std::vector<double> solution_values_p_fluid_total(n_q_points);
4278 *   fe_values_ref[p_fluid_fe].get_function_values(solution_total,
4279 *   solution_values_p_fluid_total);
4280 *  
4281 *   std::vector<Tensor<1,dim>> solution_grads_p_fluid_AD (n_q_points);
4282 *   fe_values_ref[p_fluid_fe].get_function_gradients(solution_total,
4283 *   solution_grads_p_fluid_AD);
4284 *  
4285 * @endcode
4286 *
4287 * start gauss point loop
4288 *
4289 * @code
4290 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4291 *   {
4292 *   const Tensor<2,dim,ADNumberType>
4293 *   F_AD = Physics::Elasticity::Kinematics::F(solution_grads_u[q_point]);
4294 *   ADNumberType det_F_AD = determinant(F_AD);
4295 *   const double det_F = Tensor<0,dim,double>(det_F_AD);
4296 *  
4297 *   const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4298 *   lqph = quadrature_point_history.get_data(cell);
4299 *   Assert(lqph.size() == n_q_points, ExcInternalError());
4300 *  
4301 *   const double p_fluid = solution_values_p_fluid_total[q_point];
4302 *  
4303 * @endcode
4304 *
4305 * Cauchy stress
4306 *
4307 * @code
4308 *   static const SymmetricTensor<2,dim,double>
4309 *   I (Physics::Elasticity::StandardTensors<dim>::I);
4310 *   SymmetricTensor<2,dim> sigma_E;
4311 *   const SymmetricTensor<2,dim,ADNumberType> sigma_E_AD =
4312 *   lqph[q_point]->get_Cauchy_E(F_AD);
4313 *  
4314 *   for (unsigned int i=0; i<dim; ++i)
4315 *   for (unsigned int j=0; j<dim; ++j)
4316 *   sigma_E[i][j] = Tensor<0,dim,double>(sigma_E_AD[i][j]);
4317 *  
4318 *   SymmetricTensor<2,dim> sigma_fluid_vol (I);
4319 *   sigma_fluid_vol *= -p_fluid;
4320 *   const SymmetricTensor<2,dim> sigma = sigma_E + sigma_fluid_vol;
4321 *  
4322 * @endcode
4323 *
4324 * Volumes
4325 *
4326 * @code
4327 *   const double solid_vol_fraction = (parameters.solid_vol_frac)/det_F;
4328 *  
4329 * @endcode
4330 *
4331 * Green-Lagrange strain
4332 *
4333 * @code
4334 *   const Tensor<2,dim> E_strain = 0.5*(transpose(F_AD)*F_AD - I);
4335 *  
4336 * @endcode
4337 *
4338 * Seepage velocity
4339 *
4340 * @code
4341 *   const Tensor<2,dim,ADNumberType> F_inv = invert(F_AD);
4342 *   const Tensor<1,dim,ADNumberType> grad_p_fluid_AD =
4343 *   solution_grads_p_fluid_AD[q_point]*F_inv;
4344 *   const Tensor<1,dim,ADNumberType> seepage_vel_AD =
4345 *   lqph[q_point]->get_seepage_velocity_current(F_AD, grad_p_fluid_AD);
4346 *  
4347 * @endcode
4348 *
4349 * Dissipations
4350 *
4351 * @code
4352 *   const double porous_dissipation =
4353 *   lqph[q_point]->get_porous_dissipation(F_AD, grad_p_fluid_AD);
4354 *   const double viscous_dissipation =
4355 *   lqph[q_point]->get_viscous_dissipation();
4356 *  
4357 * @endcode
4358 *
4359 * OUTPUT AVERAGED ON ELEMENTS -------------------------------------------
4360 * Both average on elements and on nodes is NOT weighted with the
4361 * integration point volume, i.e., we assume equal contribution of each
4362 * integration point to the average. Ideally, it should be weighted,
4363 * but I haven't invested time in getting it to work properly.
4364 *
4365 * @code
4366 *   if (parameters.outtype == "elements")
4367 *   {
4368 *   for (unsigned int j=0; j<dim; ++j)
4369 *   {
4370 *   cauchy_stresses_total_elements[j](cell->active_cell_index())
4371 *   += ((sigma*basis_vectors[j])*basis_vectors[j])/n_q_points;
4372 *   cauchy_stresses_E_elements[j](cell->active_cell_index())
4373 *   += ((sigma_E*basis_vectors[j])*basis_vectors[j])/n_q_points;
4374 *   stretches_elements[j](cell->active_cell_index())
4375 *   += std::sqrt(1.0+2.0*Tensor<0,dim,double>(E_strain[j][j]))
4376 *   /n_q_points;
4377 *   seepage_velocity_elements[j](cell->active_cell_index())
4378 *   += Tensor<0,dim,double>(seepage_vel_AD[j])/n_q_points;
4379 *   }
4380 *  
4381 *   porous_dissipation_elements(cell->active_cell_index())
4382 *   += porous_dissipation/n_q_points;
4383 *   viscous_dissipation_elements(cell->active_cell_index())
4384 *   += viscous_dissipation/n_q_points;
4385 *   solid_vol_fraction_elements(cell->active_cell_index())
4386 *   += solid_vol_fraction/n_q_points;
4387 *  
4388 *   cauchy_stresses_total_elements[3](cell->active_cell_index())
4389 *   += ((sigma*basis_vectors[0])*basis_vectors[1])/n_q_points; //sig_xy
4390 *   cauchy_stresses_total_elements[4](cell->active_cell_index())
4391 *   += ((sigma*basis_vectors[0])*basis_vectors[2])/n_q_points;//sig_xz
4392 *   cauchy_stresses_total_elements[5](cell->active_cell_index())
4393 *   += ((sigma*basis_vectors[1])*basis_vectors[2])/n_q_points;//sig_yz
4394 *  
4395 *   cauchy_stresses_E_elements[3](cell->active_cell_index())
4396 *   += ((sigma_E*basis_vectors[0])* basis_vectors[1])/n_q_points; //sig_xy
4397 *   cauchy_stresses_E_elements[4](cell->active_cell_index())
4398 *   += ((sigma_E*basis_vectors[0])* basis_vectors[2])/n_q_points;//sig_xz
4399 *   cauchy_stresses_E_elements[5](cell->active_cell_index())
4400 *   += ((sigma_E*basis_vectors[1])* basis_vectors[2])/n_q_points;//sig_yz
4401 *  
4402 *   }
4403 * @endcode
4404 *
4405 * OUTPUT AVERAGED ON NODES -------------------------------------------
4406 *
4407 * @code
4408 *   else if (parameters.outtype == "nodes")
4409 *   {
4410 *   for (unsigned int v=0; v<(GeometryInfo<dim>::vertices_per_cell); ++v)
4411 *   {
4412 *   types::global_dof_index local_vertex_indices =
4413 *   cell_v->vertex_dof_index(v, 0);
4414 *   counter_on_vertices_mpi(local_vertex_indices) += 1;
4415 *   for (unsigned int k=0; k<dim; ++k)
4416 *   {
4417 *   cauchy_stresses_total_vertex_mpi[k](local_vertex_indices)
4418 *   += (sigma*basis_vectors[k])*basis_vectors[k];
4419 *   cauchy_stresses_E_vertex_mpi[k](local_vertex_indices)
4420 *   += (sigma_E*basis_vectors[k])*basis_vectors[k];
4421 *   stretches_vertex_mpi[k](local_vertex_indices)
4422 *   += std::sqrt(1.0+2.0*Tensor<0,dim,double>(E_strain[k][k]));
4423 *  
4424 *   types::global_dof_index local_vertex_vec_indices =
4425 *   cell_v_vec->vertex_dof_index(v, k);
4426 *   counter_on_vertices_vec_mpi(local_vertex_vec_indices) += 1;
4427 *   seepage_velocity_vertex_vec_mpi(local_vertex_vec_indices)
4428 *   += Tensor<0,dim,double>(seepage_vel_AD[k]);
4429 *   }
4430 *  
4431 *   porous_dissipation_vertex_mpi(local_vertex_indices)
4432 *   += porous_dissipation;
4433 *   viscous_dissipation_vertex_mpi(local_vertex_indices)
4434 *   += viscous_dissipation;
4435 *   solid_vol_fraction_vertex_mpi(local_vertex_indices)
4436 *   += solid_vol_fraction;
4437 *  
4438 *   cauchy_stresses_total_vertex_mpi[3](local_vertex_indices)
4439 *   += (sigma*basis_vectors[0])*basis_vectors[1]; //sig_xy
4440 *   cauchy_stresses_total_vertex_mpi[4](local_vertex_indices)
4441 *   += (sigma*basis_vectors[0])*basis_vectors[2];//sig_xz
4442 *   cauchy_stresses_total_vertex_mpi[5](local_vertex_indices)
4443 *   += (sigma*basis_vectors[1])*basis_vectors[2]; //sig_yz
4444 *  
4445 *   cauchy_stresses_E_vertex_mpi[3](local_vertex_indices)
4446 *   += (sigma_E*basis_vectors[0])*basis_vectors[1]; //sig_xy
4447 *   cauchy_stresses_E_vertex_mpi[4](local_vertex_indices)
4448 *   += (sigma_E*basis_vectors[0])*basis_vectors[2];//sig_xz
4449 *   cauchy_stresses_E_vertex_mpi[5](local_vertex_indices)
4450 *   += (sigma_E*basis_vectors[1])*basis_vectors[2]; //sig_yz
4451 *   }
4452 *   }
4453 * @endcode
4454 *
4455 * ---------------------------------------------------------------
4456 *
4457 * @code
4458 *   } //end gauss point loop
4459 *   }//end cell loop
4460 *  
4461 * @endcode
4462 *
4463 * Different nodes might have different amount of contributions, e.g.,
4464 * corner nodes have less integration points contributing to the averaged.
4465 * This is why we need a counter and divide at the end, outside the cell loop.
4466 *
4467 * @code
4468 *   if (parameters.outtype == "nodes")
4469 *   {
4470 *   for (unsigned int d=0; d<(vertex_handler_ref.n_dofs()); ++d)
4471 *   {
4472 *   sum_counter_on_vertices[d] =
4473 *   Utilities::MPI::sum(counter_on_vertices_mpi[d],
4474 *   mpi_communicator);
4475 *   sum_porous_dissipation_vertex[d] =
4476 *   Utilities::MPI::sum(porous_dissipation_vertex_mpi[d],
4477 *   mpi_communicator);
4478 *   sum_viscous_dissipation_vertex[d] =
4479 *   Utilities::MPI::sum(viscous_dissipation_vertex_mpi[d],
4480 *   mpi_communicator);
4481 *   sum_solid_vol_fraction_vertex[d] =
4482 *   Utilities::MPI::sum(solid_vol_fraction_vertex_mpi[d],
4483 *   mpi_communicator);
4484 *  
4485 *   for (unsigned int k=0; k<num_comp_symm_tensor; ++k)
4486 *   {
4487 *   sum_cauchy_stresses_total_vertex[k][d] =
4488 *   Utilities::MPI::sum(cauchy_stresses_total_vertex_mpi[k][d],
4489 *   mpi_communicator);
4490 *   sum_cauchy_stresses_E_vertex[k][d] =
4491 *   Utilities::MPI::sum(cauchy_stresses_E_vertex_mpi[k][d],
4492 *   mpi_communicator);
4493 *   }
4494 *   for (unsigned int k=0; k<dim; ++k)
4495 *   {
4496 *   sum_stretches_vertex[k][d] =
4497 *   Utilities::MPI::sum(stretches_vertex_mpi[k][d],
4498 *   mpi_communicator);
4499 *   }
4500 *   }
4501 *  
4502 *   for (unsigned int d=0; d<(vertex_vec_handler_ref.n_dofs()); ++d)
4503 *   {
4504 *   sum_counter_on_vertices_vec[d] =
4505 *   Utilities::MPI::sum(counter_on_vertices_vec_mpi[d],
4506 *   mpi_communicator);
4507 *   sum_seepage_velocity_vertex_vec[d] =
4508 *   Utilities::MPI::sum(seepage_velocity_vertex_vec_mpi[d],
4509 *   mpi_communicator);
4510 *   }
4511 *  
4512 *   for (unsigned int d=0; d<(vertex_handler_ref.n_dofs()); ++d)
4513 *   {
4514 *   if (sum_counter_on_vertices[d]>0)
4515 *   {
4516 *   for (unsigned int i=0; i<num_comp_symm_tensor; ++i)
4517 *   {
4518 *   sum_cauchy_stresses_total_vertex[i][d] /= sum_counter_on_vertices[d];
4519 *   sum_cauchy_stresses_E_vertex[i][d] /= sum_counter_on_vertices[d];
4520 *   }
4521 *   for (unsigned int i=0; i<dim; ++i)
4522 *   {
4523 *   sum_stretches_vertex[i][d] /= sum_counter_on_vertices[d];
4524 *   }
4525 *   sum_porous_dissipation_vertex[d] /= sum_counter_on_vertices[d];
4526 *   sum_viscous_dissipation_vertex[d] /= sum_counter_on_vertices[d];
4527 *   sum_solid_vol_fraction_vertex[d] /= sum_counter_on_vertices[d];
4528 *   }
4529 *   }
4530 *  
4531 *   for (unsigned int d=0; d<(vertex_vec_handler_ref.n_dofs()); ++d)
4532 *   {
4533 *   if (sum_counter_on_vertices_vec[d]>0)
4534 *   {
4535 *   sum_seepage_velocity_vertex_vec[d] /= sum_counter_on_vertices_vec[d];
4536 *   }
4537 *   }
4538 *  
4539 *   }
4540 *  
4541 * @endcode
4542 *
4543 * Add the results to the solution to create the output file for Paraview
4544 *
4545 * @code
4546 *   DataOut<dim> data_out;
4547 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
4548 *   comp_type(dim,
4549 *   DataComponentInterpretation::component_is_part_of_vector);
4550 *   comp_type.push_back(DataComponentInterpretation::component_is_scalar);
4551 *  
4552 *   GridTools::get_subdomain_association(triangulation, partition_int);
4553 *  
4554 *   std::vector<std::string> solution_name(dim, "displacement");
4555 *   solution_name.push_back("pore_pressure");
4556 *  
4557 *   data_out.attach_dof_handler(dof_handler_ref);
4558 *   data_out.add_data_vector(solution_total,
4559 *   solution_name,
4560 *   DataOut<dim>::type_dof_data,
4561 *   comp_type);
4562 *  
4563 *   data_out.add_data_vector(solution_total,
4564 *   gradient_postprocessor);
4565 *  
4566 *   const Vector<double> partitioning(partition_int.begin(),
4567 *   partition_int.end());
4568 *  
4569 *   data_out.add_data_vector(partitioning, "partitioning");
4570 *   data_out.add_data_vector(material_id, "material_id");
4571 *  
4572 * @endcode
4573 *
4574 * Integration point results -----------------------------------------------------------
4575 *
4576 * @code
4577 *   if (parameters.outtype == "elements")
4578 *   {
4579 *   data_out.add_data_vector(cauchy_stresses_total_elements[0], "cauchy_xx");
4580 *   data_out.add_data_vector(cauchy_stresses_total_elements[1], "cauchy_yy");
4581 *   data_out.add_data_vector(cauchy_stresses_total_elements[2], "cauchy_zz");
4582 *   data_out.add_data_vector(cauchy_stresses_total_elements[3], "cauchy_xy");
4583 *   data_out.add_data_vector(cauchy_stresses_total_elements[4], "cauchy_xz");
4584 *   data_out.add_data_vector(cauchy_stresses_total_elements[5], "cauchy_yz");
4585 *  
4586 *   data_out.add_data_vector(cauchy_stresses_E_elements[0], "cauchy_E_xx");
4587 *   data_out.add_data_vector(cauchy_stresses_E_elements[1], "cauchy_E_yy");
4588 *   data_out.add_data_vector(cauchy_stresses_E_elements[2], "cauchy_E_zz");
4589 *   data_out.add_data_vector(cauchy_stresses_E_elements[3], "cauchy_E_xy");
4590 *   data_out.add_data_vector(cauchy_stresses_E_elements[4], "cauchy_E_xz");
4591 *   data_out.add_data_vector(cauchy_stresses_E_elements[5], "cauchy_E_yz");
4592 *  
4593 *   data_out.add_data_vector(stretches_elements[0], "stretch_xx");
4594 *   data_out.add_data_vector(stretches_elements[1], "stretch_yy");
4595 *   data_out.add_data_vector(stretches_elements[2], "stretch_zz");
4596 *  
4597 *   data_out.add_data_vector(seepage_velocity_elements[0], "seepage_vel_x");
4598 *   data_out.add_data_vector(seepage_velocity_elements[1], "seepage_vel_y");
4599 *   data_out.add_data_vector(seepage_velocity_elements[2], "seepage_vel_z");
4600 *  
4601 *   data_out.add_data_vector(porous_dissipation_elements, "dissipation_porous");
4602 *   data_out.add_data_vector(viscous_dissipation_elements, "dissipation_viscous");
4603 *   data_out.add_data_vector(solid_vol_fraction_elements, "solid_vol_fraction");
4604 *   }
4605 *   else if (parameters.outtype == "nodes")
4606 *   {
4607 *   data_out.add_data_vector(vertex_handler_ref,
4608 *   sum_cauchy_stresses_total_vertex[0],
4609 *   "cauchy_xx");
4610 *   data_out.add_data_vector(vertex_handler_ref,
4611 *   sum_cauchy_stresses_total_vertex[1],
4612 *   "cauchy_yy");
4613 *   data_out.add_data_vector(vertex_handler_ref,
4614 *   sum_cauchy_stresses_total_vertex[2],
4615 *   "cauchy_zz");
4616 *   data_out.add_data_vector(vertex_handler_ref,
4617 *   sum_cauchy_stresses_total_vertex[3],
4618 *   "cauchy_xy");
4619 *   data_out.add_data_vector(vertex_handler_ref,
4620 *   sum_cauchy_stresses_total_vertex[4],
4621 *   "cauchy_xz");
4622 *   data_out.add_data_vector(vertex_handler_ref,
4623 *   sum_cauchy_stresses_total_vertex[5],
4624 *   "cauchy_yz");
4625 *  
4626 *   data_out.add_data_vector(vertex_handler_ref,
4627 *   sum_cauchy_stresses_E_vertex[0],
4628 *   "cauchy_E_xx");
4629 *   data_out.add_data_vector(vertex_handler_ref,
4630 *   sum_cauchy_stresses_E_vertex[1],
4631 *   "cauchy_E_yy");
4632 *   data_out.add_data_vector(vertex_handler_ref,
4633 *   sum_cauchy_stresses_E_vertex[2],
4634 *   "cauchy_E_zz");
4635 *   data_out.add_data_vector(vertex_handler_ref,
4636 *   sum_cauchy_stresses_E_vertex[3],
4637 *   "cauchy_E_xy");
4638 *   data_out.add_data_vector(vertex_handler_ref,
4639 *   sum_cauchy_stresses_E_vertex[4],
4640 *   "cauchy_E_xz");
4641 *   data_out.add_data_vector(vertex_handler_ref,
4642 *   sum_cauchy_stresses_E_vertex[5],
4643 *   "cauchy_E_yz");
4644 *  
4645 *   data_out.add_data_vector(vertex_handler_ref,
4646 *   sum_stretches_vertex[0],
4647 *   "stretch_xx");
4648 *   data_out.add_data_vector(vertex_handler_ref,
4649 *   sum_stretches_vertex[1],
4650 *   "stretch_yy");
4651 *   data_out.add_data_vector(vertex_handler_ref,
4652 *   sum_stretches_vertex[2],
4653 *   "stretch_zz");
4654 *  
4655 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
4656 *   comp_type_vec(dim,
4657 *   DataComponentInterpretation::component_is_part_of_vector);
4658 *   std::vector<std::string> solution_name_vec(dim,"seepage_velocity");
4659 *  
4660 *   data_out.add_data_vector(vertex_vec_handler_ref,
4661 *   sum_seepage_velocity_vertex_vec,
4662 *   solution_name_vec,
4663 *   comp_type_vec);
4664 *  
4665 *   data_out.add_data_vector(vertex_handler_ref,
4666 *   sum_porous_dissipation_vertex,
4667 *   "dissipation_porous");
4668 *   data_out.add_data_vector(vertex_handler_ref,
4669 *   sum_viscous_dissipation_vertex,
4670 *   "dissipation_viscous");
4671 *   data_out.add_data_vector(vertex_handler_ref,
4672 *   sum_solid_vol_fraction_vertex,
4673 *   "solid_vol_fraction");
4674 *   }
4675 * @endcode
4676 *
4677 * ---------------------------------------------------------------------
4678 *
4679
4680 *
4681 *
4682 * @code
4683 *   data_out.build_patches(degree_displ);
4684 *  
4685 *   struct Filename
4686 *   {
4687 *   static std::string get_filename_vtu(unsigned int process,
4688 *   unsigned int timestep,
4689 *   const unsigned int n_digits = 5)
4690 *   {
4691 *   std::ostringstream filename_vtu;
4692 *   filename_vtu
4693 *   << "solution."
4694 *   << Utilities::int_to_string(process, n_digits)
4695 *   << "."
4696 *   << Utilities::int_to_string(timestep, n_digits)
4697 *   << ".vtu";
4698 *   return filename_vtu.str();
4699 *   }
4700 *  
4701 *   static std::string get_filename_pvtu(unsigned int timestep,
4702 *   const unsigned int n_digits = 5)
4703 *   {
4704 *   std::ostringstream filename_vtu;
4705 *   filename_vtu
4706 *   << "solution."
4707 *   << Utilities::int_to_string(timestep, n_digits)
4708 *   << ".pvtu";
4709 *   return filename_vtu.str();
4710 *   }
4711 *  
4712 *   static std::string get_filename_pvd (void)
4713 *   {
4714 *   std::ostringstream filename_vtu;
4715 *   filename_vtu
4716 *   << "solution.pvd";
4717 *   return filename_vtu.str();
4718 *   }
4719 *   };
4720 *  
4721 *   const std::string filename_vtu = Filename::get_filename_vtu(this_mpi_process,
4722 *   timestep);
4723 *   std::ofstream output(filename_vtu.c_str());
4724 *   data_out.write_vtu(output);
4725 *  
4726 * @endcode
4727 *
4728 * We have a collection of files written in parallel
4729 * This next set of steps should only be performed by master process
4730 *
4731 * @code
4732 *   if (this_mpi_process == 0)
4733 *   {
4734 * @endcode
4735 *
4736 * List of all files written out at this timestep by all processors
4737 *
4738 * @code
4739 *   std::vector<std::string> parallel_filenames_vtu;
4740 *   for (unsigned int p=0; p<n_mpi_processes; ++p)
4741 *   {
4742 *   parallel_filenames_vtu.push_back(Filename::get_filename_vtu(p, timestep));
4743 *   }
4744 *  
4745 *   const std::string filename_pvtu(Filename::get_filename_pvtu(timestep));
4746 *   std::ofstream pvtu_master(filename_pvtu.c_str());
4747 *   data_out.write_pvtu_record(pvtu_master,
4748 *   parallel_filenames_vtu);
4749 *  
4750 * @endcode
4751 *
4752 * Time dependent data master file
4753 *
4754 * @code
4755 *   static std::vector<std::pair<double,std::string>> time_and_name_history;
4756 *   time_and_name_history.push_back(std::make_pair(current_time,
4757 *   filename_pvtu));
4758 *   const std::string filename_pvd(Filename::get_filename_pvd());
4759 *   std::ofstream pvd_output(filename_pvd.c_str());
4760 *   DataOutBase::write_pvd_record(pvd_output, time_and_name_history);
4761 *   }
4762 *   }
4763 *  
4764 *  
4765 * @endcode
4766 *
4767 * Print results to plotting file
4768 *
4769 * @code
4770 *   template <int dim>
4771 *   void Solid<dim>::output_results_to_plot(
4772 *   const unsigned int timestep,
4773 *   const double current_time,
4774 *   TrilinosWrappers::MPI::BlockVector solution_IN,
4775 *   std::vector<Point<dim> > &tracked_vertices_IN,
4776 *   std::ofstream &plotpointfile) const
4777 *   {
4778 *   TrilinosWrappers::MPI::BlockVector solution_total(locally_owned_partitioning,
4779 *   locally_relevant_partitioning,
4780 *   mpi_communicator,
4781 *   false);
4782 *  
4783 *   (void) timestep;
4784 *   solution_total = solution_IN;
4785 *  
4786 * @endcode
4787 *
4788 * Variables needed to print the solution file for plotting
4789 *
4790 * @code
4791 *   Point<dim> reaction_force;
4792 *   Point<dim> reaction_force_pressure;
4793 *   Point<dim> reaction_force_extra;
4794 *   double total_fluid_flow = 0.0;
4795 *   double total_porous_dissipation = 0.0;
4796 *   double total_viscous_dissipation = 0.0;
4797 *   double total_solid_vol = 0.0;
4798 *   double total_vol_current = 0.0;
4799 *   double total_vol_reference = 0.0;
4800 *   std::vector<Point<dim+1>> solution_vertices(tracked_vertices_IN.size());
4801 *  
4802 * @endcode
4803 *
4804 * Auxiliar variables needed for mpi processing
4805 *
4806 * @code
4807 *   Tensor<1,dim> sum_reaction_mpi;
4808 *   Tensor<1,dim> sum_reaction_pressure_mpi;
4809 *   Tensor<1,dim> sum_reaction_extra_mpi;
4810 *   sum_reaction_mpi = 0.0;
4811 *   sum_reaction_pressure_mpi = 0.0;
4812 *   sum_reaction_extra_mpi = 0.0;
4813 *   double sum_total_flow_mpi = 0.0;
4814 *   double sum_porous_dissipation_mpi = 0.0;
4815 *   double sum_viscous_dissipation_mpi = 0.0;
4816 *   double sum_solid_vol_mpi = 0.0;
4817 *   double sum_vol_current_mpi = 0.0;
4818 *   double sum_vol_reference_mpi = 0.0;
4819 *  
4820 * @endcode
4821 *
4822 * Declare an instance of the material class object
4823 *
4824 * @code
4825 *   if (parameters.mat_type == "Neo-Hooke")
4826 *   NeoHooke<dim,ADNumberType> material(parameters,time);
4827 *   else if (parameters.mat_type == "Ogden")
4828 *   Ogden<dim,ADNumberType> material(parameters, time);
4829 *   else if (parameters.mat_type == "visco-Ogden")
4830 *   visco_Ogden <dim,ADNumberType>material(parameters,time);
4831 *   else
4832 *   Assert (false, ExcMessage("Material type not implemented"));
4833 *  
4834 * @endcode
4835 *
4836 * Define a local instance of FEValues to compute updated values required
4837 * to calculate stresses
4838 *
4839 * @code
4840 *   const UpdateFlags uf_cell(update_values | update_gradients |
4841 *   update_JxW_values);
4842 *   FEValues<dim> fe_values_ref (fe, qf_cell, uf_cell);
4843 *  
4844 * @endcode
4845 *
4846 * Iterate through elements (cells) and Gauss Points
4847 *
4848 * @code
4849 *   FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
4850 *   cell(IteratorFilters::LocallyOwnedCell(),
4851 *   dof_handler_ref.begin_active()),
4852 *   endc(IteratorFilters::LocallyOwnedCell(),
4853 *   dof_handler_ref.end());
4854 * @endcode
4855 *
4856 * start cell loop
4857 *
4858 * @code
4859 *   for (; cell!=endc; ++cell)
4860 *   {
4861 *   Assert(cell->is_locally_owned(), ExcInternalError());
4862 *   Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
4863 *  
4864 *   fe_values_ref.reinit(cell);
4865 *  
4866 *   std::vector<Tensor<2,dim>> solution_grads_u(n_q_points);
4867 *   fe_values_ref[u_fe].get_function_gradients(solution_total,
4868 *   solution_grads_u);
4869 *  
4870 *   std::vector<double> solution_values_p_fluid_total(n_q_points);
4871 *   fe_values_ref[p_fluid_fe].get_function_values(solution_total,
4872 *   solution_values_p_fluid_total);
4873 *  
4874 *   std::vector<Tensor<1,dim >> solution_grads_p_fluid_AD(n_q_points);
4875 *   fe_values_ref[p_fluid_fe].get_function_gradients(solution_total,
4876 *   solution_grads_p_fluid_AD);
4877 *  
4878 * @endcode
4879 *
4880 * start gauss point loop
4881 *
4882 * @code
4883 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4884 *   {
4885 *   const Tensor<2,dim,ADNumberType>
4886 *   F_AD = Physics::Elasticity::Kinematics::F(solution_grads_u[q_point]);
4887 *   ADNumberType det_F_AD = determinant(F_AD);
4888 *   const double det_F = Tensor<0,dim,double>(det_F_AD);
4889 *  
4890 *   const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4891 *   lqph = quadrature_point_history.get_data(cell);
4892 *   Assert(lqph.size() == n_q_points, ExcInternalError());
4893 *  
4894 *   double JxW = fe_values_ref.JxW(q_point);
4895 *  
4896 * @endcode
4897 *
4898 * Volumes
4899 *
4900 * @code
4901 *   sum_vol_current_mpi += det_F * JxW;
4902 *   sum_vol_reference_mpi += JxW;
4903 *   sum_solid_vol_mpi += parameters.solid_vol_frac * JxW * det_F;
4904 *  
4905 * @endcode
4906 *
4907 * Seepage velocity
4908 *
4909 * @code
4910 *   const Tensor<2,dim,ADNumberType> F_inv = invert(F_AD);
4911 *   const Tensor<1,dim,ADNumberType>
4912 *   grad_p_fluid_AD = solution_grads_p_fluid_AD[q_point]*F_inv;
4913 *   const Tensor<1,dim,ADNumberType> seepage_vel_AD
4914 *   = lqph[q_point]->get_seepage_velocity_current(F_AD, grad_p_fluid_AD);
4915 *  
4916 * @endcode
4917 *
4918 * Dissipations
4919 *
4920 * @code
4921 *   const double porous_dissipation =
4922 *   lqph[q_point]->get_porous_dissipation(F_AD, grad_p_fluid_AD);
4923 *   sum_porous_dissipation_mpi += porous_dissipation * det_F * JxW;
4924 *  
4925 *   const double viscous_dissipation = lqph[q_point]->get_viscous_dissipation();
4926 *   sum_viscous_dissipation_mpi += viscous_dissipation * det_F * JxW;
4927 *  
4928 * @endcode
4929 *
4930 * ---------------------------------------------------------------
4931 *
4932 * @code
4933 *   } //end gauss point loop
4934 *  
4935 * @endcode
4936 *
4937 * Compute reaction force on load boundary & total fluid flow across
4938 * drained boundary.
4939 * Define a local instance of FEFaceValues to compute values required
4940 * to calculate reaction force
4941 *
4942 * @code
4943 *   const UpdateFlags uf_face( update_values | update_gradients |
4944 *   update_normal_vectors | update_JxW_values );
4945 *   FEFaceValues<dim> fe_face_values_ref(fe, qf_face, uf_face);
4946 *  
4947 * @endcode
4948 *
4949 * start face loop
4950 *
4951 * @code
4952 *   for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
4953 *   {
4954 * @endcode
4955 *
4956 * Reaction force
4957 *
4958 * @code
4959 *   if (cell->face(face)->at_boundary() == true &&
4960 *   cell->face(face)->boundary_id() == get_reaction_boundary_id_for_output() )
4961 *   {
4962 *   fe_face_values_ref.reinit(cell, face);
4963 *  
4964 * @endcode
4965 *
4966 * Get displacement gradients for current face
4967 *
4968 * @code
4969 *   std::vector<Tensor<2,dim> > solution_grads_u_f(n_q_points_f);
4970 *   fe_face_values_ref[u_fe].get_function_gradients
4971 *   (solution_total,
4972 *   solution_grads_u_f);
4973 *  
4974 * @endcode
4975 *
4976 * Get pressure for current element
4977 *
4978 * @code
4979 *   std::vector< double > solution_values_p_fluid_total_f(n_q_points_f);
4980 *   fe_face_values_ref[p_fluid_fe].get_function_values
4981 *   (solution_total,
4982 *   solution_values_p_fluid_total_f);
4983 *  
4984 * @endcode
4985 *
4986 * start gauss points on faces loop
4987 *
4988 * @code
4989 *   for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point)
4990 *   {
4991 *   const Tensor<1,dim> &N = fe_face_values_ref.normal_vector(f_q_point);
4992 *   const double JxW_f = fe_face_values_ref.JxW(f_q_point);
4993 *  
4994 * @endcode
4995 *
4996 * Compute deformation gradient from displacements gradient
4997 * (present configuration)
4998 *
4999 * @code
5000 *   const Tensor<2,dim,ADNumberType> F_AD =
5001 *   Physics::Elasticity::Kinematics::F(solution_grads_u_f[f_q_point]);
5002 *  
5003 *   const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
5004 *   lqph = quadrature_point_history.get_data(cell);
5005 *   Assert(lqph.size() == n_q_points, ExcInternalError());
5006 *  
5007 *   const double p_fluid = solution_values_p_fluid_total[f_q_point];
5008 *  
5009 * @endcode
5010 *
5011 * Cauchy stress
5012 *
5013 * @code
5014 *   static const SymmetricTensor<2,dim,double>
5015 *   I (Physics::Elasticity::StandardTensors<dim>::I);
5016 *   SymmetricTensor<2,dim> sigma_E;
5017 *   const SymmetricTensor<2,dim,ADNumberType> sigma_E_AD =
5018 *   lqph[f_q_point]->get_Cauchy_E(F_AD);
5019 *  
5020 *   for (unsigned int i=0; i<dim; ++i)
5021 *   for (unsigned int j=0; j<dim; ++j)
5022 *   sigma_E[i][j] = Tensor<0,dim,double>(sigma_E_AD[i][j]);
5023 *  
5024 *   SymmetricTensor<2,dim> sigma_fluid_vol(I);
5025 *   sigma_fluid_vol *= -1.0*p_fluid;
5026 *   const SymmetricTensor<2,dim> sigma = sigma_E+sigma_fluid_vol;
5027 *   sum_reaction_mpi += sigma * N * JxW_f;
5028 *   sum_reaction_pressure_mpi += sigma_fluid_vol * N * JxW_f;
5029 *   sum_reaction_extra_mpi += sigma_E * N * JxW_f;
5030 *   }//end gauss points on faces loop
5031 *   }
5032 *  
5033 * @endcode
5034 *
5035 * Fluid flow
5036 *
5037 * @code
5038 *   if (cell->face(face)->at_boundary() == true &&
5039 *   (cell->face(face)->boundary_id() ==
5040 *   get_drained_boundary_id_for_output().first ||
5041 *   cell->face(face)->boundary_id() ==
5042 *   get_drained_boundary_id_for_output().second ) )
5043 *   {
5044 *   fe_face_values_ref.reinit(cell, face);
5045 *  
5046 * @endcode
5047 *
5048 * Get displacement gradients for current face
5049 *
5050 * @code
5051 *   std::vector<Tensor<2,dim>> solution_grads_u_f(n_q_points_f);
5052 *   fe_face_values_ref[u_fe].get_function_gradients
5053 *   (solution_total,
5054 *   solution_grads_u_f);
5055 *  
5056 * @endcode
5057 *
5058 * Get pressure gradients for current face
5059 *
5060 * @code
5061 *   std::vector<Tensor<1,dim>> solution_grads_p_f(n_q_points_f);
5062 *   fe_face_values_ref[p_fluid_fe].get_function_gradients
5063 *   (solution_total,
5064 *   solution_grads_p_f);
5065 *  
5066 * @endcode
5067 *
5068 * start gauss points on faces loop
5069 *
5070 * @code
5071 *   for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point)
5072 *   {
5073 *   const Tensor<1,dim> &N =
5074 *   fe_face_values_ref.normal_vector(f_q_point);
5075 *   const double JxW_f = fe_face_values_ref.JxW(f_q_point);
5076 *  
5077 * @endcode
5078 *
5079 * Deformation gradient and inverse from displacements gradient
5080 * (present configuration)
5081 *
5082 * @code
5083 *   const Tensor<2,dim,ADNumberType> F_AD
5084 *   = Physics::Elasticity::Kinematics::F(solution_grads_u_f[f_q_point]);
5085 *  
5086 *   const Tensor<2,dim,ADNumberType> F_inv_AD = invert(F_AD);
5087 *   ADNumberType det_F_AD = determinant(F_AD);
5088 *  
5089 *   const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
5090 *   lqph = quadrature_point_history.get_data(cell);
5091 *   Assert(lqph.size() == n_q_points, ExcInternalError());
5092 *  
5093 * @endcode
5094 *
5095 * Seepage velocity
5096 *
5097 * @code
5098 *   Tensor<1,dim> seepage;
5099 *   double det_F = Tensor<0,dim,double>(det_F_AD);
5100 *   const Tensor<1,dim,ADNumberType> grad_p
5101 *   = solution_grads_p_f[f_q_point]*F_inv_AD;
5102 *   const Tensor<1,dim,ADNumberType> seepage_AD
5103 *   = lqph[f_q_point]->get_seepage_velocity_current(F_AD, grad_p);
5104 *  
5105 *   for (unsigned int i=0; i<dim; ++i)
5106 *   seepage[i] = Tensor<0,dim,double>(seepage_AD[i]);
5107 *  
5108 *   sum_total_flow_mpi += (seepage/det_F) * N * JxW_f;
5109 *   }//end gauss points on faces loop
5110 *   }
5111 *   }//end face loop
5112 *   }//end cell loop
5113 *  
5114 * @endcode
5115 *
5116 * Sum the results from different MPI process and then add to the reaction_force vector
5117 * In theory, the solution on each surface (each cell) only exists in one MPI process
5118 * so, we add all MPI process, one will have the solution and the others will be zero
5119 *
5120 * @code
5121 *   for (unsigned int d=0; d<dim; ++d)
5122 *   {
5123 *   reaction_force[d] = Utilities::MPI::sum(sum_reaction_mpi[d],
5124 *   mpi_communicator);
5125 *   reaction_force_pressure[d] = Utilities::MPI::sum(sum_reaction_pressure_mpi[d],
5126 *   mpi_communicator);
5127 *   reaction_force_extra[d] = Utilities::MPI::sum(sum_reaction_extra_mpi[d],
5128 *   mpi_communicator);
5129 *   }
5130 *  
5131 * @endcode
5132 *
5133 * Same for total fluid flow, and for porous and viscous dissipations
5134 *
5135 * @code
5136 *   total_fluid_flow = Utilities::MPI::sum(sum_total_flow_mpi,
5137 *   mpi_communicator);
5138 *   total_porous_dissipation = Utilities::MPI::sum(sum_porous_dissipation_mpi,
5139 *   mpi_communicator);
5140 *   total_viscous_dissipation = Utilities::MPI::sum(sum_viscous_dissipation_mpi,
5141 *   mpi_communicator);
5142 *   total_solid_vol = Utilities::MPI::sum(sum_solid_vol_mpi,
5143 *   mpi_communicator);
5144 *   total_vol_current = Utilities::MPI::sum(sum_vol_current_mpi,
5145 *   mpi_communicator);
5146 *   total_vol_reference = Utilities::MPI::sum(sum_vol_reference_mpi,
5147 *   mpi_communicator);
5148 *  
5149 * @endcode
5150 *
5151 * Extract solution for tracked vectors
5152 * Copying an MPI::BlockVector into MPI::Vector is not possible,
5153 * so we copy each block of MPI::BlockVector into an MPI::Vector
5154 * And then we copy the MPI::Vector into "normal" Vectors
5155 *
5156 * @code
5157 *   TrilinosWrappers::MPI::Vector solution_vector_u_MPI(solution_total.block(u_block));
5158 *   TrilinosWrappers::MPI::Vector solution_vector_p_MPI(solution_total.block(p_fluid_block));
5159 *   Vector<double> solution_u_vector(solution_vector_u_MPI);
5160 *   Vector<double> solution_p_vector(solution_vector_p_MPI);
5161 *  
5162 *   if (this_mpi_process == 0)
5163 *   {
5164 * @endcode
5165 *
5166 * Append the pressure solution vector to the displacement solution vector,
5167 * creating a single solution vector equivalent to the original BlockVector
5168 * so FEFieldFunction will work with the dof_handler_ref.
5169 *
5170 * @code
5171 *   Vector<double> solution_vector(solution_p_vector.size()
5172 *   +solution_u_vector.size());
5173 *  
5174 *   for (unsigned int d=0; d<(solution_u_vector.size()); ++d)
5175 *   solution_vector[d] = solution_u_vector[d];
5176 *  
5177 *   for (unsigned int d=0; d<(solution_p_vector.size()); ++d)
5178 *   solution_vector[solution_u_vector.size()+d] = solution_p_vector[d];
5179 *  
5180 *   Functions::FEFieldFunction<dim,Vector<double>>
5181 *   find_solution(dof_handler_ref, solution_vector);
5182 *  
5183 *   for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
5184 *   {
5185 *   Vector<double> update(dim+1);
5186 *   Point<dim> pt_ref;
5187 *  
5188 *   pt_ref[0]= tracked_vertices_IN[p][0];
5189 *   pt_ref[1]= tracked_vertices_IN[p][1];
5190 *   pt_ref[2]= tracked_vertices_IN[p][2];
5191 *  
5192 *   find_solution.vector_value(pt_ref, update);
5193 *  
5194 *   for (unsigned int d=0; d<(dim+1); ++d)
5195 *   {
5196 * @endcode
5197 *
5198 * For values close to zero, set to 0.0
5199 *
5200 * @code
5201 *   if (abs(update[d])<1.5*parameters.tol_u)
5202 *   update[d] = 0.0;
5203 *   solution_vertices[p][d] = update[d];
5204 *   }
5205 *   }
5206 * @endcode
5207 *
5208 * Write the results to the plotting file.
5209 * Add two blank lines between cycles in the cyclic loading examples so GNUPLOT can detect each cycle as a different block
5210 *
5211 * @code
5212 *   if (( (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")||
5213 *   (parameters.geom_type == "Budday_cube_tension_compression")||
5214 *   (parameters.geom_type == "Budday_cube_shear_fully_fixed") ) &&
5215 *   ( (abs(current_time - parameters.end_time/3.) <0.9*parameters.delta_t)||
5216 *   (abs(current_time - 2.*parameters.end_time/3.)<0.9*parameters.delta_t) ) &&
5217 *   parameters.num_cycle_sets == 1 )
5218 *   {
5219 *   plotpointfile << std::endl<< std::endl;
5220 *   }
5221 *   if (( (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")||
5222 *   (parameters.geom_type == "Budday_cube_tension_compression")||
5223 *   (parameters.geom_type == "Budday_cube_shear_fully_fixed") ) &&
5224 *   ( (abs(current_time - parameters.end_time/9.) <0.9*parameters.delta_t)||
5225 *   (abs(current_time - 2.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5226 *   (abs(current_time - 3.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5227 *   (abs(current_time - 5.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5228 *   (abs(current_time - 7.*parameters.end_time/9.)<0.9*parameters.delta_t) ) &&
5229 *   parameters.num_cycle_sets == 2 )
5230 *   {
5231 *   plotpointfile << std::endl<< std::endl;
5232 *   }
5233 *  
5234 *   plotpointfile << std::setprecision(6) << std::scientific;
5235 *   plotpointfile << std::setw(16) << current_time << ","
5236 *   << std::setw(15) << total_vol_reference << ","
5237 *   << std::setw(15) << total_vol_current << ","
5238 *   << std::setw(15) << total_solid_vol << ",";
5239 *  
5240 *   if (current_time == 0.0)
5241 *   {
5242 *   for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
5243 *   {
5244 *   for (unsigned int d=0; d<dim; ++d)
5245 *   plotpointfile << std::setw(15) << 0.0 << ",";
5246 *  
5247 *   plotpointfile << std::setw(15) << parameters.drained_pressure << ",";
5248 *   }
5249 *   for (unsigned int d=0; d<(3*dim+2); ++d)
5250 *   plotpointfile << std::setw(15) << 0.0 << ",";
5251 *  
5252 *   plotpointfile << std::setw(15) << 0.0;
5253 *   }
5254 *   else
5255 *   {
5256 *   for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
5257 *   for (unsigned int d=0; d<(dim+1); ++d)
5258 *   plotpointfile << std::setw(15) << solution_vertices[p][d]<< ",";
5259 *  
5260 *   for (unsigned int d=0; d<dim; ++d)
5261 *   plotpointfile << std::setw(15) << reaction_force[d] << ",";
5262 *  
5263 *   for (unsigned int d=0; d<dim; ++d)
5264 *   plotpointfile << std::setw(15) << reaction_force_pressure[d] << ",";
5265 *  
5266 *   for (unsigned int d=0; d<dim; ++d)
5267 *   plotpointfile << std::setw(15) << reaction_force_extra[d] << ",";
5268 *  
5269 *   plotpointfile << std::setw(15) << total_fluid_flow << ","
5270 *   << std::setw(15) << total_porous_dissipation<< ","
5271 *   << std::setw(15) << total_viscous_dissipation;
5272 *   }
5273 *   plotpointfile << std::endl;
5274 *   }
5275 *   }
5276 *  
5277 * @endcode
5278 *
5279 * Header for console output file
5280 *
5281 * @code
5282 *   template <int dim>
5283 *   void Solid<dim>::print_console_file_header(std::ofstream &outputfile) const
5284 *   {
5285 *   outputfile << "/*-----------------------------------------------------------------------------------------";
5286 *   outputfile << "\n\n Poro-viscoelastic formulation to solve nonlinear solid mechanics problems using deal.ii";
5287 *   outputfile << "\n\n Problem setup by E Comellas and J-P Pelteret, University of Erlangen-Nuremberg, 2018";
5288 *   outputfile << "\n\n/*-----------------------------------------------------------------------------------------";
5289 *   outputfile << "\n\nCONSOLE OUTPUT: \n\n";
5290 *   }
5291 *  
5292 * @endcode
5293 *
5294 * Header for plotting output file
5295 *
5296 * @code
5297 *   template <int dim>
5298 *   void Solid<dim>::print_plot_file_header(std::vector<Point<dim> > &tracked_vertices,
5299 *   std::ofstream &plotpointfile) const
5300 *   {
5301 *   plotpointfile << "#\n# *** Solution history for tracked vertices -- DOF: 0 = Ux, 1 = Uy, 2 = Uz, 3 = P ***"
5302 *   << std::endl;
5303 *  
5304 *   for (unsigned int p=0; p<tracked_vertices.size(); ++p)
5305 *   {
5306 *   plotpointfile << "# Point " << p << " coordinates: ";
5307 *   for (unsigned int d=0; d<dim; ++d)
5308 *   {
5309 *   plotpointfile << tracked_vertices[p][d];
5310 *   if (!( (p == tracked_vertices.size()-1) && (d == dim-1) ))
5311 *   plotpointfile << ", ";
5312 *   }
5313 *   plotpointfile << std::endl;
5314 *   }
5315 *   plotpointfile << "# The reaction force is the integral over the loaded surfaces in the "
5316 *   << "undeformed configuration of the Cauchy stress times the normal surface unit vector.\n"
5317 *   << "# reac(p) corresponds to the volumetric part of the Cauchy stress due to the pore fluid pressure"
5318 *   << " and reac(E) corresponds to the extra part of the Cauchy stress due to the solid contribution."
5319 *   << std::endl
5320 *   << "# The fluid flow is the integral over the drained surfaces in the "
5321 *   << "undeformed configuration of the seepage velocity times the normal surface unit vector."
5322 *   << std::endl
5323 *   << "# Column number:"
5324 *   << std::endl
5325 *   << "#";
5326 *  
5327 *   unsigned int columns = 24;
5328 *   for (unsigned int d=1; d<columns; ++d)
5329 *   plotpointfile << std::setw(15)<< d <<",";
5330 *  
5331 *   plotpointfile << std::setw(15)<< columns
5332 *   << std::endl
5333 *   << "#"
5334 *   << std::right << std::setw(16) << "Time,"
5335 *   << std::right << std::setw(16) << "ref vol,"
5336 *   << std::right << std::setw(16) << "def vol,"
5337 *   << std::right << std::setw(16) << "solid vol,";
5338 *   for (unsigned int p=0; p<tracked_vertices.size(); ++p)
5339 *   for (unsigned int d=0; d<(dim+1); ++d)
5340 *   plotpointfile << std::right<< std::setw(11)
5341 *   <<"P" << p << "[" << d << "],";
5342 *  
5343 *   for (unsigned int d=0; d<dim; ++d)
5344 *   plotpointfile << std::right<< std::setw(13)
5345 *   << "reaction [" << d << "],";
5346 *  
5347 *   for (unsigned int d=0; d<dim; ++d)
5348 *   plotpointfile << std::right<< std::setw(13)
5349 *   << "reac(p) [" << d << "],";
5350 *  
5351 *   for (unsigned int d=0; d<dim; ++d)
5352 *   plotpointfile << std::right<< std::setw(13)
5353 *   << "reac(E) [" << d << "],";
5354 *  
5355 *   plotpointfile << std::right<< std::setw(16)<< "fluid flow,"
5356 *   << std::right<< std::setw(16)<< "porous dissip,"
5357 *   << std::right<< std::setw(15)<< "viscous dissip"
5358 *   << std::endl;
5359 *   }
5360 *  
5361 * @endcode
5362 *
5363 * Footer for console output file
5364 *
5365 * @code
5366 *   template <int dim>
5367 *   void Solid<dim>::print_console_file_footer(std::ofstream &outputfile) const
5368 *   {
5369 * @endcode
5370 *
5371 * Copy "parameters" file at end of output file.
5372 *
5373 * @code
5374 *   std::ifstream infile("parameters.prm");
5375 *   std::string content = "";
5376 *   int i;
5377 *  
5378 *   for(i=0 ; infile.eof()!=true ; i++)
5379 *   {
5380 *   char aux = infile.get();
5381 *   content += aux;
5382 *   if(aux=='\n') content += '#';
5383 *   }
5384 *  
5385 *   i--;
5386 *   content.erase(content.end()-1);
5387 *   infile.close();
5388 *  
5389 *   outputfile << "\n\n\n\n PARAMETERS FILE USED IN THIS COMPUTATION: \n#"
5390 *   << std::endl
5391 *   << content;
5392 *   }
5393 *  
5394 * @endcode
5395 *
5396 * Footer for plotting output file
5397 *
5398 * @code
5399 *   template <int dim>
5400 *   void Solid<dim>::print_plot_file_footer(std::ofstream &plotpointfile) const
5401 *   {
5402 * @endcode
5403 *
5404 * Copy "parameters" file at end of output file.
5405 *
5406 * @code
5407 *   std::ifstream infile("parameters.prm");
5408 *   std::string content = "";
5409 *   int i;
5410 *  
5411 *   for(i=0 ; infile.eof()!=true ; i++)
5412 *   {
5413 *   char aux = infile.get();
5414 *   content += aux;
5415 *   if(aux=='\n') content += '#';
5416 *   }
5417 *  
5418 *   i--;
5419 *   content.erase(content.end()-1);
5420 *   infile.close();
5421 *  
5422 *   plotpointfile << "#"<< std::endl
5423 *   << "#"<< std::endl
5424 *   << "# PARAMETERS FILE USED IN THIS COMPUTATION:" << std::endl
5425 *   << "#"<< std::endl
5426 *   << content;
5427 *   }
5428 *  
5429 *  
5430 * @endcode
5431 *
5432 *
5433 * <a name="VerificationexamplesfromEhlersandEipper1999"></a>
5434 * <h3>Verification examples from Ehlers and Eipper 1999</h3>
5435 * We group the definition of the geometry, boundary and loading conditions specific to
5436 * the verification examples from Ehlers and Eipper 1999 into specific classes.
5437 *
5438
5439 *
5440 *
5441 * <a name="BaseclassTubegeometryandboundaryconditions"></a>
5442 * <h4>Base class: Tube geometry and boundary conditions</h4>
5443 *
5444 * @code
5445 *   template <int dim>
5446 *   class VerificationEhlers1999TubeBase
5447 *   : public Solid<dim>
5448 *   {
5449 *   public:
5450 *   VerificationEhlers1999TubeBase (const Parameters::AllParameters &parameters)
5451 *   : Solid<dim> (parameters)
5452 *   {}
5453 *  
5454 *   virtual ~VerificationEhlers1999TubeBase () {}
5455 *  
5456 *   private:
5457 *   virtual void make_grid() override
5458 *   {
5459 *   GridGenerator::cylinder( this->triangulation,
5460 *   0.1,
5461 *   0.5);
5462 *  
5463 *   const double rot_angle = 3.0*numbers::PI/2.0;
5464 *   GridTools::rotate( Point<3>::unit_vector(1), rot_angle, this->triangulation);
5465 *  
5466 *   this->triangulation.reset_manifold(0);
5467 *   static const CylindricalManifold<dim> manifold_description_3d(2);
5468 *   this->triangulation.set_manifold (0, manifold_description_3d);
5469 *   GridTools::scale(this->parameters.scale, this->triangulation);
5470 *   this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5471 *   this->triangulation.reset_manifold(0);
5472 *   }
5473 *  
5474 *   virtual void define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
5475 *   {
5476 *   tracked_vertices[0][0] = 0.0*this->parameters.scale;
5477 *   tracked_vertices[0][1] = 0.0*this->parameters.scale;
5478 *   tracked_vertices[0][2] = 0.5*this->parameters.scale;
5479 *  
5480 *   tracked_vertices[1][0] = 0.0*this->parameters.scale;
5481 *   tracked_vertices[1][1] = 0.0*this->parameters.scale;
5482 *   tracked_vertices[1][2] = -0.5*this->parameters.scale;
5483 *   }
5484 *  
5485 *   virtual void make_dirichlet_constraints(AffineConstraints<double> &constraints) override
5486 *   {
5487 *   if (this->time.get_timestep() < 2)
5488 *   {
5489 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5490 *   2,
5491 *   Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5492 *   constraints,
5493 *   (this->fe.component_mask(this->pressure)));
5494 *   }
5495 *   else
5496 *   {
5497 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5498 *   2,
5499 *   Functions::ZeroFunction<dim>(this->n_components),
5500 *   constraints,
5501 *   (this->fe.component_mask(this->pressure)));
5502 *   }
5503 *  
5504 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5505 *   0,
5506 *   Functions::ZeroFunction<dim>(this->n_components),
5507 *   constraints,
5508 *   (this->fe.component_mask(this->x_displacement)|
5509 *   this->fe.component_mask(this->y_displacement) ) );
5510 *  
5511 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5512 *   1,
5513 *   Functions::ZeroFunction<dim>(this->n_components),
5514 *   constraints,
5515 *   (this->fe.component_mask(this->x_displacement) |
5516 *   this->fe.component_mask(this->y_displacement) |
5517 *   this->fe.component_mask(this->z_displacement) ));
5518 *   }
5519 *  
5520 *   virtual double
5521 *   get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5522 *   const Point<dim> &pt) const override
5523 *   {
5524 *   (void)pt;
5525 *   (void)boundary_id;
5526 *   return 0.0;
5527 *   }
5528 *  
5529 *   virtual types::boundary_id
5530 *   get_reaction_boundary_id_for_output() const override
5531 *   {
5532 *   return 2;
5533 *   }
5534 *  
5535 *   virtual std::pair<types::boundary_id,types::boundary_id>
5536 *   get_drained_boundary_id_for_output() const override
5537 *   {
5538 *   return std::make_pair(2,2);
5539 *   }
5540 *  
5541 *   virtual std::vector<double>
5542 *   get_dirichlet_load(const types::boundary_id &boundary_id,
5543 *   const int &direction) const override
5544 *   {
5545 *   std::vector<double> displ_incr(dim, 0.0);
5546 *   (void)boundary_id;
5547 *   (void)direction;
5548 *   AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));
5549 *  
5550 *   return displ_incr;
5551 *   }
5552 *   };
5553 *  
5554 * @endcode
5555 *
5556 *
5557 * <a name="DerivedclassSteploadexample"></a>
5558 * <h4>Derived class: Step load example</h4>
5559 *
5560 * @code
5561 *   template <int dim>
5562 *   class VerificationEhlers1999StepLoad
5563 *   : public VerificationEhlers1999TubeBase<dim>
5564 *   {
5565 *   public:
5566 *   VerificationEhlers1999StepLoad (const Parameters::AllParameters &parameters)
5567 *   : VerificationEhlers1999TubeBase<dim> (parameters)
5568 *   {}
5569 *  
5570 *   virtual ~VerificationEhlers1999StepLoad () {}
5571 *  
5572 *   private:
5573 *   virtual Tensor<1,dim>
5574 *   get_neumann_traction (const types::boundary_id &boundary_id,
5575 *   const Point<dim> &pt,
5576 *   const Tensor<1,dim> &N) const override
5577 *   {
5578 *   if (this->parameters.load_type == "pressure")
5579 *   {
5580 *   if (boundary_id == 2)
5581 *   {
5582 *   return this->parameters.load * N;
5583 *   }
5584 *   }
5585 *  
5586 *   (void)pt;
5587 *  
5588 *   return Tensor<1,dim>();
5589 *   }
5590 *   };
5591 *  
5592 * @endcode
5593 *
5594 *
5595 * <a name="DerivedclassLoadincreasingexample"></a>
5596 * <h4>Derived class: Load increasing example</h4>
5597 *
5598 * @code
5599 *   template <int dim>
5600 *   class VerificationEhlers1999IncreaseLoad
5601 *   : public VerificationEhlers1999TubeBase<dim>
5602 *   {
5603 *   public:
5604 *   VerificationEhlers1999IncreaseLoad (const Parameters::AllParameters &parameters)
5605 *   : VerificationEhlers1999TubeBase<dim> (parameters)
5606 *   {}
5607 *  
5608 *   virtual ~VerificationEhlers1999IncreaseLoad () {}
5609 *  
5610 *   private:
5611 *   virtual Tensor<1,dim>
5612 *   get_neumann_traction (const types::boundary_id &boundary_id,
5613 *   const Point<dim> &pt,
5614 *   const Tensor<1,dim> &N) const override
5615 *   {
5616 *   if (this->parameters.load_type == "pressure")
5617 *   {
5618 *   if (boundary_id == 2)
5619 *   {
5620 *   const double initial_load = this->parameters.load;
5621 *   const double final_load = 20.0*initial_load;
5622 *   const double initial_time = this->time.get_delta_t();
5623 *   const double final_time = this->time.get_end();
5624 *   const double current_time = this->time.get_current();
5625 *   const double load = initial_load + (final_load-initial_load)*(current_time-initial_time)/(final_time-initial_time);
5626 *   return load * N;
5627 *   }
5628 *   }
5629 *  
5630 *   (void)pt;
5631 *  
5632 *   return Tensor<1,dim>();
5633 *   }
5634 *   };
5635 *  
5636 * @endcode
5637 *
5638 *
5639 * <a name="ClassConsolidationcube"></a>
5640 * <h4>Class: Consolidation cube</h4>
5641 *
5642 * @code
5643 *   template <int dim>
5644 *   class VerificationEhlers1999CubeConsolidation
5645 *   : public Solid<dim>
5646 *   {
5647 *   public:
5648 *   VerificationEhlers1999CubeConsolidation (const Parameters::AllParameters &parameters)
5649 *   : Solid<dim> (parameters)
5650 *   {}
5651 *  
5652 *   virtual ~VerificationEhlers1999CubeConsolidation () {}
5653 *  
5654 *   private:
5655 *   virtual void
5656 *   make_grid() override
5657 *   {
5658 *   GridGenerator::hyper_rectangle(this->triangulation,
5659 *   Point<dim>(0.0, 0.0, 0.0),
5660 *   Point<dim>(1.0, 1.0, 1.0),
5661 *   true);
5662 *  
5663 *   GridTools::scale(this->parameters.scale, this->triangulation);
5664 *   this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5665 *  
5666 *   typename Triangulation<dim>::active_cell_iterator cell =
5667 *   this->triangulation.begin_active(), endc = this->triangulation.end();
5668 *   for (; cell != endc; ++cell)
5669 *   {
5670 *   for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5671 *   if (cell->face(face)->at_boundary() == true &&
5672 *   cell->face(face)->center()[2] == 1.0 * this->parameters.scale)
5673 *   {
5674 *   if (cell->face(face)->center()[0] < 0.5 * this->parameters.scale &&
5675 *   cell->face(face)->center()[1] < 0.5 * this->parameters.scale)
5676 *   cell->face(face)->set_boundary_id(100);
5677 *   else
5678 *   cell->face(face)->set_boundary_id(101);
5679 *   }
5680 *   }
5681 *   }
5682 *  
5683 *   virtual void
5684 *   define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
5685 *   {
5686 *   tracked_vertices[0][0] = 0.0*this->parameters.scale;
5687 *   tracked_vertices[0][1] = 0.0*this->parameters.scale;
5688 *   tracked_vertices[0][2] = 1.0*this->parameters.scale;
5689 *  
5690 *   tracked_vertices[1][0] = 0.0*this->parameters.scale;
5691 *   tracked_vertices[1][1] = 0.0*this->parameters.scale;
5692 *   tracked_vertices[1][2] = 0.0*this->parameters.scale;
5693 *   }
5694 *  
5695 *   virtual void
5696 *   make_dirichlet_constraints(AffineConstraints<double> &constraints) override
5697 *   {
5698 *   if (this->time.get_timestep() < 2)
5699 *   {
5700 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5701 *   101,
5702 *   Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5703 *   constraints,
5704 *   (this->fe.component_mask(this->pressure)));
5705 *   }
5706 *   else
5707 *   {
5708 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5709 *   101,
5710 *   Functions::ZeroFunction<dim>(this->n_components),
5711 *   constraints,
5712 *   (this->fe.component_mask(this->pressure)));
5713 *   }
5714 *  
5715 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5716 *   0,
5717 *   Functions::ZeroFunction<dim>(this->n_components),
5718 *   constraints,
5719 *   this->fe.component_mask(this->x_displacement));
5720 *  
5721 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5722 *   1,
5723 *   Functions::ZeroFunction<dim>(this->n_components),
5724 *   constraints,
5725 *   this->fe.component_mask(this->x_displacement));
5726 *  
5727 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5728 *   2,
5729 *   Functions::ZeroFunction<dim>(this->n_components),
5730 *   constraints,
5731 *   this->fe.component_mask(this->y_displacement));
5732 *  
5733 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5734 *   3,
5735 *   Functions::ZeroFunction<dim>(this->n_components),
5736 *   constraints,
5737 *   this->fe.component_mask(this->y_displacement));
5738 *  
5739 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5740 *   4,
5741 *   Functions::ZeroFunction<dim>(this->n_components),
5742 *   constraints,
5743 *   ( this->fe.component_mask(this->x_displacement) |
5744 *   this->fe.component_mask(this->y_displacement) |
5745 *   this->fe.component_mask(this->z_displacement) ));
5746 *   }
5747 *  
5748 *   virtual Tensor<1,dim>
5749 *   get_neumann_traction (const types::boundary_id &boundary_id,
5750 *   const Point<dim> &pt,
5751 *   const Tensor<1,dim> &N) const override
5752 *   {
5753 *   if (this->parameters.load_type == "pressure")
5754 *   {
5755 *   if (boundary_id == 100)
5756 *   {
5757 *   return this->parameters.load * N;
5758 *   }
5759 *   }
5760 *  
5761 *   (void)pt;
5762 *  
5763 *   return Tensor<1,dim>();
5764 *   }
5765 *  
5766 *   virtual double
5767 *   get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5768 *   const Point<dim> &pt) const override
5769 *   {
5770 *   (void)pt;
5771 *   (void)boundary_id;
5772 *   return 0.0;
5773 *   }
5774 *  
5775 *   virtual types::boundary_id
5776 *   get_reaction_boundary_id_for_output() const override
5777 *   {
5778 *   return 100;
5779 *   }
5780 *  
5781 *   virtual std::pair<types::boundary_id,types::boundary_id>
5782 *   get_drained_boundary_id_for_output() const override
5783 *   {
5784 *   return std::make_pair(101,101);
5785 *   }
5786 *  
5787 *   virtual std::vector<double>
5788 *   get_dirichlet_load(const types::boundary_id &boundary_id,
5789 *   const int &direction) const override
5790 *   {
5791 *   std::vector<double> displ_incr(dim, 0.0);
5792 *   (void)boundary_id;
5793 *   (void)direction;
5794 *   AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));
5795 *  
5796 *   return displ_incr;
5797 *   }
5798 *   };
5799 *  
5800 * @endcode
5801 *
5802 *
5803 * <a name="Franceschiniexperiments"></a>
5804 * <h4>Franceschini experiments</h4>
5805 *
5806 * @code
5807 *   template <int dim>
5808 *   class Franceschini2006Consolidation
5809 *   : public Solid<dim>
5810 *   {
5811 *   public:
5812 *   Franceschini2006Consolidation (const Parameters::AllParameters &parameters)
5813 *   : Solid<dim> (parameters)
5814 *   {}
5815 *  
5816 *   virtual ~Franceschini2006Consolidation () {}
5817 *  
5818 *   private:
5819 *   virtual void make_grid() override
5820 *   {
5821 *   const Point<dim-1> mesh_center(0.0, 0.0);
5822 *   const double radius = 0.5;
5823 * @endcode
5824 *
5825 * const double height = 0.27; //8.1 mm for 30 mm radius
5826 *
5827 * @code
5828 *   const double height = 0.23; //6.9 mm for 30 mm radius
5829 *   Triangulation<dim-1> triangulation_in;
5830 *   GridGenerator::hyper_ball( triangulation_in,
5831 *   mesh_center,
5832 *   radius);
5833 *  
5834 *   GridGenerator::extrude_triangulation(triangulation_in,
5835 *   2,
5836 *   height,
5837 *   this->triangulation);
5838 *  
5839 *   const CylindricalManifold<dim> cylinder_3d(2);
5840 *   const types::manifold_id cylinder_id = 0;
5841 *  
5842 *  
5843 *   this->triangulation.set_manifold(cylinder_id, cylinder_3d);
5844 *  
5845 *   for (auto cell : this->triangulation.active_cell_iterators())
5846 *   {
5847 *   for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5848 *   {
5849 *   if (cell->face(face)->at_boundary() == true)
5850 *   {
5851 *   if (cell->face(face)->center()[2] == 0.0)
5852 *   cell->face(face)->set_boundary_id(1);
5853 *  
5854 *   else if (cell->face(face)->center()[2] == height)
5855 *   cell->face(face)->set_boundary_id(2);
5856 *  
5857 *   else
5858 *   {
5859 *   cell->face(face)->set_boundary_id(0);
5860 *   cell->face(face)->set_all_manifold_ids(cylinder_id);
5861 *   }
5862 *   }
5863 *   }
5864 *   }
5865 *  
5866 *   GridTools::scale(this->parameters.scale, this->triangulation);
5867 *   this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5868 *   }
5869 *  
5870 *   virtual void define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
5871 *   {
5872 *   tracked_vertices[0][0] = 0.0*this->parameters.scale;
5873 *   tracked_vertices[0][1] = 0.0*this->parameters.scale;
5874 * @endcode
5875 *
5876 * tracked_vertices[0][2] = 0.27*this->parameters.scale;
5877 *
5878 * @code
5879 *   tracked_vertices[0][2] = 0.23*this->parameters.scale;
5880 *  
5881 *   tracked_vertices[1][0] = 0.0*this->parameters.scale;
5882 *   tracked_vertices[1][1] = 0.0*this->parameters.scale;
5883 *   tracked_vertices[1][2] = 0.0*this->parameters.scale;
5884 *   }
5885 *  
5886 *   virtual void make_dirichlet_constraints(AffineConstraints<double> &constraints) override
5887 *   {
5888 *   if (this->time.get_timestep() < 2)
5889 *   {
5890 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5891 *   1,
5892 *   Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5893 *   constraints,
5894 *   (this->fe.component_mask(this->pressure)));
5895 *  
5896 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5897 *   2,
5898 *   Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5899 *   constraints,
5900 *   (this->fe.component_mask(this->pressure)));
5901 *   }
5902 *   else
5903 *   {
5904 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5905 *   1,
5906 *   Functions::ZeroFunction<dim>(this->n_components),
5907 *   constraints,
5908 *   (this->fe.component_mask(this->pressure)));
5909 *  
5910 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5911 *   2,
5912 *   Functions::ZeroFunction<dim>(this->n_components),
5913 *   constraints,
5914 *   (this->fe.component_mask(this->pressure)));
5915 *   }
5916 *  
5917 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5918 *   0,
5919 *   Functions::ZeroFunction<dim>(this->n_components),
5920 *   constraints,
5921 *   (this->fe.component_mask(this->x_displacement)|
5922 *   this->fe.component_mask(this->y_displacement) ) );
5923 *  
5924 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5925 *   1,
5926 *   Functions::ZeroFunction<dim>(this->n_components),
5927 *   constraints,
5928 *   (this->fe.component_mask(this->x_displacement) |
5929 *   this->fe.component_mask(this->y_displacement) |
5930 *   this->fe.component_mask(this->z_displacement) ));
5931 *  
5932 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5933 *   2,
5934 *   Functions::ZeroFunction<dim>(this->n_components),
5935 *   constraints,
5936 *   (this->fe.component_mask(this->x_displacement) |
5937 *   this->fe.component_mask(this->y_displacement) ));
5938 *   }
5939 *  
5940 *   virtual double
5941 *   get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5942 *   const Point<dim> &pt) const override
5943 *   {
5944 *   (void)pt;
5945 *   (void)boundary_id;
5946 *   return 0.0;
5947 *   }
5948 *  
5949 *   virtual types::boundary_id
5950 *   get_reaction_boundary_id_for_output() const override
5951 *   {
5952 *   return 2;
5953 *   }
5954 *  
5955 *   virtual std::pair<types::boundary_id,types::boundary_id>
5956 *   get_drained_boundary_id_for_output() const override
5957 *   {
5958 *   return std::make_pair(1,2);
5959 *   }
5960 *  
5961 *   virtual std::vector<double>
5962 *   get_dirichlet_load(const types::boundary_id &boundary_id,
5963 *   const int &direction) const override
5964 *   {
5965 *   std::vector<double> displ_incr(dim, 0.0);
5966 *   (void)boundary_id;
5967 *   (void)direction;
5968 *   AssertThrow(false, ExcMessage("Displacement loading not implemented for Franceschini examples."));
5969 *  
5970 *   return displ_incr;
5971 *   }
5972 *  
5973 *   virtual Tensor<1,dim>
5974 *   get_neumann_traction (const types::boundary_id &boundary_id,
5975 *   const Point<dim> &pt,
5976 *   const Tensor<1,dim> &N) const override
5977 *   {
5978 *   if (this->parameters.load_type == "pressure")
5979 *   {
5980 *   if (boundary_id == 2)
5981 *   {
5982 *   return (this->parameters.load * N);
5983 *   /*
5984 *   const double final_load = this->parameters.load;
5985 *   const double final_load_time = 10 * this->time.get_delta_t();
5986 *   const double current_time = this->time.get_current();
5987 *  
5988 *  
5989 *   const double c = final_load_time / 2.0;
5990 *   const double r = 200.0 * 0.03 / c;
5991 *  
5992 *   const double load = final_load * std::exp(r * current_time)
5993 *   / ( std::exp(c * current_time) + std::exp(r * current_time));
5994 *   return load * N;
5995 *   */
5996 *   }
5997 *   }
5998 *  
5999 *   (void)pt;
6000 *  
6001 *   return Tensor<1,dim>();
6002 *   }
6003 *   };
6004 *  
6005 * @endcode
6006 *
6007 *
6008 * <a name="ExamplestoreproduceexperimentsbyBuddayetal2017"></a>
6009 * <h3>Examples to reproduce experiments by Budday et al. 2017</h3>
6010 * We group the definition of the geometry, boundary and loading conditions specific to
6011 * the examples to reproduce experiments by Budday et al. 2017 into specific classes.
6012 *
6013
6014 *
6015 *
6016 * <a name="BaseclassCubegeometryandloadingpattern"></a>
6017 * <h4>Base class: Cube geometry and loading pattern</h4>
6018 *
6019 * @code
6020 *   template <int dim>
6021 *   class BrainBudday2017BaseCube
6022 *   : public Solid<dim>
6023 *   {
6024 *   public:
6025 *   BrainBudday2017BaseCube (const Parameters::AllParameters &parameters)
6026 *   : Solid<dim> (parameters)
6027 *   {}
6028 *  
6029 *   virtual ~BrainBudday2017BaseCube () {}
6030 *  
6031 *   private:
6032 *   virtual void
6033 *   make_grid() override
6034 *   {
6035 *   GridGenerator::hyper_cube(this->triangulation,
6036 *   0.0,
6037 *   1.0,
6038 *   true);
6039 *  
6040 *   typename Triangulation<dim>::active_cell_iterator cell =
6041 *   this->triangulation.begin_active(), endc = this->triangulation.end();
6042 *   for (; cell != endc; ++cell)
6043 *   {
6044 *   for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
6045 *   if (cell->face(face)->at_boundary() == true &&
6046 *   ( cell->face(face)->boundary_id() == 0 ||
6047 *   cell->face(face)->boundary_id() == 1 ||
6048 *   cell->face(face)->boundary_id() == 2 ||
6049 *   cell->face(face)->boundary_id() == 3 ) )
6050 *  
6051 *   cell->face(face)->set_boundary_id(100);
6052 *  
6053 *   }
6054 *  
6055 *   GridTools::scale(this->parameters.scale, this->triangulation);
6056 *   this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
6057 *   }
6058 *  
6059 *   virtual double
6060 *   get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
6061 *   const Point<dim> &pt) const override
6062 *   {
6063 *   (void)pt;
6064 *   (void)boundary_id;
6065 *   return 0.0;
6066 *   }
6067 *  
6068 *   virtual std::pair<types::boundary_id,types::boundary_id>
6069 *   get_drained_boundary_id_for_output() const override
6070 *   {
6071 *   return std::make_pair(100,100);
6072 *   }
6073 *   };
6074 *  
6075 * @endcode
6076 *
6077 *
6078 * <a name="DerivedclassUniaxialboundaryconditions"></a>
6079 * <h4>Derived class: Uniaxial boundary conditions</h4>
6080 *
6081 * @code
6082 *   template <int dim>
6083 *   class BrainBudday2017CubeTensionCompression
6084 *   : public BrainBudday2017BaseCube<dim>
6085 *   {
6086 *   public:
6087 *   BrainBudday2017CubeTensionCompression (const Parameters::AllParameters &parameters)
6088 *   : BrainBudday2017BaseCube<dim> (parameters)
6089 *   {}
6090 *  
6091 *   virtual ~BrainBudday2017CubeTensionCompression () {}
6092 *  
6093 *   private:
6094 *   virtual void
6095 *   define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
6096 *   {
6097 *   tracked_vertices[0][0] = 0.5*this->parameters.scale;
6098 *   tracked_vertices[0][1] = 0.5*this->parameters.scale;
6099 *   tracked_vertices[0][2] = 1.0*this->parameters.scale;
6100 *  
6101 *   tracked_vertices[1][0] = 0.5*this->parameters.scale;
6102 *   tracked_vertices[1][1] = 0.5*this->parameters.scale;
6103 *   tracked_vertices[1][2] = 0.5*this->parameters.scale;
6104 *   }
6105 *  
6106 *   virtual void
6107 *   make_dirichlet_constraints(AffineConstraints<double> &constraints) override
6108 *   {
6109 *   if (this->time.get_timestep() < 2)
6110 *   {
6111 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
6112 *   100,
6113 *   Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
6114 *   constraints,
6115 *   (this->fe.component_mask(this->pressure)));
6116 *   }
6117 *   else
6118 *   {
6119 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6120 *   100,
6121 *   Functions::ZeroFunction<dim>(this->n_components),
6122 *   constraints,
6123 *   (this->fe.component_mask(this->pressure)));
6124 *   }
6125 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6126 *   4,
6127 *   Functions::ZeroFunction<dim>(this->n_components),
6128 *   constraints,
6129 *   this->fe.component_mask(this->z_displacement) );
6130 *  
6131 *   Point<dim> fix_node(0.5*this->parameters.scale, 0.5*this->parameters.scale, 0.0);
6132 *   typename DoFHandler<dim>::active_cell_iterator
6133 *   cell = this->dof_handler_ref.begin_active(), endc = this->dof_handler_ref.end();
6134 *   for (; cell != endc; ++cell)
6135 *   for (unsigned int node = 0; node < GeometryInfo<dim>::vertices_per_cell; ++node)
6136 *   {
6137 *   if ( (abs(cell->vertex(node)[2]-fix_node[2]) < (1e-6 * this->parameters.scale))
6138 *   && (abs(cell->vertex(node)[0]-fix_node[0]) < (1e-6 * this->parameters.scale)))
6139 *   constraints.add_line(cell->vertex_dof_index(node, 0));
6140 *  
6141 *   if ( (abs(cell->vertex(node)[2]-fix_node[2]) < (1e-6 * this->parameters.scale))
6142 *   && (abs(cell->vertex(node)[1]-fix_node[1]) < (1e-6 * this->parameters.scale)))
6143 *   constraints.add_line(cell->vertex_dof_index(node, 1));
6144 *   }
6145 *  
6146 *   if (this->parameters.load_type == "displacement")
6147 *   {
6148 *   const std::vector<double> value = get_dirichlet_load(5,2);
6149 *   FEValuesExtractors::Scalar direction;
6150 *   direction = this->z_displacement;
6151 *  
6152 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6153 *   5,
6154 *   Functions::ConstantFunction<dim>(value[2],this->n_components),
6155 *   constraints,
6156 *   this->fe.component_mask(direction));
6157 *   }
6158 *   }
6159 *  
6160 *   virtual Tensor<1,dim>
6161 *   get_neumann_traction (const types::boundary_id &boundary_id,
6162 *   const Point<dim> &pt,
6163 *   const Tensor<1,dim> &N) const override
6164 *   {
6165 *   if (this->parameters.load_type == "pressure")
6166 *   {
6167 *   if (boundary_id == 5)
6168 *   {
6169 *   const double final_load = this->parameters.load;
6170 *   const double current_time = this->time.get_current();
6171 *   const double final_time = this->time.get_end();
6172 *   const double num_cycles = 3.0;
6173 *  
6174 *   return final_load/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5))) * N;
6175 *   }
6176 *   }
6177 *  
6178 *   (void)pt;
6179 *  
6180 *   return Tensor<1,dim>();
6181 *   }
6182 *  
6183 *   virtual types::boundary_id
6184 *   get_reaction_boundary_id_for_output() const override
6185 *   {
6186 *   return 5;
6187 *   }
6188 *  
6189 *   virtual std::vector<double>
6190 *   get_dirichlet_load(const types::boundary_id &boundary_id,
6191 *   const int &direction) const override
6192 *   {
6193 *   std::vector<double> displ_incr(dim,0.0);
6194 *  
6195 *   if ( (boundary_id == 5) && (direction == 2) )
6196 *   {
6197 *   const double final_displ = this->parameters.load;
6198 *   const double current_time = this->time.get_current();
6199 *   const double final_time = this->time.get_end();
6200 *   const double delta_time = this->time.get_delta_t();
6201 *   const double num_cycles = 3.0;
6202 *   double current_displ = 0.0;
6203 *   double previous_displ = 0.0;
6204 *  
6205 *   if (this->parameters.num_cycle_sets == 1)
6206 *   {
6207 *   current_displ = final_displ/2.0 * (1.0
6208 *   - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5)));
6209 *   previous_displ = final_displ/2.0 * (1.0
6210 *   - std::sin(numbers::PI * (2.0*num_cycles*(current_time-delta_time)/final_time + 0.5)));
6211 *   }
6212 *   else
6213 *   {
6214 *   if ( current_time <= (final_time*1.0/3.0) )
6215 *   {
6216 *   current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6217 *   (2.0*num_cycles*current_time/(final_time*1.0/3.0) + 0.5)));
6218 *   previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6219 *   (2.0*num_cycles*(current_time-delta_time)/(final_time*1.0/3.0) + 0.5)));
6220 *   }
6221 *   else
6222 *   {
6223 *   current_displ = final_displ * (1.0 - std::sin(numbers::PI *
6224 *   (2.0*num_cycles*current_time / (final_time*2.0/3.0)
6225 *   - (num_cycles - 0.5) )));
6226 *   previous_displ = final_displ * (1.0 - std::sin(numbers::PI *
6227 *   (2.0*num_cycles*(current_time-delta_time) / (final_time*2.0/3.0)
6228 *   - (num_cycles - 0.5))));
6229 *   }
6230 *   }
6231 *   displ_incr[2] = current_displ - previous_displ;
6232 *   }
6233 *   return displ_incr;
6234 *   }
6235 *   };
6236 *  
6237 * @endcode
6238 *
6239 *
6240 * <a name="DerivedclassNolateraldisplacementinloadingsurfaces"></a>
6241 * <h4>Derived class: No lateral displacement in loading surfaces</h4>
6242 *
6243 * @code
6244 *   template <int dim>
6245 *   class BrainBudday2017CubeTensionCompressionFullyFixed
6246 *   : public BrainBudday2017BaseCube<dim>
6247 *   {
6248 *   public:
6249 *   BrainBudday2017CubeTensionCompressionFullyFixed (const Parameters::AllParameters &parameters)
6250 *   : BrainBudday2017BaseCube<dim> (parameters)
6251 *   {}
6252 *  
6253 *   virtual ~BrainBudday2017CubeTensionCompressionFullyFixed () {}
6254 *  
6255 *   private:
6256 *   virtual void
6257 *   define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
6258 *   {
6259 *   tracked_vertices[0][0] = 0.5*this->parameters.scale;
6260 *   tracked_vertices[0][1] = 0.5*this->parameters.scale;
6261 *   tracked_vertices[0][2] = 1.0*this->parameters.scale;
6262 *  
6263 *   tracked_vertices[1][0] = 0.5*this->parameters.scale;
6264 *   tracked_vertices[1][1] = 0.5*this->parameters.scale;
6265 *   tracked_vertices[1][2] = 0.5*this->parameters.scale;
6266 *   }
6267 *  
6268 *   virtual void
6269 *   make_dirichlet_constraints(AffineConstraints<double> &constraints) override
6270 *   {
6271 *   if (this->time.get_timestep() < 2)
6272 *   {
6273 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
6274 *   100,
6275 *   Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
6276 *   constraints,
6277 *   (this->fe.component_mask(this->pressure)));
6278 *   }
6279 *   else
6280 *   {
6281 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6282 *   100,
6283 *   Functions::ZeroFunction<dim>(this->n_components),
6284 *   constraints,
6285 *   (this->fe.component_mask(this->pressure)));
6286 *   }
6287 *  
6288 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6289 *   4,
6290 *   Functions::ZeroFunction<dim>(this->n_components),
6291 *   constraints,
6292 *   (this->fe.component_mask(this->x_displacement) |
6293 *   this->fe.component_mask(this->y_displacement) |
6294 *   this->fe.component_mask(this->z_displacement) ));
6295 *  
6296 *  
6297 *   if (this->parameters.load_type == "displacement")
6298 *   {
6299 *   const std::vector<double> value = get_dirichlet_load(5,2);
6300 *   FEValuesExtractors::Scalar direction;
6301 *   direction = this->z_displacement;
6302 *  
6303 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6304 *   5,
6305 *   Functions::ConstantFunction<dim>(value[2],this->n_components),
6306 *   constraints,
6307 *   this->fe.component_mask(direction) );
6308 *  
6309 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6310 *   5,
6311 *   Functions::ZeroFunction<dim>(this->n_components),
6312 *   constraints,
6313 *   (this->fe.component_mask(this->x_displacement) |
6314 *   this->fe.component_mask(this->y_displacement) ));
6315 *   }
6316 *   }
6317 *  
6318 *   virtual Tensor<1,dim>
6319 *   get_neumann_traction (const types::boundary_id &boundary_id,
6320 *   const Point<dim> &pt,
6321 *   const Tensor<1,dim> &N) const override
6322 *   {
6323 *   if (this->parameters.load_type == "pressure")
6324 *   {
6325 *   if (boundary_id == 5)
6326 *   {
6327 *   const double final_load = this->parameters.load;
6328 *   const double current_time = this->time.get_current();
6329 *   const double final_time = this->time.get_end();
6330 *   const double num_cycles = 3.0;
6331 *  
6332 *   return final_load/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5))) * N;
6333 *   }
6334 *   }
6335 *  
6336 *   (void)pt;
6337 *  
6338 *   return Tensor<1,dim>();
6339 *   }
6340 *  
6341 *   virtual types::boundary_id
6342 *   get_reaction_boundary_id_for_output() const override
6343 *   {
6344 *   return 5;
6345 *   }
6346 *  
6347 *   virtual std::vector<double>
6348 *   get_dirichlet_load(const types::boundary_id &boundary_id,
6349 *   const int &direction) const override
6350 *   {
6351 *   std::vector<double> displ_incr(dim,0.0);
6352 *  
6353 *   if ( (boundary_id == 5) && (direction == 2) )
6354 *   {
6355 *   const double final_displ = this->parameters.load;
6356 *   const double current_time = this->time.get_current();
6357 *   const double final_time = this->time.get_end();
6358 *   const double delta_time = this->time.get_delta_t();
6359 *   const double num_cycles = 3.0;
6360 *   double current_displ = 0.0;
6361 *   double previous_displ = 0.0;
6362 *  
6363 *   if (this->parameters.num_cycle_sets == 1)
6364 *   {
6365 *   current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5)));
6366 *   previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*(current_time-delta_time)/final_time + 0.5)));
6367 *   }
6368 *   else
6369 *   {
6370 *   if ( current_time <= (final_time*1.0/3.0) )
6371 *   {
6372 *   current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6373 *   (2.0*num_cycles*current_time/(final_time*1.0/3.0) + 0.5)));
6374 *   previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6375 *   (2.0*num_cycles*(current_time-delta_time)/(final_time*1.0/3.0) + 0.5)));
6376 *   }
6377 *   else
6378 *   {
6379 *   current_displ = final_displ * (1.0 - std::sin(numbers::PI *
6380 *   (2.0*num_cycles*current_time / (final_time*2.0/3.0)
6381 *   - (num_cycles - 0.5) )));
6382 *   previous_displ = final_displ * (1.0 - std::sin(numbers::PI *
6383 *   (2.0*num_cycles*(current_time-delta_time) / (final_time*2.0/3.0)
6384 *   - (num_cycles - 0.5))));
6385 *   }
6386 *   }
6387 *   displ_incr[2] = current_displ - previous_displ;
6388 *   }
6389 *   return displ_incr;
6390 *   }
6391 *   };
6392 *  
6393 * @endcode
6394 *
6395 *
6396 * <a name="DerivedclassNolateralorverticaldisplacementinloadingsurface"></a>
6397 * <h4>Derived class: No lateral or vertical displacement in loading surface</h4>
6398 *
6399 * @code
6400 *   template <int dim>
6401 *   class BrainBudday2017CubeShearFullyFixed
6402 *   : public BrainBudday2017BaseCube<dim>
6403 *   {
6404 *   public:
6405 *   BrainBudday2017CubeShearFullyFixed (const Parameters::AllParameters &parameters)
6406 *   : BrainBudday2017BaseCube<dim> (parameters)
6407 *   {}
6408 *  
6409 *   virtual ~BrainBudday2017CubeShearFullyFixed () {}
6410 *  
6411 *   private:
6412 *   virtual void
6413 *   define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
6414 *   {
6415 *   tracked_vertices[0][0] = 0.75*this->parameters.scale;
6416 *   tracked_vertices[0][1] = 0.5*this->parameters.scale;
6417 *   tracked_vertices[0][2] = 0.0*this->parameters.scale;
6418 *  
6419 *   tracked_vertices[1][0] = 0.25*this->parameters.scale;
6420 *   tracked_vertices[1][1] = 0.5*this->parameters.scale;
6421 *   tracked_vertices[1][2] = 0.0*this->parameters.scale;
6422 *   }
6423 *  
6424 *   virtual void
6425 *   make_dirichlet_constraints(AffineConstraints<double> &constraints) override
6426 *   {
6427 *   if (this->time.get_timestep() < 2)
6428 *   {
6429 *   VectorTools::interpolate_boundary_values(this->dof_handler_ref,
6430 *   100,
6431 *   Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
6432 *   constraints,
6433 *   (this->fe.component_mask(this->pressure)));
6434 *   }
6435 *   else
6436 *   {
6437 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6438 *   100,
6439 *   Functions::ZeroFunction<dim>(this->n_components),
6440 *   constraints,
6441 *   (this->fe.component_mask(this->pressure)));
6442 *   }
6443 *  
6444 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6445 *   5,
6446 *   Functions::ZeroFunction<dim>(this->n_components),
6447 *   constraints,
6448 *   (this->fe.component_mask(this->x_displacement) |
6449 *   this->fe.component_mask(this->y_displacement) |
6450 *   this->fe.component_mask(this->z_displacement) ));
6451 *  
6452 *  
6453 *   if (this->parameters.load_type == "displacement")
6454 *   {
6455 *   const std::vector<double> value = get_dirichlet_load(4,0);
6456 *   FEValuesExtractors::Scalar direction;
6457 *   direction = this->x_displacement;
6458 *  
6459 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6460 *   4,
6461 *   Functions::ConstantFunction<dim>(value[0],this->n_components),
6462 *   constraints,
6463 *   this->fe.component_mask(direction));
6464 *  
6465 *   VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6466 *   4,
6467 *   Functions::ZeroFunction<dim>(this->n_components),
6468 *   constraints,
6469 *   (this->fe.component_mask(this->y_displacement) |
6470 *   this->fe.component_mask(this->z_displacement) ));
6471 *   }
6472 *   }
6473 *  
6474 *   virtual Tensor<1,dim>
6475 *   get_neumann_traction (const types::boundary_id &boundary_id,
6476 *   const Point<dim> &pt,
6477 *   const Tensor<1,dim> &N) const override
6478 *   {
6479 *   if (this->parameters.load_type == "pressure")
6480 *   {
6481 *   if (boundary_id == 4)
6482 *   {
6483 *   const double final_load = this->parameters.load;
6484 *   const double current_time = this->time.get_current();
6485 *   const double final_time = this->time.get_end();
6486 *   const double num_cycles = 3.0;
6487 *   const Tensor<1,3> axis ({0.0,1.0,0.0});
6488 *   const double angle = numbers::PI;
6489 *   static const Tensor< 2, dim, double> R(Physics::Transformations::Rotations::rotation_matrix_3d(axis,angle));
6490 *  
6491 *   return (final_load * (std::sin(2.0*(numbers::PI)*num_cycles*current_time/final_time)) * (R * N));
6492 *   }
6493 *   }
6494 *  
6495 *   (void)pt;
6496 *  
6497 *   return Tensor<1,dim>();
6498 *   }
6499 *  
6500 *   virtual types::boundary_id
6501 *   get_reaction_boundary_id_for_output() const override
6502 *   {
6503 *   return 4;
6504 *   }
6505 *  
6506 *   virtual std::vector<double>
6507 *   get_dirichlet_load(const types::boundary_id &boundary_id,
6508 *   const int &direction) const override
6509 *   {
6510 *   std::vector<double> displ_incr (dim, 0.0);
6511 *  
6512 *   if ( (boundary_id == 4) && (direction == 0) )
6513 *   {
6514 *   const double final_displ = this->parameters.load;
6515 *   const double current_time = this->time.get_current();
6516 *   const double final_time = this->time.get_end();
6517 *   const double delta_time = this->time.get_delta_t();
6518 *   const double num_cycles = 3.0;
6519 *   double current_displ = 0.0;
6520 *   double previous_displ = 0.0;
6521 *  
6522 *   if (this->parameters.num_cycle_sets == 1)
6523 *   {
6524 *   current_displ = final_displ * (std::sin(2.0*(numbers::PI)*num_cycles*current_time/final_time));
6525 *   previous_displ = final_displ * (std::sin(2.0*(numbers::PI)*num_cycles*(current_time-delta_time)/final_time));
6526 *   }
6527 *   else
6528 *   {
6529 *   AssertThrow(false, ExcMessage("Problem type not defined. Budday shear experiments implemented only for one set of cycles."));
6530 *   }
6531 *   displ_incr[0] = current_displ - previous_displ;
6532 *   }
6533 *   return displ_incr;
6534 *   }
6535 *   };
6536 *  
6537 *   }
6538 *  
6539 * @endcode
6540 *
6541 *
6542 * <a name="Mainfunction"></a>
6543 * <h3>Main function</h3>
6544 * Lastly we provide the main driver function which is similar to the other tutorials.
6545 *
6546 * @code
6547 *   int main (int argc, char *argv[])
6548 *   {
6549 *   using namespace dealii;
6550 *   using namespace NonLinearPoroViscoElasticity;
6551 *  
6552 *   const unsigned int n_tbb_processes = 1;
6553 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, n_tbb_processes);
6554 *  
6555 *   try
6556 *   {
6557 *   Parameters::AllParameters parameters ("parameters.prm");
6558 *   if (parameters.geom_type == "Ehlers_tube_step_load")
6559 *   {
6560 *   VerificationEhlers1999StepLoad<3> solid_3d(parameters);
6561 *   solid_3d.run();
6562 *   }
6563 *   else if (parameters.geom_type == "Ehlers_tube_increase_load")
6564 *   {
6565 *   VerificationEhlers1999IncreaseLoad<3> solid_3d(parameters);
6566 *   solid_3d.run();
6567 *   }
6568 *   else if (parameters.geom_type == "Ehlers_cube_consolidation")
6569 *   {
6570 *   VerificationEhlers1999CubeConsolidation<3> solid_3d(parameters);
6571 *   solid_3d.run();
6572 *   }
6573 *   else if (parameters.geom_type == "Franceschini_consolidation")
6574 *   {
6575 *   Franceschini2006Consolidation<3> solid_3d(parameters);
6576 *   solid_3d.run();
6577 *   }
6578 *   else if (parameters.geom_type == "Budday_cube_tension_compression")
6579 *   {
6580 *   BrainBudday2017CubeTensionCompression<3> solid_3d(parameters);
6581 *   solid_3d.run();
6582 *   }
6583 *   else if (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")
6584 *   {
6585 *   BrainBudday2017CubeTensionCompressionFullyFixed<3> solid_3d(parameters);
6586 *   solid_3d.run();
6587 *   }
6588 *   else if (parameters.geom_type == "Budday_cube_shear_fully_fixed")
6589 *   {
6590 *   BrainBudday2017CubeShearFullyFixed<3> solid_3d(parameters);
6591 *   solid_3d.run();
6592 *   }
6593 *   else
6594 *   {
6595 *   AssertThrow(false, ExcMessage("Problem type not defined. Current setting: " + parameters.geom_type));
6596 *   }
6597 *  
6598 *   }
6599 *   catch (std::exception &exc)
6600 *   {
6601 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
6602 *   {
6603 *   std::cerr << std::endl << std::endl
6604 *   << "----------------------------------------------------"
6605 *   << std::endl;
6606 *   std::cerr << "Exception on processing: " << std::endl << exc.what()
6607 *   << std::endl << "Aborting!" << std::endl
6608 *   << "----------------------------------------------------"
6609 *   << std::endl;
6610 *  
6611 *   return 1;
6612 *   }
6613 *   }
6614 *   catch (...)
6615 *   {
6616 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
6617 *   {
6618 *   std::cerr << std::endl << std::endl
6619 *   << "----------------------------------------------------"
6620 *   << std::endl;
6621 *   std::cerr << "Unknown exception!" << std::endl << "Aborting!"
6622 *   << std::endl
6623 *   << "----------------------------------------------------"
6624 *   << std::endl;
6625 *   return 1;
6626 *   }
6627 *   }
6628 *   return 0;
6629 *   }
6630 * @endcode
6631
6632
6633*/
FE_Q
Definition fe_q.h:551
MPI_Comm
Point
Definition point.h:112
TimerOutput::wall_times
@ wall_times
Definition timer.h:652
TrilinosWrappers::MPI::Vector::reinit
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
Definition trilinos_vector.cc:196
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:586
first
Point< 2 > first
Definition grid_out.cc:4615
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1703
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition dof_handler.h:441
linear_operator
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
Definition linear_operator.h:1408
MeshWorker::loop
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:439
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition dof_tools_sparsity.cc:56
UpdateFlags
UpdateFlags
Definition fe_update_flags.h:65
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:76
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition fe_update_flags.h:142
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition fe_update_flags.h:135
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:82
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition fe_update_flags.h:128
Algorithms::Events::initial
const Event initial
Definition event.cc:65
DerivativeApproximation::internal::approximate
void approximate(SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Definition derivative_approximation.cc:921
DoFRenumbering::component_wise
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
Definition dof_renumbering.cc:682
DoFRenumbering::Cuthill_McKee
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 >())
Definition dof_renumbering.cc:367
DoFTools::locally_owned_dofs_per_subdomain
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_tools.cc:1376
DoFTools::count_dofs_per_fe_block
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 >())
Definition dof_tools.cc:1973
DoFTools::locally_relevant_dofs_per_subdomain
std::vector< IndexSet > locally_relevant_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_tools.cc:1471
DoFTools::count_dofs_with_subdomain_association
unsigned int count_dofs_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
Definition dof_tools.cc:1637
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:2202
GridTools::count_cells_with_subdomain_association
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition grid_tools.cc:4411
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria)
Definition grid_tools.cc:144
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition tria_iterator_base.h:38
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:58
LAPACKSupport::diagonal
@ diagonal
Matrix is diagonal.
Definition lapack_support.h:119
LocalIntegrators::Advection::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:75
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:472
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:50
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
TensorAccessors::extract
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
Definition tensor_accessors.h:219
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition trilinos_sparse_matrix.cc:65
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:150
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161
Utilities::MPI::reduce
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
WorkStream::internal::tbb_no_coloring::run
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)
Definition work_stream.h:472
deal_II_exceptions::internals::abort
void abort(const ExceptionBase &exc) noexcept
Definition exceptions.cc:460
internal::MatrixFreeFunctions::check
bool check(const ConstraintKinds kind_in, const unsigned int dim)
Definition hanging_nodes_internal.h:93
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition vector_operations_internal.h:66
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition p4est_wrappers.cc:339
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:71
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6008
types::material_id
unsigned int material_id
Definition types.h:164
triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition p4est_wrappers.cc:69
determinant
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition symmetric_tensor.h:2800
symmetrize
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition symmetric_tensor.h:3536
invert
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Definition symmetric_tensor.h:3458
SymmetricTensorEigenvectorMethod
SymmetricTensorEigenvectorMethod
Definition symmetric_tensor.h:3259
SymmetricTensorEigenvectorMethod::jacobi
@ jacobi
SymmetricTensorEigenvectorMethod::ql_implicit_shifts
@ ql_implicit_shifts
Generated by doxygen 1.9.7