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