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