678 *
and import all the deal.II function and class names into it:
681 * namespace NonLinearPoroViscoElasticity
683 * using namespace dealii;
688 * <a name=
"nonlinear-poro-viscoelasticity.cc-Runtimeparameters"></a>
702 * <a name=
"nonlinear-poro-viscoelasticity.cc-FiniteElementsystem"></a>
724 *
prm.enter_subsection(
"Finite element system");
726 *
prm.declare_entry(
"Polynomial degree displ",
"2",
728 *
"Displacement system polynomial order");
730 *
prm.declare_entry(
"Polynomial degree pore",
"1",
732 *
"Pore pressure system polynomial order");
734 *
prm.declare_entry(
"Quadrature order",
"3",
736 *
"Gauss quadrature order");
738 *
prm.leave_subsection();
743 *
prm.enter_subsection(
"Finite element system");
747 *
quad_order = prm.get_integer(
"Quadrature order");
749 *
prm.leave_subsection();
755 * <a name=
"nonlinear-poro-viscoelasticity.cc-Geometry"></a>
781 *
prm.enter_subsection(
"Geometry");
783 *
prm.declare_entry(
"Geometry type",
"Ehlers_tube_step_load",
785 *
"|Ehlers_tube_increase_load"
786 *
"|Ehlers_cube_consolidation"
787 *
"|Franceschini_consolidation"
788 *
"|Budday_cube_tension_compression"
789 *
"|Budday_cube_tension_compression_fully_fixed"
790 *
"|Budday_cube_shear_fully_fixed"),
791 *
"Type of geometry used. "
792 *
"For Ehlers verification examples see Ehlers and Eipper (1999). "
793 *
"For Franceschini brain consolidation see Franceschini et al. (2006)"
794 *
"For Budday brain examples see Budday et al. (2017)");
796 *
prm.declare_entry(
"Global refinement",
"1",
798 *
"Global refinement level");
800 *
prm.declare_entry(
"Grid scale",
"1.0",
802 *
"Global grid scaling factor");
804 *
prm.declare_entry(
"Load type",
"pressure",
806 *
"Type of loading");
808 *
prm.declare_entry(
"Load value",
"-7.5e+6",
812 *
prm.declare_entry(
"Number of cycle sets",
"1",
814 *
"Number of times each set of 3 cycles is repeated, only for "
815 *
"Budday_cube_tension_compression and Budday_cube_tension_compression_fully_fixed. "
816 *
"Load value is doubled in second set, load rate is kept constant."
817 *
"Final time indicates end of second cycle set.");
819 *
prm.declare_entry(
"Fluid flow value",
"0.0",
821 *
"Prescribed fluid flow. Not implemented in any example yet.");
823 *
prm.declare_entry(
"Drained pressure",
"0.0",
825 *
"Increase of pressure value at drained boundary w.r.t the atmospheric pressure.");
827 *
prm.leave_subsection();
832 *
prm.enter_subsection(
"Geometry");
836 *
scale = prm.get_double(
"Grid scale");
838 *
load = prm.get_double(
"Load value");
840 *
fluid_flow = prm.get_double(
"Fluid flow value");
843 *
prm.leave_subsection();
849 * <a name=
"nonlinear-poro-viscoelasticity.cc-Materials"></a>
901 *
prm.enter_subsection(
"Material properties");
903 *
prm.declare_entry(
"material",
"Neo-Hooke",
905 *
"Type of material used in the problem");
907 *
prm.declare_entry(
"lambda",
"8.375e6",
909 *
"First Lamé parameter for extension function related to compactation point in solid material [Pa].");
911 *
prm.declare_entry(
"shear modulus",
"5.583e6",
913 *
"shear modulus for Neo-Hooke materials [Pa].");
915 *
prm.declare_entry(
"eigen solver",
"QL Implicit Shifts",
917 *
"The type of eigen solver to be used for Ogden and visco-Ogden models.");
919 *
prm.declare_entry(
"mu1",
"0.0",
921 *
"Shear material parameter 'mu1' for Ogden material [Pa].");
923 *
prm.declare_entry(
"mu2",
"0.0",
925 *
"Shear material parameter 'mu2' for Ogden material [Pa].");
927 *
prm.declare_entry(
"mu3",
"0.0",
929 *
"Shear material parameter 'mu1' for Ogden material [Pa].");
931 *
prm.declare_entry(
"alpha1",
"1.0",
933 *
"Stiffness material parameter 'alpha1' for Ogden material [-].");
935 *
prm.declare_entry(
"alpha2",
"1.0",
937 *
"Stiffness material parameter 'alpha2' for Ogden material [-].");
939 *
prm.declare_entry(
"alpha3",
"1.0",
941 *
"Stiffness material parameter 'alpha3' for Ogden material [-].");
943 *
prm.declare_entry(
"mu1_1",
"0.0",
945 *
"Shear material parameter 'mu1' for first viscous mode in Ogden material [Pa].");
947 *
prm.declare_entry(
"mu2_1",
"0.0",
949 *
"Shear material parameter 'mu2' for first viscous mode in Ogden material [Pa].");
951 *
prm.declare_entry(
"mu3_1",
"0.0",
953 *
"Shear material parameter 'mu1' for first viscous mode in Ogden material [Pa].");
955 *
prm.declare_entry(
"alpha1_1",
"1.0",
957 *
"Stiffness material parameter 'alpha1' for first viscous mode in Ogden material [-].");
959 *
prm.declare_entry(
"alpha2_1",
"1.0",
961 *
"Stiffness material parameter 'alpha2' for first viscous mode in Ogden material [-].");
963 *
prm.declare_entry(
"alpha3_1",
"1.0",
965 *
"Stiffness material parameter 'alpha3' for first viscous mode in Ogden material [-].");
967 *
prm.declare_entry(
"viscosity_1",
"1e-10",
969 *
"Deformation-independent viscosity parameter 'eta_1' for first viscous mode in Ogden material [-].");
971 *
prm.declare_entry(
"seepage definition",
"Ehlers",
973 *
"Type of formulation used to define the seepage velocity in the problem. "
974 *
"Choose between Markert formulation of deformation-dependent intrinsic permeability "
975 *
"and Ehlers formulation of deformation-dependent Darcy flow coefficient.");
977 *
prm.declare_entry(
"initial solid volume fraction",
"0.67",
979 *
"Initial porosity (solid volume fraction, 0 < n_0s < 1)");
981 *
prm.declare_entry(
"kappa",
"0.0",
983 *
"Deformation-dependency control parameter for specific permeability (kappa >= 0)");
985 *
prm.declare_entry(
"initial intrinsic permeability",
"0.0",
987 *
"Initial intrinsic permeability parameter [m^2] (isotropic permeability). To be used with Markert formulation.");
989 *
prm.declare_entry(
"fluid viscosity",
"0.0",
991 *
"Effective shear viscosity parameter of the fluid [Pa·s, (N·s)/m^2]. To be used with Markert formulation.");
993 *
prm.declare_entry(
"initial Darcy coefficient",
"1.0e-4",
995 *
"Initial Darcy flow coefficient [m/s] (isotropic permeability). To be used with Ehlers formulation.");
997 *
prm.declare_entry(
"fluid weight",
"1.0e4",
999 *
"Effective weight of the fluid [N/m^3]. To be used with Ehlers formulation.");
1001 *
prm.declare_entry(
"gravity term",
"false",
1003 *
"Gravity term considered (true) or neglected (false)");
1005 *
prm.declare_entry(
"fluid density",
"1.0",
1007 *
"Real (or effective) density of the fluid");
1009 *
prm.declare_entry(
"solid density",
"1.0",
1011 *
"Real (or effective) density of the solid");
1013 *
prm.declare_entry(
"gravity direction",
"2",
1015 *
"Direction of gravity (unit vector 0 for x, 1 for y, 2 for z)");
1017 *
prm.declare_entry(
"gravity value",
"-9.81",
1019 *
"Value of gravity (be careful to have consistent units!)");
1021 *
prm.leave_subsection();
1026 *
prm.enter_subsection(
"Material properties");
1034 *
lambda = prm.get_double(
"lambda");
1035 *
mu = prm.get_double(
"shear modulus");
1054 *
fluid_type = prm.get(
"seepage definition");
1055 *
solid_vol_frac = prm.get_double(
"initial solid volume fraction");
1060 *
weight_FR = prm.get_double(
"fluid weight");
1067 *
density_FR = prm.get_double(
"fluid density");
1068 *
density_SR = prm.get_double(
"solid density");
1074 *
"'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0."));
1078 *
"'initial Darcy coefficient' and 'fluid weight' greater than 0.0."));
1090 *
prm.leave_subsection();
1096 * <a name=
"nonlinear-poro-viscoelasticity.cc-Nonlinearsolver"></a>
1105 *
struct NonlinearSolver
1121 *
prm.enter_subsection(
"Nonlinear solver");
1123 *
prm.declare_entry(
"Max iterations Newton-Raphson",
"15",
1125 *
"Number of Newton-Raphson iterations allowed");
1127 *
prm.declare_entry(
"Tolerance force",
"1.0e-8",
1129 *
"Force residual tolerance");
1131 *
prm.declare_entry(
"Tolerance displacement",
"1.0e-6",
1133 *
"Displacement error tolerance");
1135 *
prm.declare_entry(
"Tolerance pore pressure",
"1.0e-6",
1137 *
"Pore pressure error tolerance");
1139 *
prm.leave_subsection();
1144 *
prm.enter_subsection(
"Nonlinear solver");
1147 *
tol_f = prm.get_double(
"Tolerance force");
1148 *
tol_u = prm.get_double(
"Tolerance displacement");
1149 *
tol_p_fluid = prm.get_double(
"Tolerance pore pressure");
1151 *
prm.leave_subsection();
1157 * <a name=
"nonlinear-poro-viscoelasticity.cc-Time"></a>
1175 *
prm.enter_subsection(
"Time");
1177 *
prm.declare_entry(
"End time",
"10.0",
1181 *
prm.declare_entry(
"Time step size",
"0.002",
1183 *
"Time step size. The value must be larger than the displacement error tolerance defined.");
1185 *
prm.leave_subsection();
1190 *
prm.enter_subsection(
"Time");
1192 *
end_time = prm.get_double(
"End time");
1193 *
delta_t = prm.get_double(
"Time step size");
1195 *
prm.leave_subsection();
1202 * <a name=
"nonlinear-poro-viscoelasticity.cc-Output"></a>
1223 *
prm.enter_subsection(
"Output parameters");
1225 *
prm.declare_entry(
"Output files",
"true",
1227 *
"Paraview output files to generate.");
1228 *
prm.declare_entry(
"Time step number output",
"1",
1230 *
"Output data for time steps multiple of the given "
1231 *
"integer value.");
1232 *
prm.declare_entry(
"Averaged results",
"nodes",
1234 *
"Output data associated with integration point values"
1235 *
" averaged on elements or on nodes.");
1237 *
prm.leave_subsection();
1242 *
prm.enter_subsection(
"Output parameters");
1246 *
outtype = prm.get(
"Averaged results");
1248 *
prm.leave_subsection();
1254 * <a name=
"nonlinear-poro-viscoelasticity.cc-Allparameters"></a>
1262 *
public NonlinearSolver,
1278 *
declare_parameters(prm);
1280 *
parse_parameters(prm);
1285 *
FESystem::declare_parameters(prm);
1286 *
Geometry::declare_parameters(prm);
1287 *
Materials::declare_parameters(prm);
1288 *
NonlinearSolver::declare_parameters(prm);
1289 *
Time::declare_parameters(prm);
1290 *
OutputParam::declare_parameters(prm);
1295 *
FESystem::parse_parameters(prm);
1296 *
Geometry::parse_parameters(prm);
1297 *
Materials::parse_parameters(prm);
1298 *
NonlinearSolver::parse_parameters(prm);
1299 *
Time::parse_parameters(prm);
1300 *
OutputParam::parse_parameters(prm);
1307 * <a name=
"nonlinear-poro-viscoelasticity.cc-Timeclass"></a>
1317 *
const double delta_t)
1354 *
const double delta_t;
1360 * <a name=
"nonlinear-poro-viscoelasticity.cc-Constitutiveequationforthesolidcomponentofthebiphasicmaterial"></a>
1366 * <a name=
"nonlinear-poro-viscoelasticity.cc-Baseclassgenerichyperelasticmaterial"></a>
1383 *
template <
int dim,
typename NumberType = Sacado::Fad::DFad<
double> >
1390 *
n_OS (parameters.solid_vol_frac),
1391 *
lambda (parameters.lambda),
1435 *
const double n_OS;
1451 *
return ( NumberType(lambda * (1.0-
n_OS)*(1.0-
n_OS)
1462 * <a name=
"nonlinear-poro-viscoelasticity.cc-DerivedclassNeoHookeanhyperelasticmaterial"></a>
1466 *
template <
int dim,
typename NumberType = Sacado::Fad::DFad<
double> >
1514 *
const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1518 *
for (
unsigned int d=0;
d<dim; ++
d)
1521 *
return (
mu*(
B_ev-I) );
1529 * <a name=
"nonlinear-poro-viscoelasticity.cc-DerivedclassOgdenhyperelasticmaterial"></a>
1533 *
template <
int dim,
typename NumberType = Sacado::Fad::DFad<
double> >
1541 *
mu({parameters.mu1_infty,
1542 *
parameters.mu2_infty,
1543 *
parameters.mu3_infty}),
1544 *
alpha({parameters.alpha1_infty,
1545 *
parameters.alpha2_infty,
1546 *
parameters.alpha3_infty})
1558 *
std::vector<double>
mu;
1559 *
std::vector<double>
alpha;
1567 *
const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1574 *
for (
unsigned int i = 0; i < 3; ++i)
1576 *
for (
unsigned int A = 0;
A < dim; ++
A)
1594 * <a name=
"nonlinear-poro-viscoelasticity.cc-DerivedclassSinglemodeOgdenviscoelasticmaterial"></a>
1611 *
parameters.mu2_infty,
1612 *
parameters.mu3_infty}),
1614 *
parameters.alpha2_infty,
1615 *
parameters.alpha3_infty}),
1617 *
parameters.mu2_mode_1,
1618 *
parameters.mu3_mode_1}),
1620 *
parameters.alpha2_mode_1,
1621 *
parameters.alpha3_mode_1}),
1637 *
const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1642 *
for (
int a = 0; a < dim; ++a)
1648 *
const double tolerance = 1
e-8;
1666 *
for (
unsigned int a = 0; a < dim; ++a)
1674 *
for (
unsigned int a = 0; a < dim; ++a)
1677 *
for (
unsigned int a = 0; a < dim; ++a)
1684 *
for (
unsigned int b = 0;
b < dim; ++
b)
1695 *
for (
unsigned int a = 0; a < dim; ++a)
1702 *
AssertThrow(
false, ExcMessage(
"No convergence in local Newton iteration for the "
1703 *
"viscoelastic exponential time integration algorithm."));
1707 *
for (
unsigned int a = 0; a < dim; ++a)
1714 *
for (
unsigned int a = 0; a < dim; ++a)
1717 *
for (
unsigned int a = 0; a < dim; ++a)
1728 *
for (
unsigned int a = 0; a < dim; ++a)
1741 *
for (
unsigned int a = 0; a < dim; ++a)
1742 *
for (
unsigned int b = 0;
b < dim; ++
b)
1781 *
std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
eigen_B;
1788 *
for (
unsigned int i = 0; i < 3; ++i)
1790 *
for (
unsigned int A = 0;
A < dim; ++
A)
1811 *
get_beta_mode_1(std::vector< NumberType > &lambda,
const int &A)
const
1813 *
NumberType
beta = 0.0;
1815 *
for (
unsigned int i = 0; i < 3; ++i)
1818 *
NumberType aux = 0.0;
1819 *
for (
int p = 0; p < dim; ++p)
1834 *
const int &
B )
const
1836 *
NumberType
gamma = 0.0;
1840 *
for (
unsigned int i = 0; i < 3; ++i)
1842 *
NumberType aux = 0.0;
1843 *
for (
int p = 0; p < dim; ++p)
1846 *
aux *= 1.0/(dim*dim);
1855 *
for (
unsigned int i = 0; i < 3; ++i)
1857 *
NumberType aux = 0.0;
1858 *
for (
int p = 0; p < dim; ++p)
1861 *
aux *= 1.0/(dim*dim);
1878 * <a name=
"nonlinear-poro-viscoelasticity.cc-Constitutiveequationforthefluidcomponentofthebiphasicmaterial"></a>
1886 *
template <
int dim,
typename NumberType = Sacado::Fad::DFad<
double> >
1893 *
n_OS(parameters.solid_vol_frac),
1926 *
"Material_Darcy_Fluid --> Only Markert "
1927 *
"and Ehlers formulations have been implemented."));
1955 *
"Material_Darcy_Fluid --> Only Markert and Ehlers "
1956 *
"formulations have been implemented."));
1965 *
const double n_OS;
2003 *
const NumberType fraction = (1.0 - (
n_OS /
det_F) )/(1.0 -
n_OS);
2026 * <a name=
"nonlinear-poro-viscoelasticity.cc-Quadraturepointhistory"></a>
2036 *
template <
int dim,
typename NumberType = Sacado::Fad::DFad<
double> >
2049 *
if (parameters.mat_type ==
"Neo-Hooke")
2051 *
else if (parameters.mat_type ==
"Ogden")
2053 *
else if (parameters.mat_type ==
"visco-Ogden")
2117 *
if (parameters.gravity_term ==
true)
2123 *
= parameters.density_SR * parameters.solid_vol_frac
2124 *
+ parameters.density_FR
2125 *
* (
det_F_AD - parameters.solid_vol_frac);
2128 *
gravity_vector[parameters.gravity_direction] = parameters.gravity_value;
2135 *
std::shared_ptr< Material_Hyperelastic<dim, NumberType> >
solid_material;
2136 *
std::shared_ptr< Material_Darcy_Fluid<dim, NumberType> >
fluid_material;
2142 * <a name=
"nonlinear-poro-viscoelasticity.cc-Nonlinearporoviscoelasticsolid"></a>
2222 *
virtual std::pair<types::boundary_id,types::boundary_id>
2226 *
const int &direction)
const = 0;
2270 *
const double current_time,
2273 *
const double current_time,
2356 *
const unsigned int dofs_per_cell;
2375 *
static const unsigned int n_blocks = 2;
2376 *
static const unsigned int n_components = dim+1;
2436 *
const unsigned int n_q_points;
2504 *
if (
rhs.p_fluid != 0.0)
2548 * <a name=
"nonlinear-poro-viscoelasticity.cc-ImplementationofthecodeSolidcodeclass"></a>
2551 * <a name=
"nonlinear-poro-viscoelasticity.cc-Publicinterface"></a>
2556 *
template <
int dim>
2562 *
pcout(std::cout, this_mpi_process == 0),
2563 *
parameters(parameters),
2565 *
time(parameters.end_time, parameters.delta_t),
2576 *
fe(
FE_Q<dim>(parameters.poly_degree_displ), dim,
2577 *
FE_Q<dim>(parameters.poly_degree_pore), 1 ),
2579 *
dofs_per_cell (fe.dofs_per_cell),
2587 *
qf_cell(parameters.quad_order),
2588 *
qf_face(parameters.quad_order),
2589 *
n_q_points (
qf_cell.size()),
2592 *
Assert(dim==3,
ExcMessage(
"This problem only works in 3 space dimensions."));
2601 *
template <
int dim>
2612 *
template <
int dim>
2628 *
if (this_mpi_process == 0)
2630 *
outfile.open(
"console-output.sol");
2657 *
if (this_mpi_process == 0)
2659 *
pointfile.open(
"data-for-gnuplot.sol");
2665 *
Print results
to output file
2668 *
if (parameters.outfiles_requested ==
"true")
2671 *
time.get_current(),
2676 *
time.get_current(),
2688 *
time.increment_time();
2695 *
pcout <<
"\nSolver:";
2696 *
pcout <<
"\n CST = make constraints";
2697 *
pcout <<
"\n ASM_SYS = assemble system";
2698 *
pcout <<
"\n SLV = linear solver \n";
2706 *
outfile <<
"\n CST = make constraints";
2707 *
outfile <<
"\n ASM_SYS = assemble system";
2708 *
outfile <<
"\n SLV = linear solver \n";
2710 *
while ( (time.get_end() - time.get_current()) > -1.0*parameters.tol_u )
2745 *
if (( (time.get_timestep()%parameters.timestep_output) == 0 )
2746 *
&& (parameters.outfiles_requested ==
"true") )
2749 *
time.get_current(),
2754 *
time.get_current(),
2764 *
time.increment_time();
2772 *
if (this_mpi_process == 0)
2790 * <a name=
"nonlinear-poro-viscoelasticity.cc-Privateinterface"></a>
2797 *
template <
int dim>
2802 *
std::vector<types::global_dof_index> local_dof_indices;
2808 *
local_dof_indices(dofs_per_cell)
2823 *
template <
int dim>
2824 *
template <
typename NumberType>
2842 *
std::vector<NumberType> local_dof_values;
2853 *
std::vector<std::vector<Tensor<1,dim>>>
Nx;
2854 *
std::vector<std::vector<double>>
Nx_p_fluid;
2860 *
std::vector<std::vector<Tensor<2,dim, NumberType>>>
grad_Nx;
2861 *
std::vector<std::vector<SymmetricTensor<2,dim, NumberType>>>
symm_grad_Nx;
2862 *
std::vector<std::vector<Tensor<1,dim, NumberType>>>
grad_Nx_p_fluid;
2872 *
local_dof_values(
fe_cell.dofs_per_cell),
2888 *
rhs.fe_values_ref.get_quadrature(),
2889 *
rhs.fe_values_ref.get_update_flags()),
2891 *
rhs.fe_face_values_ref.get_quadrature(),
2892 *
rhs.fe_face_values_ref.get_update_flags()),
2893 *
local_dof_values(
rhs.local_dof_values),
2907 *
const unsigned int n_q_points =
Nx_p_fluid.size();
2908 *
const unsigned int n_dofs_per_cell =
Nx_p_fluid[0].size();
2912 *
for (
unsigned int k = 0;
k < n_dofs_per_cell; ++
k)
2914 *
local_dof_values[
k] = 0.0;
2935 *
for (
unsigned int k = 0;
k < n_dofs_per_cell; ++
k)
2958 *
template <
int dim>
2961 *
pcout <<
" CST " << std::flush;
2962 *
outfile <<
" CST " << std::flush;
2970 *
constraints.
clear();
2976 *
if (constraints.is_inhomogeneously_constrained(i) ==
true)
2977 *
constraints.set_inhomogeneity(i,0.0);
2979 *
constraints.close();
2987 *
template <
int dim>
2991 *
timerfile.enter_subsection(
"Setup system");
3026 *
locally_owned_dofs.
clear();
3056 *
pcout <<
"\nTriangulation:\n"
3057 *
<<
" Number of active cells: "
3059 *
<<
" (by partition:";
3061 *
pcout << (p==0 ?
' ' :
'+')
3065 *
pcout <<
" Number of degrees of freedom: "
3067 *
<<
" (by partition:";
3069 *
pcout << (p==0 ?
' ' :
'+')
3073 *
pcout <<
" Number of degrees of freedom per block: "
3074 *
<<
"[n_u, n_p_fluid] = ["
3086 *
outfile <<
"\nTriangulation:\n"
3087 *
<<
" Number of active cells: "
3089 *
<<
" (by partition:";
3091 *
outfile << (p==0 ?
' ' :
'+')
3095 *
outfile <<
" Number of degrees of freedom: "
3097 *
<<
" (by partition:";
3099 *
outfile << (p==0 ?
' ' :
'+')
3103 *
outfile <<
" Number of degrees of freedom per block: "
3104 *
<<
"[n_u, n_p_fluid] = ["
3118 *
for (
unsigned int ii = 0;
ii < n_components; ++
ii)
3119 *
for (
unsigned int jj = 0;
jj < n_components; ++
jj)
3139 *
mpi_communicator);
3142 *
false, this_mpi_process);
3167 *
mpi_communicator);
3169 *
false, this_mpi_process);
3172 *
system_rhs_nb.reinit(locally_owned_dofs, mpi_communicator);
3191 *
template <
int dim>
3197 *
for (
unsigned int k = 0;
k < fe.dofs_per_cell; ++
k)
3199 *
const unsigned int k_group = fe.system_to_base_index(
k).first.first;
3216 *
template <
int dim>
3219 *
pcout <<
"\nSetting up quadrature point data..." << std::endl;
3220 *
outfile <<
"\nSetting up quadrature point data..." << std::endl;
3240 *
for (; cell!=
endc; ++cell)
3242 *
Assert(cell->is_locally_owned(), ExcInternalError());
3243 *
Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
3245 *
const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >
3247 *
Assert(
lqph.size() == n_q_points, ExcInternalError());
3259 *
template <
int dim>
3267 *
pcout << std::endl
3269 *
<< time.get_timestep()
3271 *
<< time.get_current()
3276 *
<< time.get_timestep()
3278 *
<< time.get_current()
3377 *
&&
system_rhs.l2_norm() <= parameters.tol_f) )
3379 *
pcout <<
"\n ***** CONVERGED! ***** "
3387 *
<<
" " << std::endl;
3388 *
outfile <<
"\n ***** CONVERGED! ***** "
3396 *
<<
" " << std::endl;
3446 *
pcout <<
" | " << std::fixed << std::setprecision(3)
3447 *
<< std::setw(7) << std::scientific
3455 *
<<
" " << std::endl;
3457 *
outfile <<
" | " << std::fixed << std::setprecision(3)
3458 *
<< std::setw(7) << std::scientific
3466 *
<<
" " << std::endl;
3491 *
template <
int dim>
3494 *
static const unsigned int l_width = 120;
3496 *
for (
unsigned int i = 0; i <
l_width; ++i)
3502 *
pcout << std::endl;
3505 *
pcout <<
"\n SOLVER STEP | SYS_RES "
3506 *
<<
"RES_NORM RES_U RES_P "
3507 *
<<
"NU_NORM NU_U NU_P " << std::endl;
3508 *
outfile <<
"\n SOLVER STEP | SYS_RES "
3509 *
<<
"RES_NORM RES_U RES_P "
3510 *
<<
"NU_NORM NU_U NU_P " << std::endl;
3512 *
for (
unsigned int i = 0; i <
l_width; ++i)
3517 *
pcout << std::endl << std::endl;
3518 *
outfile << std::endl << std::endl;
3526 *
template <
int dim>
3529 *
static const unsigned int l_width = 120;
3531 *
for (
unsigned int i = 0; i <
l_width; ++i)
3536 *
pcout << std::endl << std::endl;
3537 *
outfile << std::endl << std::endl;
3539 *
pcout <<
"Relative errors:" << std::endl
3540 *
<<
"Displacement: "
3542 *
<<
"Force (displ): "
3544 *
<<
"Pore pressure: "
3546 *
<<
"Force (pore): "
3548 *
outfile <<
"Relative errors:" << std::endl
3549 *
<<
"Displacement: "
3551 *
<<
"Force (displ): "
3553 *
<<
"Pore pressure: "
3555 *
<<
"Force (pore): "
3564 *
template <
int dim>
3580 *
template <
int dim>
3599 *
template <
int dim>
3625 *
template <
int dim>
3629 *
timerfile.enter_subsection(
"Assemble system");
3630 *
pcout <<
" ASM_SYS " << std::flush;
3631 *
outfile <<
" ASM_SYS " << std::flush;
3665 *
for (; cell !=
endc; ++cell)
3667 *
Assert(cell->is_locally_owned(), ExcInternalError());
3668 *
Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
3689 *
template <
int dim>
3692 *
constraints.distribute_local_to_global(
data.cell_matrix,
3694 *
data.local_dof_indices,
3698 *
constraints.distribute_local_to_global(
data.cell_matrix,
3700 *
data.local_dof_indices,
3710 *
template <
int dim>
3716 *
Assert(cell->is_locally_owned(), ExcInternalError());
3720 *
scratch.fe_values_ref.reinit(cell);
3721 *
cell->get_dof_indices(
data.local_dof_indices);
3728 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
3735 *
scratch.local_dof_values[
k] = scratch.solution_total[
data.local_dof_indices[
k]];
3741 *
scratch.local_dof_values[
k].diff(
k, dofs_per_cell);
3750 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
3752 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
3754 *
const unsigned int k_group = fe.system_to_base_index(
k).first.first;
3758 *
scratch.fe_values_ref[
u_fe].gradient(
k,
q);
3759 *
for (
unsigned int dd = 0;
dd < dim; ++
dd)
3761 *
for (
unsigned int ee = 0;
ee < dim; ++
ee)
3763 *
scratch.solution_grads_u_total[
q][
dd][
ee]
3774 *
scratch.solution_values_p_fluid_total[
q]
3775 *
+= scratch.local_dof_values[
k] *
Nx_p;
3776 *
for (
unsigned int dd = 0;
dd < dim; ++
dd)
3778 *
scratch.solution_grads_p_fluid_total[
q][
dd]
3792 *
const std::vector<std::shared_ptr<const PointHistory<dim, ADNumberType> > >
3794 *
Assert(
lqph.size() == n_q_points, ExcInternalError());
3809 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
3811 *
const unsigned int i_group = fe.system_to_base_index(i).first.first;
3817 *
scratch.grad_Nx[
q_point][i] =
3819 *
scratch.symm_grad_Nx[
q_point][i] =
3824 *
scratch.Nx_p_fluid[
q_point][i] =
3826 *
scratch.grad_Nx_p_fluid[
q_point][i] =
3884 *
const std::vector<Tensor<1,dim>> &
Nu = scratch.Nx[
q_point];
3885 *
const std::vector<SymmetricTensor<2, dim, ADNumberType>>
3887 *
const std::vector<double> &
Np = scratch.Nx_p_fluid[
q_point];
3888 *
const std::vector<Tensor<1, dim, ADNumberType> > &
grad_Np
3889 *
= scratch.grad_Nx_p_fluid[
q_point];
3892 *
const double JxW = scratch.fe_values_ref.JxW(
q_point);
3894 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
3896 *
const unsigned int i_group = fe.system_to_base_index(i).first.first;
3923 *
if (cell->face(face)->at_boundary() ==
true)
3925 *
scratch.fe_face_values_ref.reinit(cell, face);
3930 *
= scratch.fe_face_values_ref.normal_vector(
f_q_point);
3932 *
= scratch.fe_face_values_ref.quadrature_point(
f_q_point);
3940 *
const double JxW_f = scratch.fe_face_values_ref.JxW(
f_q_point);
3942 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
3944 *
const unsigned int i_group = fe.system_to_base_index(i).first.first;
3949 *
= fe.system_to_component_index(i).first;
3951 *
= scratch.fe_face_values_ref.shape_value(i,
f_q_point);
3957 *
= scratch.fe_face_values_ref.shape_value(i,
f_q_point);
3970 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
3974 *
data.cell_rhs(i) -=
R_i.val();
3975 *
for (
unsigned int j=0;
j<dofs_per_cell; ++
j)
3976 *
data.cell_matrix(i,
j) +=
R_i.fastAccessDx(
j);
3985 *
template <
int dim>
3993 *
for (; cell!=
endc; ++cell)
3995 *
Assert(cell->is_locally_owned(), ExcInternalError());
3996 *
Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
3998 *
const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >
4000 *
Assert(
lqph.size() == n_q_points, ExcInternalError());
4012 *
template <
int dim>
4017 *
timerfile.enter_subsection(
"Linear solver");
4018 *
pcout <<
" SLV " << std::flush;
4019 *
outfile <<
" SLV " << std::flush;
4034 *
for (
unsigned int i=0; i<locally_owned_dofs.n_elements(); ++i)
4037 *
= locally_owned_dofs.nth_index_in_set(i);
4051 *
template <
int dim>
4065 *
evaluate_vector_field
4071 *
for (
unsigned int p=0; p<
input_data.solution_gradients.size(); ++p)
4074 *
for (
unsigned int d=0;
d<dim; ++
d)
4092 *
const double current_time,
4201 * -----------------------------------------------------------------------
4209 *
for (
unsigned int i=0; i<dim; ++i)
4217 *
if (parameters.mat_type ==
"Neo-Hooke")
4219 *
else if (parameters.mat_type ==
"Ogden")
4221 *
else if (parameters.mat_type ==
"visco-Ogden")
4224 *
Assert (
false, ExcMessage(
"Material type not implemented"));
4257 *
Assert(cell->is_locally_owned(), ExcInternalError());
4258 *
Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
4261 *
static_cast<int>(cell->material_id());
4289 *
const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4291 *
Assert(
lqph.size() == n_q_points, ExcInternalError());
4306 *
for (
unsigned int i=0; i<dim; ++i)
4307 *
for (
unsigned int j=0;
j<dim; ++
j)
4355 *
but I
haven't invested time in getting it to work properly.
4358 * if (parameters.outtype == "elements")
4360 * for (unsigned int j=0; j<dim; ++j)
4362 * cauchy_stresses_total_elements[j](cell->active_cell_index())
4363 * += ((sigma*basis_vectors[j])*basis_vectors[j])/n_q_points;
4364 * cauchy_stresses_E_elements[j](cell->active_cell_index())
4365 * += ((sigma_E*basis_vectors[j])*basis_vectors[j])/n_q_points;
4366 * stretches_elements[j](cell->active_cell_index())
4367 * += std::sqrt(1.0+2.0*Tensor<0,dim,double>(E_strain[j][j]))
4369 * seepage_velocity_elements[j](cell->active_cell_index())
4370 * += Tensor<0,dim,double>(seepage_vel_AD[j])/n_q_points;
4373 * porous_dissipation_elements(cell->active_cell_index())
4374 * += porous_dissipation/n_q_points;
4375 * viscous_dissipation_elements(cell->active_cell_index())
4376 * += viscous_dissipation/n_q_points;
4377 * solid_vol_fraction_elements(cell->active_cell_index())
4378 * += solid_vol_fraction/n_q_points;
4380 * cauchy_stresses_total_elements[3](cell->active_cell_index())
4381 * += ((sigma*basis_vectors[0])*basis_vectors[1])/n_q_points; //sig_xy
4382 * cauchy_stresses_total_elements[4](cell->active_cell_index())
4383 * += ((sigma*basis_vectors[0])*basis_vectors[2])/n_q_points;//sig_xz
4384 * cauchy_stresses_total_elements[5](cell->active_cell_index())
4385 * += ((sigma*basis_vectors[1])*basis_vectors[2])/n_q_points;//sig_yz
4387 * cauchy_stresses_E_elements[3](cell->active_cell_index())
4388 * += ((sigma_E*basis_vectors[0])* basis_vectors[1])/n_q_points; //sig_xy
4389 * cauchy_stresses_E_elements[4](cell->active_cell_index())
4390 * += ((sigma_E*basis_vectors[0])* basis_vectors[2])/n_q_points;//sig_xz
4391 * cauchy_stresses_E_elements[5](cell->active_cell_index())
4392 * += ((sigma_E*basis_vectors[1])* basis_vectors[2])/n_q_points;//sig_yz
4397 * OUTPUT AVERAGED ON NODES -------------------------------------------
4400 * else if (parameters.outtype == "nodes")
4402 * for (unsigned int v=0; v<(GeometryInfo<dim>::vertices_per_cell); ++v)
4404 * types::global_dof_index local_vertex_indices =
4405 * cell_v->vertex_dof_index(v, 0);
4406 * counter_on_vertices_mpi(local_vertex_indices) += 1;
4407 * for (unsigned int k=0; k<dim; ++k)
4409 * cauchy_stresses_total_vertex_mpi[k](local_vertex_indices)
4410 * += (sigma*basis_vectors[k])*basis_vectors[k];
4411 * cauchy_stresses_E_vertex_mpi[k](local_vertex_indices)
4412 * += (sigma_E*basis_vectors[k])*basis_vectors[k];
4413 * stretches_vertex_mpi[k](local_vertex_indices)
4414 * += std::sqrt(1.0+2.0*Tensor<0,dim,double>(E_strain[k][k]));
4416 * types::global_dof_index local_vertex_vec_indices =
4417 * cell_v_vec->vertex_dof_index(v, k);
4418 * counter_on_vertices_vec_mpi(local_vertex_vec_indices) += 1;
4419 * seepage_velocity_vertex_vec_mpi(local_vertex_vec_indices)
4420 * += Tensor<0,dim,double>(seepage_vel_AD[k]);
4423 * porous_dissipation_vertex_mpi(local_vertex_indices)
4424 * += porous_dissipation;
4425 * viscous_dissipation_vertex_mpi(local_vertex_indices)
4426 * += viscous_dissipation;
4427 * solid_vol_fraction_vertex_mpi(local_vertex_indices)
4428 * += solid_vol_fraction;
4430 * cauchy_stresses_total_vertex_mpi[3](local_vertex_indices)
4431 * += (sigma*basis_vectors[0])*basis_vectors[1]; //sig_xy
4432 * cauchy_stresses_total_vertex_mpi[4](local_vertex_indices)
4433 * += (sigma*basis_vectors[0])*basis_vectors[2];//sig_xz
4434 * cauchy_stresses_total_vertex_mpi[5](local_vertex_indices)
4435 * += (sigma*basis_vectors[1])*basis_vectors[2]; //sig_yz
4437 * cauchy_stresses_E_vertex_mpi[3](local_vertex_indices)
4438 * += (sigma_E*basis_vectors[0])*basis_vectors[1]; //sig_xy
4439 * cauchy_stresses_E_vertex_mpi[4](local_vertex_indices)
4440 * += (sigma_E*basis_vectors[0])*basis_vectors[2];//sig_xz
4441 * cauchy_stresses_E_vertex_mpi[5](local_vertex_indices)
4442 * += (sigma_E*basis_vectors[1])*basis_vectors[2]; //sig_yz
4447 * ---------------------------------------------------------------
4450 * } //end gauss point loop
4455 * Different nodes might have different amount of contributions, e.g.,
4456 * corner nodes have less integration points contributing to the averaged.
4457 * This is why we need a counter and divide at the end, outside the cell loop.
4460 * if (parameters.outtype == "nodes")
4462 * for (unsigned int d=0; d<(vertex_handler_ref.n_dofs()); ++d)
4464 * sum_counter_on_vertices[d] =
4465 * Utilities::MPI::sum(counter_on_vertices_mpi[d],
4466 * mpi_communicator);
4467 * sum_porous_dissipation_vertex[d] =
4468 * Utilities::MPI::sum(porous_dissipation_vertex_mpi[d],
4469 * mpi_communicator);
4470 * sum_viscous_dissipation_vertex[d] =
4471 * Utilities::MPI::sum(viscous_dissipation_vertex_mpi[d],
4472 * mpi_communicator);
4473 * sum_solid_vol_fraction_vertex[d] =
4474 * Utilities::MPI::sum(solid_vol_fraction_vertex_mpi[d],
4475 * mpi_communicator);
4477 * for (unsigned int k=0; k<num_comp_symm_tensor; ++k)
4479 * sum_cauchy_stresses_total_vertex[k][d] =
4480 * Utilities::MPI::sum(cauchy_stresses_total_vertex_mpi[k][d],
4481 * mpi_communicator);
4482 * sum_cauchy_stresses_E_vertex[k][d] =
4483 * Utilities::MPI::sum(cauchy_stresses_E_vertex_mpi[k][d],
4484 * mpi_communicator);
4486 * for (unsigned int k=0; k<dim; ++k)
4488 * sum_stretches_vertex[k][d] =
4489 * Utilities::MPI::sum(stretches_vertex_mpi[k][d],
4490 * mpi_communicator);
4494 * for (unsigned int d=0; d<(vertex_vec_handler_ref.n_dofs()); ++d)
4496 * sum_counter_on_vertices_vec[d] =
4497 * Utilities::MPI::sum(counter_on_vertices_vec_mpi[d],
4498 * mpi_communicator);
4499 * sum_seepage_velocity_vertex_vec[d] =
4500 * Utilities::MPI::sum(seepage_velocity_vertex_vec_mpi[d],
4501 * mpi_communicator);
4504 * for (unsigned int d=0; d<(vertex_handler_ref.n_dofs()); ++d)
4506 * if (sum_counter_on_vertices[d]>0)
4508 * for (unsigned int i=0; i<num_comp_symm_tensor; ++i)
4510 * sum_cauchy_stresses_total_vertex[i][d] /= sum_counter_on_vertices[d];
4511 * sum_cauchy_stresses_E_vertex[i][d] /= sum_counter_on_vertices[d];
4513 * for (unsigned int i=0; i<dim; ++i)
4515 * sum_stretches_vertex[i][d] /= sum_counter_on_vertices[d];
4517 * sum_porous_dissipation_vertex[d] /= sum_counter_on_vertices[d];
4518 * sum_viscous_dissipation_vertex[d] /= sum_counter_on_vertices[d];
4519 * sum_solid_vol_fraction_vertex[d] /= sum_counter_on_vertices[d];
4523 * for (unsigned int d=0; d<(vertex_vec_handler_ref.n_dofs()); ++d)
4525 * if (sum_counter_on_vertices_vec[d]>0)
4527 * sum_seepage_velocity_vertex_vec[d] /= sum_counter_on_vertices_vec[d];
4535 * Add the results to the solution to create the output file for Paraview
4538 * DataOut<dim> data_out;
4539 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
4541 * DataComponentInterpretation::component_is_part_of_vector);
4542 * comp_type.push_back(DataComponentInterpretation::component_is_scalar);
4544 * GridTools::get_subdomain_association(triangulation, partition_int);
4546 * std::vector<std::string> solution_name(dim, "displacement");
4547 * solution_name.push_back("pore_pressure");
4549 * data_out.attach_dof_handler(dof_handler_ref);
4550 * data_out.add_data_vector(solution_total,
4552 * DataOut<dim>::type_dof_data,
4555 * data_out.add_data_vector(solution_total,
4556 * gradient_postprocessor);
4558 * const Vector<double> partitioning(partition_int.begin(),
4559 * partition_int.end());
4561 * data_out.add_data_vector(partitioning, "partitioning");
4562 * data_out.add_data_vector(material_id, "material_id");
4566 * Integration point results -----------------------------------------------------------
4569 * if (parameters.outtype == "elements")
4571 * data_out.add_data_vector(cauchy_stresses_total_elements[0], "cauchy_xx");
4572 * data_out.add_data_vector(cauchy_stresses_total_elements[1], "cauchy_yy");
4573 * data_out.add_data_vector(cauchy_stresses_total_elements[2], "cauchy_zz");
4574 * data_out.add_data_vector(cauchy_stresses_total_elements[3], "cauchy_xy");
4575 * data_out.add_data_vector(cauchy_stresses_total_elements[4], "cauchy_xz");
4576 * data_out.add_data_vector(cauchy_stresses_total_elements[5], "cauchy_yz");
4578 * data_out.add_data_vector(cauchy_stresses_E_elements[0], "cauchy_E_xx");
4579 * data_out.add_data_vector(cauchy_stresses_E_elements[1], "cauchy_E_yy");
4580 * data_out.add_data_vector(cauchy_stresses_E_elements[2], "cauchy_E_zz");
4581 * data_out.add_data_vector(cauchy_stresses_E_elements[3], "cauchy_E_xy");
4582 * data_out.add_data_vector(cauchy_stresses_E_elements[4], "cauchy_E_xz");
4583 * data_out.add_data_vector(cauchy_stresses_E_elements[5], "cauchy_E_yz");
4585 * data_out.add_data_vector(stretches_elements[0], "stretch_xx");
4586 * data_out.add_data_vector(stretches_elements[1], "stretch_yy");
4587 * data_out.add_data_vector(stretches_elements[2], "stretch_zz");
4589 * data_out.add_data_vector(seepage_velocity_elements[0], "seepage_vel_x");
4590 * data_out.add_data_vector(seepage_velocity_elements[1], "seepage_vel_y");
4591 * data_out.add_data_vector(seepage_velocity_elements[2], "seepage_vel_z");
4593 * data_out.add_data_vector(porous_dissipation_elements, "dissipation_porous");
4594 * data_out.add_data_vector(viscous_dissipation_elements, "dissipation_viscous");
4595 * data_out.add_data_vector(solid_vol_fraction_elements, "solid_vol_fraction");
4597 * else if (parameters.outtype == "nodes")
4599 * data_out.add_data_vector(vertex_handler_ref,
4600 * sum_cauchy_stresses_total_vertex[0],
4602 * data_out.add_data_vector(vertex_handler_ref,
4603 * sum_cauchy_stresses_total_vertex[1],
4605 * data_out.add_data_vector(vertex_handler_ref,
4606 * sum_cauchy_stresses_total_vertex[2],
4608 * data_out.add_data_vector(vertex_handler_ref,
4609 * sum_cauchy_stresses_total_vertex[3],
4611 * data_out.add_data_vector(vertex_handler_ref,
4612 * sum_cauchy_stresses_total_vertex[4],
4614 * data_out.add_data_vector(vertex_handler_ref,
4615 * sum_cauchy_stresses_total_vertex[5],
4618 * data_out.add_data_vector(vertex_handler_ref,
4619 * sum_cauchy_stresses_E_vertex[0],
4621 * data_out.add_data_vector(vertex_handler_ref,
4622 * sum_cauchy_stresses_E_vertex[1],
4624 * data_out.add_data_vector(vertex_handler_ref,
4625 * sum_cauchy_stresses_E_vertex[2],
4627 * data_out.add_data_vector(vertex_handler_ref,
4628 * sum_cauchy_stresses_E_vertex[3],
4630 * data_out.add_data_vector(vertex_handler_ref,
4631 * sum_cauchy_stresses_E_vertex[4],
4633 * data_out.add_data_vector(vertex_handler_ref,
4634 * sum_cauchy_stresses_E_vertex[5],
4637 * data_out.add_data_vector(vertex_handler_ref,
4638 * sum_stretches_vertex[0],
4640 * data_out.add_data_vector(vertex_handler_ref,
4641 * sum_stretches_vertex[1],
4643 * data_out.add_data_vector(vertex_handler_ref,
4644 * sum_stretches_vertex[2],
4647 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
4648 * comp_type_vec(dim,
4649 * DataComponentInterpretation::component_is_part_of_vector);
4650 * std::vector<std::string> solution_name_vec(dim,"seepage_velocity");
4652 * data_out.add_data_vector(vertex_vec_handler_ref,
4653 * sum_seepage_velocity_vertex_vec,
4654 * solution_name_vec,
4657 * data_out.add_data_vector(vertex_handler_ref,
4658 * sum_porous_dissipation_vertex,
4659 * "dissipation_porous");
4660 * data_out.add_data_vector(vertex_handler_ref,
4661 * sum_viscous_dissipation_vertex,
4662 * "dissipation_viscous");
4663 * data_out.add_data_vector(vertex_handler_ref,
4664 * sum_solid_vol_fraction_vertex,
4665 * "solid_vol_fraction");
4669 * ---------------------------------------------------------------------
4675 * data_out.build_patches(degree_displ);
4679 * static std::string get_filename_vtu(unsigned int process,
4680 * unsigned int timestep,
4681 * const unsigned int n_digits = 5)
4683 * std::ostringstream filename_vtu;
4686 * << Utilities::int_to_string(process, n_digits)
4688 * << Utilities::int_to_string(timestep, n_digits)
4690 * return filename_vtu.str();
4693 * static std::string get_filename_pvtu(unsigned int timestep,
4694 * const unsigned int n_digits = 5)
4696 * std::ostringstream filename_vtu;
4699 * << Utilities::int_to_string(timestep, n_digits)
4701 * return filename_vtu.str();
4704 * static std::string get_filename_pvd (void)
4706 * std::ostringstream filename_vtu;
4708 * << "solution.pvd";
4709 * return filename_vtu.str();
4713 * const std::string filename_vtu = Filename::get_filename_vtu(this_mpi_process,
4715 * std::ofstream output(filename_vtu.c_str());
4716 * data_out.write_vtu(output);
4720 * We have a collection of files written in parallel
4721 * This next set of steps should only be performed by master process
4724 * if (this_mpi_process == 0)
4728 * List of all files written out at this timestep by all processors
4731 * std::vector<std::string> parallel_filenames_vtu;
4732 * for (unsigned int p=0; p<n_mpi_processes; ++p)
4734 * parallel_filenames_vtu.push_back(Filename::get_filename_vtu(p, timestep));
4737 * const std::string filename_pvtu(Filename::get_filename_pvtu(timestep));
4738 * std::ofstream pvtu_master(filename_pvtu.c_str());
4739 * data_out.write_pvtu_record(pvtu_master,
4740 * parallel_filenames_vtu);
4744 * Time dependent data master file
4747 * static std::vector<std::pair<double,std::string>> time_and_name_history;
4748 * time_and_name_history.push_back(std::make_pair(current_time,
4750 * const std::string filename_pvd(Filename::get_filename_pvd());
4751 * std::ofstream pvd_output(filename_pvd.c_str());
4752 * DataOutBase::write_pvd_record(pvd_output, time_and_name_history);
4759 * Print results to plotting file
4762 * template <int dim>
4763 * void Solid<dim>::output_results_to_plot(
4764 * const unsigned int timestep,
4765 * const double current_time,
4766 * TrilinosWrappers::MPI::BlockVector solution_IN,
4767 * std::vector<Point<dim> > &tracked_vertices_IN,
4768 * std::ofstream &plotpointfile) const
4770 * TrilinosWrappers::MPI::BlockVector solution_total(locally_owned_partitioning,
4771 * locally_relevant_partitioning,
4776 * solution_total = solution_IN;
4780 * Variables needed to print the solution file for plotting
4783 * Point<dim> reaction_force;
4784 * Point<dim> reaction_force_pressure;
4785 * Point<dim> reaction_force_extra;
4786 * double total_fluid_flow = 0.0;
4787 * double total_porous_dissipation = 0.0;
4788 * double total_viscous_dissipation = 0.0;
4789 * double total_solid_vol = 0.0;
4790 * double total_vol_current = 0.0;
4791 * double total_vol_reference = 0.0;
4792 * std::vector<Point<dim+1>> solution_vertices(tracked_vertices_IN.size());
4796 * Auxiliary variables needed for mpi processing
4799 * Tensor<1,dim> sum_reaction_mpi;
4800 * Tensor<1,dim> sum_reaction_pressure_mpi;
4801 * Tensor<1,dim> sum_reaction_extra_mpi;
4802 * sum_reaction_mpi = 0.0;
4803 * sum_reaction_pressure_mpi = 0.0;
4804 * sum_reaction_extra_mpi = 0.0;
4805 * double sum_total_flow_mpi = 0.0;
4806 * double sum_porous_dissipation_mpi = 0.0;
4807 * double sum_viscous_dissipation_mpi = 0.0;
4808 * double sum_solid_vol_mpi = 0.0;
4809 * double sum_vol_current_mpi = 0.0;
4810 * double sum_vol_reference_mpi = 0.0;
4814 * Declare an instance of the material class object
4817 * if (parameters.mat_type == "Neo-Hooke")
4818 * NeoHooke<dim,ADNumberType> material(parameters,time);
4819 * else if (parameters.mat_type == "Ogden")
4820 * Ogden<dim,ADNumberType> material(parameters, time);
4821 * else if (parameters.mat_type == "visco-Ogden")
4822 * visco_Ogden <dim,ADNumberType>material(parameters,time);
4824 * Assert (false, ExcMessage("Material type not implemented"));
4828 * Define a local instance of FEValues to compute updated values required
4829 * to calculate stresses
4832 * const UpdateFlags uf_cell(update_values | update_gradients |
4833 * update_JxW_values);
4834 * FEValues<dim> fe_values_ref (fe, qf_cell, uf_cell);
4838 * Iterate through elements (cells) and Gauss Points
4841 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
4842 * cell(IteratorFilters::LocallyOwnedCell(),
4843 * dof_handler_ref.begin_active()),
4844 * endc(IteratorFilters::LocallyOwnedCell(),
4845 * dof_handler_ref.end());
4851 * for (; cell!=endc; ++cell)
4853 * Assert(cell->is_locally_owned(), ExcInternalError());
4854 * Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
4856 * fe_values_ref.reinit(cell);
4858 * std::vector<Tensor<2,dim>> solution_grads_u(n_q_points);
4859 * fe_values_ref[u_fe].get_function_gradients(solution_total,
4860 * solution_grads_u);
4862 * std::vector<double> solution_values_p_fluid_total(n_q_points);
4863 * fe_values_ref[p_fluid_fe].get_function_values(solution_total,
4864 * solution_values_p_fluid_total);
4866 * std::vector<Tensor<1,dim >> solution_grads_p_fluid_AD(n_q_points);
4867 * fe_values_ref[p_fluid_fe].get_function_gradients(solution_total,
4868 * solution_grads_p_fluid_AD);
4872 * start gauss point loop
4875 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4877 * const Tensor<2,dim,ADNumberType>
4878 * F_AD = Physics::Elasticity::Kinematics::F(solution_grads_u[q_point]);
4879 * ADNumberType det_F_AD = determinant(F_AD);
4880 * const double det_F = Tensor<0,dim,double>(det_F_AD);
4882 * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4883 * lqph = quadrature_point_history.get_data(cell);
4884 * Assert(lqph.size() == n_q_points, ExcInternalError());
4886 * double JxW = fe_values_ref.JxW(q_point);
4893 * sum_vol_current_mpi += det_F * JxW;
4894 * sum_vol_reference_mpi += JxW;
4895 * sum_solid_vol_mpi += parameters.solid_vol_frac * JxW * det_F;
4902 * const Tensor<2,dim,ADNumberType> F_inv = invert(F_AD);
4903 * const Tensor<1,dim,ADNumberType>
4904 * grad_p_fluid_AD = solution_grads_p_fluid_AD[q_point]*F_inv;
4905 * const Tensor<1,dim,ADNumberType> seepage_vel_AD
4906 * = lqph[q_point]->get_seepage_velocity_current(F_AD, grad_p_fluid_AD);
4913 * const double porous_dissipation =
4914 * lqph[q_point]->get_porous_dissipation(F_AD, grad_p_fluid_AD);
4915 * sum_porous_dissipation_mpi += porous_dissipation * det_F * JxW;
4917 * const double viscous_dissipation = lqph[q_point]->get_viscous_dissipation();
4918 * sum_viscous_dissipation_mpi += viscous_dissipation * det_F * JxW;
4922 * ---------------------------------------------------------------
4925 * } //end gauss point loop
4929 * Compute reaction force on load boundary & total fluid flow across
4931 * Define a local instance of FEFaceValues to compute values required
4932 * to calculate reaction force
4935 * const UpdateFlags uf_face( update_values | update_gradients |
4936 * update_normal_vectors | update_JxW_values );
4937 * FEFaceValues<dim> fe_face_values_ref(fe, qf_face, uf_face);
4944 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
4951 * if (cell->face(face)->at_boundary() == true &&
4952 * cell->face(face)->boundary_id() == get_reaction_boundary_id_for_output() )
4954 * fe_face_values_ref.reinit(cell, face);
4958 * Get displacement gradients for current face
4961 * std::vector<Tensor<2,dim> > solution_grads_u_f(n_q_points_f);
4962 * fe_face_values_ref[u_fe].get_function_gradients
4964 * solution_grads_u_f);
4968 * Get pressure for current element
4971 * std::vector< double > solution_values_p_fluid_total_f(n_q_points_f);
4972 * fe_face_values_ref[p_fluid_fe].get_function_values
4974 * solution_values_p_fluid_total_f);
4978 * start gauss points on faces loop
4981 * for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point)
4983 * const Tensor<1,dim> &N = fe_face_values_ref.normal_vector(f_q_point);
4984 * const double JxW_f = fe_face_values_ref.JxW(f_q_point);
4988 * Compute deformation gradient from displacements gradient
4989 * (present configuration)
4992 * const Tensor<2,dim,ADNumberType> F_AD =
4993 * Physics::Elasticity::Kinematics::F(solution_grads_u_f[f_q_point]);
4995 * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4996 * lqph = quadrature_point_history.get_data(cell);
4997 * Assert(lqph.size() == n_q_points, ExcInternalError());
4999 * const double p_fluid = solution_values_p_fluid_total[f_q_point];
5006 * static const SymmetricTensor<2,dim,double>
5007 * I (Physics::Elasticity::StandardTensors<dim>::I);
5008 * SymmetricTensor<2,dim> sigma_E;
5009 * const SymmetricTensor<2,dim,ADNumberType> sigma_E_AD =
5010 * lqph[f_q_point]->get_Cauchy_E(F_AD);
5012 * for (unsigned int i=0; i<dim; ++i)
5013 * for (unsigned int j=0; j<dim; ++j)
5014 * sigma_E[i][j] = Tensor<0,dim,double>(sigma_E_AD[i][j]);
5016 * SymmetricTensor<2,dim> sigma_fluid_vol(I);
5017 * sigma_fluid_vol *= -1.0*p_fluid;
5018 * const SymmetricTensor<2,dim> sigma = sigma_E+sigma_fluid_vol;
5019 * sum_reaction_mpi += sigma * N * JxW_f;
5020 * sum_reaction_pressure_mpi += sigma_fluid_vol * N * JxW_f;
5021 * sum_reaction_extra_mpi += sigma_E * N * JxW_f;
5022 * }//end gauss points on faces loop
5030 * if (cell->face(face)->at_boundary() == true &&
5031 * (cell->face(face)->boundary_id() ==
5032 * get_drained_boundary_id_for_output().first ||
5033 * cell->face(face)->boundary_id() ==
5034 * get_drained_boundary_id_for_output().second ) )
5036 * fe_face_values_ref.reinit(cell, face);
5040 * Get displacement gradients for current face
5043 * std::vector<Tensor<2,dim>> solution_grads_u_f(n_q_points_f);
5044 * fe_face_values_ref[u_fe].get_function_gradients
5046 * solution_grads_u_f);
5050 * Get pressure gradients for current face
5053 * std::vector<Tensor<1,dim>> solution_grads_p_f(n_q_points_f);
5054 * fe_face_values_ref[p_fluid_fe].get_function_gradients
5056 * solution_grads_p_f);
5060 * start gauss points on faces loop
5063 * for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point)
5065 * const Tensor<1,dim> &N =
5066 * fe_face_values_ref.normal_vector(f_q_point);
5067 * const double JxW_f = fe_face_values_ref.JxW(f_q_point);
5071 * Deformation gradient and inverse from displacements gradient
5072 * (present configuration)
5075 * const Tensor<2,dim,ADNumberType> F_AD
5076 * = Physics::Elasticity::Kinematics::F(solution_grads_u_f[f_q_point]);
5078 * const Tensor<2,dim,ADNumberType> F_inv_AD = invert(F_AD);
5079 * ADNumberType det_F_AD = determinant(F_AD);
5081 * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
5082 * lqph = quadrature_point_history.get_data(cell);
5083 * Assert(lqph.size() == n_q_points, ExcInternalError());
5090 * Tensor<1,dim> seepage;
5091 * double det_F = Tensor<0,dim,double>(det_F_AD);
5092 * const Tensor<1,dim,ADNumberType> grad_p
5093 * = solution_grads_p_f[f_q_point]*F_inv_AD;
5094 * const Tensor<1,dim,ADNumberType> seepage_AD
5095 * = lqph[f_q_point]->get_seepage_velocity_current(F_AD, grad_p);
5097 * for (unsigned int i=0; i<dim; ++i)
5098 * seepage[i] = Tensor<0,dim,double>(seepage_AD[i]);
5100 * sum_total_flow_mpi += (seepage/det_F) * N * JxW_f;
5101 * }//end gauss points on faces loop
5108 * Sum the results from different MPI process and then add to the reaction_force vector
5109 * In theory, the solution on each surface (each cell) only exists in one MPI process
5110 * so, we add all MPI process, one will have the solution and the others will be zero
5113 * for (unsigned int d=0; d<dim; ++d)
5115 * reaction_force[d] = Utilities::MPI::sum(sum_reaction_mpi[d],
5116 * mpi_communicator);
5117 * reaction_force_pressure[d] = Utilities::MPI::sum(sum_reaction_pressure_mpi[d],
5118 * mpi_communicator);
5119 * reaction_force_extra[d] = Utilities::MPI::sum(sum_reaction_extra_mpi[d],
5120 * mpi_communicator);
5125 * Same for total fluid flow, and for porous and viscous dissipations
5128 * total_fluid_flow = Utilities::MPI::sum(sum_total_flow_mpi,
5129 * mpi_communicator);
5130 * total_porous_dissipation = Utilities::MPI::sum(sum_porous_dissipation_mpi,
5131 * mpi_communicator);
5132 * total_viscous_dissipation = Utilities::MPI::sum(sum_viscous_dissipation_mpi,
5133 * mpi_communicator);
5134 * total_solid_vol = Utilities::MPI::sum(sum_solid_vol_mpi,
5135 * mpi_communicator);
5136 * total_vol_current = Utilities::MPI::sum(sum_vol_current_mpi,
5137 * mpi_communicator);
5138 * total_vol_reference = Utilities::MPI::sum(sum_vol_reference_mpi,
5139 * mpi_communicator);
5143 * Extract solution for tracked vectors
5144 * Copying an MPI::BlockVector into MPI::Vector is not possible,
5145 * so we copy each block of MPI::BlockVector into an MPI::Vector
5146 * And then we copy the MPI::Vector into "normal" Vectors
5149 * TrilinosWrappers::MPI::Vector solution_vector_u_MPI(solution_total.block(u_block));
5150 * TrilinosWrappers::MPI::Vector solution_vector_p_MPI(solution_total.block(p_fluid_block));
5151 * Vector<double> solution_u_vector(solution_vector_u_MPI);
5152 * Vector<double> solution_p_vector(solution_vector_p_MPI);
5154 * if (this_mpi_process == 0)
5158 * Append the pressure solution vector to the displacement solution vector,
5159 * creating a single solution vector equivalent to the original BlockVector
5160 * so FEFieldFunction will work with the dof_handler_ref.
5163 * Vector<double> solution_vector(solution_p_vector.size()
5164 * +solution_u_vector.size());
5166 * for (unsigned int d=0; d<(solution_u_vector.size()); ++d)
5167 * solution_vector[d] = solution_u_vector[d];
5169 * for (unsigned int d=0; d<(solution_p_vector.size()); ++d)
5170 * solution_vector[solution_u_vector.size()+d] = solution_p_vector[d];
5172 * Functions::FEFieldFunction<dim,Vector<double>>
5173 * find_solution(dof_handler_ref, solution_vector);
5175 * for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
5177 * Vector<double> update(dim+1);
5178 * Point<dim> pt_ref;
5180 * pt_ref[0]= tracked_vertices_IN[p][0];
5181 * pt_ref[1]= tracked_vertices_IN[p][1];
5182 * pt_ref[2]= tracked_vertices_IN[p][2];
5184 * find_solution.vector_value(pt_ref, update);
5186 * for (unsigned int d=0; d<(dim+1); ++d)
5190 * For values close to zero, set to 0.0
5193 * if (abs(update[d])<1.5*parameters.tol_u)
5195 * solution_vertices[p][d] = update[d];
5200 * Write the results to the plotting file.
5201 * Add two blank lines between cycles in the cyclic loading examples so GNUPLOT can detect each cycle as a different block
5204 * if (( (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")||
5205 * (parameters.geom_type == "Budday_cube_tension_compression")||
5206 * (parameters.geom_type == "Budday_cube_shear_fully_fixed") ) &&
5207 * ( (abs(current_time - parameters.end_time/3.) <0.9*parameters.delta_t)||
5208 * (abs(current_time - 2.*parameters.end_time/3.)<0.9*parameters.delta_t) ) &&
5209 * parameters.num_cycle_sets == 1 )
5211 * plotpointfile << std::endl<< std::endl;
5213 * if (( (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")||
5214 * (parameters.geom_type == "Budday_cube_tension_compression")||
5215 * (parameters.geom_type == "Budday_cube_shear_fully_fixed") ) &&
5216 * ( (abs(current_time - parameters.end_time/9.) <0.9*parameters.delta_t)||
5217 * (abs(current_time - 2.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5218 * (abs(current_time - 3.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5219 * (abs(current_time - 5.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5220 * (abs(current_time - 7.*parameters.end_time/9.)<0.9*parameters.delta_t) ) &&
5221 * parameters.num_cycle_sets == 2 )
5223 * plotpointfile << std::endl<< std::endl;
5226 * plotpointfile << std::setprecision(6) << std::scientific;
5227 * plotpointfile << std::setw(16) << current_time << ","
5228 * << std::setw(15) << total_vol_reference << ","
5229 * << std::setw(15) << total_vol_current << ","
5230 * << std::setw(15) << total_solid_vol << ",";
5232 * if (current_time == 0.0)
5234 * for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
5236 * for (unsigned int d=0; d<dim; ++d)
5237 * plotpointfile << std::setw(15) << 0.0 << ",";
5239 * plotpointfile << std::setw(15) << parameters.drained_pressure << ",";
5241 * for (unsigned int d=0; d<(3*dim+2); ++d)
5242 * plotpointfile << std::setw(15) << 0.0 << ",";
5244 * plotpointfile << std::setw(15) << 0.0;
5248 * for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
5249 * for (unsigned int d=0; d<(dim+1); ++d)
5250 * plotpointfile << std::setw(15) << solution_vertices[p][d]<< ",";
5252 * for (unsigned int d=0; d<dim; ++d)
5253 * plotpointfile << std::setw(15) << reaction_force[d] << ",";
5255 * for (unsigned int d=0; d<dim; ++d)
5256 * plotpointfile << std::setw(15) << reaction_force_pressure[d] << ",";
5258 * for (unsigned int d=0; d<dim; ++d)
5259 * plotpointfile << std::setw(15) << reaction_force_extra[d] << ",";
5261 * plotpointfile << std::setw(15) << total_fluid_flow << ","
5262 * << std::setw(15) << total_porous_dissipation<< ","
5263 * << std::setw(15) << total_viscous_dissipation;
5265 * plotpointfile << std::endl;
5271 * Header for console output file
5274 * template <int dim>
5275 * void Solid<dim>::print_console_file_header(std::ofstream &outputfile) const
5277 * outputfile << "/*-----------------------------------------------------------------------------------------";
5278 * outputfile << "\n\n Poro-viscoelastic formulation to solve nonlinear solid mechanics problems using deal.ii";
5279 * outputfile << "\n\n Problem setup by E Comellas and J-P Pelteret, University of Erlangen-Nuremberg, 2018";
5280 * outputfile << "\n\n/*-----------------------------------------------------------------------------------------";
5281 * outputfile << "\n\nCONSOLE OUTPUT: \n\n";
5286 * Header for plotting output file
5289 * template <int dim>
5290 * void Solid<dim>::print_plot_file_header(std::vector<Point<dim> > &tracked_vertices,
5291 * std::ofstream &plotpointfile) const
5293 * plotpointfile << "#\n# *** Solution history for tracked vertices -- DOF: 0 = Ux, 1 = Uy, 2 = Uz, 3 = P ***"
5296 * for (unsigned int p=0; p<tracked_vertices.size(); ++p)
5298 * plotpointfile << "# Point " << p << " coordinates: ";
5299 * for (unsigned int d=0; d<dim; ++d)
5301 * plotpointfile << tracked_vertices[p][d];
5302 * if (!( (p == tracked_vertices.size()-1) && (d == dim-1) ))
5303 * plotpointfile << ", ";
5305 * plotpointfile << std::endl;
5307 * plotpointfile << "# The reaction force is the integral over the loaded surfaces in the "
5308 * << "undeformed configuration of the Cauchy stress times the normal surface unit vector.\n"
5309 * << "# reac(p) corresponds to the volumetric part of the Cauchy stress due to the pore fluid pressure"
5310 * << " and reac(E) corresponds to the extra part of the Cauchy stress due to the solid contribution."
5312 * << "# The fluid flow is the integral over the drained surfaces in the "
5313 * << "undeformed configuration of the seepage velocity times the normal surface unit vector."
5315 * << "# Column number:"
5319 * unsigned int columns = 24;
5320 * for (unsigned int d=1; d<columns; ++d)
5321 * plotpointfile << std::setw(15)<< d <<",";
5323 * plotpointfile << std::setw(15)<< columns
5326 * << std::right << std::setw(16) << "Time,"
5327 * << std::right << std::setw(16) << "ref vol,"
5328 * << std::right << std::setw(16) << "def vol,"
5329 * << std::right << std::setw(16) << "solid vol,";
5330 * for (unsigned int p=0; p<tracked_vertices.size(); ++p)
5331 * for (unsigned int d=0; d<(dim+1); ++d)
5332 * plotpointfile << std::right<< std::setw(11)
5333 * <<"P" << p << "[" << d << "],";
5335 * for (unsigned int d=0; d<dim; ++d)
5336 * plotpointfile << std::right<< std::setw(13)
5337 * << "reaction [" << d << "],";
5339 * for (unsigned int d=0; d<dim; ++d)
5340 * plotpointfile << std::right<< std::setw(13)
5341 * << "reac(p) [" << d << "],";
5343 * for (unsigned int d=0; d<dim; ++d)
5344 * plotpointfile << std::right<< std::setw(13)
5345 * << "reac(E) [" << d << "],";
5347 * plotpointfile << std::right<< std::setw(16)<< "fluid flow,"
5348 * << std::right<< std::setw(16)<< "porous dissip,"
5349 * << std::right<< std::setw(15)<< "viscous dissip"
5355 * Footer for console output file
5358 * template <int dim>
5359 * void Solid<dim>::print_console_file_footer(std::ofstream &outputfile) const
5363 * Copy "parameters" file at end of output file.
5366 * std::ifstream infile("parameters.prm");
5367 * std::string content = "";
5370 * for(i=0 ; infile.eof()!=true ; i++)
5372 * char aux = infile.get();
5374 * if(aux=='\n
') content += '#
';
5378 * content.erase(content.end()-1);
5381 * outputfile << "\n\n\n\n PARAMETERS FILE USED IN THIS COMPUTATION: \n#"
5388 * Footer for plotting output file
5391 * template <int dim>
5392 * void Solid<dim>::print_plot_file_footer(std::ofstream &plotpointfile) const
5396 * Copy "parameters" file at end of output file.
5399 * std::ifstream infile("parameters.prm");
5400 * std::string content = "";
5403 * for(i=0 ; infile.eof()!=true ; i++)
5405 * char aux = infile.get();
5407 * if(aux=='\n
') content += '#
';
5411 * content.erase(content.end()-1);
5414 * plotpointfile << "#"<< std::endl
5415 * << "#"<< std::endl
5416 * << "# PARAMETERS FILE USED IN THIS COMPUTATION:" << std::endl
5417 * << "#"<< std::endl
5425 * <a name="nonlinear-poro-viscoelasticity.cc-VerificationexamplesfromEhlersandEipper1999"></a>
5426 * <h3>Verification examples from Ehlers and Eipper 1999</h3>
5427 * We group the definition of the geometry, boundary and loading conditions specific to
5428 * the verification examples from Ehlers and Eipper 1999 into specific classes.
5433 * <a name="nonlinear-poro-viscoelasticity.cc-BaseclassTubegeometryandboundaryconditions"></a>
5434 * <h4>Base class: Tube geometry and boundary conditions</h4>
5437 * template <int dim>
5438 * class VerificationEhlers1999TubeBase
5439 * : public Solid<dim>
5442 * VerificationEhlers1999TubeBase (const Parameters::AllParameters ¶meters)
5443 * : Solid<dim> (parameters)
5446 * virtual ~VerificationEhlers1999TubeBase () {}
5449 * virtual void make_grid() override
5451 * GridGenerator::cylinder( this->triangulation,
5455 * const double rot_angle = 3.0*numbers::PI/2.0;
5456 * GridTools::rotate( Point<3>::unit_vector(1), rot_angle, this->triangulation);
5458 * this->triangulation.reset_manifold(0);
5459 * static const CylindricalManifold<dim> manifold_description_3d(2);
5460 * this->triangulation.set_manifold (0, manifold_description_3d);
5461 * GridTools::scale(this->parameters.scale, this->triangulation);
5462 * this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5463 * this->triangulation.reset_manifold(0);
5466 * virtual void define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
5468 * tracked_vertices[0][0] = 0.0*this->parameters.scale;
5469 * tracked_vertices[0][1] = 0.0*this->parameters.scale;
5470 * tracked_vertices[0][2] = 0.5*this->parameters.scale;
5472 * tracked_vertices[1][0] = 0.0*this->parameters.scale;
5473 * tracked_vertices[1][1] = 0.0*this->parameters.scale;
5474 * tracked_vertices[1][2] = -0.5*this->parameters.scale;
5477 * virtual void make_dirichlet_constraints(AffineConstraints<double> &constraints) override
5479 * if (this->time.get_timestep() < 2)
5481 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5483 * Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5485 * (this->fe.component_mask(this->pressure)));
5489 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5491 * Functions::ZeroFunction<dim>(this->n_components),
5493 * (this->fe.component_mask(this->pressure)));
5496 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5498 * Functions::ZeroFunction<dim>(this->n_components),
5500 * (this->fe.component_mask(this->x_displacement)|
5501 * this->fe.component_mask(this->y_displacement) ) );
5503 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5505 * Functions::ZeroFunction<dim>(this->n_components),
5507 * (this->fe.component_mask(this->x_displacement) |
5508 * this->fe.component_mask(this->y_displacement) |
5509 * this->fe.component_mask(this->z_displacement) ));
5513 * get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5514 * const Point<dim> &pt) const override
5517 * (void)boundary_id;
5521 * virtual types::boundary_id
5522 * get_reaction_boundary_id_for_output() const override
5527 * virtual std::pair<types::boundary_id,types::boundary_id>
5528 * get_drained_boundary_id_for_output() const override
5530 * return std::make_pair(2,2);
5533 * virtual std::vector<double>
5534 * get_dirichlet_load(const types::boundary_id &boundary_id,
5535 * const int &direction) const override
5537 * std::vector<double> displ_incr(dim, 0.0);
5538 * (void)boundary_id;
5540 * AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));
5542 * return displ_incr;
5549 * <a name="nonlinear-poro-viscoelasticity.cc-DerivedclassSteploadexample"></a>
5550 * <h4>Derived class: Step load example</h4>
5553 * template <int dim>
5554 * class VerificationEhlers1999StepLoad
5555 * : public VerificationEhlers1999TubeBase<dim>
5558 * VerificationEhlers1999StepLoad (const Parameters::AllParameters ¶meters)
5559 * : VerificationEhlers1999TubeBase<dim> (parameters)
5562 * virtual ~VerificationEhlers1999StepLoad () {}
5565 * virtual Tensor<1,dim>
5566 * get_neumann_traction (const types::boundary_id &boundary_id,
5567 * const Point<dim> &pt,
5568 * const Tensor<1,dim> &N) const override
5570 * if (this->parameters.load_type == "pressure")
5572 * if (boundary_id == 2)
5574 * return this->parameters.load * N;
5580 * return Tensor<1,dim>();
5587 * <a name="nonlinear-poro-viscoelasticity.cc-DerivedclassLoadincreasingexample"></a>
5588 * <h4>Derived class: Load increasing example</h4>
5591 * template <int dim>
5592 * class VerificationEhlers1999IncreaseLoad
5593 * : public VerificationEhlers1999TubeBase<dim>
5596 * VerificationEhlers1999IncreaseLoad (const Parameters::AllParameters ¶meters)
5597 * : VerificationEhlers1999TubeBase<dim> (parameters)
5600 * virtual ~VerificationEhlers1999IncreaseLoad () {}
5603 * virtual Tensor<1,dim>
5604 * get_neumann_traction (const types::boundary_id &boundary_id,
5605 * const Point<dim> &pt,
5606 * const Tensor<1,dim> &N) const override
5608 * if (this->parameters.load_type == "pressure")
5610 * if (boundary_id == 2)
5612 * const double initial_load = this->parameters.load;
5613 * const double final_load = 20.0*initial_load;
5614 * const double initial_time = this->time.get_delta_t();
5615 * const double final_time = this->time.get_end();
5616 * const double current_time = this->time.get_current();
5617 * const double load = initial_load + (final_load-initial_load)*(current_time-initial_time)/(final_time-initial_time);
5624 * return Tensor<1,dim>();
5631 * <a name="nonlinear-poro-viscoelasticity.cc-ClassConsolidationcube"></a>
5632 * <h4>Class: Consolidation cube</h4>
5635 * template <int dim>
5636 * class VerificationEhlers1999CubeConsolidation
5637 * : public Solid<dim>
5640 * VerificationEhlers1999CubeConsolidation (const Parameters::AllParameters ¶meters)
5641 * : Solid<dim> (parameters)
5644 * virtual ~VerificationEhlers1999CubeConsolidation () {}
5648 * make_grid() override
5650 * GridGenerator::hyper_rectangle(this->triangulation,
5651 * Point<dim>(0.0, 0.0, 0.0),
5652 * Point<dim>(1.0, 1.0, 1.0),
5655 * GridTools::scale(this->parameters.scale, this->triangulation);
5656 * this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5658 * typename Triangulation<dim>::active_cell_iterator cell =
5659 * this->triangulation.begin_active(), endc = this->triangulation.end();
5660 * for (; cell != endc; ++cell)
5662 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5663 * if (cell->face(face)->at_boundary() == true &&
5664 * cell->face(face)->center()[2] == 1.0 * this->parameters.scale)
5666 * if (cell->face(face)->center()[0] < 0.5 * this->parameters.scale &&
5667 * cell->face(face)->center()[1] < 0.5 * this->parameters.scale)
5668 * cell->face(face)->set_boundary_id(100);
5670 * cell->face(face)->set_boundary_id(101);
5676 * define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
5678 * tracked_vertices[0][0] = 0.0*this->parameters.scale;
5679 * tracked_vertices[0][1] = 0.0*this->parameters.scale;
5680 * tracked_vertices[0][2] = 1.0*this->parameters.scale;
5682 * tracked_vertices[1][0] = 0.0*this->parameters.scale;
5683 * tracked_vertices[1][1] = 0.0*this->parameters.scale;
5684 * tracked_vertices[1][2] = 0.0*this->parameters.scale;
5688 * make_dirichlet_constraints(AffineConstraints<double> &constraints) override
5690 * if (this->time.get_timestep() < 2)
5692 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5694 * Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5696 * (this->fe.component_mask(this->pressure)));
5700 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5702 * Functions::ZeroFunction<dim>(this->n_components),
5704 * (this->fe.component_mask(this->pressure)));
5707 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5709 * Functions::ZeroFunction<dim>(this->n_components),
5711 * this->fe.component_mask(this->x_displacement));
5713 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5715 * Functions::ZeroFunction<dim>(this->n_components),
5717 * this->fe.component_mask(this->x_displacement));
5719 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5721 * Functions::ZeroFunction<dim>(this->n_components),
5723 * this->fe.component_mask(this->y_displacement));
5725 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5727 * Functions::ZeroFunction<dim>(this->n_components),
5729 * this->fe.component_mask(this->y_displacement));
5731 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5733 * Functions::ZeroFunction<dim>(this->n_components),
5735 * ( this->fe.component_mask(this->x_displacement) |
5736 * this->fe.component_mask(this->y_displacement) |
5737 * this->fe.component_mask(this->z_displacement) ));
5740 * virtual Tensor<1,dim>
5741 * get_neumann_traction (const types::boundary_id &boundary_id,
5742 * const Point<dim> &pt,
5743 * const Tensor<1,dim> &N) const override
5745 * if (this->parameters.load_type == "pressure")
5747 * if (boundary_id == 100)
5749 * return this->parameters.load * N;
5755 * return Tensor<1,dim>();
5759 * get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5760 * const Point<dim> &pt) const override
5763 * (void)boundary_id;
5767 * virtual types::boundary_id
5768 * get_reaction_boundary_id_for_output() const override
5773 * virtual std::pair<types::boundary_id,types::boundary_id>
5774 * get_drained_boundary_id_for_output() const override
5776 * return std::make_pair(101,101);
5779 * virtual std::vector<double>
5780 * get_dirichlet_load(const types::boundary_id &boundary_id,
5781 * const int &direction) const override
5783 * std::vector<double> displ_incr(dim, 0.0);
5784 * (void)boundary_id;
5786 * AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));
5788 * return displ_incr;
5795 * <a name="nonlinear-poro-viscoelasticity.cc-Franceschiniexperiments"></a>
5796 * <h4>Franceschini experiments</h4>
5799 * template <int dim>
5800 * class Franceschini2006Consolidation
5801 * : public Solid<dim>
5804 * Franceschini2006Consolidation (const Parameters::AllParameters ¶meters)
5805 * : Solid<dim> (parameters)
5808 * virtual ~Franceschini2006Consolidation () {}
5811 * virtual void make_grid() override
5813 * const Point<dim-1> mesh_center(0.0, 0.0);
5814 * const double radius = 0.5;
5817 * const double height = 0.27; //8.1 mm for 30 mm radius
5820 * const double height = 0.23; //6.9 mm for 30 mm radius
5821 * Triangulation<dim-1> triangulation_in;
5822 * GridGenerator::hyper_ball( triangulation_in,
5826 * GridGenerator::extrude_triangulation(triangulation_in,
5829 * this->triangulation);
5831 * const CylindricalManifold<dim> cylinder_3d(2);
5832 * const types::manifold_id cylinder_id = 0;
5835 * this->triangulation.set_manifold(cylinder_id, cylinder_3d);
5837 * for (auto cell : this->triangulation.active_cell_iterators())
5839 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5841 * if (cell->face(face)->at_boundary() == true)
5843 * if (cell->face(face)->center()[2] == 0.0)
5844 * cell->face(face)->set_boundary_id(1);
5846 * else if (cell->face(face)->center()[2] == height)
5847 * cell->face(face)->set_boundary_id(2);
5851 * cell->face(face)->set_boundary_id(0);
5852 * cell->face(face)->set_all_manifold_ids(cylinder_id);
5858 * GridTools::scale(this->parameters.scale, this->triangulation);
5859 * this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5862 * virtual void define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
5864 * tracked_vertices[0][0] = 0.0*this->parameters.scale;
5865 * tracked_vertices[0][1] = 0.0*this->parameters.scale;
5868 * tracked_vertices[0][2] = 0.27*this->parameters.scale;
5871 * tracked_vertices[0][2] = 0.23*this->parameters.scale;
5873 * tracked_vertices[1][0] = 0.0*this->parameters.scale;
5874 * tracked_vertices[1][1] = 0.0*this->parameters.scale;
5875 * tracked_vertices[1][2] = 0.0*this->parameters.scale;
5878 * virtual void make_dirichlet_constraints(AffineConstraints<double> &constraints) override
5880 * if (this->time.get_timestep() < 2)
5882 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5884 * Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5886 * (this->fe.component_mask(this->pressure)));
5888 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5890 * Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5892 * (this->fe.component_mask(this->pressure)));
5896 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5898 * Functions::ZeroFunction<dim>(this->n_components),
5900 * (this->fe.component_mask(this->pressure)));
5902 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5904 * Functions::ZeroFunction<dim>(this->n_components),
5906 * (this->fe.component_mask(this->pressure)));
5909 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5911 * Functions::ZeroFunction<dim>(this->n_components),
5913 * (this->fe.component_mask(this->x_displacement)|
5914 * this->fe.component_mask(this->y_displacement) ) );
5916 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5918 * Functions::ZeroFunction<dim>(this->n_components),
5920 * (this->fe.component_mask(this->x_displacement) |
5921 * this->fe.component_mask(this->y_displacement) |
5922 * this->fe.component_mask(this->z_displacement) ));
5924 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5926 * Functions::ZeroFunction<dim>(this->n_components),
5928 * (this->fe.component_mask(this->x_displacement) |
5929 * this->fe.component_mask(this->y_displacement) ));
5933 * get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5934 * const Point<dim> &pt) const override
5937 * (void)boundary_id;
5941 * virtual types::boundary_id
5942 * get_reaction_boundary_id_for_output() const override
5947 * virtual std::pair<types::boundary_id,types::boundary_id>
5948 * get_drained_boundary_id_for_output() const override
5950 * return std::make_pair(1,2);
5953 * virtual std::vector<double>
5954 * get_dirichlet_load(const types::boundary_id &boundary_id,
5955 * const int &direction) const override
5957 * std::vector<double> displ_incr(dim, 0.0);
5958 * (void)boundary_id;
5960 * AssertThrow(false, ExcMessage("Displacement loading not implemented for Franceschini examples."));
5962 * return displ_incr;
5965 * virtual Tensor<1,dim>
5966 * get_neumann_traction (const types::boundary_id &boundary_id,
5967 * const Point<dim> &pt,
5968 * const Tensor<1,dim> &N) const override
5970 * if (this->parameters.load_type == "pressure")
5972 * if (boundary_id == 2)
5974 * return (this->parameters.load * N);
5976 * const double final_load = this->parameters.load;
5977 * const double final_load_time = 10 * this->time.get_delta_t();
5978 * const double current_time = this->time.get_current();
5981 * const double c = final_load_time / 2.0;
5982 * const double r = 200.0 * 0.03 / c;
5984 * const double load = final_load * std::exp(r * current_time)
5985 * / ( std::exp(c * current_time) + std::exp(r * current_time));
5993 * return Tensor<1,dim>();
6000 * <a name="nonlinear-poro-viscoelasticity.cc-ExamplestoreproduceexperimentsbyBuddayetal2017"></a>
6001 * <h3>Examples to reproduce experiments by Budday et al. 2017</h3>
6002 * We group the definition of the geometry, boundary and loading conditions specific to
6003 * the examples to reproduce experiments by Budday et al. 2017 into specific classes.
6008 * <a name="nonlinear-poro-viscoelasticity.cc-BaseclassCubegeometryandloadingpattern"></a>
6009 * <h4>Base class: Cube geometry and loading pattern</h4>
6012 * template <int dim>
6013 * class BrainBudday2017BaseCube
6014 * : public Solid<dim>
6017 * BrainBudday2017BaseCube (const Parameters::AllParameters ¶meters)
6018 * : Solid<dim> (parameters)
6021 * virtual ~BrainBudday2017BaseCube () {}
6025 * make_grid() override
6027 * GridGenerator::hyper_cube(this->triangulation,
6032 * typename Triangulation<dim>::active_cell_iterator cell =
6033 * this->triangulation.begin_active(), endc = this->triangulation.end();
6034 * for (; cell != endc; ++cell)
6036 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
6037 * if (cell->face(face)->at_boundary() == true &&
6038 * ( cell->face(face)->boundary_id() == 0 ||
6039 * cell->face(face)->boundary_id() == 1 ||
6040 * cell->face(face)->boundary_id() == 2 ||
6041 * cell->face(face)->boundary_id() == 3 ) )
6043 * cell->face(face)->set_boundary_id(100);
6047 * GridTools::scale(this->parameters.scale, this->triangulation);
6048 * this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
6052 * get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
6053 * const Point<dim> &pt) const override
6056 * (void)boundary_id;
6060 * virtual std::pair<types::boundary_id,types::boundary_id>
6061 * get_drained_boundary_id_for_output() const override
6063 * return std::make_pair(100,100);
6070 * <a name="nonlinear-poro-viscoelasticity.cc-DerivedclassUniaxialboundaryconditions"></a>
6071 * <h4>Derived class: Uniaxial boundary conditions</h4>
6074 * template <int dim>
6075 * class BrainBudday2017CubeTensionCompression
6076 * : public BrainBudday2017BaseCube<dim>
6079 * BrainBudday2017CubeTensionCompression (const Parameters::AllParameters ¶meters)
6080 * : BrainBudday2017BaseCube<dim> (parameters)
6083 * virtual ~BrainBudday2017CubeTensionCompression () {}
6087 * define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
6089 * tracked_vertices[0][0] = 0.5*this->parameters.scale;
6090 * tracked_vertices[0][1] = 0.5*this->parameters.scale;
6091 * tracked_vertices[0][2] = 1.0*this->parameters.scale;
6093 * tracked_vertices[1][0] = 0.5*this->parameters.scale;
6094 * tracked_vertices[1][1] = 0.5*this->parameters.scale;
6095 * tracked_vertices[1][2] = 0.5*this->parameters.scale;
6099 * make_dirichlet_constraints(AffineConstraints<double> &constraints) override
6101 * if (this->time.get_timestep() < 2)
6103 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
6105 * Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
6107 * (this->fe.component_mask(this->pressure)));
6111 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6113 * Functions::ZeroFunction<dim>(this->n_components),
6115 * (this->fe.component_mask(this->pressure)));
6117 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6119 * Functions::ZeroFunction<dim>(this->n_components),
6121 * this->fe.component_mask(this->z_displacement) );
6123 * Point<dim> fix_node(0.5*this->parameters.scale, 0.5*this->parameters.scale, 0.0);
6124 * typename DoFHandler<dim>::active_cell_iterator
6125 * cell = this->dof_handler_ref.begin_active(), endc = this->dof_handler_ref.end();
6126 * for (; cell != endc; ++cell)
6127 * for (unsigned int node = 0; node < GeometryInfo<dim>::vertices_per_cell; ++node)
6129 * if ( (abs(cell->vertex(node)[2]-fix_node[2]) < (1e-6 * this->parameters.scale))
6130 * && (abs(cell->vertex(node)[0]-fix_node[0]) < (1e-6 * this->parameters.scale)))
6131 * constraints.add_line(cell->vertex_dof_index(node, 0));
6133 * if ( (abs(cell->vertex(node)[2]-fix_node[2]) < (1e-6 * this->parameters.scale))
6134 * && (abs(cell->vertex(node)[1]-fix_node[1]) < (1e-6 * this->parameters.scale)))
6135 * constraints.add_line(cell->vertex_dof_index(node, 1));
6138 * if (this->parameters.load_type == "displacement")
6140 * const std::vector<double> value = get_dirichlet_load(5,2);
6141 * const FEValuesExtractors::Scalar direction(this->z_displacement);
6143 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6145 * Functions::ConstantFunction<dim>(value[2],this->n_components),
6147 * this->fe.component_mask(direction));
6151 * virtual Tensor<1,dim>
6152 * get_neumann_traction (const types::boundary_id &boundary_id,
6153 * const Point<dim> &pt,
6154 * const Tensor<1,dim> &N) const override
6156 * if (this->parameters.load_type == "pressure")
6158 * if (boundary_id == 5)
6160 * const double final_load = this->parameters.load;
6161 * const double current_time = this->time.get_current();
6162 * const double final_time = this->time.get_end();
6163 * const double num_cycles = 3.0;
6165 * return final_load/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5))) * N;
6171 * return Tensor<1,dim>();
6174 * virtual types::boundary_id
6175 * get_reaction_boundary_id_for_output() const override
6180 * virtual std::vector<double>
6181 * get_dirichlet_load(const types::boundary_id &boundary_id,
6182 * const int &direction) const override
6184 * std::vector<double> displ_incr(dim,0.0);
6186 * if ( (boundary_id == 5) && (direction == 2) )
6188 * const double final_displ = this->parameters.load;
6189 * const double current_time = this->time.get_current();
6190 * const double final_time = this->time.get_end();
6191 * const double delta_time = this->time.get_delta_t();
6192 * const double num_cycles = 3.0;
6193 * double current_displ = 0.0;
6194 * double previous_displ = 0.0;
6196 * if (this->parameters.num_cycle_sets == 1)
6198 * current_displ = final_displ/2.0 * (1.0
6199 * - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5)));
6200 * previous_displ = final_displ/2.0 * (1.0
6201 * - std::sin(numbers::PI * (2.0*num_cycles*(current_time-delta_time)/final_time + 0.5)));
6205 * if ( current_time <= (final_time*1.0/3.0) )
6207 * current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6208 * (2.0*num_cycles*current_time/(final_time*1.0/3.0) + 0.5)));
6209 * previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6210 * (2.0*num_cycles*(current_time-delta_time)/(final_time*1.0/3.0) + 0.5)));
6214 * current_displ = final_displ * (1.0 - std::sin(numbers::PI *
6215 * (2.0*num_cycles*current_time / (final_time*2.0/3.0)
6216 * - (num_cycles - 0.5) )));
6217 * previous_displ = final_displ * (1.0 - std::sin(numbers::PI *
6218 * (2.0*num_cycles*(current_time-delta_time) / (final_time*2.0/3.0)
6219 * - (num_cycles - 0.5))));
6222 * displ_incr[2] = current_displ - previous_displ;
6224 * return displ_incr;
6231 * <a name="nonlinear-poro-viscoelasticity.cc-DerivedclassNolateraldisplacementinloadingsurfaces"></a>
6232 * <h4>Derived class: No lateral displacement in loading surfaces</h4>
6235 * template <int dim>
6236 * class BrainBudday2017CubeTensionCompressionFullyFixed
6237 * : public BrainBudday2017BaseCube<dim>
6240 * BrainBudday2017CubeTensionCompressionFullyFixed (const Parameters::AllParameters ¶meters)
6241 * : BrainBudday2017BaseCube<dim> (parameters)
6244 * virtual ~BrainBudday2017CubeTensionCompressionFullyFixed () {}
6248 * define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
6250 * tracked_vertices[0][0] = 0.5*this->parameters.scale;
6251 * tracked_vertices[0][1] = 0.5*this->parameters.scale;
6252 * tracked_vertices[0][2] = 1.0*this->parameters.scale;
6254 * tracked_vertices[1][0] = 0.5*this->parameters.scale;
6255 * tracked_vertices[1][1] = 0.5*this->parameters.scale;
6256 * tracked_vertices[1][2] = 0.5*this->parameters.scale;
6260 * make_dirichlet_constraints(AffineConstraints<double> &constraints) override
6262 * if (this->time.get_timestep() < 2)
6264 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
6266 * Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
6268 * (this->fe.component_mask(this->pressure)));
6272 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6274 * Functions::ZeroFunction<dim>(this->n_components),
6276 * (this->fe.component_mask(this->pressure)));
6279 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6281 * Functions::ZeroFunction<dim>(this->n_components),
6283 * (this->fe.component_mask(this->x_displacement) |
6284 * this->fe.component_mask(this->y_displacement) |
6285 * this->fe.component_mask(this->z_displacement) ));
6288 * if (this->parameters.load_type == "displacement")
6290 * const std::vector<double> value = get_dirichlet_load(5,2);
6291 * const FEValuesExtractors::Scalar direction(this->z_displacement);
6293 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6295 * Functions::ConstantFunction<dim>(value[2],this->n_components),
6297 * this->fe.component_mask(direction) );
6299 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6301 * Functions::ZeroFunction<dim>(this->n_components),
6303 * (this->fe.component_mask(this->x_displacement) |
6304 * this->fe.component_mask(this->y_displacement) ));
6308 * virtual Tensor<1,dim>
6309 * get_neumann_traction (const types::boundary_id &boundary_id,
6310 * const Point<dim> &pt,
6311 * const Tensor<1,dim> &N) const override
6313 * if (this->parameters.load_type == "pressure")
6315 * if (boundary_id == 5)
6317 * const double final_load = this->parameters.load;
6318 * const double current_time = this->time.get_current();
6319 * const double final_time = this->time.get_end();
6320 * const double num_cycles = 3.0;
6322 * return final_load/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5))) * N;
6328 * return Tensor<1,dim>();
6331 * virtual types::boundary_id
6332 * get_reaction_boundary_id_for_output() const override
6337 * virtual std::vector<double>
6338 * get_dirichlet_load(const types::boundary_id &boundary_id,
6339 * const int &direction) const override
6341 * std::vector<double> displ_incr(dim,0.0);
6343 * if ( (boundary_id == 5) && (direction == 2) )
6345 * const double final_displ = this->parameters.load;
6346 * const double current_time = this->time.get_current();
6347 * const double final_time = this->time.get_end();
6348 * const double delta_time = this->time.get_delta_t();
6349 * const double num_cycles = 3.0;
6350 * double current_displ = 0.0;
6351 * double previous_displ = 0.0;
6353 * if (this->parameters.num_cycle_sets == 1)
6355 * current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5)));
6356 * previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*(current_time-delta_time)/final_time + 0.5)));
6360 * if ( current_time <= (final_time*1.0/3.0) )
6362 * current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6363 * (2.0*num_cycles*current_time/(final_time*1.0/3.0) + 0.5)));
6364 * previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6365 * (2.0*num_cycles*(current_time-delta_time)/(final_time*1.0/3.0) + 0.5)));
6369 * current_displ = final_displ * (1.0 - std::sin(numbers::PI *
6370 * (2.0*num_cycles*current_time / (final_time*2.0/3.0)
6371 * - (num_cycles - 0.5) )));
6372 * previous_displ = final_displ * (1.0 - std::sin(numbers::PI *
6373 * (2.0*num_cycles*(current_time-delta_time) / (final_time*2.0/3.0)
6374 * - (num_cycles - 0.5))));
6377 * displ_incr[2] = current_displ - previous_displ;
6379 * return displ_incr;
6386 * <a name="nonlinear-poro-viscoelasticity.cc-DerivedclassNolateralorverticaldisplacementinloadingsurface"></a>
6387 * <h4>Derived class: No lateral or vertical displacement in loading surface</h4>
6390 * template <int dim>
6391 * class BrainBudday2017CubeShearFullyFixed
6392 * : public BrainBudday2017BaseCube<dim>
6395 * BrainBudday2017CubeShearFullyFixed (const Parameters::AllParameters ¶meters)
6396 * : BrainBudday2017BaseCube<dim> (parameters)
6399 * virtual ~BrainBudday2017CubeShearFullyFixed () {}
6403 * define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices) override
6405 * tracked_vertices[0][0] = 0.75*this->parameters.scale;
6406 * tracked_vertices[0][1] = 0.5*this->parameters.scale;
6407 * tracked_vertices[0][2] = 0.0*this->parameters.scale;
6409 * tracked_vertices[1][0] = 0.25*this->parameters.scale;
6410 * tracked_vertices[1][1] = 0.5*this->parameters.scale;
6411 * tracked_vertices[1][2] = 0.0*this->parameters.scale;
6415 * make_dirichlet_constraints(AffineConstraints<double> &constraints) override
6417 * if (this->time.get_timestep() < 2)
6419 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
6421 * Functions::ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
6423 * (this->fe.component_mask(this->pressure)));
6427 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6429 * Functions::ZeroFunction<dim>(this->n_components),
6431 * (this->fe.component_mask(this->pressure)));
6434 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6436 * Functions::ZeroFunction<dim>(this->n_components),
6438 * (this->fe.component_mask(this->x_displacement) |
6439 * this->fe.component_mask(this->y_displacement) |
6440 * this->fe.component_mask(this->z_displacement) ));
6443 * if (this->parameters.load_type == "displacement")
6445 * const std::vector<double> value = get_dirichlet_load(4,0);
6446 * const FEValuesExtractors::Scalar direction(this->x_displacement);
6448 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6450 * Functions::ConstantFunction<dim>(value[0],this->n_components),
6452 * this->fe.component_mask(direction));
6454 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6456 * Functions::ZeroFunction<dim>(this->n_components),
6458 * (this->fe.component_mask(this->y_displacement) |
6459 * this->fe.component_mask(this->z_displacement) ));
6463 * virtual Tensor<1,dim>
6464 * get_neumann_traction (const types::boundary_id &boundary_id,
6465 * const Point<dim> &pt,
6466 * const Tensor<1,dim> &N) const override
6468 * if (this->parameters.load_type == "pressure")
6470 * if (boundary_id == 4)
6472 * const double final_load = this->parameters.load;
6473 * const double current_time = this->time.get_current();
6474 * const double final_time = this->time.get_end();
6475 * const double num_cycles = 3.0;
6476 * const Tensor<1,3> axis ({0.0,1.0,0.0});
6477 * const double angle = numbers::PI;
6478 * static const Tensor< 2, dim, double> R(Physics::Transformations::Rotations::rotation_matrix_3d(axis,angle));
6480 * return (final_load * (std::sin(2.0*(numbers::PI)*num_cycles*current_time/final_time)) * (R * N));
6486 * return Tensor<1,dim>();
6489 * virtual types::boundary_id
6490 * get_reaction_boundary_id_for_output() const override
6495 * virtual std::vector<double>
6496 * get_dirichlet_load(const types::boundary_id &boundary_id,
6497 * const int &direction) const override
6499 * std::vector<double> displ_incr (dim, 0.0);
6501 * if ( (boundary_id == 4) && (direction == 0) )
6503 * const double final_displ = this->parameters.load;
6504 * const double current_time = this->time.get_current();
6505 * const double final_time = this->time.get_end();
6506 * const double delta_time = this->time.get_delta_t();
6507 * const double num_cycles = 3.0;
6508 * double current_displ = 0.0;
6509 * double previous_displ = 0.0;
6511 * if (this->parameters.num_cycle_sets == 1)
6513 * current_displ = final_displ * (std::sin(2.0*(numbers::PI)*num_cycles*current_time/final_time));
6514 * previous_displ = final_displ * (std::sin(2.0*(numbers::PI)*num_cycles*(current_time-delta_time)/final_time));
6518 * AssertThrow(false, ExcMessage("Problem type not defined. Budday shear experiments implemented only for one set of cycles."));
6520 * displ_incr[0] = current_displ - previous_displ;
6522 * return displ_incr;
6531 * <a name="nonlinear-poro-viscoelasticity.cc-Mainfunction"></a>
6532 * <h3>Main function</h3>
6533 * Lastly we provide the main driver function which is similar to the other tutorials.
6536 * int main (int argc, char *argv[])
6538 * using namespace dealii;
6539 * using namespace NonLinearPoroViscoElasticity;
6541 * const unsigned int n_tbb_processes = 1;
6542 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, n_tbb_processes);
6546 * Parameters::AllParameters parameters ("parameters.prm");
6547 * if (parameters.geom_type == "Ehlers_tube_step_load")
6549 * VerificationEhlers1999StepLoad<3> solid_3d(parameters);
6552 * else if (parameters.geom_type == "Ehlers_tube_increase_load")
6554 * VerificationEhlers1999IncreaseLoad<3> solid_3d(parameters);
6557 * else if (parameters.geom_type == "Ehlers_cube_consolidation")
6559 * VerificationEhlers1999CubeConsolidation<3> solid_3d(parameters);
6562 * else if (parameters.geom_type == "Franceschini_consolidation")
6564 * Franceschini2006Consolidation<3> solid_3d(parameters);
6567 * else if (parameters.geom_type == "Budday_cube_tension_compression")
6569 * BrainBudday2017CubeTensionCompression<3> solid_3d(parameters);
6572 * else if (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")
6574 * BrainBudday2017CubeTensionCompressionFullyFixed<3> solid_3d(parameters);
6577 * else if (parameters.geom_type == "Budday_cube_shear_fully_fixed")
6579 * BrainBudday2017CubeShearFullyFixed<3> solid_3d(parameters);
6584 * AssertThrow(false, ExcMessage("Problem type not defined. Current setting: " + parameters.geom_type));
6588 * catch (std::exception &exc)
6590 * if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
6592 * std::cerr << std::endl << std::endl
6593 * << "----------------------------------------------------"
6595 * std::cerr << "Exception on processing: " << std::endl << exc.what()
6596 * << std::endl << "Aborting!" << std::endl
6597 * << "----------------------------------------------------"
6605 * if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
6607 * std::cerr << std::endl << std::endl
6608 * << "----------------------------------------------------"
6610 * std::cerr << "Unknown exception!" << std::endl << "Aborting!"
6612 * << "----------------------------------------------------"
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
constexpr types::blas_int zero
constexpr types::blas_int one
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.)
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void abort(const ExceptionBase &exc) noexcept
bool check(const ConstraintKinds kind_in, const unsigned int dim)
long double gamma(const unsigned int n)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
SymmetricTensorEigenvectorMethod