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