323 * We start by including all the necessary deal.II header files and some
C++
324 * related ones. They have been discussed in detail in previous tutorial
325 * programs, so you need only refer to past tutorials
for details.
331 * #include <deal.II/base/function.h>
332 * #include <deal.II/base/parameter_handler.h>
333 * #include <deal.II/base/
point.h>
334 * #include <deal.II/base/quadrature_lib.h>
335 * #include <deal.II/base/symmetric_tensor.h>
336 * #include <deal.II/base/tensor.h>
337 * #include <deal.II/base/timer.h>
338 * #include <deal.II/base/work_stream.h>
339 * #include <deal.II/base/mpi.h>
340 * #include <deal.II/base/quadrature_point_data.h>
342 * #include <deal.II/differentiation/ad.h>
344 * #include <deal.II/distributed/shared_tria.h>
346 * #include <deal.II/dofs/dof_renumbering.h>
347 * #include <deal.II/dofs/dof_tools.h>
348 * #include <deal.II/dofs/dof_accessor.h>
350 * #include <deal.II/grid/filtered_iterator.h>
351 * #include <deal.II/grid/grid_generator.h>
352 * #include <deal.II/grid/grid_tools.h>
353 * #include <deal.II/grid/grid_in.h>
354 * #include <deal.II/grid/grid_out.h>
355 * #include <deal.II/grid/manifold_lib.h>
356 * #include <deal.II/grid/tria_accessor.h>
357 * #include <deal.II/grid/tria_boundary_lib.h>
358 * #include <deal.II/grid/tria_iterator.h>
360 * #include <deal.II/fe/fe_dgp_monomial.h>
361 * #include <deal.II/fe/fe_q.h>
362 * #include <deal.II/fe/fe_system.h>
363 * #include <deal.II/fe/fe_tools.h>
364 * #include <deal.II/fe/fe_values.h>
366 * #include <deal.II/lac/block_sparsity_pattern.h>
367 * #include <deal.II/lac/affine_constraints.h>
368 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
369 * #include <deal.II/lac/full_matrix.h>
371 * #include <deal.II/lac/packaged_operation.h>
373 * #include <deal.II/lac/trilinos_block_sparse_matrix.h>
374 * #include <deal.II/lac/trilinos_linear_operator.h>
375 * #include <deal.II/lac/trilinos_parallel_block_vector.h>
376 * #include <deal.II/lac/trilinos_precondition.h>
377 * #include <deal.II/lac/trilinos_sparse_matrix.h>
378 * #include <deal.II/lac/trilinos_sparsity_pattern.h>
379 * #include <deal.II/lac/trilinos_solver.h>
380 * #include <deal.II/lac/trilinos_vector.h>
382 * #include <deal.II/lac/block_vector.h>
383 * #include <deal.II/lac/vector.h>
385 * #include <deal.II/numerics/data_postprocessor.h>
386 * #include <deal.II/numerics/data_out.h>
387 * #include <deal.II/numerics/data_out_faces.h>
388 * #include <deal.II/numerics/fe_field_function.h>
389 * #include <deal.II/numerics/vector_tools.h>
391 * #include <deal.II/physics/transformations.h>
392 * #include <deal.II/physics/elasticity/kinematics.h>
393 * #include <deal.II/physics/elasticity/standard_tensors.h>
395 * #include <iostream>
403 * We create a
namespace for everything that relates to
404 * the nonlinear poro-viscoelastic formulation,
405 * and
import all the deal.II function and
class names into it:
408 * namespace NonLinearPoroViscoElasticity
415 * <a name=
"Runtimeparameters"></a>
416 * <h3>Run-time parameters</h3>
420 * introduced by the user through the file
"parameters.prm"
423 *
namespace Parameters
428 * <a name=
"FiniteElementsystem"></a>
429 * <h4>Finite Element system</h4>
430 * Here we specify the polynomial order used to
approximate the solution,
431 * both
for the displacements and pressure unknowns.
432 * The quadrature order should be adjusted accordingly.
437 *
unsigned int poly_degree_displ;
438 *
unsigned int poly_degree_pore;
439 *
unsigned int quad_order;
450 * prm.enter_subsection(
"Finite element system");
452 * prm.declare_entry(
"Polynomial degree displ",
"2",
454 *
"Displacement system polynomial order");
456 * prm.declare_entry(
"Polynomial degree pore",
"1",
458 *
"Pore pressure system polynomial order");
460 * prm.declare_entry(
"Quadrature order",
"3",
462 *
"Gauss quadrature order");
464 * prm.leave_subsection();
469 * prm.enter_subsection(
"Finite element system");
471 * poly_degree_displ = prm.get_integer(
"Polynomial degree displ");
472 * poly_degree_pore = prm.get_integer(
"Polynomial degree pore");
473 * quad_order = prm.get_integer(
"Quadrature order");
475 * prm.leave_subsection();
481 * <a name=
"Geometry"></a>
483 * These parameters are related to the geometry definition and mesh generation.
484 * We select the type of problem to solve and introduce the desired load
values.
489 * std::string geom_type;
490 *
unsigned int global_refinement;
492 * std::string load_type;
494 *
unsigned int num_cycle_sets;
496 *
double drained_pressure;
507 * prm.enter_subsection(
"Geometry");
509 * prm.declare_entry(
"Geometry type",
"Ehlers_tube_step_load",
511 *
"|Ehlers_tube_increase_load"
512 *
"|Ehlers_cube_consolidation"
513 *
"|Franceschini_consolidation"
514 *
"|Budday_cube_tension_compression"
515 *
"|Budday_cube_tension_compression_fully_fixed"
516 *
"|Budday_cube_shear_fully_fixed"),
517 *
"Type of geometry used. "
518 *
"For Ehlers verification examples see Ehlers and Eipper (1999). "
519 *
"For Franceschini brain consolidation see Franceschini et al. (2006)"
520 *
"For Budday brain examples see Budday et al. (2017)");
522 * prm.declare_entry(
"Global refinement",
"1",
524 *
"Global refinement level");
526 * prm.declare_entry(
"Grid scale",
"1.0",
528 *
"Global grid scaling factor");
530 * prm.declare_entry(
"Load type",
"pressure",
532 *
"Type of loading");
534 * prm.declare_entry(
"Load value",
"-7.5e+6",
538 * prm.declare_entry(
"Number of cycle sets",
"1",
540 *
"Number of times each set of 3 cycles is repeated, only for "
541 *
"Budday_cube_tension_compression and Budday_cube_tension_compression_fully_fixed. "
542 *
"Load value is doubled in second set, load rate is kept constant."
543 *
"Final time indicates end of second cycle set.");
545 * prm.declare_entry(
"Fluid flow value",
"0.0",
547 *
"Prescribed fluid flow. Not implemented in any example yet.");
549 * prm.declare_entry(
"Drained pressure",
"0.0",
551 *
"Increase of pressure value at drained boundary w.r.t the atmospheric pressure.");
553 * prm.leave_subsection();
558 * prm.enter_subsection(
"Geometry");
560 * geom_type = prm.get(
"Geometry type");
561 * global_refinement = prm.get_integer(
"Global refinement");
562 *
scale = prm.get_double(
"Grid scale");
563 * load_type = prm.get(
"Load type");
564 * load = prm.get_double(
"Load value");
565 * num_cycle_sets = prm.get_integer(
"Number of cycle sets");
566 * fluid_flow = prm.get_double(
"Fluid flow value");
567 * drained_pressure = prm.get_double(
"Drained pressure");
569 * prm.leave_subsection();
575 * <a name=
"Materials"></a>
580 * Here we select the type of material
for the solid component
581 * and define the corresponding material parameters.
582 * Then we define he fluid data, including the type of
583 * seepage velocity definition to use.
588 * std::string mat_type;
594 *
double alpha1_infty;
595 *
double alpha2_infty;
596 *
double alpha3_infty;
600 *
double alpha1_mode_1;
601 *
double alpha2_mode_1;
602 *
double alpha3_mode_1;
603 *
double viscosity_mode_1;
604 * std::string fluid_type;
605 *
double solid_vol_frac;
606 *
double kappa_darcy;
607 *
double init_intrinsic_perm;
608 *
double viscosity_FR;
609 *
double init_darcy_coef;
612 *
int gravity_direction;
613 *
double gravity_value;
627 * prm.enter_subsection(
"Material properties");
629 * prm.declare_entry(
"material",
"Neo-Hooke",
631 *
"Type of material used in the problem");
633 * prm.declare_entry(
"lambda",
"8.375e6",
635 *
"First Lamé parameter for extension function related to compactation point in solid material [Pa].");
637 * prm.declare_entry(
"shear modulus",
"5.583e6",
639 *
"shear modulus for Neo-Hooke materials [Pa].");
641 * prm.declare_entry(
"eigen solver",
"QL Implicit Shifts",
643 *
"The type of eigen solver to be used for Ogden and visco-Ogden models.");
645 * prm.declare_entry(
"mu1",
"0.0",
647 *
"Shear material parameter 'mu1' for Ogden material [Pa].");
649 * prm.declare_entry(
"mu2",
"0.0",
651 *
"Shear material parameter 'mu2' for Ogden material [Pa].");
653 * prm.declare_entry(
"mu3",
"0.0",
655 *
"Shear material parameter 'mu1' for Ogden material [Pa].");
657 * prm.declare_entry(
"alpha1",
"1.0",
659 *
"Stiffness material parameter 'alpha1' for Ogden material [-].");
661 * prm.declare_entry(
"alpha2",
"1.0",
663 *
"Stiffness material parameter 'alpha2' for Ogden material [-].");
665 * prm.declare_entry(
"alpha3",
"1.0",
667 *
"Stiffness material parameter 'alpha3' for Ogden material [-].");
669 * prm.declare_entry(
"mu1_1",
"0.0",
671 *
"Shear material parameter 'mu1' for first viscous mode in Ogden material [Pa].");
673 * prm.declare_entry(
"mu2_1",
"0.0",
675 *
"Shear material parameter 'mu2' for first viscous mode in Ogden material [Pa].");
677 * prm.declare_entry(
"mu3_1",
"0.0",
679 *
"Shear material parameter 'mu1' for first viscous mode in Ogden material [Pa].");
681 * prm.declare_entry(
"alpha1_1",
"1.0",
683 *
"Stiffness material parameter 'alpha1' for first viscous mode in Ogden material [-].");
685 * prm.declare_entry(
"alpha2_1",
"1.0",
687 *
"Stiffness material parameter 'alpha2' for first viscous mode in Ogden material [-].");
689 * prm.declare_entry(
"alpha3_1",
"1.0",
691 *
"Stiffness material parameter 'alpha3' for first viscous mode in Ogden material [-].");
693 * prm.declare_entry(
"viscosity_1",
"1e-10",
695 *
"Deformation-independent viscosity parameter 'eta_1' for first viscous mode in Ogden material [-].");
697 * prm.declare_entry(
"seepage definition",
"Ehlers",
699 *
"Type of formulation used to define the seepage velocity in the problem. "
700 *
"Choose between Markert formulation of deformation-dependent intrinsic permeability "
701 *
"and Ehlers formulation of deformation-dependent Darcy flow coefficient.");
703 * prm.declare_entry(
"initial solid volume fraction",
"0.67",
705 *
"Initial porosity (solid volume fraction, 0 < n_0s < 1)");
707 * prm.declare_entry(
"kappa",
"0.0",
709 *
"Deformation-dependency control parameter for specific permeability (kappa >= 0)");
711 * prm.declare_entry(
"initial intrinsic permeability",
"0.0",
713 *
"Initial intrinsic permeability parameter [m^2] (isotropic permeability). To be used with Markert formulation.");
715 * prm.declare_entry(
"fluid viscosity",
"0.0",
717 *
"Effective shear viscosity parameter of the fluid [Pa·s, (N·s)/m^2]. To be used with Markert formulation.");
719 * prm.declare_entry(
"initial Darcy coefficient",
"1.0e-4",
721 *
"Initial Darcy flow coefficient [m/s] (isotropic permeability). To be used with Ehlers formulation.");
723 * prm.declare_entry(
"fluid weight",
"1.0e4",
725 *
"Effective weight of the fluid [N/m^3]. To be used with Ehlers formulation.");
727 * prm.declare_entry(
"gravity term",
"false",
729 *
"Gravity term considered (true) or neglected (false)");
731 * prm.declare_entry(
"fluid density",
"1.0",
733 *
"Real (or effective) density of the fluid");
735 * prm.declare_entry(
"solid density",
"1.0",
737 *
"Real (or effective) density of the solid");
739 * prm.declare_entry(
"gravity direction",
"2",
741 *
"Direction of gravity (unit vector 0 for x, 1 for y, 2 for z)");
743 * prm.declare_entry(
"gravity value",
"-9.81",
745 *
"Value of gravity (be careful to have consistent units!)");
747 * prm.leave_subsection();
752 * prm.enter_subsection(
"Material properties");
759 * mat_type = prm.get(
"material");
760 *
lambda = prm.get_double(
"lambda");
761 *
mu = prm.get_double(
"shear modulus");
762 * mu1_infty = prm.get_double(
"mu1");
763 * mu2_infty = prm.get_double(
"mu2");
764 * mu3_infty = prm.get_double(
"mu3");
765 * alpha1_infty = prm.get_double(
"alpha1");
766 * alpha2_infty = prm.get_double(
"alpha2");
767 * alpha3_infty = prm.get_double(
"alpha3");
768 * mu1_mode_1 = prm.get_double(
"mu1_1");
769 * mu2_mode_1 = prm.get_double(
"mu2_1");
770 * mu3_mode_1 = prm.get_double(
"mu3_1");
771 * alpha1_mode_1 = prm.get_double(
"alpha1_1");
772 * alpha2_mode_1 = prm.get_double(
"alpha2_1");
773 * alpha3_mode_1 = prm.get_double(
"alpha3_1");
774 * viscosity_mode_1 = prm.get_double(
"viscosity_1");
780 * fluid_type = prm.get(
"seepage definition");
781 * solid_vol_frac = prm.get_double(
"initial solid volume fraction");
782 * kappa_darcy = prm.get_double(
"kappa");
783 * init_intrinsic_perm = prm.get_double(
"initial intrinsic permeability");
784 * viscosity_FR = prm.get_double(
"fluid viscosity");
785 * init_darcy_coef = prm.get_double(
"initial Darcy coefficient");
786 * weight_FR = prm.get_double(
"fluid weight");
792 * gravity_term = prm.get_bool(
"gravity term");
793 * density_FR = prm.get_double(
"fluid density");
794 * density_SR = prm.get_double(
"solid density");
795 * gravity_direction = prm.get_integer(
"gravity direction");
796 * gravity_value = prm.get_double(
"gravity value");
798 *
if ( (fluid_type ==
"Markert") && ((init_intrinsic_perm == 0.0) || (viscosity_FR == 0.0)) )
800 *
"'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0."));
802 *
if ( (fluid_type ==
"Ehlers") && ((init_darcy_coef == 0.0) || (weight_FR == 0.0)) )
804 *
"'initial Darcy coefficient' and 'fluid weight' greater than 0.0."));
806 *
const std::string eigen_solver_type = prm.get(
"eigen solver");
807 *
if (eigen_solver_type ==
"QL Implicit Shifts")
809 *
else if (eigen_solver_type ==
"Jacobi")
816 * prm.leave_subsection();
822 * <a name=
"Nonlinearsolver"></a>
823 * <h4>Nonlinear solver</h4>
827 * We now define the tolerances and the maximum number of iterations
for the
828 * Newton-Raphson scheme used to solve the nonlinear system of governing equations.
831 *
struct NonlinearSolver
833 *
unsigned int max_iterations_NR;
836 *
double tol_p_fluid;
847 * prm.enter_subsection(
"Nonlinear solver");
849 * prm.declare_entry(
"Max iterations Newton-Raphson",
"15",
851 *
"Number of Newton-Raphson iterations allowed");
853 * prm.declare_entry(
"Tolerance force",
"1.0e-8",
855 *
"Force residual tolerance");
857 * prm.declare_entry(
"Tolerance displacement",
"1.0e-6",
859 *
"Displacement error tolerance");
861 * prm.declare_entry(
"Tolerance pore pressure",
"1.0e-6",
863 *
"Pore pressure error tolerance");
865 * prm.leave_subsection();
870 * prm.enter_subsection(
"Nonlinear solver");
872 * max_iterations_NR = prm.get_integer(
"Max iterations Newton-Raphson");
873 * tol_f = prm.get_double(
"Tolerance force");
874 * tol_u = prm.get_double(
"Tolerance displacement");
875 * tol_p_fluid = prm.get_double(
"Tolerance pore pressure");
877 * prm.leave_subsection();
883 * <a name=
"Time"></a>
885 * Here we
set the timestep size @f$ \varDelta t @f$ and the simulation
end-time.
901 * prm.enter_subsection(
"Time");
903 * prm.declare_entry(
"End time",
"10.0",
907 * prm.declare_entry(
"Time step size",
"0.002",
909 *
"Time step size. The value must be larger than the displacement error tolerance defined.");
911 * prm.leave_subsection();
916 * prm.enter_subsection(
"Time");
918 * end_time = prm.get_double(
"End time");
919 * delta_t = prm.get_double(
"Time step size");
921 * prm.leave_subsection();
928 * <a name=
"Output"></a>
930 * We can choose the frequency of the data
for the output files.
936 * std::string outfiles_requested;
937 *
unsigned int timestep_output;
938 * std::string outtype;
949 * prm.enter_subsection(
"Output parameters");
951 * prm.declare_entry(
"Output files",
"true",
953 *
"Paraview output files to generate.");
954 * prm.declare_entry(
"Time step number output",
"1",
956 *
"Output data for time steps multiple of the given "
958 * prm.declare_entry(
"Averaged results",
"nodes",
960 *
"Output data associated with integration point values"
961 *
" averaged on elements or on nodes.");
963 * prm.leave_subsection();
968 * prm.enter_subsection(
"Output parameters");
970 * outfiles_requested = prm.get(
"Output files");
971 * timestep_output = prm.get_integer(
"Time step number output");
972 * outtype = prm.get(
"Averaged results");
974 * prm.leave_subsection();
980 * <a name=
"Allparameters"></a>
981 * <h4>All parameters</h4>
982 * We
finally consolidate all of the above structures into a single container that holds all the
run-time selections.
985 *
struct AllParameters :
public FESystem,
988 *
public NonlinearSolver,
992 * AllParameters(
const std::string &input_file);
1001 * AllParameters::AllParameters(
const std::string &input_file)
1004 * declare_parameters(prm);
1006 * parse_parameters(prm);
1011 * FESystem::declare_parameters(prm);
1012 * Geometry::declare_parameters(prm);
1013 * Materials::declare_parameters(prm);
1014 * NonlinearSolver::declare_parameters(prm);
1015 * Time::declare_parameters(prm);
1016 * OutputParam::declare_parameters(prm);
1021 * FESystem::parse_parameters(prm);
1022 * Geometry::parse_parameters(prm);
1023 * Materials::parse_parameters(prm);
1024 * NonlinearSolver::parse_parameters(prm);
1025 * Time::parse_parameters(prm);
1026 * OutputParam::parse_parameters(prm);
1033 * <a name=
"Timeclass"></a>
1034 * <h3>Time
class</h3>
1035 *
A simple
class to store time data.
1036 * For simplicity we assume a constant time step size.
1042 * Time (
const double time_end,
1043 *
const double delta_t)
1046 * time_current(0.0),
1047 * time_end(time_end),
1054 *
double get_current() const
1056 *
return time_current;
1058 *
double get_end() const
1062 *
double get_delta_t() const
1066 *
unsigned int get_timestep() const
1070 *
void increment_time ()
1072 * time_current += delta_t;
1077 *
unsigned int timestep;
1078 *
double time_current;
1080 *
const double delta_t;
1086 * <a name=
"Constitutiveequationforthesolidcomponentofthebiphasicmaterial"></a>
1087 * <h3>Constitutive equation
for the solid component of the biphasic material</h3>
1092 * <a name=
"Baseclassgenerichyperelasticmaterial"></a>
1093 * <h4>Base
class:
generic hyperelastic material</h4>
1094 * The ``extra
" Kirchhoff stress in the solid component is the sum of isochoric
1095 * and a volumetric part.
1096 * @f$\mathbf{\tau} = \mathbf{\tau}_E^{(\bullet)} + \mathbf{\tau}^{\textrm{vol}}@f$
1097 * The deviatoric part changes depending on the type of material model selected:
1098 * Neo-Hooken hyperelasticity, Ogden hyperelasticiy,
1099 * or a single-mode finite viscoelasticity based on the Ogden hyperelastic model.
1100 * In this base class we declare it as a virtual function,
1101 * and it will be defined for each model type in the corresponding derived class.
1102 * We define here the volumetric component, which depends on the
1103 * extension function @f$U(J_S)@f$ selected, and in this case is the same for all models.
1104 * We use the function proposed by
1105 * Ehlers & Eipper 1999 doi:10.1023/A:1006565509095
1106 * We also define some public functions to access and update the internal variables.
1109 * template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1110 * class Material_Hyperelastic
1113 * Material_Hyperelastic(const Parameters::AllParameters ¶meters,
1116 * n_OS (parameters.solid_vol_frac),
1117 * lambda (parameters.lambda),
1120 * det_F_converged (1.0),
1121 * eigen_solver (parameters.eigen_solver)
1123 * ~Material_Hyperelastic()
1126 * SymmetricTensor<2, dim, NumberType>
1127 * get_tau_E(const Tensor<2,dim, NumberType> &F) const
1129 * return ( get_tau_E_base(F) + get_tau_E_ext_func(F) );
1132 * SymmetricTensor<2, dim, NumberType>
1133 * get_Cauchy_E(const Tensor<2, dim, NumberType> &F) const
1135 * const NumberType det_F = determinant(F);
1136 * Assert(det_F > 0, ExcInternalError());
1137 * return get_tau_E(F)*NumberType(1/det_F);
1141 * get_converged_det_F() const
1143 * return det_F_converged;
1147 * update_end_timestep()
1149 * det_F_converged = det_F;
1153 * update_internal_equilibrium( const Tensor<2, dim, NumberType> &F )
1155 * det_F = Tensor<0,dim,double>(determinant(F));
1159 * get_viscous_dissipation( ) const = 0;
1161 * const double n_OS;
1162 * const double lambda;
1165 * double det_F_converged;
1166 * const enum SymmetricTensorEigenvectorMethod eigen_solver;
1169 * SymmetricTensor<2, dim, NumberType>
1170 * get_tau_E_ext_func(const Tensor<2,dim, NumberType> &F) const
1172 * const NumberType det_F = determinant(F);
1173 * Assert(det_F > 0, ExcInternalError());
1175 * static const SymmetricTensor< 2, dim, double>
1176 * I (Physics::Elasticity::StandardTensors<dim>::I);
1177 * return ( NumberType(lambda * (1.0-n_OS)*(1.0-n_OS)
1178 * * (det_F/(1.0-n_OS) - det_F/(det_F-n_OS))) * I );
1181 * virtual SymmetricTensor<2, dim, NumberType>
1182 * get_tau_E_base(const Tensor<2,dim, NumberType> &F) const = 0;
1188 * <a name="DerivedclassNeoHookeanhyperelasticmaterial
"></a>
1189 * <h4>Derived class: Neo-Hookean hyperelastic material</h4>
1192 * template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1193 * class NeoHooke : public Material_Hyperelastic < dim, NumberType >
1196 * NeoHooke(const Parameters::AllParameters ¶meters,
1199 * Material_Hyperelastic< dim, NumberType > (parameters,time),
1202 * virtual ~NeoHooke()
1206 * get_viscous_dissipation() const
1214 * SymmetricTensor<2, dim, NumberType>
1215 * get_tau_E_base(const Tensor<2,dim, NumberType> &F) const
1217 * static const SymmetricTensor< 2, dim, double>
1218 * I (Physics::Elasticity::StandardTensors<dim>::I);
1220 * const bool use_standard_model = true;
1222 * if (use_standard_model)
1226 * Standard Neo-Hooke
1229 * return ( mu * ( symmetrize(F * transpose(F)) - I ) );
1235 * Neo-Hooke in terms of principal stretches
1238 * const SymmetricTensor<2, dim, NumberType>
1239 * B = symmetrize(F * transpose(F));
1240 * const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1241 * eigen_B = eigenvectors(B, this->eigen_solver);
1243 * SymmetricTensor<2, dim, NumberType> B_ev;
1244 * for (unsigned int d=0; d<dim; ++d)
1245 * B_ev += eigen_B[d].first*symmetrize(outer_product(eigen_B[d].second,eigen_B[d].second));
1247 * return ( mu*(B_ev-I) );
1255 * <a name="DerivedclassOgdenhyperelasticmaterial
"></a>
1256 * <h4>Derived class: Ogden hyperelastic material</h4>
1259 * template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1260 * class Ogden : public Material_Hyperelastic < dim, NumberType >
1263 * Ogden(const Parameters::AllParameters ¶meters,
1266 * Material_Hyperelastic< dim, NumberType > (parameters,time),
1267 * mu({parameters.mu1_infty,
1268 * parameters.mu2_infty,
1269 * parameters.mu3_infty}),
1270 * alpha({parameters.alpha1_infty,
1271 * parameters.alpha2_infty,
1272 * parameters.alpha3_infty})
1278 * get_viscous_dissipation() const
1284 * std::vector<double> mu;
1285 * std::vector<double> alpha;
1287 * SymmetricTensor<2, dim, NumberType>
1288 * get_tau_E_base(const Tensor<2,dim, NumberType> &F) const
1290 * const SymmetricTensor<2, dim, NumberType>
1291 * B = symmetrize(F * transpose(F));
1293 * const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1294 * eigen_B = eigenvectors(B, this->eigen_solver);
1296 * SymmetricTensor<2, dim, NumberType> tau;
1297 * static const SymmetricTensor< 2, dim, double>
1298 * I (Physics::Elasticity::StandardTensors<dim>::I);
1300 * for (unsigned int i = 0; i < 3; ++i)
1302 * for (unsigned int A = 0; A < dim; ++A)
1304 * SymmetricTensor<2, dim, NumberType> tau_aux1 = symmetrize(
1305 * outer_product(eigen_B[A].second,eigen_B[A].second));
1306 * tau_aux1 *= mu[i]*std::pow(eigen_B[A].first, (alpha[i]/2.) );
1309 * SymmetricTensor<2, dim, NumberType> tau_aux2 (I);
1310 * tau_aux2 *= mu[i];
1320 * <a name="DerivedclassSinglemodeOgdenviscoelasticmaterial
"></a>
1321 * <h4>Derived class: Single-mode Ogden viscoelastic material</h4>
1322 * We use the finite viscoelastic model described in
1323 * Reese & Govindjee (1998) doi:10.1016/S0020-7683(97)00217-5
1324 * The algorithm for the implicit exponential time integration is given in
1325 * Budday et al. (2017) doi: 10.1016/j.actbio.2017.06.024
1328 * template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1329 * class visco_Ogden : public Material_Hyperelastic < dim, NumberType >
1332 * visco_Ogden(const Parameters::AllParameters ¶meters,
1335 * Material_Hyperelastic< dim, NumberType > (parameters,time),
1336 * mu_infty({parameters.mu1_infty,
1337 * parameters.mu2_infty,
1338 * parameters.mu3_infty}),
1339 * alpha_infty({parameters.alpha1_infty,
1340 * parameters.alpha2_infty,
1341 * parameters.alpha3_infty}),
1342 * mu_mode_1({parameters.mu1_mode_1,
1343 * parameters.mu2_mode_1,
1344 * parameters.mu3_mode_1}),
1345 * alpha_mode_1({parameters.alpha1_mode_1,
1346 * parameters.alpha2_mode_1,
1347 * parameters.alpha3_mode_1}),
1348 * viscosity_mode_1(parameters.viscosity_mode_1),
1349 * Cinv_v_1(Physics::Elasticity::StandardTensors<dim>::I),
1350 * Cinv_v_1_converged(Physics::Elasticity::StandardTensors<dim>::I)
1352 * virtual ~visco_Ogden()
1356 * update_internal_equilibrium( const Tensor<2, dim, NumberType> &F )
1358 * Material_Hyperelastic < dim, NumberType >::update_internal_equilibrium(F);
1360 * this->Cinv_v_1 = this->Cinv_v_1_converged;
1361 * SymmetricTensor<2, dim, NumberType> B_e_1_tr = symmetrize(F * this->Cinv_v_1 * transpose(F));
1363 * const std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim >
1364 * eigen_B_e_1_tr = eigenvectors(B_e_1_tr, this->eigen_solver);
1366 * Tensor< 1, dim, NumberType > lambdas_e_1_tr;
1367 * Tensor< 1, dim, NumberType > epsilon_e_1_tr;
1368 * for (int a = 0; a < dim; ++a)
1370 * lambdas_e_1_tr[a] = std::sqrt(eigen_B_e_1_tr[a].first);
1371 * epsilon_e_1_tr[a] = std::log(lambdas_e_1_tr[a]);
1374 * const double tolerance = 1e-8;
1375 * double residual_check = tolerance*10.0;
1376 * Tensor< 1, dim, NumberType > residual;
1377 * Tensor< 2, dim, NumberType > tangent;
1378 * static const SymmetricTensor< 2, dim, double> I(Physics::Elasticity::StandardTensors<dim>::I);
1379 * NumberType J_e_1 = std::sqrt(determinant(B_e_1_tr));
1381 * std::vector<NumberType> lambdas_e_1_iso(dim);
1382 * SymmetricTensor<2, dim, NumberType> B_e_1;
1383 * int iteration = 0;
1385 * Tensor< 1, dim, NumberType > lambdas_e_1;
1386 * Tensor< 1, dim, NumberType > epsilon_e_1;
1387 * epsilon_e_1 = epsilon_e_1_tr;
1389 * while(residual_check > tolerance)
1391 * NumberType aux_J_e_1 = 1.0;
1392 * for (unsigned int a = 0; a < dim; ++a)
1394 * lambdas_e_1[a] = std::exp(epsilon_e_1[a]);
1395 * aux_J_e_1 *= lambdas_e_1[a];
1398 * J_e_1 = aux_J_e_1;
1400 * for (unsigned int a = 0; a < dim; ++a)
1401 * lambdas_e_1_iso[a] = lambdas_e_1[a]*std::pow(J_e_1,-1.0/dim);
1403 * for (unsigned int a = 0; a < dim; ++a)
1405 * residual[a] = get_beta_mode_1(lambdas_e_1_iso, a);
1406 * residual[a] *= this->time.get_delta_t()/(2.0*viscosity_mode_1);
1407 * residual[a] += epsilon_e_1[a];
1408 * residual[a] -= epsilon_e_1_tr[a];
1410 * for (unsigned int b = 0; b < dim; ++b)
1412 * tangent[a][b] = get_gamma_mode_1(lambdas_e_1_iso, a, b);
1413 * tangent[a][b] *= this->time.get_delta_t()/(2.0*viscosity_mode_1);
1414 * tangent[a][b] += I[a][b];
1418 * epsilon_e_1 -= invert(tangent)*residual;
1420 * residual_check = 0.0;
1421 * for (unsigned int a = 0; a < dim; ++a)
1423 * if ( std::abs(residual[a]) > residual_check)
1424 * residual_check = std::abs(Tensor<0,dim,double>(residual[a]));
1427 * if (iteration > 15 )
1428 * AssertThrow(false, ExcMessage("No convergence in local Newton iteration
for the
"
1429 * "viscoelastic exponential time integration algorithm.
"));
1432 * NumberType aux_J_e_1 = 1.0;
1433 * for (unsigned int a = 0; a < dim; ++a)
1435 * lambdas_e_1[a] = std::exp(epsilon_e_1[a]);
1436 * aux_J_e_1 *= lambdas_e_1[a];
1438 * J_e_1 = aux_J_e_1;
1440 * for (unsigned int a = 0; a < dim; ++a)
1441 * lambdas_e_1_iso[a] = lambdas_e_1[a]*std::pow(J_e_1,-1.0/dim);
1443 * for (unsigned int a = 0; a < dim; ++a)
1445 * SymmetricTensor<2, dim, NumberType>
1446 * B_e_1_aux = symmetrize(outer_product(eigen_B_e_1_tr[a].second,eigen_B_e_1_tr[a].second));
1447 * B_e_1_aux *= lambdas_e_1[a] * lambdas_e_1[a];
1448 * B_e_1 += B_e_1_aux;
1451 * Tensor<2, dim, NumberType>Cinv_v_1_AD = symmetrize(invert(F) * B_e_1 * invert(transpose(F)));
1453 * this->tau_neq_1 = 0;
1454 * for (unsigned int a = 0; a < dim; ++a)
1456 * SymmetricTensor<2, dim, NumberType>
1457 * tau_neq_1_aux = symmetrize(outer_product(eigen_B_e_1_tr[a].second,eigen_B_e_1_tr[a].second));
1458 * tau_neq_1_aux *= get_beta_mode_1(lambdas_e_1_iso, a);
1459 * this->tau_neq_1 += tau_neq_1_aux;
1467 * for (unsigned int a = 0; a < dim; ++a)
1468 * for (unsigned int b = 0; b < dim; ++b)
1469 * this->Cinv_v_1[a][b]= Tensor<0,dim,double>(Cinv_v_1_AD[a][b]);
1472 * void update_end_timestep()
1474 * Material_Hyperelastic < dim, NumberType >::update_end_timestep();
1475 * this->Cinv_v_1_converged = this->Cinv_v_1;
1478 * double get_viscous_dissipation() const
1480 * NumberType dissipation_term = get_tau_E_neq() * get_tau_E_neq(); //Double contract the two SymmetricTensor
1481 * dissipation_term /= (2*viscosity_mode_1);
1483 * return dissipation_term.val();
1487 * std::vector<double> mu_infty;
1488 * std::vector<double> alpha_infty;
1489 * std::vector<double> mu_mode_1;
1490 * std::vector<double> alpha_mode_1;
1491 * double viscosity_mode_1;
1492 * SymmetricTensor<2, dim, double> Cinv_v_1;
1493 * SymmetricTensor<2, dim, double> Cinv_v_1_converged;
1494 * SymmetricTensor<2, dim, NumberType> tau_neq_1;
1496 * SymmetricTensor<2, dim, NumberType>
1497 * get_tau_E_base(const Tensor<2,dim, NumberType> &F) const
1499 * return ( get_tau_E_neq() + get_tau_E_eq(F) );
1502 * SymmetricTensor<2, dim, NumberType>
1503 * get_tau_E_eq(const Tensor<2,dim, NumberType> &F) const
1505 * const SymmetricTensor<2, dim, NumberType> B = symmetrize(F * transpose(F));
1507 * std::array< std::pair< NumberType, Tensor< 1, dim, NumberType > >, dim > eigen_B;
1508 * eigen_B = eigenvectors(B, this->eigen_solver);
1510 * SymmetricTensor<2, dim, NumberType> tau;
1511 * static const SymmetricTensor< 2, dim, double>
1512 * I (Physics::Elasticity::StandardTensors<dim>::I);
1514 * for (unsigned int i = 0; i < 3; ++i)
1516 * for (unsigned int A = 0; A < dim; ++A)
1518 * SymmetricTensor<2, dim, NumberType> tau_aux1 = symmetrize(
1519 * outer_product(eigen_B[A].second,eigen_B[A].second));
1520 * tau_aux1 *= mu_infty[i]*std::pow(eigen_B[A].first, (alpha_infty[i]/2.) );
1523 * SymmetricTensor<2, dim, NumberType> tau_aux2 (I);
1524 * tau_aux2 *= mu_infty[i];
1530 * SymmetricTensor<2, dim, NumberType>
1531 * get_tau_E_neq() const
1537 * get_beta_mode_1(std::vector< NumberType > &lambda, const int &A) const
1539 * NumberType beta = 0.0;
1541 * for (unsigned int i = 0; i < 3; ++i) //3rd-order Ogden model
1544 * NumberType aux = 0.0;
1545 * for (int p = 0; p < dim; ++p)
1546 * aux += std::pow(lambda[p],alpha_mode_1[i]);
1549 * aux += std::pow(lambda[A], alpha_mode_1[i]);
1550 * aux *= mu_mode_1[i];
1558 * get_gamma_mode_1(std::vector< NumberType > &lambda,
1560 * const int &B ) const
1562 * NumberType gamma = 0.0;
1566 * for (unsigned int i = 0; i < 3; ++i)
1568 * NumberType aux = 0.0;
1569 * for (int p = 0; p < dim; ++p)
1570 * aux += std::pow(lambda[p],alpha_mode_1[i]);
1572 * aux *= 1.0/(dim*dim);
1573 * aux += 1.0/dim * std::pow(lambda[A], alpha_mode_1[i]);
1574 * aux *= mu_mode_1[i]*alpha_mode_1[i];
1581 * for (unsigned int i = 0; i < 3; ++i)
1583 * NumberType aux = 0.0;
1584 * for (int p = 0; p < dim; ++p)
1585 * aux += std::pow(lambda[p],alpha_mode_1[i]);
1587 * aux *= 1.0/(dim*dim);
1588 * aux -= 1.0/dim * std::pow(lambda[A], alpha_mode_1[i]);
1589 * aux -= 1.0/dim * std::pow(lambda[B], alpha_mode_1[i]);
1590 * aux *= mu_mode_1[i]*alpha_mode_1[i];
1604 * <a name="Constitutiveequationforthefluidcomponentofthebiphasicmaterial
"></a>
1605 * <h3>Constitutive equation for the fluid component of the biphasic material</h3>
1606 * We consider two slightly different definitions to define the seepage velocity with a Darcy-like law.
1607 * Ehlers & Eipper 1999, doi:10.1023/A:1006565509095
1608 * Markert 2007, doi:10.1007/s11242-007-9107-6
1609 * The selection of one or another is made by the user via the parameters file.
1612 * template <int dim, typename NumberType = Sacado::Fad::DFad<double> >
1613 * class Material_Darcy_Fluid
1616 * Material_Darcy_Fluid(const Parameters::AllParameters ¶meters)
1618 * fluid_type(parameters.fluid_type),
1619 * n_OS(parameters.solid_vol_frac),
1620 * initial_intrinsic_permeability(parameters.init_intrinsic_perm),
1621 * viscosity_FR(parameters.viscosity_FR),
1622 * initial_darcy_coefficient(parameters.init_darcy_coef),
1623 * weight_FR(parameters.weight_FR),
1624 * kappa_darcy(parameters.kappa_darcy),
1625 * gravity_term(parameters.gravity_term),
1626 * density_FR(parameters.density_FR),
1627 * gravity_direction(parameters.gravity_direction),
1628 * gravity_value(parameters.gravity_value)
1630 * Assert(kappa_darcy >= 0, ExcInternalError());
1632 * ~Material_Darcy_Fluid()
1635 * Tensor<1, dim, NumberType> get_seepage_velocity_current
1636 * (const Tensor<2,dim, NumberType> &F,
1637 * const Tensor<1,dim, NumberType> &grad_p_fluid) const
1639 * const NumberType det_F = determinant(F);
1640 * Assert(det_F > 0.0, ExcInternalError());
1642 * Tensor<2, dim, NumberType> permeability_term;
1644 * if (fluid_type == "Markert
")
1645 * permeability_term = get_instrinsic_permeability_current(F) / viscosity_FR;
1647 * else if (fluid_type == "Ehlers
")
1648 * permeability_term = get_darcy_flow_current(F) / weight_FR;
1651 * AssertThrow(false, ExcMessage(
1652 * "Material_Darcy_Fluid --> Only Markert
"
1653 * "and Ehlers formulations have been implemented.
"));
1655 * return ( -1.0 * permeability_term * det_F
1656 * * (grad_p_fluid - get_body_force_FR_current()) );
1659 * double get_porous_dissipation(const Tensor<2,dim, NumberType> &F,
1660 * const Tensor<1,dim, NumberType> &grad_p_fluid) const
1662 * NumberType dissipation_term;
1663 * Tensor<1, dim, NumberType> seepage_velocity;
1664 * Tensor<2, dim, NumberType> permeability_term;
1666 * const NumberType det_F = determinant(F);
1667 * Assert(det_F > 0.0, ExcInternalError());
1669 * if (fluid_type == "Markert
")
1671 * permeability_term = get_instrinsic_permeability_current(F) / viscosity_FR;
1672 * seepage_velocity = get_seepage_velocity_current(F,grad_p_fluid);
1674 * else if (fluid_type == "Ehlers
")
1676 * permeability_term = get_darcy_flow_current(F) / weight_FR;
1677 * seepage_velocity = get_seepage_velocity_current(F,grad_p_fluid);
1680 * AssertThrow(false, ExcMessage(
1681 * "Material_Darcy_Fluid --> Only Markert and Ehlers
"
1682 * "formulations have been implemented.
"));
1684 * dissipation_term = ( invert(permeability_term) * seepage_velocity ) * seepage_velocity;
1685 * dissipation_term *= 1.0/(det_F*det_F);
1686 * return Tensor<0,dim,double>(dissipation_term);
1690 * const std::string fluid_type;
1691 * const double n_OS;
1692 * const double initial_intrinsic_permeability;
1693 * const double viscosity_FR;
1694 * const double initial_darcy_coefficient;
1695 * const double weight_FR;
1696 * const double kappa_darcy;
1697 * const bool gravity_term;
1698 * const double density_FR;
1699 * const int gravity_direction;
1700 * const double gravity_value;
1702 * Tensor<2, dim, NumberType>
1703 * get_instrinsic_permeability_current(const Tensor<2,dim, NumberType> &F) const
1705 * static const SymmetricTensor< 2, dim, double>
1706 * I (Physics::Elasticity::StandardTensors<dim>::I);
1707 * const Tensor<2, dim, NumberType> initial_instrinsic_permeability_tensor
1708 * = Tensor<2, dim, double>(initial_intrinsic_permeability * I);
1710 * const NumberType det_F = determinant(F);
1711 * Assert(det_F > 0.0, ExcInternalError());
1713 * const NumberType fraction = (det_F - n_OS)/(1 - n_OS);
1714 * return ( NumberType (std::pow(fraction, kappa_darcy))
1715 * * initial_instrinsic_permeability_tensor );
1718 * Tensor<2, dim, NumberType>
1719 * get_darcy_flow_current(const Tensor<2,dim, NumberType> &F) const
1721 * static const SymmetricTensor< 2, dim, double>
1722 * I (Physics::Elasticity::StandardTensors<dim>::I);
1723 * const Tensor<2, dim, NumberType> initial_darcy_flow_tensor
1724 * = Tensor<2, dim, double>(initial_darcy_coefficient * I);
1726 * const NumberType det_F = determinant(F);
1727 * Assert(det_F > 0.0, ExcInternalError());
1729 * const NumberType fraction = (1.0 - (n_OS / det_F) )/(1.0 - n_OS);
1730 * return ( NumberType (std::pow(fraction, kappa_darcy))
1731 * * initial_darcy_flow_tensor);
1734 * Tensor<1, dim, NumberType>
1735 * get_body_force_FR_current() const
1737 * Tensor<1, dim, NumberType> body_force_FR_current;
1739 * if (gravity_term == true)
1741 * Tensor<1, dim, NumberType> gravity_vector;
1742 * gravity_vector[gravity_direction] = gravity_value;
1743 * body_force_FR_current = density_FR * gravity_vector;
1745 * return body_force_FR_current;
1752 * <a name="Quadraturepointhistory
"></a>
1753 * <h3>Quadrature point history</h3>
1754 * As seen in @ref step_18 "step-18
", the <code> PointHistory </code> class offers a method
1755 * for storing data at the quadrature points. Here each quadrature point
1756 * holds a pointer to a material description. Thus, different material models
1757 * can be used in different regions of the domain. Among other data, we
1758 * choose to store the ``extra" Kirchhoff stress @f$\boldsymbol{\tau}_E@f$ and
1759 * the dissipation
values @f$\mathcal{D}_p@f$ and @f$\mathcal{D}_v@f$.
1762 *
template <
int dim,
typename NumberType = Sacado::Fad::DFad<
double> >
1763 *
class PointHistory
1769 *
virtual ~PointHistory()
1772 *
void setup_lqp (
const Parameters::AllParameters ¶meters,
1775 *
if (parameters.mat_type ==
"Neo-Hooke")
1776 * solid_material.reset(
new NeoHooke<dim,NumberType>(parameters,time));
1777 *
else if (parameters.mat_type ==
"Ogden")
1778 * solid_material.reset(
new Ogden<dim,NumberType>(parameters,time));
1779 *
else if (parameters.mat_type ==
"visco-Ogden")
1780 * solid_material.reset(
new visco_Ogden<dim,NumberType>(parameters,time));
1784 * fluid_material.reset(
new Material_Darcy_Fluid<dim,NumberType>(parameters));
1790 *
return solid_material->get_tau_E(
F);
1796 *
return solid_material->get_Cauchy_E(
F);
1800 * get_converged_det_F() const
1802 *
return solid_material->get_converged_det_F();
1806 * update_end_timestep()
1808 * solid_material->update_end_timestep();
1814 * solid_material->update_internal_equilibrium(
F);
1818 * get_viscous_dissipation() const
1820 *
return solid_material->get_viscous_dissipation();
1827 *
return fluid_material->get_seepage_velocity_current(
F, grad_p_fluid);
1834 *
return fluid_material->get_porous_dissipation(
F, grad_p_fluid);
1839 *
const Parameters::AllParameters ¶meters)
const
1843 *
if (parameters.gravity_term ==
true)
1848 *
const NumberType overall_density_ref
1849 * = parameters.density_SR * parameters.solid_vol_frac
1850 * + parameters.density_FR
1851 * * (det_F_AD - parameters.solid_vol_frac);
1854 * gravity_vector[parameters.gravity_direction] = parameters.gravity_value;
1855 * body_force = overall_density_ref * gravity_vector;
1858 *
return body_force;
1861 * std::shared_ptr< Material_Hyperelastic<dim, NumberType> > solid_material;
1862 * std::shared_ptr< Material_Darcy_Fluid<dim, NumberType> > fluid_material;
1868 * <a name=
"Nonlinearporoviscoelasticsolid"></a>
1869 * <h3>Nonlinear poro-viscoelastic solid</h3>
1870 * The Solid
class is the central class as it represents the problem at hand:
1871 * the nonlinear poro-viscoelastic solid
1874 * template <int dim>
1878 * Solid(
const Parameters::AllParameters ¶meters);
1883 *
using ADNumberType = Sacado::Fad::DFad<double>;
1885 * std::ofstream outfile;
1886 * std::ofstream pointfile;
1888 *
struct PerTaskData_ASM;
1889 *
template<
typename NumberType =
double>
struct ScratchData_ASM;
1896 *
virtual void make_grid() = 0;
1900 * Define points
for post-processing
1903 *
virtual void define_tracked_vertices(std::vector<
Point<dim> > &tracked_vertices) = 0;
1907 * Set up the finite element system to be solved:
1914 * Extract sub-blocks from the global
matrix
1917 *
void determine_component_extractors();
1921 * Several
functions to
assemble the system and right hand side matrices
using multithreading.
1924 *
void assemble_system
1926 *
void assemble_system_one_cell
1928 * ScratchData_ASM<ADNumberType> &scratch,
1929 * PerTaskData_ASM &data)
const;
1930 *
void copy_local_to_global_system(
const PerTaskData_ASM &data);
1934 * Define boundary conditions
1937 *
virtual void make_constraints(
const int &it_nr);
1943 *
virtual double get_prescribed_fluid_flow
1947 * get_reaction_boundary_id_for_output ()
const = 0;
1948 *
virtual std::pair<types::boundary_id,types::boundary_id>
1949 * get_drained_boundary_id_for_output ()
const = 0;
1950 *
virtual std::vector<double> get_dirichlet_load
1952 *
const int &direction)
const = 0;
1956 * Create and update the quadrature points.
1963 * Solve non-linear system
using a Newton-Raphson scheme
1970 * Solve the linearized equations
using a direct solver
1977 * Retrieve the solution
1985 * Store the converged
values of the
internal variables at the
end of each timestep
1988 *
void update_end_timestep();
1992 * Post-processing and writing data to files
1995 *
void output_results_to_vtu(
const unsigned int timestep,
1996 *
const double current_time,
1998 *
void output_results_to_plot(
const unsigned int timestep,
1999 *
const double current_time,
2001 * std::vector<
Point<dim> > &tracked_vertices,
2002 * std::ofstream &pointfile)
const;
2006 * Headers and footer
for the output files
2009 *
void print_console_file_header( std::ofstream &outfile)
const;
2010 *
void print_plot_file_header(std::vector<
Point<dim> > &tracked_vertices,
2011 * std::ofstream &pointfile)
const;
2012 *
void print_console_file_footer(std::ofstream &outfile)
const;
2013 *
void print_plot_file_footer( std::ofstream &pointfile)
const;
2027 *
A collection of the parameters used to describe the problem setup
2030 *
const Parameters::AllParameters ¶meters;
2041 * Keep track of the current time and the time spent evaluating certain
functions
2050 *
A storage
object for quadrature
point information.
2057 * Integers to store polynomial degree (needed
for output)
2060 *
const unsigned int degree_displ;
2061 *
const unsigned int degree_pore;
2065 * Declare an instance of
dealii FESystem class (finite element definition)
2079 * Integer to store DoFs per element (
this value will be used often)
2082 *
const unsigned int dofs_per_cell;
2086 * Declare an instance of
dealii Extractor objects used to retrieve information from the solution vectors
2087 * We will use
"u_fe" and
"p_fluid_fe"as subscript in
operator [] expressions on
FEValues and
FEFaceValues
2088 * objects to
extract the components of the displacement vector and fluid pressure, respectively.
2096 * Description of how the block-system is arranged. There are 3 blocks:
2097 * 0 - vector DOF displacements u
2098 * 1 - scalar DOF fluid pressure p_fluid
2101 *
static const unsigned int n_blocks = 2;
2102 *
static const unsigned int n_components = dim+1;
2103 *
static const unsigned int first_u_component = 0;
2104 *
static const unsigned int p_fluid_component = dim;
2127 * std::vector<unsigned int> block_component;
2134 * std::vector<IndexSet> all_locally_owned_dofs;
2137 * std::vector<IndexSet> locally_owned_partitioning;
2138 * std::vector<IndexSet> locally_relevant_partitioning;
2140 * std::vector<types::global_dof_index> dofs_per_block;
2141 * std::vector<types::global_dof_index> element_indices_u;
2142 * std::vector<types::global_dof_index> element_indices_p_fluid;
2146 * Declare an instance of
dealii QGauss class (The Gauss-Legendre family of quadrature rules
for numerical integration)
2147 * Gauss Points in element, with n quadrature points (in each space direction <dim> )
2153 * Gauss Points on element faces (used
for definition of BCs)
2156 *
const QGauss<dim - 1> qf_face;
2159 * Integer to store num GPs per element (
this value will be used often)
2162 *
const unsigned int n_q_points;
2165 * Integer to store num GPs per face (
this value will be used often)
2168 *
const unsigned int n_q_points_f;
2172 * Declare an instance of
dealii AffineConstraints class (linear constraints on DoFs due to hanging nodes or BCs)
2179 * Declare an instance of
dealii classes necessary
for FE system
set-up and assembly
2187 * Right hand side vector of forces
2193 * Total displacement
values + pressure (accumulated solution to FE system)
2200 * Non-block system
for the direct solver. We will
copy the block system into these to solve the linearized system of equations.
2208 * We define variables to store norms and update norms and normalisation factors.
2215 *
norm(1.0), u(1.0), p_fluid(1.0)
2224 *
void normalise(
const Errors &rhs)
2226 *
if (rhs.norm != 0.0)
2230 *
if (rhs.p_fluid != 0.0)
2231 * p_fluid /= rhs.p_fluid;
2234 *
double norm, u, p_fluid;
2239 * Declare several instances of the
"Error" structure
2242 * Errors error_residual, error_residual_0, error_residual_norm, error_update,
2243 * error_update_0, error_update_norm;
2247 * Methods to calculate error measures
2250 *
void get_error_residual(Errors &error_residual_OUT);
2251 *
void get_error_update
2253 * Errors &error_update_OUT);
2257 * Print information to screen
2260 *
void print_conv_header();
2261 *
void print_conv_footer();
2266 * modifying the input variables inside the
functions will change them outside the function.
2274 * <a name=
"ImplementationofthecodeSolidcodeclass"></a>
2275 * <h3>Implementation of the <code>Solid</code>
class</h3>
2277 * <a name=
"Publicinterface"></a>
2278 * <h4>Public interface</h4>
2279 * We initialise the Solid
class using data extracted from the parameter file.
2282 *
template <
int dim>
2283 * Solid<dim>::Solid(
const Parameters::AllParameters ¶meters)
2285 * mpi_communicator(MPI_COMM_WORLD),
2289 * parameters(parameters),
2291 * time(parameters.end_time, parameters.delta_t),
2292 * timerconsole( mpi_communicator,
2296 * timerfile( mpi_communicator,
2300 * degree_displ(parameters.poly_degree_displ),
2301 * degree_pore(parameters.poly_degree_pore),
2302 * fe(
FE_Q<dim>(parameters.poly_degree_displ), dim,
2303 *
FE_Q<dim>(parameters.poly_degree_pore), 1 ),
2305 * dofs_per_cell (fe.dofs_per_cell),
2306 * u_fe(first_u_component),
2307 * p_fluid_fe(p_fluid_component),
2308 * x_displacement(first_u_component),
2309 * y_displacement(first_u_component+1),
2310 * z_displacement(first_u_component+2),
2311 * pressure(p_fluid_component),
2313 * qf_cell(parameters.quad_order),
2314 * qf_face(parameters.quad_order),
2315 * n_q_points (qf_cell.size()),
2316 * n_q_points_f (qf_face.size())
2318 *
Assert(dim==3,
ExcMessage(
"This problem only works in 3 space dimensions."));
2319 * determine_component_extractors();
2324 * The
class destructor simply clears the data held by the DOFHandler
2327 *
template <
int dim>
2328 * Solid<dim>::~Solid()
2330 * dof_handler_ref.clear();
2335 * Runs the 3D solid problem
2338 *
template <
int dim>
2343 * The current solution increment is defined as a block vector to reflect the structure
2344 * of the PDE system, with multiple solution components
2356 * outfile.open(
"console-output.sol");
2357 * print_console_file_header(outfile);
2369 * Assign DOFs and create the stiffness and right-hand-side force vector
2372 * system_setup(solution_delta);
2376 * Define points
for post-processing
2379 * std::vector<Point<dim> > tracked_vertices (2);
2380 * define_tracked_vertices(tracked_vertices);
2381 * std::vector<Point<dim>> reaction_force;
2385 * pointfile.open(
"data-for-gnuplot.sol");
2386 * print_plot_file_header(tracked_vertices, pointfile);
2391 * Print results to output file
2394 *
if (parameters.outfiles_requested ==
"true")
2396 * output_results_to_vtu(time.get_timestep(),
2397 * time.get_current(),
2401 * output_results_to_plot(time.get_timestep(),
2402 * time.get_current(),
2409 * Increment time step (=load step)
2410 * NOTE: In solving the quasi-
static problem, the time becomes a loading parameter,
2411 * i.e. we increase the loading linearly with time, making the two concepts interchangeable.
2414 * time.increment_time();
2418 * Print information on screen
2421 * pcout <<
"\nSolver:";
2422 * pcout <<
"\n CST = make constraints";
2423 * pcout <<
"\n ASM_SYS = assemble system";
2424 * pcout <<
"\n SLV = linear solver \n";
2428 * Print information on file
2431 * outfile <<
"\nSolver:";
2432 * outfile <<
"\n CST = make constraints";
2433 * outfile <<
"\n ASM_SYS = assemble system";
2434 * outfile <<
"\n SLV = linear solver \n";
2436 *
while ( (time.get_end() - time.get_current()) > -1.0*parameters.tol_u )
2440 * Initialize the current solution increment to
zero
2443 * solution_delta = 0.0;
2447 * Solve the non-linear system
using a Newton-Rapshon scheme
2450 * solve_nonlinear_timestep(solution_delta);
2454 * Add the computed solution increment to total solution
2457 * solution_n += solution_delta;
2464 * update_end_timestep();
2471 *
if (( (time.get_timestep()%parameters.timestep_output) == 0 )
2472 * && (parameters.outfiles_requested ==
"true") )
2474 * output_results_to_vtu(time.get_timestep(),
2475 * time.get_current(),
2479 * output_results_to_plot(time.get_timestep(),
2480 * time.get_current(),
2487 * Increment the time step (=load step)
2490 * time.increment_time();
2495 * Print the footers and close files
2500 * print_plot_file_footer(pointfile);
2501 * pointfile.close ();
2502 * print_console_file_footer(outfile);
2506 * NOTE: ideally, we should close the outfile here [ >> outfile.close (); ]
2507 * But
if we
do, then the timer output will not be printed. That is why we leave it open.
2516 * <a name=
"Privateinterface"></a>
2517 * <h4>Private interface</h4>
2518 * We define the structures needed
for parallelization with Threading Building Blocks (TBB)
2519 * Tangent
matrix and right-hand side force vector assembly structures.
2520 * PerTaskData_ASM stores local contributions
2523 *
template <
int dim>
2524 *
struct Solid<dim>::PerTaskData_ASM
2528 * std::vector<types::global_dof_index> local_dof_indices;
2530 * PerTaskData_ASM(
const unsigned int dofs_per_cell)
2533 * cell_rhs(dofs_per_cell),
2534 * local_dof_indices(dofs_per_cell)
2546 * ScratchData_ASM stores larger objects used during the assembly
2549 *
template <
int dim>
2550 *
template <
typename NumberType>
2551 *
struct Solid<dim>::ScratchData_ASM
2557 * Integration helper
2568 * std::vector<NumberType> local_dof_values;
2569 * std::vector<Tensor<2, dim, NumberType> > solution_grads_u_total;
2570 * std::vector<NumberType> solution_values_p_fluid_total;
2571 * std::vector<Tensor<1, dim, NumberType> > solution_grads_p_fluid_total;
2572 * std::vector<Tensor<1, dim, NumberType> > solution_grads_face_p_fluid_total;
2579 * std::vector<std::vector<Tensor<1,dim>>> Nx;
2580 * std::vector<std::vector<double>> Nx_p_fluid;
2586 * std::vector<std::vector<Tensor<2,dim, NumberType>>> grad_Nx;
2587 * std::vector<std::vector<SymmetricTensor<2,dim, NumberType>>> symm_grad_Nx;
2588 * std::vector<std::vector<Tensor<1,dim, NumberType>>> grad_Nx_p_fluid;
2595 * solution_total (solution_total),
2596 * fe_values_ref(fe_cell, qf_cell, uf_cell),
2597 * fe_face_values_ref(fe_cell, qf_face, uf_face),
2598 * local_dof_values(fe_cell.dofs_per_cell),
2599 * solution_grads_u_total(qf_cell.size()),
2600 * solution_values_p_fluid_total(qf_cell.size()),
2601 * solution_grads_p_fluid_total(qf_cell.size()),
2602 * solution_grads_face_p_fluid_total(qf_face.size()),
2603 * Nx(qf_cell.size(), std::vector<
Tensor<1,dim>>(fe_cell.dofs_per_cell)),
2604 * Nx_p_fluid(qf_cell.size(), std::vector<double>(fe_cell.dofs_per_cell)),
2610 * ScratchData_ASM(
const ScratchData_ASM &rhs)
2612 * solution_total (rhs.solution_total),
2613 * fe_values_ref(rhs.fe_values_ref.get_fe(),
2614 * rhs.fe_values_ref.get_quadrature(),
2615 * rhs.fe_values_ref.get_update_flags()),
2616 * fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
2617 * rhs.fe_face_values_ref.get_quadrature(),
2618 * rhs.fe_face_values_ref.get_update_flags()),
2619 * local_dof_values(rhs.local_dof_values),
2620 * solution_grads_u_total(rhs.solution_grads_u_total),
2621 * solution_values_p_fluid_total(rhs.solution_values_p_fluid_total),
2622 * solution_grads_p_fluid_total(rhs.solution_grads_p_fluid_total),
2623 * solution_grads_face_p_fluid_total(rhs.solution_grads_face_p_fluid_total),
2625 * Nx_p_fluid(rhs.Nx_p_fluid),
2626 * grad_Nx(rhs.grad_Nx),
2627 * symm_grad_Nx(rhs.symm_grad_Nx),
2628 * grad_Nx_p_fluid(rhs.grad_Nx_p_fluid)
2633 *
const unsigned int n_q_points = Nx_p_fluid.size();
2634 *
const unsigned int n_dofs_per_cell = Nx_p_fluid[0].size();
2638 *
for (
unsigned int k = 0; k < n_dofs_per_cell; ++k)
2640 * local_dof_values[k] = 0.0;
2651 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2657 * solution_grads_u_total[q_point] = 0.0;
2658 * solution_values_p_fluid_total[q_point] = 0.0;
2659 * solution_grads_p_fluid_total[q_point] = 0.0;
2661 *
for (
unsigned int k = 0; k < n_dofs_per_cell; ++k)
2663 * Nx[q_point][k] = 0.0;
2664 * Nx_p_fluid[q_point][k] = 0.0;
2665 * grad_Nx[q_point][k] = 0.0;
2666 * symm_grad_Nx[q_point][k] = 0.0;
2667 * grad_Nx_p_fluid[q_point][k] = 0.0;
2671 *
const unsigned int n_f_q_points = solution_grads_face_p_fluid_total.size();
2674 *
for (
unsigned int f_q_point = 0; f_q_point < n_f_q_points; ++f_q_point)
2675 * solution_grads_face_p_fluid_total[f_q_point] = 0.0;
2681 * Define the boundary conditions on the mesh
2684 *
template <
int dim>
2685 *
void Solid<dim>::make_constraints(
const int &it_nr_IN)
2687 * pcout <<
" CST " << std::flush;
2688 * outfile <<
" CST " << std::flush;
2690 *
if (it_nr_IN > 1)
return;
2692 *
const bool apply_dirichlet_bc = (it_nr_IN == 0);
2694 *
if (apply_dirichlet_bc)
2696 * constraints.clear();
2697 * make_dirichlet_constraints(constraints);
2701 *
for (
unsigned int i=0; i<dof_handler_ref.n_dofs(); ++i)
2702 *
if (constraints.is_inhomogeneously_constrained(i) ==
true)
2703 * constraints.set_inhomogeneity(i,0.0);
2705 * constraints.close();
2710 * Set-up the FE system
2713 *
template <
int dim>
2721 * Determine number of components per block
2724 * std::vector<unsigned int> block_component(n_components, u_block);
2725 * block_component[p_fluid_component] = p_fluid_block;
2729 * The DOF handler is initialised and we renumber the grid in an efficient manner.
2732 * dof_handler_ref.distribute_dofs(fe);
2738 * Count the number of DoFs in each block
2741 * dofs_per_block.clear();
2747 * Setup the sparsity pattern and tangent
matrix
2751 * std::vector<IndexSet> all_locally_relevant_dofs
2754 * locally_owned_dofs.
clear();
2755 * locally_owned_partitioning.clear();
2759 * locally_relevant_dofs.
clear();
2760 * locally_relevant_partitioning.clear();
2764 * locally_owned_partitioning.reserve(
n_blocks);
2765 * locally_relevant_partitioning.reserve(
n_blocks);
2770 * = std::accumulate(dofs_per_block.begin(),
2771 * std::next(dofs_per_block.begin(),
b), 0);
2773 * = std::accumulate(dofs_per_block.begin(),
2774 * std::next(dofs_per_block.begin(),
b+1), 0);
2775 * locally_owned_partitioning.push_back(locally_owned_dofs.
get_view(idx_begin, idx_end));
2776 * locally_relevant_partitioning.push_back(locally_relevant_dofs.
get_view(idx_begin, idx_end));
2781 * Print information on screen
2784 * pcout <<
"\nTriangulation:\n"
2785 * <<
" Number of active cells: "
2787 * <<
" (by partition:";
2789 * pcout << (p==0 ?
' ' :
'+')
2793 * pcout <<
" Number of degrees of freedom: "
2794 * << dof_handler_ref.n_dofs()
2795 * <<
" (by partition:";
2797 * pcout << (p==0 ?
' ' :
'+')
2801 * pcout <<
" Number of degrees of freedom per block: "
2802 * <<
"[n_u, n_p_fluid] = ["
2803 * << dofs_per_block[u_block]
2805 * << dofs_per_block[p_fluid_block]
2811 * Print information to file
2814 * outfile <<
"\nTriangulation:\n"
2815 * <<
" Number of active cells: "
2817 * <<
" (by partition:";
2819 * outfile << (p==0 ?
' ' :
'+')
2823 * outfile <<
" Number of degrees of freedom: "
2824 * << dof_handler_ref.n_dofs()
2825 * <<
" (by partition:";
2827 * outfile << (p==0 ?
' ' :
'+')
2831 * outfile <<
" Number of degrees of freedom per block: "
2832 * <<
"[n_u, n_p_fluid] = ["
2833 * << dofs_per_block[u_block]
2835 * << dofs_per_block[p_fluid_block]
2841 * We optimise the sparsity pattern to reflect
this structure and prevent
2842 * unnecessary data creation
for the right-
diagonal block components.
2846 *
for (
unsigned int ii = 0; ii < n_components; ++ii)
2847 *
for (
unsigned int jj = 0; jj < n_components; ++jj)
2851 * Identify
"zero" matrix components of FE-system (The two components
do not couple)
2854 *
if (((ii == p_fluid_component) && (jj < p_fluid_component))
2855 * || ((ii < p_fluid_component) && (jj == p_fluid_component)) )
2860 * The rest of components
always couple
2867 * mpi_communicator);
2875 * Reinitialize the (sparse) tangent
matrix with the given sparsity pattern.
2878 * tangent_matrix.reinit (bsp);
2882 * Initialize the right hand side and solution vectors with number of DoFs
2885 * system_rhs.reinit(locally_owned_partitioning, mpi_communicator);
2886 * solution_n.reinit(locally_owned_partitioning, mpi_communicator);
2887 * solution_delta_OUT.reinit(locally_owned_partitioning, mpi_communicator);
2895 * mpi_communicator);
2899 * tangent_matrix_nb.reinit (sp);
2900 * system_rhs_nb.
reinit(locally_owned_dofs, mpi_communicator);
2904 * Set up the quadrature
point history
2915 * Component extractors: used to
extract sub-blocks from the global
matrix
2916 * Description of which local element DOFs are attached to which block component
2919 *
template <
int dim>
2920 *
void Solid<dim>::determine_component_extractors()
2922 * element_indices_u.clear();
2923 * element_indices_p_fluid.clear();
2925 *
for (
unsigned int k = 0; k < fe.dofs_per_cell; ++k)
2927 *
const unsigned int k_group = fe.system_to_base_index(k).first.first;
2928 *
if (k_group == u_block)
2929 * element_indices_u.push_back(k);
2930 *
else if (k_group == p_fluid_block)
2931 * element_indices_p_fluid.push_back(k);
2941 * Set-up quadrature
point history (QPH) data objects
2944 *
template <
int dim>
2945 *
void Solid<dim>::setup_qph()
2947 * pcout <<
"\nSetting up quadrature point data..." << std::endl;
2948 * outfile <<
"\nSetting up quadrature point data..." << std::endl;
2952 * Create QPH data objects.
2955 * quadrature_point_history.initialize(
triangulation.begin_active(),
2960 * Setup the
initial quadrature
point data
using the info stored in parameters
2965 * dof_handler_ref.begin_active()),
2967 * dof_handler_ref.end());
2968 *
for (; cell!=endc; ++cell)
2973 *
const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >
2974 * lqph = quadrature_point_history.get_data(cell);
2977 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2978 * lqph[q_point]->setup_lqp(parameters, time);
2984 * Solve the non-linear system
using a Newton-Raphson scheme
2987 *
template <
int dim>
2992 * Print the load step
2995 * pcout << std::endl
2997 * << time.get_timestep()
2999 * << time.get_current()
3002 * outfile << std::endl
3004 * << time.get_timestep()
3006 * << time.get_current()
3012 * Declare newton_update vector (solution of a Newton iteration),
3013 * which must have as many positions as global DoFs.
3017 * (locally_owned_partitioning, mpi_communicator);
3021 * Reset the error storage objects
3024 * error_residual.reset();
3025 * error_residual_0.reset();
3026 * error_residual_norm.reset();
3027 * error_update.reset();
3028 * error_update_0.reset();
3029 * error_update_norm.reset();
3031 * print_conv_header();
3035 * Declare and initialize iterator
for the Newton-Raphson algorithm steps
3038 *
unsigned int newton_iteration = 0;
3042 * Iterate until error is below tolerance or
max number iterations are reached
3045 *
while(newton_iteration < parameters.max_iterations_NR)
3047 * pcout <<
" " << std::setw(2) << newton_iteration <<
" " << std::flush;
3048 * outfile <<
" " << std::setw(2) << newton_iteration <<
" " << std::flush;
3052 * Initialize global stiffness
matrix and global force vector to
zero
3055 * tangent_matrix = 0.0;
3058 * tangent_matrix_nb = 0.0;
3059 * system_rhs_nb = 0.0;
3063 * Apply boundary conditions
3066 * make_constraints(newton_iteration);
3067 * assemble_system(solution_delta_OUT);
3071 * Compute the rhs residual (error between external and
internal forces in FE system)
3074 * get_error_residual(error_residual);
3078 * error_residual in
first iteration is stored to normalize posterior error measures
3081 *
if (newton_iteration == 0)
3082 * error_residual_0 = error_residual;
3086 * Determine the normalised residual error
3089 * error_residual_norm = error_residual;
3090 * error_residual_norm.normalise(error_residual_0);
3094 * If both errors are below the tolerances, exit the
loop.
3095 * We need to check the residual vector directly
for convergence
3096 * in the load steps where no external forces or displacements are imposed.
3099 *
if ( ((newton_iteration > 0)
3100 * && (error_update_norm.u <= parameters.tol_u)
3101 * && (error_update_norm.p_fluid <= parameters.tol_p_fluid)
3102 * && (error_residual_norm.u <= parameters.tol_f)
3103 * && (error_residual_norm.p_fluid <= parameters.tol_f))
3104 * || ( (newton_iteration > 0)
3105 * && system_rhs.l2_norm() <= parameters.tol_f) )
3107 * pcout <<
"\n ***** CONVERGED! ***** "
3108 * << system_rhs.l2_norm() <<
" "
3109 * <<
" " << error_residual_norm.norm
3110 * <<
" " << error_residual_norm.u
3111 * <<
" " << error_residual_norm.p_fluid
3112 * <<
" " << error_update_norm.norm
3113 * <<
" " << error_update_norm.u
3114 * <<
" " << error_update_norm.p_fluid
3115 * <<
" " << std::endl;
3116 * outfile <<
"\n ***** CONVERGED! ***** "
3117 * << system_rhs.l2_norm() <<
" "
3118 * <<
" " << error_residual_norm.norm
3119 * <<
" " << error_residual_norm.u
3120 * <<
" " << error_residual_norm.p_fluid
3121 * <<
" " << error_update_norm.norm
3122 * <<
" " << error_update_norm.u
3123 * <<
" " << error_update_norm.p_fluid
3124 * <<
" " << std::endl;
3125 * print_conv_footer();
3132 * Solve the linearized system
3135 * solve_linear_system(newton_update);
3136 * constraints.distribute(newton_update);
3140 * Compute the displacement error
3143 * get_error_update(newton_update, error_update);
3147 * error_update in
first iteration is stored to normalize posterior error measures
3150 *
if (newton_iteration == 0)
3151 * error_update_0 = error_update;
3155 * Determine the normalised Newton update error
3158 * error_update_norm = error_update;
3159 * error_update_norm.normalise(error_update_0);
3163 * Determine the normalised residual error
3166 * error_residual_norm = error_residual;
3167 * error_residual_norm.normalise(error_residual_0);
3174 * pcout <<
" | " << std::fixed << std::setprecision(3)
3175 * << std::setw(7) << std::scientific
3176 * << system_rhs.l2_norm()
3177 * <<
" " << error_residual_norm.norm
3178 * <<
" " << error_residual_norm.u
3179 * <<
" " << error_residual_norm.p_fluid
3180 * <<
" " << error_update_norm.norm
3181 * <<
" " << error_update_norm.u
3182 * <<
" " << error_update_norm.p_fluid
3183 * <<
" " << std::endl;
3185 * outfile <<
" | " << std::fixed << std::setprecision(3)
3186 * << std::setw(7) << std::scientific
3187 * << system_rhs.l2_norm()
3188 * <<
" " << error_residual_norm.norm
3189 * <<
" " << error_residual_norm.u
3190 * <<
" " << error_residual_norm.p_fluid
3191 * <<
" " << error_update_norm.norm
3192 * <<
" " << error_update_norm.u
3193 * <<
" " << error_update_norm.p_fluid
3194 * <<
" " << std::endl;
3201 * solution_delta_OUT += newton_update;
3202 * newton_update = 0.0;
3203 * newton_iteration++;
3208 * If maximum allowed number of iterations
for Newton algorithm are reached, print non-convergence message and
abort program
3211 *
AssertThrow (newton_iteration < parameters.max_iterations_NR,
ExcMessage(
"No convergence in nonlinear solver!"));
3216 * Prints the header
for convergence info on console
3219 *
template <
int dim>
3220 *
void Solid<dim>::print_conv_header()
3222 *
static const unsigned int l_width = 120;
3224 *
for (
unsigned int i = 0; i < l_width; ++i)
3230 * pcout << std::endl;
3231 * outfile << std::endl;
3233 * pcout <<
"\n SOLVER STEP | SYS_RES "
3234 * <<
"RES_NORM RES_U RES_P "
3235 * <<
"NU_NORM NU_U NU_P " << std::endl;
3236 * outfile <<
"\n SOLVER STEP | SYS_RES "
3237 * <<
"RES_NORM RES_U RES_P "
3238 * <<
"NU_NORM NU_U NU_P " << std::endl;
3240 *
for (
unsigned int i = 0; i < l_width; ++i)
3245 * pcout << std::endl << std::endl;
3246 * outfile << std::endl << std::endl;
3251 * Prints the footer
for convergence info on console
3254 *
template <
int dim>
3255 *
void Solid<dim>::print_conv_footer()
3257 *
static const unsigned int l_width = 120;
3259 *
for (
unsigned int i = 0; i < l_width; ++i)
3264 * pcout << std::endl << std::endl;
3265 * outfile << std::endl << std::endl;
3267 * pcout <<
"Relative errors:" << std::endl
3268 * <<
"Displacement: "
3269 * << error_update.u / error_update_0.u << std::endl
3270 * <<
"Force (displ): "
3271 * << error_residual.u / error_residual_0.u << std::endl
3272 * <<
"Pore pressure: "
3273 * << error_update.p_fluid / error_update_0.p_fluid << std::endl
3274 * <<
"Force (pore): "
3275 * << error_residual.p_fluid / error_residual_0.p_fluid << std::endl;
3276 * outfile <<
"Relative errors:" << std::endl
3277 * <<
"Displacement: "
3278 * << error_update.u / error_update_0.u << std::endl
3279 * <<
"Force (displ): "
3280 * << error_residual.u / error_residual_0.u << std::endl
3281 * <<
"Pore pressure: "
3282 * << error_update.p_fluid / error_update_0.p_fluid << std::endl
3283 * <<
"Force (pore): "
3284 * << error_residual.p_fluid / error_residual_0.p_fluid << std::endl;
3289 * Determine the
true residual error
for the problem
3292 *
template <
int dim>
3293 *
void Solid<dim>::get_error_residual(Errors &error_residual_OUT)
3296 * constraints.set_zero(error_res);
3298 * error_residual_OUT.norm = error_res.l2_norm();
3299 * error_residual_OUT.u = error_res.block(u_block).l2_norm();
3300 * error_residual_OUT.p_fluid = error_res.block(p_fluid_block).l2_norm();
3305 * Determine the
true Newton update error
for the problem
3308 *
template <
int dim>
3309 *
void Solid<dim>::get_error_update
3311 * Errors &error_update_OUT)
3314 * constraints.set_zero(error_ud);
3316 * error_update_OUT.norm = error_ud.l2_norm();
3317 * error_update_OUT.u = error_ud.block(u_block).l2_norm();
3318 * error_update_OUT.p_fluid = error_ud.block(p_fluid_block).l2_norm();
3323 * Compute the total solution, which is
valid at any Newton step. This is required as, to reduce
3324 * computational error, the total solution is only updated at the
end of the timestep.
3327 *
template <
int dim>
3333 * Cell interpolation -> Ghosted vector
3337 * solution_total (locally_owned_partitioning,
3338 * locally_relevant_partitioning,
3342 * solution_total = solution_n;
3343 * tmp = solution_delta_IN;
3344 * solution_total += tmp;
3345 *
return solution_total;
3350 * Compute elemental stiffness tensor and right-hand side force vector, and
assemble into global ones
3353 *
template <
int dim>
3358 * pcout <<
" ASM_SYS " << std::flush;
3359 * outfile <<
" ASM_SYS " << std::flush;
3365 * Info given to
FEValues and
FEFaceValues constructors, to indicate which data will be needed at each element.
3379 * Setup a
copy of the data structures required
for the process and pass them, along with the
3383 * PerTaskData_ASM per_task_data(dofs_per_cell);
3384 * ScratchData_ASM<ADNumberType> scratch_data(fe, qf_cell, uf_cell,
3390 * dof_handler_ref.begin_active()),
3392 * dof_handler_ref.end());
3393 *
for (; cell != endc; ++cell)
3398 * assemble_system_one_cell(cell, scratch_data, per_task_data);
3399 * copy_local_to_global_system(per_task_data);
3413 * Add the local elemental contribution to the global stiffness tensor
3414 * We
do it twice,
for the block and the non-block systems
3417 *
template <
int dim>
3418 *
void Solid<dim>::copy_local_to_global_system (
const PerTaskData_ASM &data)
3420 * constraints.distribute_local_to_global(data.cell_matrix,
3422 * data.local_dof_indices,
3426 * constraints.distribute_local_to_global(data.cell_matrix,
3428 * data.local_dof_indices,
3429 * tangent_matrix_nb,
3435 * Compute stiffness
matrix and corresponding rhs
for one element
3438 *
template <
int dim>
3439 *
void Solid<dim>::assemble_system_one_cell
3441 * ScratchData_ASM<ADNumberType> &scratch,
3442 * PerTaskData_ASM &data)
const
3448 * scratch.fe_values_ref.reinit(cell);
3449 * cell->get_dof_indices(data.local_dof_indices);
3453 * Setup automatic differentiation
3456 *
for (
unsigned int k = 0; k < dofs_per_cell; ++k)
3460 * Initialise the dofs
for the cell
using the current solution.
3463 * scratch.local_dof_values[k] = scratch.solution_total[data.local_dof_indices[k]];
3466 * Mark
this cell DoF as an independent variable
3469 * scratch.local_dof_values[k].diff(k, dofs_per_cell);
3474 * Update the quadrature
point solution
3475 * Compute the
values and
gradients of the solution in terms of the AD variables
3478 *
for (
unsigned int q = 0; q < n_q_points; ++q)
3480 *
for (
unsigned int k = 0; k < dofs_per_cell; ++k)
3482 *
const unsigned int k_group = fe.system_to_base_index(k).first.first;
3483 *
if (k_group == u_block)
3486 * scratch.fe_values_ref[u_fe].gradient(k, q);
3487 *
for (
unsigned int dd = 0; dd < dim; ++dd)
3489 *
for (
unsigned int ee = 0; ee < dim; ++ee)
3491 * scratch.solution_grads_u_total[q][dd][ee]
3492 * += scratch.local_dof_values[k] * Grad_Nx_u[dd][ee];
3496 *
else if (k_group == p_fluid_block)
3498 *
const double Nx_p = scratch.fe_values_ref[p_fluid_fe].value(k, q);
3500 * scratch.fe_values_ref[p_fluid_fe].gradient(k, q);
3502 * scratch.solution_values_p_fluid_total[q]
3503 * += scratch.local_dof_values[k] * Nx_p;
3504 *
for (
unsigned int dd = 0; dd < dim; ++dd)
3506 * scratch.solution_grads_p_fluid_total[q][dd]
3507 * += scratch.local_dof_values[k] * Grad_Nx_p[dd];
3517 * Set up pointer
"lgph" to the PointHistory
object of
this element
3520 *
const std::vector<std::shared_ptr<const PointHistory<dim, ADNumberType> > >
3521 * lqph = quadrature_point_history.get_data(cell);
3530 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3537 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
3539 *
const unsigned int i_group = fe.system_to_base_index(i).first.first;
3541 *
if (i_group == u_block)
3543 * scratch.Nx[q_point][i] =
3544 * scratch.fe_values_ref[u_fe].value(i, q_point);
3545 * scratch.grad_Nx[q_point][i] =
3546 * scratch.fe_values_ref[u_fe].gradient(i, q_point)*F_inv_AD;
3547 * scratch.symm_grad_Nx[q_point][i] =
3550 *
else if (i_group == p_fluid_block)
3552 * scratch.Nx_p_fluid[q_point][i] =
3553 * scratch.fe_values_ref[p_fluid_fe].value(i, q_point);
3554 * scratch.grad_Nx_p_fluid[q_point][i] =
3555 * scratch.fe_values_ref[p_fluid_fe].gradient(i, q_point)*F_inv_AD;
3564 * Assemble the stiffness
matrix and rhs vector
3567 * std::vector<ADNumberType> residual_ad (dofs_per_cell, ADNumberType(0.0));
3568 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3572 *
const ADNumberType det_F_AD =
determinant(F_AD);
3577 *
const ADNumberType p_fluid = scratch.solution_values_p_fluid_total[q_point];
3580 * PointHistory<dim, ADNumberType> *lqph_q_point_nc =
3581 *
const_cast<PointHistory<dim, ADNumberType>*
>(lqph[q_point].get());
3582 * lqph_q_point_nc->update_internal_equilibrium(F_AD);
3587 * Get some info from constitutive model of solid
3593 * tau_E = lqph[q_point]->get_tau_E(F_AD);
3595 * tau_fluid_vol *= -1.0 * p_fluid * det_F_AD;
3599 * Get some info from constitutive model of fluid
3602 *
const ADNumberType det_F_aux = lqph[q_point]->get_converged_det_F();
3605 * = lqph[q_point]->get_overall_body_force(F_AD, parameters);
3609 * Define some aliases to make the assembly process easier to follow
3612 *
const std::vector<Tensor<1,dim>> &Nu = scratch.Nx[q_point];
3613 *
const std::vector<SymmetricTensor<2, dim, ADNumberType>>
3614 * &symm_grad_Nu = scratch.symm_grad_Nx[q_point];
3615 *
const std::vector<double> &Np = scratch.Nx_p_fluid[q_point];
3616 *
const std::vector<Tensor<1, dim, ADNumberType> > &grad_Np
3617 * = scratch.grad_Nx_p_fluid[q_point];
3619 * = scratch.solution_grads_p_fluid_total[q_point]*F_inv_AD;
3620 *
const double JxW = scratch.fe_values_ref.JxW(q_point);
3622 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
3624 *
const unsigned int i_group = fe.system_to_base_index(i).first.first;
3626 *
if (i_group == u_block)
3628 * residual_ad[i] += symm_grad_Nu[i] * ( tau_E + tau_fluid_vol ) * JxW;
3629 * residual_ad[i] -= Nu[i] * overall_body_force * JxW;
3631 *
else if (i_group == p_fluid_block)
3634 * = lqph[q_point]->get_seepage_velocity_current(F_AD, grad_p);
3635 * residual_ad[i] += Np[i] * (det_F_AD - det_F_converged) * JxW;
3636 * residual_ad[i] -= time.get_delta_t() * grad_Np[i]
3637 * * seepage_vel_current * JxW;
3646 * Assemble the Neumann contribution (external force contribution).
3649 *
for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3651 *
if (cell->face(face)->at_boundary() ==
true)
3653 * scratch.fe_face_values_ref.reinit(cell, face);
3655 *
for (
unsigned int f_q_point = 0; f_q_point < n_q_points_f; ++f_q_point)
3658 * = scratch.fe_face_values_ref.normal_vector(f_q_point);
3660 * = scratch.fe_face_values_ref.quadrature_point(f_q_point);
3662 * = get_neumann_traction(cell->face(face)->boundary_id(), pt,
N);
3664 * = get_prescribed_fluid_flow(cell->face(face)->boundary_id(), pt);
3666 *
if ( (traction.
norm() < 1
e-12) && (
std::abs(flow) < 1
e-12) )
continue;
3668 *
const double JxW_f = scratch.fe_face_values_ref.JxW(f_q_point);
3670 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
3672 *
const unsigned int i_group = fe.system_to_base_index(i).first.first;
3674 *
if ((i_group == u_block) && (traction.
norm() > 1
e-12))
3676 *
const unsigned int component_i
3677 * = fe.system_to_component_index(i).first;
3679 * = scratch.fe_face_values_ref.shape_value(i, f_q_point);
3680 * residual_ad[i] -= (Nu_f * traction[component_i]) * JxW_f;
3682 *
if ((i_group == p_fluid_block) && (
std::abs(flow) > 1
e-12))
3685 * = scratch.fe_face_values_ref.shape_value(i, f_q_point);
3686 * residual_ad[i] -= (Nu_p * flow) * JxW_f;
3695 * Linearise the residual
3698 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
3700 *
const ADNumberType &R_i = residual_ad[i];
3702 * data.cell_rhs(i) -= R_i.val();
3703 *
for (
unsigned int j=0; j<dofs_per_cell; ++j)
3704 * data.cell_matrix(i,j) += R_i.fastAccessDx(j);
3713 *
template <
int dim>
3714 *
void Solid<dim>::update_end_timestep()
3718 * dof_handler_ref.begin_active()),
3720 * dof_handler_ref.end());
3721 *
for (; cell!=endc; ++cell)
3726 *
const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >
3727 * lqph = quadrature_point_history.get_data(cell);
3729 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3730 * lqph[q_point]->update_end_timestep();
3737 * Solve the linearized equations
3740 *
template <
int dim>
3746 * pcout <<
" SLV " << std::flush;
3747 * outfile <<
" SLV " << std::flush;
3750 * newton_update_nb.
reinit(locally_owned_dofs, mpi_communicator);
3753 * 1.0e-6 * system_rhs_nb.
l2_norm());
3755 * solver.solve(tangent_matrix_nb, newton_update_nb, system_rhs_nb);
3759 * Copy the non-block solution back to block system
3762 *
for (
unsigned int i=0; i<locally_owned_dofs.
n_elements(); ++i)
3766 * newton_update_OUT(idx_i) = newton_update_nb(idx_i);
3776 * Class to be able to output results correctly when
using Paraview
3779 *
template<
int dim,
class DH=DoFHandler<dim> >
3780 *
class FilteredDataOut :
public DataOut<dim, DH>
3783 * FilteredDataOut ()
3786 *
virtual ~FilteredDataOut() {}
3792 * cell = this->dofs->begin_active();
3793 *
while ((cell != this->dofs->end()) &&
3794 * (!cell->is_locally_owned()))
3802 *
if (old_cell != this->dofs->end())
3807 * (predicate,old_cell));
3814 *
template<
int dim,
class DH=DoFHandler<dim> >
3815 *
class FilteredDataOutFaces :
public DataOutFaces<dim,DH>
3818 * FilteredDataOutFaces ()
3821 *
virtual ~FilteredDataOutFaces() {}
3827 * cell = this->dofs->begin_active();
3828 *
while ((cell!=this->dofs->end()) && (!cell->is_locally_owned()))
3836 *
if (old_cell!=this->dofs->end())
3841 * (predicate,old_cell));
3850 * Class to compute
gradient of the pressure
3853 *
template <
int dim>
3857 * GradientPostprocessor (
const unsigned int p_fluid_component)
3861 * p_fluid_component (p_fluid_component)
3864 *
virtual ~GradientPostprocessor(){}
3867 * evaluate_vector_field
3869 * std::vector<Vector<double> > &computed_quantities)
const
3872 * computed_quantities.size());
3873 *
for (
unsigned int p=0; p<input_data.solution_gradients.size(); ++p)
3876 *
for (
unsigned int d=0;
d<dim; ++
d)
3877 * computed_quantities[p][
d]
3878 * = input_data.solution_gradients[p][p_fluid_component][
d];
3883 *
const unsigned int p_fluid_component;
3889 * Print results to
vtu file
3892 *
template <
int dim>
void Solid<dim>::output_results_to_vtu
3893 * (
const unsigned int timestep,
3894 *
const double current_time,
3898 * locally_relevant_partitioning,
3901 * solution_total = solution_IN;
3904 * std::vector<types::subdomain_id> partition_int(
triangulation.n_active_cells());
3905 * GradientPostprocessor<dim> gradient_postprocessor(p_fluid_component);
3909 * Declare local variables with number of stress components
3913 *
unsigned int num_comp_symm_tensor = 6;
3917 * Declare local vectors to store
values
3918 * OUTPUT AVERAGED ON ELEMENTS -------------------------------------------
3921 * std::vector<Vector<double>>cauchy_stresses_total_elements
3922 * (num_comp_symm_tensor,
3924 * std::vector<Vector<double>>cauchy_stresses_E_elements
3925 * (num_comp_symm_tensor,
3927 * std::vector<Vector<double>>stretches_elements
3930 * std::vector<Vector<double>>seepage_velocity_elements
3942 * OUTPUT AVERAGED ON NODES ----------------------------------------------
3943 * We need to create a
new FE space with a single dof per node to avoid
3944 * duplication of the output on nodes
for our problem with dim+1 dofs.
3949 * vertex_handler_ref.distribute_dofs(fe_vertex);
3955 * (vertex_handler_ref.n_dofs());
3957 * (vertex_handler_ref.n_dofs());
3959 * std::vector<Vector<double>>cauchy_stresses_total_vertex_mpi
3960 * (num_comp_symm_tensor,
3962 * std::vector<Vector<double>>sum_cauchy_stresses_total_vertex
3963 * (num_comp_symm_tensor,
3965 * std::vector<Vector<double>>cauchy_stresses_E_vertex_mpi
3966 * (num_comp_symm_tensor,
3968 * std::vector<Vector<double>>sum_cauchy_stresses_E_vertex
3969 * (num_comp_symm_tensor,
3971 * std::vector<Vector<double>>stretches_vertex_mpi
3974 * std::vector<Vector<double>>sum_stretches_vertex
3977 *
Vector<double> porous_dissipation_vertex_mpi(vertex_handler_ref.n_dofs());
3978 *
Vector<double> sum_porous_dissipation_vertex(vertex_handler_ref.n_dofs());
3979 *
Vector<double> viscous_dissipation_vertex_mpi(vertex_handler_ref.n_dofs());
3980 *
Vector<double> sum_viscous_dissipation_vertex(vertex_handler_ref.n_dofs());
3981 *
Vector<double> solid_vol_fraction_vertex_mpi(vertex_handler_ref.n_dofs());
3982 *
Vector<double> sum_solid_vol_fraction_vertex(vertex_handler_ref.n_dofs());
3986 * We need to create a
new FE space with a dim dof per node to
3987 * be able to ouput data on nodes in vector form
3992 * vertex_vec_handler_ref.distribute_dofs(fe_vertex_vec);
3997 *
Vector<double> seepage_velocity_vertex_vec_mpi(vertex_vec_handler_ref.n_dofs());
3998 *
Vector<double> sum_seepage_velocity_vertex_vec(vertex_vec_handler_ref.n_dofs());
3999 *
Vector<double> counter_on_vertices_vec_mpi(vertex_vec_handler_ref.n_dofs());
4000 *
Vector<double> sum_counter_on_vertices_vec(vertex_vec_handler_ref.n_dofs());
4003 * -----------------------------------------------------------------------
4007 * Declare and initialize local unit vectors (to construct tensor basis)
4010 * std::vector<Tensor<1,dim>> basis_vectors (dim,
Tensor<1,dim>() );
4011 *
for (
unsigned int i=0; i<dim; ++i)
4012 * basis_vectors[i][i] = 1;
4016 * Declare an instance of the material
class object
4019 *
if (parameters.mat_type ==
"Neo-Hooke")
4020 * NeoHooke<dim,ADNumberType> material(parameters,time);
4021 *
else if (parameters.mat_type ==
"Ogden")
4022 * Ogden<dim,ADNumberType> material(parameters,time);
4023 *
else if (parameters.mat_type ==
"visco-Ogden")
4024 * visco_Ogden <dim,ADNumberType>material(parameters,time);
4030 * Define a local instance of
FEValues to compute updated
values required
4031 * to calculate stresses
4040 * Iterate through elements (cells) and Gauss Points
4045 * dof_handler_ref.begin_active()),
4047 * dof_handler_ref.end()),
4049 * vertex_handler_ref.begin_active()),
4051 * vertex_vec_handler_ref.begin_active());
4057 *
for (; cell!=endc; ++cell, ++cell_v, ++cell_v_vec)
4063 *
static_cast<int>(cell->material_id());
4065 * fe_values_ref.reinit(cell);
4067 * std::vector<Tensor<2,dim>> solution_grads_u(n_q_points);
4068 * fe_values_ref[u_fe].get_function_gradients(solution_total,
4069 * solution_grads_u);
4071 * std::vector<double> solution_values_p_fluid_total(n_q_points);
4072 * fe_values_ref[p_fluid_fe].get_function_values(solution_total,
4073 * solution_values_p_fluid_total);
4075 * std::vector<Tensor<1,dim>> solution_grads_p_fluid_AD (n_q_points);
4076 * fe_values_ref[p_fluid_fe].get_function_gradients(solution_total,
4077 * solution_grads_p_fluid_AD);
4084 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
4091 *
const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4092 * lqph = quadrature_point_history.get_data(cell);
4095 *
const double p_fluid = solution_values_p_fluid_total[q_point];
4106 * lqph[q_point]->get_Cauchy_E(F_AD);
4108 *
for (
unsigned int i=0; i<dim; ++i)
4109 *
for (
unsigned int j=0; j<dim; ++j)
4113 * sigma_fluid_vol *= -p_fluid;
4121 *
const double solid_vol_fraction = (parameters.solid_vol_frac)/det_F;
4125 * Green-Lagrange strain
4137 * solution_grads_p_fluid_AD[q_point]*F_inv;
4139 * lqph[q_point]->get_seepage_velocity_current(F_AD, grad_p_fluid_AD);
4146 *
const double porous_dissipation =
4147 * lqph[q_point]->get_porous_dissipation(F_AD, grad_p_fluid_AD);
4148 *
const double viscous_dissipation =
4149 * lqph[q_point]->get_viscous_dissipation();
4153 * OUTPUT AVERAGED ON ELEMENTS -------------------------------------------
4154 * Both average on elements and on nodes is NOT weighted with the
4156 * integration
point to the average. Ideally, it should be weighted,
4157 * but I haven
't invested time in getting it to work properly.
4160 * if (parameters.outtype == "elements")
4162 * for (unsigned int j=0; j<dim; ++j)
4164 * cauchy_stresses_total_elements[j](cell->active_cell_index())
4165 * += ((sigma*basis_vectors[j])*basis_vectors[j])/n_q_points;
4166 * cauchy_stresses_E_elements[j](cell->active_cell_index())
4167 * += ((sigma_E*basis_vectors[j])*basis_vectors[j])/n_q_points;
4168 * stretches_elements[j](cell->active_cell_index())
4169 * += std::sqrt(1.0+2.0*Tensor<0,dim,double>(E_strain[j][j]))
4171 * seepage_velocity_elements[j](cell->active_cell_index())
4172 * += Tensor<0,dim,double>(seepage_vel_AD[j])/n_q_points;
4175 * porous_dissipation_elements(cell->active_cell_index())
4176 * += porous_dissipation/n_q_points;
4177 * viscous_dissipation_elements(cell->active_cell_index())
4178 * += viscous_dissipation/n_q_points;
4179 * solid_vol_fraction_elements(cell->active_cell_index())
4180 * += solid_vol_fraction/n_q_points;
4182 * cauchy_stresses_total_elements[3](cell->active_cell_index())
4183 * += ((sigma*basis_vectors[0])*basis_vectors[1])/n_q_points; //sig_xy
4184 * cauchy_stresses_total_elements[4](cell->active_cell_index())
4185 * += ((sigma*basis_vectors[0])*basis_vectors[2])/n_q_points;//sig_xz
4186 * cauchy_stresses_total_elements[5](cell->active_cell_index())
4187 * += ((sigma*basis_vectors[1])*basis_vectors[2])/n_q_points;//sig_yz
4189 * cauchy_stresses_E_elements[3](cell->active_cell_index())
4190 * += ((sigma_E*basis_vectors[0])* basis_vectors[1])/n_q_points; //sig_xy
4191 * cauchy_stresses_E_elements[4](cell->active_cell_index())
4192 * += ((sigma_E*basis_vectors[0])* basis_vectors[2])/n_q_points;//sig_xz
4193 * cauchy_stresses_E_elements[5](cell->active_cell_index())
4194 * += ((sigma_E*basis_vectors[1])* basis_vectors[2])/n_q_points;//sig_yz
4199 * OUTPUT AVERAGED ON NODES -------------------------------------------
4202 * else if (parameters.outtype == "nodes")
4204 * for (unsigned int v=0; v<(GeometryInfo<dim>::vertices_per_cell); ++v)
4206 * types::global_dof_index local_vertex_indices =
4207 * cell_v->vertex_dof_index(v, 0);
4208 * counter_on_vertices_mpi(local_vertex_indices) += 1;
4209 * for (unsigned int k=0; k<dim; ++k)
4211 * cauchy_stresses_total_vertex_mpi[k](local_vertex_indices)
4212 * += (sigma*basis_vectors[k])*basis_vectors[k];
4213 * cauchy_stresses_E_vertex_mpi[k](local_vertex_indices)
4214 * += (sigma_E*basis_vectors[k])*basis_vectors[k];
4215 * stretches_vertex_mpi[k](local_vertex_indices)
4216 * += std::sqrt(1.0+2.0*Tensor<0,dim,double>(E_strain[k][k]));
4218 * types::global_dof_index local_vertex_vec_indices =
4219 * cell_v_vec->vertex_dof_index(v, k);
4220 * counter_on_vertices_vec_mpi(local_vertex_vec_indices) += 1;
4221 * seepage_velocity_vertex_vec_mpi(local_vertex_vec_indices)
4222 * += Tensor<0,dim,double>(seepage_vel_AD[k]);
4225 * porous_dissipation_vertex_mpi(local_vertex_indices)
4226 * += porous_dissipation;
4227 * viscous_dissipation_vertex_mpi(local_vertex_indices)
4228 * += viscous_dissipation;
4229 * solid_vol_fraction_vertex_mpi(local_vertex_indices)
4230 * += solid_vol_fraction;
4232 * cauchy_stresses_total_vertex_mpi[3](local_vertex_indices)
4233 * += (sigma*basis_vectors[0])*basis_vectors[1]; //sig_xy
4234 * cauchy_stresses_total_vertex_mpi[4](local_vertex_indices)
4235 * += (sigma*basis_vectors[0])*basis_vectors[2];//sig_xz
4236 * cauchy_stresses_total_vertex_mpi[5](local_vertex_indices)
4237 * += (sigma*basis_vectors[1])*basis_vectors[2]; //sig_yz
4239 * cauchy_stresses_E_vertex_mpi[3](local_vertex_indices)
4240 * += (sigma_E*basis_vectors[0])*basis_vectors[1]; //sig_xy
4241 * cauchy_stresses_E_vertex_mpi[4](local_vertex_indices)
4242 * += (sigma_E*basis_vectors[0])*basis_vectors[2];//sig_xz
4243 * cauchy_stresses_E_vertex_mpi[5](local_vertex_indices)
4244 * += (sigma_E*basis_vectors[1])*basis_vectors[2]; //sig_yz
4249 * ---------------------------------------------------------------
4252 * } //end gauss point loop
4257 * Different nodes might have different amount of contributions, e.g.,
4258 * corner nodes have less integration points contributing to the averaged.
4259 * This is why we need a counter and divide at the end, outside the cell loop.
4262 * if (parameters.outtype == "nodes")
4264 * for (unsigned int d=0; d<(vertex_handler_ref.n_dofs()); ++d)
4266 * sum_counter_on_vertices[d] =
4267 * Utilities::MPI::sum(counter_on_vertices_mpi[d],
4268 * mpi_communicator);
4269 * sum_porous_dissipation_vertex[d] =
4270 * Utilities::MPI::sum(porous_dissipation_vertex_mpi[d],
4271 * mpi_communicator);
4272 * sum_viscous_dissipation_vertex[d] =
4273 * Utilities::MPI::sum(viscous_dissipation_vertex_mpi[d],
4274 * mpi_communicator);
4275 * sum_solid_vol_fraction_vertex[d] =
4276 * Utilities::MPI::sum(solid_vol_fraction_vertex_mpi[d],
4277 * mpi_communicator);
4279 * for (unsigned int k=0; k<num_comp_symm_tensor; ++k)
4281 * sum_cauchy_stresses_total_vertex[k][d] =
4282 * Utilities::MPI::sum(cauchy_stresses_total_vertex_mpi[k][d],
4283 * mpi_communicator);
4284 * sum_cauchy_stresses_E_vertex[k][d] =
4285 * Utilities::MPI::sum(cauchy_stresses_E_vertex_mpi[k][d],
4286 * mpi_communicator);
4288 * for (unsigned int k=0; k<dim; ++k)
4290 * sum_stretches_vertex[k][d] =
4291 * Utilities::MPI::sum(stretches_vertex_mpi[k][d],
4292 * mpi_communicator);
4296 * for (unsigned int d=0; d<(vertex_vec_handler_ref.n_dofs()); ++d)
4298 * sum_counter_on_vertices_vec[d] =
4299 * Utilities::MPI::sum(counter_on_vertices_vec_mpi[d],
4300 * mpi_communicator);
4301 * sum_seepage_velocity_vertex_vec[d] =
4302 * Utilities::MPI::sum(seepage_velocity_vertex_vec_mpi[d],
4303 * mpi_communicator);
4306 * for (unsigned int d=0; d<(vertex_handler_ref.n_dofs()); ++d)
4308 * if (sum_counter_on_vertices[d]>0)
4310 * for (unsigned int i=0; i<num_comp_symm_tensor; ++i)
4312 * sum_cauchy_stresses_total_vertex[i][d] /= sum_counter_on_vertices[d];
4313 * sum_cauchy_stresses_E_vertex[i][d] /= sum_counter_on_vertices[d];
4315 * for (unsigned int i=0; i<dim; ++i)
4317 * sum_stretches_vertex[i][d] /= sum_counter_on_vertices[d];
4319 * sum_porous_dissipation_vertex[d] /= sum_counter_on_vertices[d];
4320 * sum_viscous_dissipation_vertex[d] /= sum_counter_on_vertices[d];
4321 * sum_solid_vol_fraction_vertex[d] /= sum_counter_on_vertices[d];
4325 * for (unsigned int d=0; d<(vertex_vec_handler_ref.n_dofs()); ++d)
4327 * if (sum_counter_on_vertices_vec[d]>0)
4329 * sum_seepage_velocity_vertex_vec[d] /= sum_counter_on_vertices_vec[d];
4337 * Add the results to the solution to create the output file for Paraview
4340 * FilteredDataOut<dim> data_out;
4341 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
4343 * DataComponentInterpretation::component_is_part_of_vector);
4344 * comp_type.push_back(DataComponentInterpretation::component_is_scalar);
4346 * GridTools::get_subdomain_association(triangulation, partition_int);
4348 * std::vector<std::string> solution_name(dim, "displacement");
4349 * solution_name.push_back("pore_pressure");
4351 * data_out.attach_dof_handler(dof_handler_ref);
4352 * data_out.add_data_vector(solution_total,
4354 * DataOut<dim>::type_dof_data,
4357 * data_out.add_data_vector(solution_total,
4358 * gradient_postprocessor);
4360 * const Vector<double> partitioning(partition_int.begin(),
4361 * partition_int.end());
4363 * data_out.add_data_vector(partitioning, "partitioning");
4364 * data_out.add_data_vector(material_id, "material_id");
4368 * Integration point results -----------------------------------------------------------
4371 * if (parameters.outtype == "elements")
4373 * data_out.add_data_vector(cauchy_stresses_total_elements[0], "cauchy_xx");
4374 * data_out.add_data_vector(cauchy_stresses_total_elements[1], "cauchy_yy");
4375 * data_out.add_data_vector(cauchy_stresses_total_elements[2], "cauchy_zz");
4376 * data_out.add_data_vector(cauchy_stresses_total_elements[3], "cauchy_xy");
4377 * data_out.add_data_vector(cauchy_stresses_total_elements[4], "cauchy_xz");
4378 * data_out.add_data_vector(cauchy_stresses_total_elements[5], "cauchy_yz");
4380 * data_out.add_data_vector(cauchy_stresses_E_elements[0], "cauchy_E_xx");
4381 * data_out.add_data_vector(cauchy_stresses_E_elements[1], "cauchy_E_yy");
4382 * data_out.add_data_vector(cauchy_stresses_E_elements[2], "cauchy_E_zz");
4383 * data_out.add_data_vector(cauchy_stresses_E_elements[3], "cauchy_E_xy");
4384 * data_out.add_data_vector(cauchy_stresses_E_elements[4], "cauchy_E_xz");
4385 * data_out.add_data_vector(cauchy_stresses_E_elements[5], "cauchy_E_yz");
4387 * data_out.add_data_vector(stretches_elements[0], "stretch_xx");
4388 * data_out.add_data_vector(stretches_elements[1], "stretch_yy");
4389 * data_out.add_data_vector(stretches_elements[2], "stretch_zz");
4391 * data_out.add_data_vector(seepage_velocity_elements[0], "seepage_vel_x");
4392 * data_out.add_data_vector(seepage_velocity_elements[1], "seepage_vel_y");
4393 * data_out.add_data_vector(seepage_velocity_elements[2], "seepage_vel_z");
4395 * data_out.add_data_vector(porous_dissipation_elements, "dissipation_porous");
4396 * data_out.add_data_vector(viscous_dissipation_elements, "dissipation_viscous");
4397 * data_out.add_data_vector(solid_vol_fraction_elements, "solid_vol_fraction");
4399 * else if (parameters.outtype == "nodes")
4401 * data_out.add_data_vector(vertex_handler_ref,
4402 * sum_cauchy_stresses_total_vertex[0],
4404 * data_out.add_data_vector(vertex_handler_ref,
4405 * sum_cauchy_stresses_total_vertex[1],
4407 * data_out.add_data_vector(vertex_handler_ref,
4408 * sum_cauchy_stresses_total_vertex[2],
4410 * data_out.add_data_vector(vertex_handler_ref,
4411 * sum_cauchy_stresses_total_vertex[3],
4413 * data_out.add_data_vector(vertex_handler_ref,
4414 * sum_cauchy_stresses_total_vertex[4],
4416 * data_out.add_data_vector(vertex_handler_ref,
4417 * sum_cauchy_stresses_total_vertex[5],
4420 * data_out.add_data_vector(vertex_handler_ref,
4421 * sum_cauchy_stresses_E_vertex[0],
4423 * data_out.add_data_vector(vertex_handler_ref,
4424 * sum_cauchy_stresses_E_vertex[1],
4426 * data_out.add_data_vector(vertex_handler_ref,
4427 * sum_cauchy_stresses_E_vertex[2],
4429 * data_out.add_data_vector(vertex_handler_ref,
4430 * sum_cauchy_stresses_E_vertex[3],
4432 * data_out.add_data_vector(vertex_handler_ref,
4433 * sum_cauchy_stresses_E_vertex[4],
4435 * data_out.add_data_vector(vertex_handler_ref,
4436 * sum_cauchy_stresses_E_vertex[5],
4439 * data_out.add_data_vector(vertex_handler_ref,
4440 * sum_stretches_vertex[0],
4442 * data_out.add_data_vector(vertex_handler_ref,
4443 * sum_stretches_vertex[1],
4445 * data_out.add_data_vector(vertex_handler_ref,
4446 * sum_stretches_vertex[2],
4449 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
4450 * comp_type_vec(dim,
4451 * DataComponentInterpretation::component_is_part_of_vector);
4452 * std::vector<std::string> solution_name_vec(dim,"seepage_velocity");
4454 * data_out.add_data_vector(vertex_vec_handler_ref,
4455 * sum_seepage_velocity_vertex_vec,
4456 * solution_name_vec,
4459 * data_out.add_data_vector(vertex_handler_ref,
4460 * sum_porous_dissipation_vertex,
4461 * "dissipation_porous");
4462 * data_out.add_data_vector(vertex_handler_ref,
4463 * sum_viscous_dissipation_vertex,
4464 * "dissipation_viscous");
4465 * data_out.add_data_vector(vertex_handler_ref,
4466 * sum_solid_vol_fraction_vertex,
4467 * "solid_vol_fraction");
4471 * ---------------------------------------------------------------------
4477 * data_out.build_patches(degree_displ);
4481 * static std::string get_filename_vtu(unsigned int process,
4482 * unsigned int timestep,
4483 * const unsigned int n_digits = 5)
4485 * std::ostringstream filename_vtu;
4488 * << Utilities::int_to_string(process, n_digits)
4490 * << Utilities::int_to_string(timestep, n_digits)
4492 * return filename_vtu.str();
4495 * static std::string get_filename_pvtu(unsigned int timestep,
4496 * const unsigned int n_digits = 5)
4498 * std::ostringstream filename_vtu;
4501 * << Utilities::int_to_string(timestep, n_digits)
4503 * return filename_vtu.str();
4506 * static std::string get_filename_pvd (void)
4508 * std::ostringstream filename_vtu;
4510 * << "solution.pvd";
4511 * return filename_vtu.str();
4515 * const std::string filename_vtu = Filename::get_filename_vtu(this_mpi_process,
4517 * std::ofstream output(filename_vtu.c_str());
4518 * data_out.write_vtu(output);
4522 * We have a collection of files written in parallel
4523 * This next set of steps should only be performed by master process
4526 * if (this_mpi_process == 0)
4530 * List of all files written out at this timestep by all processors
4533 * std::vector<std::string> parallel_filenames_vtu;
4534 * for (unsigned int p=0; p<n_mpi_processes; ++p)
4536 * parallel_filenames_vtu.push_back(Filename::get_filename_vtu(p, timestep));
4539 * const std::string filename_pvtu(Filename::get_filename_pvtu(timestep));
4540 * std::ofstream pvtu_master(filename_pvtu.c_str());
4541 * data_out.write_pvtu_record(pvtu_master,
4542 * parallel_filenames_vtu);
4546 * Time dependent data master file
4549 * static std::vector<std::pair<double,std::string>> time_and_name_history;
4550 * time_and_name_history.push_back(std::make_pair(current_time,
4552 * const std::string filename_pvd(Filename::get_filename_pvd());
4553 * std::ofstream pvd_output(filename_pvd.c_str());
4554 * DataOutBase::write_pvd_record(pvd_output, time_and_name_history);
4561 * Print results to plotting file
4564 * template <int dim>
4565 * void Solid<dim>::output_results_to_plot(
4566 * const unsigned int timestep,
4567 * const double current_time,
4568 * TrilinosWrappers::MPI::BlockVector solution_IN,
4569 * std::vector<Point<dim> > &tracked_vertices_IN,
4570 * std::ofstream &plotpointfile) const
4572 * TrilinosWrappers::MPI::BlockVector solution_total(locally_owned_partitioning,
4573 * locally_relevant_partitioning,
4578 * solution_total = solution_IN;
4582 * Variables needed to print the solution file for plotting
4585 * Point<dim> reaction_force;
4586 * Point<dim> reaction_force_pressure;
4587 * Point<dim> reaction_force_extra;
4588 * double total_fluid_flow = 0.0;
4589 * double total_porous_dissipation = 0.0;
4590 * double total_viscous_dissipation = 0.0;
4591 * double total_solid_vol = 0.0;
4592 * double total_vol_current = 0.0;
4593 * double total_vol_reference = 0.0;
4594 * std::vector<Point<dim+1>> solution_vertices(tracked_vertices_IN.size());
4598 * Auxiliar variables needed for mpi processing
4601 * Tensor<1,dim> sum_reaction_mpi;
4602 * Tensor<1,dim> sum_reaction_pressure_mpi;
4603 * Tensor<1,dim> sum_reaction_extra_mpi;
4604 * sum_reaction_mpi = 0.0;
4605 * sum_reaction_pressure_mpi = 0.0;
4606 * sum_reaction_extra_mpi = 0.0;
4607 * double sum_total_flow_mpi = 0.0;
4608 * double sum_porous_dissipation_mpi = 0.0;
4609 * double sum_viscous_dissipation_mpi = 0.0;
4610 * double sum_solid_vol_mpi = 0.0;
4611 * double sum_vol_current_mpi = 0.0;
4612 * double sum_vol_reference_mpi = 0.0;
4616 * Declare an instance of the material class object
4619 * if (parameters.mat_type == "Neo-Hooke")
4620 * NeoHooke<dim,ADNumberType> material(parameters,time);
4621 * else if (parameters.mat_type == "Ogden")
4622 * Ogden<dim,ADNumberType> material(parameters, time);
4623 * else if (parameters.mat_type == "visco-Ogden")
4624 * visco_Ogden <dim,ADNumberType>material(parameters,time);
4626 * Assert (false, ExcMessage("Material type not implemented"));
4630 * Define a local instance of FEValues to compute updated values required
4631 * to calculate stresses
4634 * const UpdateFlags uf_cell(update_values | update_gradients |
4635 * update_JxW_values);
4636 * FEValues<dim> fe_values_ref (fe, qf_cell, uf_cell);
4640 * Iterate through elements (cells) and Gauss Points
4643 * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
4644 * cell(IteratorFilters::LocallyOwnedCell(),
4645 * dof_handler_ref.begin_active()),
4646 * endc(IteratorFilters::LocallyOwnedCell(),
4647 * dof_handler_ref.end());
4653 * for (; cell!=endc; ++cell)
4655 * Assert(cell->is_locally_owned(), ExcInternalError());
4656 * Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());
4658 * fe_values_ref.reinit(cell);
4660 * std::vector<Tensor<2,dim>> solution_grads_u(n_q_points);
4661 * fe_values_ref[u_fe].get_function_gradients(solution_total,
4662 * solution_grads_u);
4664 * std::vector<double> solution_values_p_fluid_total(n_q_points);
4665 * fe_values_ref[p_fluid_fe].get_function_values(solution_total,
4666 * solution_values_p_fluid_total);
4668 * std::vector<Tensor<1,dim >> solution_grads_p_fluid_AD(n_q_points);
4669 * fe_values_ref[p_fluid_fe].get_function_gradients(solution_total,
4670 * solution_grads_p_fluid_AD);
4674 * start gauss point loop
4677 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4679 * const Tensor<2,dim,ADNumberType>
4680 * F_AD = Physics::Elasticity::Kinematics::F(solution_grads_u[q_point]);
4681 * ADNumberType det_F_AD = determinant(F_AD);
4682 * const double det_F = Tensor<0,dim,double>(det_F_AD);
4684 * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4685 * lqph = quadrature_point_history.get_data(cell);
4686 * Assert(lqph.size() == n_q_points, ExcInternalError());
4688 * double JxW = fe_values_ref.JxW(q_point);
4695 * sum_vol_current_mpi += det_F * JxW;
4696 * sum_vol_reference_mpi += JxW;
4697 * sum_solid_vol_mpi += parameters.solid_vol_frac * JxW * det_F;
4704 * const Tensor<2,dim,ADNumberType> F_inv = invert(F_AD);
4705 * const Tensor<1,dim,ADNumberType>
4706 * grad_p_fluid_AD = solution_grads_p_fluid_AD[q_point]*F_inv;
4707 * const Tensor<1,dim,ADNumberType> seepage_vel_AD
4708 * = lqph[q_point]->get_seepage_velocity_current(F_AD, grad_p_fluid_AD);
4715 * const double porous_dissipation =
4716 * lqph[q_point]->get_porous_dissipation(F_AD, grad_p_fluid_AD);
4717 * sum_porous_dissipation_mpi += porous_dissipation * det_F * JxW;
4719 * const double viscous_dissipation = lqph[q_point]->get_viscous_dissipation();
4720 * sum_viscous_dissipation_mpi += viscous_dissipation * det_F * JxW;
4724 * ---------------------------------------------------------------
4727 * } //end gauss point loop
4731 * Compute reaction force on load boundary & total fluid flow across
4733 * Define a local instance of FEFaceValues to compute values required
4734 * to calculate reaction force
4737 * const UpdateFlags uf_face( update_values | update_gradients |
4738 * update_normal_vectors | update_JxW_values );
4739 * FEFaceValues<dim> fe_face_values_ref(fe, qf_face, uf_face);
4746 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
4753 * if (cell->face(face)->at_boundary() == true &&
4754 * cell->face(face)->boundary_id() == get_reaction_boundary_id_for_output() )
4756 * fe_face_values_ref.reinit(cell, face);
4760 * Get displacement gradients for current face
4763 * std::vector<Tensor<2,dim> > solution_grads_u_f(n_q_points_f);
4764 * fe_face_values_ref[u_fe].get_function_gradients
4766 * solution_grads_u_f);
4770 * Get pressure for current element
4773 * std::vector< double > solution_values_p_fluid_total_f(n_q_points_f);
4774 * fe_face_values_ref[p_fluid_fe].get_function_values
4776 * solution_values_p_fluid_total_f);
4780 * start gauss points on faces loop
4783 * for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point)
4785 * const Tensor<1,dim> &N = fe_face_values_ref.normal_vector(f_q_point);
4786 * const double JxW_f = fe_face_values_ref.JxW(f_q_point);
4790 * Compute deformation gradient from displacements gradient
4791 * (present configuration)
4794 * const Tensor<2,dim,ADNumberType> F_AD =
4795 * Physics::Elasticity::Kinematics::F(solution_grads_u_f[f_q_point]);
4797 * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4798 * lqph = quadrature_point_history.get_data(cell);
4799 * Assert(lqph.size() == n_q_points, ExcInternalError());
4801 * const double p_fluid = solution_values_p_fluid_total[f_q_point];
4808 * static const SymmetricTensor<2,dim,double>
4809 * I (Physics::Elasticity::StandardTensors<dim>::I);
4810 * SymmetricTensor<2,dim> sigma_E;
4811 * const SymmetricTensor<2,dim,ADNumberType> sigma_E_AD =
4812 * lqph[f_q_point]->get_Cauchy_E(F_AD);
4814 * for (unsigned int i=0; i<dim; ++i)
4815 * for (unsigned int j=0; j<dim; ++j)
4816 * sigma_E[i][j] = Tensor<0,dim,double>(sigma_E_AD[i][j]);
4818 * SymmetricTensor<2,dim> sigma_fluid_vol(I);
4819 * sigma_fluid_vol *= -1.0*p_fluid;
4820 * const SymmetricTensor<2,dim> sigma = sigma_E+sigma_fluid_vol;
4821 * sum_reaction_mpi += sigma * N * JxW_f;
4822 * sum_reaction_pressure_mpi += sigma_fluid_vol * N * JxW_f;
4823 * sum_reaction_extra_mpi += sigma_E * N * JxW_f;
4824 * }//end gauss points on faces loop
4832 * if (cell->face(face)->at_boundary() == true &&
4833 * (cell->face(face)->boundary_id() ==
4834 * get_drained_boundary_id_for_output().first ||
4835 * cell->face(face)->boundary_id() ==
4836 * get_drained_boundary_id_for_output().second ) )
4838 * fe_face_values_ref.reinit(cell, face);
4842 * Get displacement gradients for current face
4845 * std::vector<Tensor<2,dim>> solution_grads_u_f(n_q_points_f);
4846 * fe_face_values_ref[u_fe].get_function_gradients
4848 * solution_grads_u_f);
4852 * Get pressure gradients for current face
4855 * std::vector<Tensor<1,dim>> solution_grads_p_f(n_q_points_f);
4856 * fe_face_values_ref[p_fluid_fe].get_function_gradients
4858 * solution_grads_p_f);
4862 * start gauss points on faces loop
4865 * for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point)
4867 * const Tensor<1,dim> &N =
4868 * fe_face_values_ref.normal_vector(f_q_point);
4869 * const double JxW_f = fe_face_values_ref.JxW(f_q_point);
4873 * Deformation gradient and inverse from displacements gradient
4874 * (present configuration)
4877 * const Tensor<2,dim,ADNumberType> F_AD
4878 * = Physics::Elasticity::Kinematics::F(solution_grads_u_f[f_q_point]);
4880 * const Tensor<2,dim,ADNumberType> F_inv_AD = invert(F_AD);
4881 * ADNumberType det_F_AD = determinant(F_AD);
4883 * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType>>>
4884 * lqph = quadrature_point_history.get_data(cell);
4885 * Assert(lqph.size() == n_q_points, ExcInternalError());
4892 * Tensor<1,dim> seepage;
4893 * double det_F = Tensor<0,dim,double>(det_F_AD);
4894 * const Tensor<1,dim,ADNumberType> grad_p
4895 * = solution_grads_p_f[f_q_point]*F_inv_AD;
4896 * const Tensor<1,dim,ADNumberType> seepage_AD
4897 * = lqph[f_q_point]->get_seepage_velocity_current(F_AD, grad_p);
4899 * for (unsigned int i=0; i<dim; ++i)
4900 * seepage[i] = Tensor<0,dim,double>(seepage_AD[i]);
4902 * sum_total_flow_mpi += (seepage/det_F) * N * JxW_f;
4903 * }//end gauss points on faces loop
4910 * Sum the results from different MPI process and then add to the reaction_force vector
4911 * In theory, the solution on each surface (each cell) only exists in one MPI process
4912 * so, we add all MPI process, one will have the solution and the others will be zero
4915 * for (unsigned int d=0; d<dim; ++d)
4917 * reaction_force[d] = Utilities::MPI::sum(sum_reaction_mpi[d],
4918 * mpi_communicator);
4919 * reaction_force_pressure[d] = Utilities::MPI::sum(sum_reaction_pressure_mpi[d],
4920 * mpi_communicator);
4921 * reaction_force_extra[d] = Utilities::MPI::sum(sum_reaction_extra_mpi[d],
4922 * mpi_communicator);
4927 * Same for total fluid flow, and for porous and viscous dissipations
4930 * total_fluid_flow = Utilities::MPI::sum(sum_total_flow_mpi,
4931 * mpi_communicator);
4932 * total_porous_dissipation = Utilities::MPI::sum(sum_porous_dissipation_mpi,
4933 * mpi_communicator);
4934 * total_viscous_dissipation = Utilities::MPI::sum(sum_viscous_dissipation_mpi,
4935 * mpi_communicator);
4936 * total_solid_vol = Utilities::MPI::sum(sum_solid_vol_mpi,
4937 * mpi_communicator);
4938 * total_vol_current = Utilities::MPI::sum(sum_vol_current_mpi,
4939 * mpi_communicator);
4940 * total_vol_reference = Utilities::MPI::sum(sum_vol_reference_mpi,
4941 * mpi_communicator);
4945 * Extract solution for tracked vectors
4946 * Copying an MPI::BlockVector into MPI::Vector is not possible,
4947 * so we copy each block of MPI::BlockVector into an MPI::Vector
4948 * And then we copy the MPI::Vector into "normal" Vectors
4951 * TrilinosWrappers::MPI::Vector solution_vector_u_MPI(solution_total.block(u_block));
4952 * TrilinosWrappers::MPI::Vector solution_vector_p_MPI(solution_total.block(p_fluid_block));
4953 * Vector<double> solution_u_vector(solution_vector_u_MPI);
4954 * Vector<double> solution_p_vector(solution_vector_p_MPI);
4956 * if (this_mpi_process == 0)
4960 * Append the pressure solution vector to the displacement solution vector,
4961 * creating a single solution vector equivalent to the original BlockVector
4962 * so FEFieldFunction will work with the dof_handler_ref.
4965 * Vector<double> solution_vector(solution_p_vector.size()
4966 * +solution_u_vector.size());
4968 * for (unsigned int d=0; d<(solution_u_vector.size()); ++d)
4969 * solution_vector[d] = solution_u_vector[d];
4971 * for (unsigned int d=0; d<(solution_p_vector.size()); ++d)
4972 * solution_vector[solution_u_vector.size()+d] = solution_p_vector[d];
4974 * Functions::FEFieldFunction<dim,DoFHandler<dim>,Vector<double>>
4975 * find_solution(dof_handler_ref, solution_vector);
4977 * for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
4979 * Vector<double> update(dim+1);
4980 * Point<dim> pt_ref;
4982 * pt_ref[0]= tracked_vertices_IN[p][0];
4983 * pt_ref[1]= tracked_vertices_IN[p][1];
4984 * pt_ref[2]= tracked_vertices_IN[p][2];
4986 * find_solution.vector_value(pt_ref, update);
4988 * for (unsigned int d=0; d<(dim+1); ++d)
4992 * For values close to zero, set to 0.0
4995 * if (abs(update[d])<1.5*parameters.tol_u)
4997 * solution_vertices[p][d] = update[d];
5002 * Write the results to the plotting file.
5003 * Add two blank lines between cycles in the cyclic loading examples so GNUPLOT can detect each cycle as a different block
5006 * if (( (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")||
5007 * (parameters.geom_type == "Budday_cube_tension_compression")||
5008 * (parameters.geom_type == "Budday_cube_shear_fully_fixed") ) &&
5009 * ( (abs(current_time - parameters.end_time/3.) <0.9*parameters.delta_t)||
5010 * (abs(current_time - 2.*parameters.end_time/3.)<0.9*parameters.delta_t) ) &&
5011 * parameters.num_cycle_sets == 1 )
5013 * plotpointfile << std::endl<< std::endl;
5015 * if (( (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")||
5016 * (parameters.geom_type == "Budday_cube_tension_compression")||
5017 * (parameters.geom_type == "Budday_cube_shear_fully_fixed") ) &&
5018 * ( (abs(current_time - parameters.end_time/9.) <0.9*parameters.delta_t)||
5019 * (abs(current_time - 2.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5020 * (abs(current_time - 3.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5021 * (abs(current_time - 5.*parameters.end_time/9.)<0.9*parameters.delta_t)||
5022 * (abs(current_time - 7.*parameters.end_time/9.)<0.9*parameters.delta_t) ) &&
5023 * parameters.num_cycle_sets == 2 )
5025 * plotpointfile << std::endl<< std::endl;
5028 * plotpointfile << std::setprecision(6) << std::scientific;
5029 * plotpointfile << std::setw(16) << current_time << ","
5030 * << std::setw(15) << total_vol_reference << ","
5031 * << std::setw(15) << total_vol_current << ","
5032 * << std::setw(15) << total_solid_vol << ",";
5034 * if (current_time == 0.0)
5036 * for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
5038 * for (unsigned int d=0; d<dim; ++d)
5039 * plotpointfile << std::setw(15) << 0.0 << ",";
5041 * plotpointfile << std::setw(15) << parameters.drained_pressure << ",";
5043 * for (unsigned int d=0; d<(3*dim+2); ++d)
5044 * plotpointfile << std::setw(15) << 0.0 << ",";
5046 * plotpointfile << std::setw(15) << 0.0;
5050 * for (unsigned int p=0; p<tracked_vertices_IN.size(); ++p)
5051 * for (unsigned int d=0; d<(dim+1); ++d)
5052 * plotpointfile << std::setw(15) << solution_vertices[p][d]<< ",";
5054 * for (unsigned int d=0; d<dim; ++d)
5055 * plotpointfile << std::setw(15) << reaction_force[d] << ",";
5057 * for (unsigned int d=0; d<dim; ++d)
5058 * plotpointfile << std::setw(15) << reaction_force_pressure[d] << ",";
5060 * for (unsigned int d=0; d<dim; ++d)
5061 * plotpointfile << std::setw(15) << reaction_force_extra[d] << ",";
5063 * plotpointfile << std::setw(15) << total_fluid_flow << ","
5064 * << std::setw(15) << total_porous_dissipation<< ","
5065 * << std::setw(15) << total_viscous_dissipation;
5067 * plotpointfile << std::endl;
5073 * Header for console output file
5076 * template <int dim>
5077 * void Solid<dim>::print_console_file_header(std::ofstream &outputfile) const
5079 * outputfile << "/*-----------------------------------------------------------------------------------------";
5080 * outputfile << "\n\n Poro-viscoelastic formulation to solve nonlinear solid mechanics problems using deal.ii";
5081 * outputfile << "\n\n Problem setup by E Comellas and J-P Pelteret, University of Erlangen-Nuremberg, 2018";
5082 * outputfile << "\n\n/*-----------------------------------------------------------------------------------------";
5083 * outputfile << "\n\nCONSOLE OUTPUT: \n\n";
5088 * Header for plotting output file
5091 * template <int dim>
5092 * void Solid<dim>::print_plot_file_header(std::vector<Point<dim> > &tracked_vertices,
5093 * std::ofstream &plotpointfile) const
5095 * plotpointfile << "#\n# *** Solution history for tracked vertices -- DOF: 0 = Ux, 1 = Uy, 2 = Uz, 3 = P ***"
5098 * for (unsigned int p=0; p<tracked_vertices.size(); ++p)
5100 * plotpointfile << "# Point " << p << " coordinates: ";
5101 * for (unsigned int d=0; d<dim; ++d)
5103 * plotpointfile << tracked_vertices[p][d];
5104 * if (!( (p == tracked_vertices.size()-1) && (d == dim-1) ))
5105 * plotpointfile << ", ";
5107 * plotpointfile << std::endl;
5109 * plotpointfile << "# The reaction force is the integral over the loaded surfaces in the "
5110 * << "undeformed configuration of the Cauchy stress times the normal surface unit vector.\n"
5111 * << "# reac(p) corresponds to the volumetric part of the Cauchy stress due to the pore fluid pressure"
5112 * << " and reac(E) corresponds to the extra part of the Cauchy stress due to the solid contribution."
5114 * << "# The fluid flow is the integral over the drained surfaces in the "
5115 * << "undeformed configuration of the seepage velocity times the normal surface unit vector."
5117 * << "# Column number:"
5121 * unsigned int columns = 24;
5122 * for (unsigned int d=1; d<columns; ++d)
5123 * plotpointfile << std::setw(15)<< d <<",";
5125 * plotpointfile << std::setw(15)<< columns
5128 * << std::right << std::setw(16) << "Time,"
5129 * << std::right << std::setw(16) << "ref vol,"
5130 * << std::right << std::setw(16) << "def vol,"
5131 * << std::right << std::setw(16) << "solid vol,";
5132 * for (unsigned int p=0; p<tracked_vertices.size(); ++p)
5133 * for (unsigned int d=0; d<(dim+1); ++d)
5134 * plotpointfile << std::right<< std::setw(11)
5135 * <<"P" << p << "[" << d << "],";
5137 * for (unsigned int d=0; d<dim; ++d)
5138 * plotpointfile << std::right<< std::setw(13)
5139 * << "reaction [" << d << "],";
5141 * for (unsigned int d=0; d<dim; ++d)
5142 * plotpointfile << std::right<< std::setw(13)
5143 * << "reac(p) [" << d << "],";
5145 * for (unsigned int d=0; d<dim; ++d)
5146 * plotpointfile << std::right<< std::setw(13)
5147 * << "reac(E) [" << d << "],";
5149 * plotpointfile << std::right<< std::setw(16)<< "fluid flow,"
5150 * << std::right<< std::setw(16)<< "porous dissip,"
5151 * << std::right<< std::setw(15)<< "viscous dissip"
5157 * Footer for console output file
5160 * template <int dim>
5161 * void Solid<dim>::print_console_file_footer(std::ofstream &outputfile) const
5165 * Copy "parameters" file at end of output file.
5168 * std::ifstream infile("parameters.prm");
5169 * std::string content = "";
5172 * for(i=0 ; infile.eof()!=true ; i++)
5174 * char aux = infile.get();
5176 * if(aux=='\n
') content += '#
';
5180 * content.erase(content.end()-1);
5183 * outputfile << "\n\n\n\n PARAMETERS FILE USED IN THIS COMPUTATION: \n#"
5190 * Footer for plotting output file
5193 * template <int dim>
5194 * void Solid<dim>::print_plot_file_footer(std::ofstream &plotpointfile) const
5198 * Copy "parameters" file at end of output file.
5201 * std::ifstream infile("parameters.prm");
5202 * std::string content = "";
5205 * for(i=0 ; infile.eof()!=true ; i++)
5207 * char aux = infile.get();
5209 * if(aux=='\n
') content += '#
';
5213 * content.erase(content.end()-1);
5216 * plotpointfile << "#"<< std::endl
5217 * << "#"<< std::endl
5218 * << "# PARAMETERS FILE USED IN THIS COMPUTATION:" << std::endl
5219 * << "#"<< std::endl
5227 * <a name="VerificationexamplesfromEhlersandEipper1999"></a>
5228 * <h3>Verification examples from Ehlers and Eipper 1999</h3>
5229 * We group the definition of the geometry, boundary and loading conditions specific to
5230 * the verification examples from Ehlers and Eipper 1999 into specific classes.
5235 * <a name="BaseclassTubegeometryandboundaryconditions"></a>
5236 * <h4>Base class: Tube geometry and boundary conditions</h4>
5239 * template <int dim>
5240 * class VerificationEhlers1999TubeBase
5241 * : public Solid<dim>
5244 * VerificationEhlers1999TubeBase (const Parameters::AllParameters ¶meters)
5245 * : Solid<dim> (parameters)
5248 * virtual ~VerificationEhlers1999TubeBase () {}
5251 * virtual void make_grid()
5253 * GridGenerator::cylinder( this->triangulation,
5257 * const double rot_angle = 3.0*numbers::PI/2.0;
5258 * GridTools::rotate( rot_angle, 1, this->triangulation);
5260 * this->triangulation.reset_manifold(0);
5261 * static const CylindricalManifold<dim> manifold_description_3d(2);
5262 * this->triangulation.set_manifold (0, manifold_description_3d);
5263 * GridTools::scale(this->parameters.scale, this->triangulation);
5264 * this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5265 * this->triangulation.reset_manifold(0);
5268 * virtual void define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices)
5270 * tracked_vertices[0][0] = 0.0*this->parameters.scale;
5271 * tracked_vertices[0][1] = 0.0*this->parameters.scale;
5272 * tracked_vertices[0][2] = 0.5*this->parameters.scale;
5274 * tracked_vertices[1][0] = 0.0*this->parameters.scale;
5275 * tracked_vertices[1][1] = 0.0*this->parameters.scale;
5276 * tracked_vertices[1][2] = -0.5*this->parameters.scale;
5279 * virtual void make_dirichlet_constraints(AffineConstraints<double> &constraints)
5281 * if (this->time.get_timestep() < 2)
5283 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5285 * ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5287 * (this->fe.component_mask(this->pressure)));
5291 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5293 * ZeroFunction<dim>(this->n_components),
5295 * (this->fe.component_mask(this->pressure)));
5298 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5300 * ZeroFunction<dim>(this->n_components),
5302 * (this->fe.component_mask(this->x_displacement)|
5303 * this->fe.component_mask(this->y_displacement) ) );
5305 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5307 * ZeroFunction<dim>(this->n_components),
5309 * (this->fe.component_mask(this->x_displacement) |
5310 * this->fe.component_mask(this->y_displacement) |
5311 * this->fe.component_mask(this->z_displacement) ));
5315 * get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5316 * const Point<dim> &pt) const
5319 * (void)boundary_id;
5323 * virtual types::boundary_id
5324 * get_reaction_boundary_id_for_output() const
5329 * virtual std::pair<types::boundary_id,types::boundary_id>
5330 * get_drained_boundary_id_for_output() const
5332 * return std::make_pair(2,2);
5335 * virtual std::vector<double>
5336 * get_dirichlet_load(const types::boundary_id &boundary_id,
5337 * const int &direction) const
5339 * std::vector<double> displ_incr(dim, 0.0);
5340 * (void)boundary_id;
5342 * AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));
5344 * return displ_incr;
5351 * <a name="DerivedclassSteploadexample"></a>
5352 * <h4>Derived class: Step load example</h4>
5355 * template <int dim>
5356 * class VerificationEhlers1999StepLoad
5357 * : public VerificationEhlers1999TubeBase<dim>
5360 * VerificationEhlers1999StepLoad (const Parameters::AllParameters ¶meters)
5361 * : VerificationEhlers1999TubeBase<dim> (parameters)
5364 * virtual ~VerificationEhlers1999StepLoad () {}
5367 * virtual Tensor<1,dim>
5368 * get_neumann_traction (const types::boundary_id &boundary_id,
5369 * const Point<dim> &pt,
5370 * const Tensor<1,dim> &N) const
5372 * if (this->parameters.load_type == "pressure")
5374 * if (boundary_id == 2)
5376 * return this->parameters.load * N;
5382 * return Tensor<1,dim>();
5389 * <a name="DerivedclassLoadincreasingexample"></a>
5390 * <h4>Derived class: Load increasing example</h4>
5393 * template <int dim>
5394 * class VerificationEhlers1999IncreaseLoad
5395 * : public VerificationEhlers1999TubeBase<dim>
5398 * VerificationEhlers1999IncreaseLoad (const Parameters::AllParameters ¶meters)
5399 * : VerificationEhlers1999TubeBase<dim> (parameters)
5402 * virtual ~VerificationEhlers1999IncreaseLoad () {}
5405 * virtual Tensor<1,dim>
5406 * get_neumann_traction (const types::boundary_id &boundary_id,
5407 * const Point<dim> &pt,
5408 * const Tensor<1,dim> &N) const
5410 * if (this->parameters.load_type == "pressure")
5412 * if (boundary_id == 2)
5414 * const double initial_load = this->parameters.load;
5415 * const double final_load = 20.0*initial_load;
5416 * const double initial_time = this->time.get_delta_t();
5417 * const double final_time = this->time.get_end();
5418 * const double current_time = this->time.get_current();
5419 * const double load = initial_load + (final_load-initial_load)*(current_time-initial_time)/(final_time-initial_time);
5426 * return Tensor<1,dim>();
5433 * <a name="ClassConsolidationcube"></a>
5434 * <h4>Class: Consolidation cube</h4>
5437 * template <int dim>
5438 * class VerificationEhlers1999CubeConsolidation
5439 * : public Solid<dim>
5442 * VerificationEhlers1999CubeConsolidation (const Parameters::AllParameters ¶meters)
5443 * : Solid<dim> (parameters)
5446 * virtual ~VerificationEhlers1999CubeConsolidation () {}
5452 * GridGenerator::hyper_rectangle(this->triangulation,
5453 * Point<dim>(0.0, 0.0, 0.0),
5454 * Point<dim>(1.0, 1.0, 1.0),
5457 * GridTools::scale(this->parameters.scale, this->triangulation);
5458 * this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5460 * typename Triangulation<dim>::active_cell_iterator cell =
5461 * this->triangulation.begin_active(), endc = this->triangulation.end();
5462 * for (; cell != endc; ++cell)
5464 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5465 * if (cell->face(face)->at_boundary() == true &&
5466 * cell->face(face)->center()[2] == 1.0 * this->parameters.scale)
5468 * if (cell->face(face)->center()[0] < 0.5 * this->parameters.scale &&
5469 * cell->face(face)->center()[1] < 0.5 * this->parameters.scale)
5470 * cell->face(face)->set_boundary_id(100);
5472 * cell->face(face)->set_boundary_id(101);
5478 * define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices)
5480 * tracked_vertices[0][0] = 0.0*this->parameters.scale;
5481 * tracked_vertices[0][1] = 0.0*this->parameters.scale;
5482 * tracked_vertices[0][2] = 1.0*this->parameters.scale;
5484 * tracked_vertices[1][0] = 0.0*this->parameters.scale;
5485 * tracked_vertices[1][1] = 0.0*this->parameters.scale;
5486 * tracked_vertices[1][2] = 0.0*this->parameters.scale;
5490 * make_dirichlet_constraints(AffineConstraints<double> &constraints)
5492 * if (this->time.get_timestep() < 2)
5494 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5496 * ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5498 * (this->fe.component_mask(this->pressure)));
5502 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5504 * ZeroFunction<dim>(this->n_components),
5506 * (this->fe.component_mask(this->pressure)));
5509 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5511 * ZeroFunction<dim>(this->n_components),
5513 * this->fe.component_mask(this->x_displacement));
5515 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5517 * ZeroFunction<dim>(this->n_components),
5519 * this->fe.component_mask(this->x_displacement));
5521 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5523 * ZeroFunction<dim>(this->n_components),
5525 * this->fe.component_mask(this->y_displacement));
5527 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5529 * ZeroFunction<dim>(this->n_components),
5531 * this->fe.component_mask(this->y_displacement));
5533 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5535 * ZeroFunction<dim>(this->n_components),
5537 * ( this->fe.component_mask(this->x_displacement) |
5538 * this->fe.component_mask(this->y_displacement) |
5539 * this->fe.component_mask(this->z_displacement) ));
5542 * virtual Tensor<1,dim>
5543 * get_neumann_traction (const types::boundary_id &boundary_id,
5544 * const Point<dim> &pt,
5545 * const Tensor<1,dim> &N) const
5547 * if (this->parameters.load_type == "pressure")
5549 * if (boundary_id == 100)
5551 * return this->parameters.load * N;
5557 * return Tensor<1,dim>();
5561 * get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5562 * const Point<dim> &pt) const
5565 * (void)boundary_id;
5569 * virtual types::boundary_id
5570 * get_reaction_boundary_id_for_output() const
5575 * virtual std::pair<types::boundary_id,types::boundary_id>
5576 * get_drained_boundary_id_for_output() const
5578 * return std::make_pair(101,101);
5581 * virtual std::vector<double>
5582 * get_dirichlet_load(const types::boundary_id &boundary_id,
5583 * const int &direction) const
5585 * std::vector<double> displ_incr(dim, 0.0);
5586 * (void)boundary_id;
5588 * AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));
5590 * return displ_incr;
5597 * <a name="Franceschiniexperiments"></a>
5598 * <h4>Franceschini experiments</h4>
5601 * template <int dim>
5602 * class Franceschini2006Consolidation
5603 * : public Solid<dim>
5606 * Franceschini2006Consolidation (const Parameters::AllParameters ¶meters)
5607 * : Solid<dim> (parameters)
5610 * virtual ~Franceschini2006Consolidation () {}
5613 * virtual void make_grid()
5615 * const Point<dim-1> mesh_center(0.0, 0.0);
5616 * const double radius = 0.5;
5619 * const double height = 0.27; //8.1 mm for 30 mm radius
5622 * const double height = 0.23; //6.9 mm for 30 mm radius
5623 * Triangulation<dim-1> triangulation_in;
5624 * GridGenerator::hyper_ball( triangulation_in,
5628 * GridGenerator::extrude_triangulation(triangulation_in,
5631 * this->triangulation);
5633 * const CylindricalManifold<dim> cylinder_3d(2);
5634 * const types::manifold_id cylinder_id = 0;
5637 * this->triangulation.set_manifold(cylinder_id, cylinder_3d);
5639 * for (auto cell : this->triangulation.active_cell_iterators())
5641 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5643 * if (cell->face(face)->at_boundary() == true)
5645 * if (cell->face(face)->center()[2] == 0.0)
5646 * cell->face(face)->set_boundary_id(1);
5648 * else if (cell->face(face)->center()[2] == height)
5649 * cell->face(face)->set_boundary_id(2);
5653 * cell->face(face)->set_boundary_id(0);
5654 * cell->face(face)->set_all_manifold_ids(cylinder_id);
5660 * GridTools::scale(this->parameters.scale, this->triangulation);
5661 * this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5664 * virtual void define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices)
5666 * tracked_vertices[0][0] = 0.0*this->parameters.scale;
5667 * tracked_vertices[0][1] = 0.0*this->parameters.scale;
5670 * tracked_vertices[0][2] = 0.27*this->parameters.scale;
5673 * tracked_vertices[0][2] = 0.23*this->parameters.scale;
5675 * tracked_vertices[1][0] = 0.0*this->parameters.scale;
5676 * tracked_vertices[1][1] = 0.0*this->parameters.scale;
5677 * tracked_vertices[1][2] = 0.0*this->parameters.scale;
5680 * virtual void make_dirichlet_constraints(AffineConstraints<double> &constraints)
5682 * if (this->time.get_timestep() < 2)
5684 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5686 * ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5688 * (this->fe.component_mask(this->pressure)));
5690 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5692 * ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5694 * (this->fe.component_mask(this->pressure)));
5698 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5700 * ZeroFunction<dim>(this->n_components),
5702 * (this->fe.component_mask(this->pressure)));
5704 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5706 * ZeroFunction<dim>(this->n_components),
5708 * (this->fe.component_mask(this->pressure)));
5711 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5713 * ZeroFunction<dim>(this->n_components),
5715 * (this->fe.component_mask(this->x_displacement)|
5716 * this->fe.component_mask(this->y_displacement) ) );
5718 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5720 * ZeroFunction<dim>(this->n_components),
5722 * (this->fe.component_mask(this->x_displacement) |
5723 * this->fe.component_mask(this->y_displacement) |
5724 * this->fe.component_mask(this->z_displacement) ));
5726 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5728 * ZeroFunction<dim>(this->n_components),
5730 * (this->fe.component_mask(this->x_displacement) |
5731 * this->fe.component_mask(this->y_displacement) ));
5735 * get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5736 * const Point<dim> &pt) const
5739 * (void)boundary_id;
5743 * virtual types::boundary_id
5744 * get_reaction_boundary_id_for_output() const
5749 * virtual std::pair<types::boundary_id,types::boundary_id>
5750 * get_drained_boundary_id_for_output() const
5752 * return std::make_pair(1,2);
5755 * virtual std::vector<double>
5756 * get_dirichlet_load(const types::boundary_id &boundary_id,
5757 * const int &direction) const
5759 * std::vector<double> displ_incr(dim, 0.0);
5760 * (void)boundary_id;
5762 * AssertThrow(false, ExcMessage("Displacement loading not implemented for Franceschini examples."));
5764 * return displ_incr;
5767 * virtual Tensor<1,dim>
5768 * get_neumann_traction (const types::boundary_id &boundary_id,
5769 * const Point<dim> &pt,
5770 * const Tensor<1,dim> &N) const
5772 * if (this->parameters.load_type == "pressure")
5774 * if (boundary_id == 2)
5776 * return (this->parameters.load * N);
5778 * const double final_load = this->parameters.load;
5779 * const double final_load_time = 10 * this->time.get_delta_t();
5780 * const double current_time = this->time.get_current();
5783 * const double c = final_load_time / 2.0;
5784 * const double r = 200.0 * 0.03 / c;
5786 * const double load = final_load * std::exp(r * current_time)
5787 * / ( std::exp(c * current_time) + std::exp(r * current_time));
5795 * return Tensor<1,dim>();
5802 * <a name="ExamplestoreproduceexperimentsbyBuddayetal2017"></a>
5803 * <h3>Examples to reproduce experiments by Budday et al. 2017</h3>
5804 * We group the definition of the geometry, boundary and loading conditions specific to
5805 * the examples to reproduce experiments by Budday et al. 2017 into specific classes.
5810 * <a name="BaseclassCubegeometryandloadingpattern"></a>
5811 * <h4>Base class: Cube geometry and loading pattern</h4>
5814 * template <int dim>
5815 * class BrainBudday2017BaseCube
5816 * : public Solid<dim>
5819 * BrainBudday2017BaseCube (const Parameters::AllParameters ¶meters)
5820 * : Solid<dim> (parameters)
5823 * virtual ~BrainBudday2017BaseCube () {}
5829 * GridGenerator::hyper_cube(this->triangulation,
5834 * typename Triangulation<dim>::active_cell_iterator cell =
5835 * this->triangulation.begin_active(), endc = this->triangulation.end();
5836 * for (; cell != endc; ++cell)
5838 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5839 * if (cell->face(face)->at_boundary() == true &&
5840 * ( cell->face(face)->boundary_id() == 0 ||
5841 * cell->face(face)->boundary_id() == 1 ||
5842 * cell->face(face)->boundary_id() == 2 ||
5843 * cell->face(face)->boundary_id() == 3 ) )
5845 * cell->face(face)->set_boundary_id(100);
5849 * GridTools::scale(this->parameters.scale, this->triangulation);
5850 * this->triangulation.refine_global(std::max (1U, this->parameters.global_refinement));
5854 * get_prescribed_fluid_flow (const types::boundary_id &boundary_id,
5855 * const Point<dim> &pt) const
5858 * (void)boundary_id;
5862 * virtual std::pair<types::boundary_id,types::boundary_id>
5863 * get_drained_boundary_id_for_output() const
5865 * return std::make_pair(100,100);
5872 * <a name="DerivedclassUniaxialboundaryconditions"></a>
5873 * <h4>Derived class: Uniaxial boundary conditions</h4>
5876 * template <int dim>
5877 * class BrainBudday2017CubeTensionCompression
5878 * : public BrainBudday2017BaseCube<dim>
5881 * BrainBudday2017CubeTensionCompression (const Parameters::AllParameters ¶meters)
5882 * : BrainBudday2017BaseCube<dim> (parameters)
5885 * virtual ~BrainBudday2017CubeTensionCompression () {}
5889 * define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices)
5891 * tracked_vertices[0][0] = 0.5*this->parameters.scale;
5892 * tracked_vertices[0][1] = 0.5*this->parameters.scale;
5893 * tracked_vertices[0][2] = 1.0*this->parameters.scale;
5895 * tracked_vertices[1][0] = 0.5*this->parameters.scale;
5896 * tracked_vertices[1][1] = 0.5*this->parameters.scale;
5897 * tracked_vertices[1][2] = 0.5*this->parameters.scale;
5901 * make_dirichlet_constraints(AffineConstraints<double> &constraints)
5903 * if (this->time.get_timestep() < 2)
5905 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
5907 * ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
5909 * (this->fe.component_mask(this->pressure)));
5913 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5915 * ZeroFunction<dim>(this->n_components),
5917 * (this->fe.component_mask(this->pressure)));
5919 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5921 * ZeroFunction<dim>(this->n_components),
5923 * this->fe.component_mask(this->z_displacement) );
5925 * Point<dim> fix_node(0.5*this->parameters.scale, 0.5*this->parameters.scale, 0.0);
5926 * typename DoFHandler<dim>::active_cell_iterator
5927 * cell = this->dof_handler_ref.begin_active(), endc = this->dof_handler_ref.end();
5928 * for (; cell != endc; ++cell)
5929 * for (unsigned int node = 0; node < GeometryInfo<dim>::vertices_per_cell; ++node)
5931 * if ( (abs(cell->vertex(node)[2]-fix_node[2]) < (1e-6 * this->parameters.scale))
5932 * && (abs(cell->vertex(node)[0]-fix_node[0]) < (1e-6 * this->parameters.scale)))
5933 * constraints.add_line(cell->vertex_dof_index(node, 0));
5935 * if ( (abs(cell->vertex(node)[2]-fix_node[2]) < (1e-6 * this->parameters.scale))
5936 * && (abs(cell->vertex(node)[1]-fix_node[1]) < (1e-6 * this->parameters.scale)))
5937 * constraints.add_line(cell->vertex_dof_index(node, 1));
5940 * if (this->parameters.load_type == "displacement")
5942 * const std::vector<double> value = get_dirichlet_load(5,2);
5943 * FEValuesExtractors::Scalar direction;
5944 * direction = this->z_displacement;
5946 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
5948 * ConstantFunction<dim>(value[2],this->n_components),
5950 * this->fe.component_mask(direction));
5954 * virtual Tensor<1,dim>
5955 * get_neumann_traction (const types::boundary_id &boundary_id,
5956 * const Point<dim> &pt,
5957 * const Tensor<1,dim> &N) const
5959 * if (this->parameters.load_type == "pressure")
5961 * if (boundary_id == 5)
5963 * const double final_load = this->parameters.load;
5964 * const double current_time = this->time.get_current();
5965 * const double final_time = this->time.get_end();
5966 * const double num_cycles = 3.0;
5968 * return final_load/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5))) * N;
5974 * return Tensor<1,dim>();
5977 * virtual types::boundary_id
5978 * get_reaction_boundary_id_for_output() const
5983 * virtual std::vector<double>
5984 * get_dirichlet_load(const types::boundary_id &boundary_id,
5985 * const int &direction) const
5987 * std::vector<double> displ_incr(dim,0.0);
5989 * if ( (boundary_id == 5) && (direction == 2) )
5991 * const double final_displ = this->parameters.load;
5992 * const double current_time = this->time.get_current();
5993 * const double final_time = this->time.get_end();
5994 * const double delta_time = this->time.get_delta_t();
5995 * const double num_cycles = 3.0;
5996 * double current_displ = 0.0;
5997 * double previous_displ = 0.0;
5999 * if (this->parameters.num_cycle_sets == 1)
6001 * current_displ = final_displ/2.0 * (1.0
6002 * - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5)));
6003 * previous_displ = final_displ/2.0 * (1.0
6004 * - std::sin(numbers::PI * (2.0*num_cycles*(current_time-delta_time)/final_time + 0.5)));
6008 * if ( current_time <= (final_time*1.0/3.0) )
6010 * current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6011 * (2.0*num_cycles*current_time/(final_time*1.0/3.0) + 0.5)));
6012 * previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6013 * (2.0*num_cycles*(current_time-delta_time)/(final_time*1.0/3.0) + 0.5)));
6017 * current_displ = final_displ * (1.0 - std::sin(numbers::PI *
6018 * (2.0*num_cycles*current_time / (final_time*2.0/3.0)
6019 * - (num_cycles - 0.5) )));
6020 * previous_displ = final_displ * (1.0 - std::sin(numbers::PI *
6021 * (2.0*num_cycles*(current_time-delta_time) / (final_time*2.0/3.0)
6022 * - (num_cycles - 0.5))));
6025 * displ_incr[2] = current_displ - previous_displ;
6027 * return displ_incr;
6034 * <a name="DerivedclassNolateraldisplacementinloadingsurfaces"></a>
6035 * <h4>Derived class: No lateral displacement in loading surfaces</h4>
6038 * template <int dim>
6039 * class BrainBudday2017CubeTensionCompressionFullyFixed
6040 * : public BrainBudday2017BaseCube<dim>
6043 * BrainBudday2017CubeTensionCompressionFullyFixed (const Parameters::AllParameters ¶meters)
6044 * : BrainBudday2017BaseCube<dim> (parameters)
6047 * virtual ~BrainBudday2017CubeTensionCompressionFullyFixed () {}
6051 * define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices)
6053 * tracked_vertices[0][0] = 0.5*this->parameters.scale;
6054 * tracked_vertices[0][1] = 0.5*this->parameters.scale;
6055 * tracked_vertices[0][2] = 1.0*this->parameters.scale;
6057 * tracked_vertices[1][0] = 0.5*this->parameters.scale;
6058 * tracked_vertices[1][1] = 0.5*this->parameters.scale;
6059 * tracked_vertices[1][2] = 0.5*this->parameters.scale;
6063 * make_dirichlet_constraints(AffineConstraints<double> &constraints)
6065 * if (this->time.get_timestep() < 2)
6067 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
6069 * ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
6071 * (this->fe.component_mask(this->pressure)));
6075 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6077 * ZeroFunction<dim>(this->n_components),
6079 * (this->fe.component_mask(this->pressure)));
6082 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6084 * ZeroFunction<dim>(this->n_components),
6086 * (this->fe.component_mask(this->x_displacement) |
6087 * this->fe.component_mask(this->y_displacement) |
6088 * this->fe.component_mask(this->z_displacement) ));
6091 * if (this->parameters.load_type == "displacement")
6093 * const std::vector<double> value = get_dirichlet_load(5,2);
6094 * FEValuesExtractors::Scalar direction;
6095 * direction = this->z_displacement;
6097 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6099 * ConstantFunction<dim>(value[2],this->n_components),
6101 * this->fe.component_mask(direction) );
6103 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6105 * ZeroFunction<dim>(this->n_components),
6107 * (this->fe.component_mask(this->x_displacement) |
6108 * this->fe.component_mask(this->y_displacement) ));
6112 * virtual Tensor<1,dim>
6113 * get_neumann_traction (const types::boundary_id &boundary_id,
6114 * const Point<dim> &pt,
6115 * const Tensor<1,dim> &N) const
6117 * if (this->parameters.load_type == "pressure")
6119 * if (boundary_id == 5)
6121 * const double final_load = this->parameters.load;
6122 * const double current_time = this->time.get_current();
6123 * const double final_time = this->time.get_end();
6124 * const double num_cycles = 3.0;
6126 * return final_load/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5))) * N;
6132 * return Tensor<1,dim>();
6135 * virtual types::boundary_id
6136 * get_reaction_boundary_id_for_output() const
6141 * virtual std::vector<double>
6142 * get_dirichlet_load(const types::boundary_id &boundary_id,
6143 * const int &direction) const
6145 * std::vector<double> displ_incr(dim,0.0);
6147 * if ( (boundary_id == 5) && (direction == 2) )
6149 * const double final_displ = this->parameters.load;
6150 * const double current_time = this->time.get_current();
6151 * const double final_time = this->time.get_end();
6152 * const double delta_time = this->time.get_delta_t();
6153 * const double num_cycles = 3.0;
6154 * double current_displ = 0.0;
6155 * double previous_displ = 0.0;
6157 * if (this->parameters.num_cycle_sets == 1)
6159 * current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*current_time/final_time + 0.5)));
6160 * previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI * (2.0*num_cycles*(current_time-delta_time)/final_time + 0.5)));
6164 * if ( current_time <= (final_time*1.0/3.0) )
6166 * current_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6167 * (2.0*num_cycles*current_time/(final_time*1.0/3.0) + 0.5)));
6168 * previous_displ = final_displ/2.0 * (1.0 - std::sin(numbers::PI *
6169 * (2.0*num_cycles*(current_time-delta_time)/(final_time*1.0/3.0) + 0.5)));
6173 * current_displ = final_displ * (1.0 - std::sin(numbers::PI *
6174 * (2.0*num_cycles*current_time / (final_time*2.0/3.0)
6175 * - (num_cycles - 0.5) )));
6176 * previous_displ = final_displ * (1.0 - std::sin(numbers::PI *
6177 * (2.0*num_cycles*(current_time-delta_time) / (final_time*2.0/3.0)
6178 * - (num_cycles - 0.5))));
6181 * displ_incr[2] = current_displ - previous_displ;
6183 * return displ_incr;
6190 * <a name="DerivedclassNolateralorverticaldisplacementinloadingsurface"></a>
6191 * <h4>Derived class: No lateral or vertical displacement in loading surface</h4>
6194 * template <int dim>
6195 * class BrainBudday2017CubeShearFullyFixed
6196 * : public BrainBudday2017BaseCube<dim>
6199 * BrainBudday2017CubeShearFullyFixed (const Parameters::AllParameters ¶meters)
6200 * : BrainBudday2017BaseCube<dim> (parameters)
6203 * virtual ~BrainBudday2017CubeShearFullyFixed () {}
6207 * define_tracked_vertices(std::vector<Point<dim> > &tracked_vertices)
6209 * tracked_vertices[0][0] = 0.75*this->parameters.scale;
6210 * tracked_vertices[0][1] = 0.5*this->parameters.scale;
6211 * tracked_vertices[0][2] = 0.0*this->parameters.scale;
6213 * tracked_vertices[1][0] = 0.25*this->parameters.scale;
6214 * tracked_vertices[1][1] = 0.5*this->parameters.scale;
6215 * tracked_vertices[1][2] = 0.0*this->parameters.scale;
6219 * make_dirichlet_constraints(AffineConstraints<double> &constraints)
6221 * if (this->time.get_timestep() < 2)
6223 * VectorTools::interpolate_boundary_values(this->dof_handler_ref,
6225 * ConstantFunction<dim>(this->parameters.drained_pressure,this->n_components),
6227 * (this->fe.component_mask(this->pressure)));
6231 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6233 * ZeroFunction<dim>(this->n_components),
6235 * (this->fe.component_mask(this->pressure)));
6238 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6240 * ZeroFunction<dim>(this->n_components),
6242 * (this->fe.component_mask(this->x_displacement) |
6243 * this->fe.component_mask(this->y_displacement) |
6244 * this->fe.component_mask(this->z_displacement) ));
6247 * if (this->parameters.load_type == "displacement")
6249 * const std::vector<double> value = get_dirichlet_load(4,0);
6250 * FEValuesExtractors::Scalar direction;
6251 * direction = this->x_displacement;
6253 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6255 * ConstantFunction<dim>(value[0],this->n_components),
6257 * this->fe.component_mask(direction));
6259 * VectorTools::interpolate_boundary_values( this->dof_handler_ref,
6261 * ZeroFunction<dim>(this->n_components),
6263 * (this->fe.component_mask(this->y_displacement) |
6264 * this->fe.component_mask(this->z_displacement) ));
6268 * virtual Tensor<1,dim>
6269 * get_neumann_traction (const types::boundary_id &boundary_id,
6270 * const Point<dim> &pt,
6271 * const Tensor<1,dim> &N) const
6273 * if (this->parameters.load_type == "pressure")
6275 * if (boundary_id == 4)
6277 * const double final_load = this->parameters.load;
6278 * const double current_time = this->time.get_current();
6279 * const double final_time = this->time.get_end();
6280 * const double num_cycles = 3.0;
6281 * const Point< 3, double> axis (0.0,1.0,0.0);
6282 * const double angle = numbers::PI;
6283 * static const Tensor< 2, dim, double> R(Physics::Transformations::Rotations::rotation_matrix_3d(axis,angle));
6285 * return (final_load * (std::sin(2.0*(numbers::PI)*num_cycles*current_time/final_time)) * (R * N));
6291 * return Tensor<1,dim>();
6294 * virtual types::boundary_id
6295 * get_reaction_boundary_id_for_output() const
6300 * virtual std::vector<double>
6301 * get_dirichlet_load(const types::boundary_id &boundary_id,
6302 * const int &direction) const
6304 * std::vector<double> displ_incr (dim, 0.0);
6306 * if ( (boundary_id == 4) && (direction == 0) )
6308 * const double final_displ = this->parameters.load;
6309 * const double current_time = this->time.get_current();
6310 * const double final_time = this->time.get_end();
6311 * const double delta_time = this->time.get_delta_t();
6312 * const double num_cycles = 3.0;
6313 * double current_displ = 0.0;
6314 * double previous_displ = 0.0;
6316 * if (this->parameters.num_cycle_sets == 1)
6318 * current_displ = final_displ * (std::sin(2.0*(numbers::PI)*num_cycles*current_time/final_time));
6319 * previous_displ = final_displ * (std::sin(2.0*(numbers::PI)*num_cycles*(current_time-delta_time)/final_time));
6323 * AssertThrow(false, ExcMessage("Problem type not defined. Budday shear experiments implemented only for one set of cycles."));
6325 * displ_incr[0] = current_displ - previous_displ;
6327 * return displ_incr;
6336 * <a name="Mainfunction"></a>
6337 * <h3>Main function</h3>
6338 * Lastly we provide the main driver function which is similar to the other tutorials.
6341 * int main (int argc, char *argv[])
6343 * using namespace dealii;
6344 * using namespace NonLinearPoroViscoElasticity;
6346 * const unsigned int n_tbb_processes = 1;
6347 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, n_tbb_processes);
6351 * Parameters::AllParameters parameters ("parameters.prm");
6352 * if (parameters.geom_type == "Ehlers_tube_step_load")
6354 * VerificationEhlers1999StepLoad<3> solid_3d(parameters);
6357 * else if (parameters.geom_type == "Ehlers_tube_increase_load")
6359 * VerificationEhlers1999IncreaseLoad<3> solid_3d(parameters);
6362 * else if (parameters.geom_type == "Ehlers_cube_consolidation")
6364 * VerificationEhlers1999CubeConsolidation<3> solid_3d(parameters);
6367 * else if (parameters.geom_type == "Franceschini_consolidation")
6369 * Franceschini2006Consolidation<3> solid_3d(parameters);
6372 * else if (parameters.geom_type == "Budday_cube_tension_compression")
6374 * BrainBudday2017CubeTensionCompression<3> solid_3d(parameters);
6377 * else if (parameters.geom_type == "Budday_cube_tension_compression_fully_fixed")
6379 * BrainBudday2017CubeTensionCompressionFullyFixed<3> solid_3d(parameters);
6382 * else if (parameters.geom_type == "Budday_cube_shear_fully_fixed")
6384 * BrainBudday2017CubeShearFullyFixed<3> solid_3d(parameters);
6389 * AssertThrow(false, ExcMessage("Problem type not defined. Current setting: " + parameters.geom_type));
6393 * catch (std::exception &exc)
6395 * if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
6397 * std::cerr << std::endl << std::endl
6398 * << "----------------------------------------------------"
6400 * std::cerr << "Exception on processing: " << std::endl << exc.what()
6401 * << std::endl << "Aborting!" << std::endl
6402 * << "----------------------------------------------------"
6410 * if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
6412 * std::cerr << std::endl << std::endl
6413 * << "----------------------------------------------------"
6415 * std::cerr << "Unknown exception!" << std::endl << "Aborting!"
6417 * << "----------------------------------------------------"
typename DataOut_DoFData< DoFHandlerType, dimension - 1, dimension >::cell_iterator cell_iterator
typename DataOut_DoFData< DoFHandlerType, DoFHandlerType::dimension, DoFHandlerType::space_dimension >::cell_iterator cell_iterator
size_type n_elements() const
IndexSet get_view(const size_type begin, const size_type end) const
size_type nth_index_in_set(const size_type local_index) const
virtual void parse_input(std::istream &input, const std::string &filename="input file", const std::string &last_line="", const bool skip_undefined=false)
numbers::NumberTraits< Number >::real_type norm() const
void leave_subsection(const std::string §ion_name="")
void enter_subsection(const std::string §ion_name)
@ 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.
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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 Matrix &matrix)
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
real_type l2_norm() const
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
void compress(::VectorOperation::values operation)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void approximate(SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
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.
static const types::blas_int zero
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
static const 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.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
std::enable_if< IsBlockVector< VectorType >::value, unsignedint >::type 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 this_mpi_process(const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void abort(const ExceptionBase &exc) noexcept
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 > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
SymmetricTensorEigenvectorMethod