273 * <a name=
"cook_membrane.cc-Runtimeparameters"></a>
287 * <a name=
"cook_membrane.cc-Assemblymethod"></a>
311 *
prm.enter_subsection(
"Assembly method");
313 *
prm.declare_entry(
"Automatic differentiation order",
"0",
315 *
"The automatic differentiation order to be used in the assembly of the linear system.\n"
316 *
"# Order = 0: Both the residual and linearisation are computed manually.\n"
317 *
"# Order = 1: The residual is computed manually but the linearisation is performed using AD.\n"
318 *
"# Order = 2: Both the residual and linearisation are computed using AD.");
320 *
prm.leave_subsection();
325 *
prm.enter_subsection(
"Assembly method");
329 *
prm.leave_subsection();
335 * <a name=
"cook_membrane.cc-FiniteElementsystem"></a>
359 *
prm.enter_subsection(
"Finite element system");
361 *
prm.declare_entry(
"Polynomial degree",
"2",
363 *
"Displacement system polynomial order");
365 *
prm.declare_entry(
"Quadrature order",
"3",
367 *
"Gauss quadrature order");
369 *
prm.leave_subsection();
374 *
prm.enter_subsection(
"Finite element system");
376 *
poly_degree = prm.get_integer(
"Polynomial degree");
377 *
quad_order = prm.get_integer(
"Quadrature order");
379 *
prm.leave_subsection();
385 * <a name=
"cook_membrane.cc-Geometry"></a>
407 *
prm.enter_subsection(
"Geometry");
409 *
prm.declare_entry(
"Elements per edge",
"32",
411 *
"Number of elements per long edge of the beam");
413 *
prm.declare_entry(
"Grid scale",
"1e-3",
415 *
"Global grid scaling factor");
417 *
prm.leave_subsection();
422 *
prm.enter_subsection(
"Geometry");
425 *
scale = prm.get_double(
"Grid scale");
427 *
prm.leave_subsection();
433 * <a name=
"cook_membrane.cc-Materials"></a>
456 *
prm.enter_subsection(
"Material properties");
458 *
prm.declare_entry(
"Poisson's ratio",
"0.3",
460 *
"Poisson's ratio");
462 *
prm.declare_entry(
"Shear modulus",
"0.4225e6",
466 *
prm.leave_subsection();
471 *
prm.enter_subsection(
"Material properties");
473 *
nu = prm.get_double(
"Poisson's ratio");
474 *
mu = prm.get_double(
"Shear modulus");
476 *
prm.leave_subsection();
482 * <a name=
"cook_membrane.cc-Linearsolver"></a>
509 *
prm.enter_subsection(
"Linear solver");
511 *
prm.declare_entry(
"Solver type",
"CG",
513 *
"Type of solver used to solve the linear system");
515 *
prm.declare_entry(
"Residual",
"1e-6",
517 *
"Linear solver residual (scaled by residual norm)");
519 *
prm.declare_entry(
"Max iteration multiplier",
"1",
521 *
"Linear solver iterations (multiples of the system matrix size)");
523 *
prm.declare_entry(
"Preconditioner type",
"ssor",
525 *
"Type of preconditioner");
527 *
prm.declare_entry(
"Preconditioner relaxation",
"0.65",
529 *
"Preconditioner relaxation value");
531 *
prm.leave_subsection();
536 *
prm.enter_subsection(
"Linear solver");
538 *
type_lin = prm.get(
"Solver type");
539 *
tol_lin = prm.get_double(
"Residual");
544 *
prm.leave_subsection();
550 * <a name=
"cook_membrane.cc-Nonlinearsolver"></a>
560 *
struct NonlinearSolver
575 *
prm.enter_subsection(
"Nonlinear solver");
577 *
prm.declare_entry(
"Max iterations Newton-Raphson",
"10",
579 *
"Number of Newton-Raphson iterations allowed");
581 *
prm.declare_entry(
"Tolerance force",
"1.0e-9",
583 *
"Force residual tolerance");
585 *
prm.declare_entry(
"Tolerance displacement",
"1.0e-6",
587 *
"Displacement error tolerance");
589 *
prm.leave_subsection();
594 *
prm.enter_subsection(
"Nonlinear solver");
597 *
tol_f = prm.get_double(
"Tolerance force");
598 *
tol_u = prm.get_double(
"Tolerance displacement");
600 *
prm.leave_subsection();
606 * <a name=
"cook_membrane.cc-Time"></a>
628 *
prm.enter_subsection(
"Time");
630 *
prm.declare_entry(
"End time",
"1",
634 *
prm.declare_entry(
"Time step size",
"0.1",
638 *
prm.leave_subsection();
643 *
prm.enter_subsection(
"Time");
645 *
end_time = prm.get_double(
"End time");
646 *
delta_t = prm.get_double(
"Time step size");
648 *
prm.leave_subsection();
654 * <a name=
"cook_membrane.cc-Allparameters"></a>
669 *
public NonlinearSolver,
685 *
declare_parameters(prm);
687 *
parse_parameters(prm);
692 *
AssemblyMethod::declare_parameters(prm);
693 *
FESystem::declare_parameters(prm);
694 *
Geometry::declare_parameters(prm);
695 *
Materials::declare_parameters(prm);
696 *
LinearSolver::declare_parameters(prm);
697 *
NonlinearSolver::declare_parameters(prm);
698 *
Time::declare_parameters(prm);
703 *
AssemblyMethod::parse_parameters(prm);
704 *
FESystem::parse_parameters(prm);
705 *
Geometry::parse_parameters(prm);
706 *
Materials::parse_parameters(prm);
707 *
LinearSolver::parse_parameters(prm);
708 *
NonlinearSolver::parse_parameters(prm);
709 *
Time::parse_parameters(prm);
717 * <a name=
"cook_membrane.cc-Timeclass"></a>
731 *
const double delta_t)
768 *
const double delta_t;
774 * <a name=
"cook_membrane.cc-CompressibleneoHookeanmaterialwithinaonefieldformulation"></a>
799 *
href=
"http://en.wikipedia.org/wiki/Lam%C3%A9_parameters">
Lame's first
804 * The following class will be used to characterize the material we work with,
805 * and provides a central point that one would need to modify if one were to
806 * implement a different material model. For it to work, we will store one
807 * object of this type per quadrature point, and in each of these objects
808 * store the current state (characterized by the values or measures of the
809 * displacement field) so that we can compute the elastic coefficients
810 * linearized around the current state.
813 * template <int dim,typename NumberType>
814 * class Material_Compressible_Neo_Hook_One_Field
817 * Material_Compressible_Neo_Hook_One_Field(const double mu,
820 * kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
823 * Assert(kappa > 0, ExcInternalError());
826 * ~Material_Compressible_Neo_Hook_One_Field()
831 * The first function is the total energy
832 * @f$\Psi = \Psi_{\textrm{iso}} + \Psi_{\textrm{vol}}@f$.
836 * get_Psi(const NumberType &det_F,
837 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
839 * return get_Psi_vol(det_F) + get_Psi_iso(b_bar);
844 * The second function determines the Kirchhoff stress @f$\boldsymbol{\tau}
845 * = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}@f$
848 * SymmetricTensor<2,dim,NumberType>
849 * get_tau(const NumberType &det_F,
850 * const SymmetricTensor<2,dim,NumberType> &b_bar)
854 * See Holzapfel p231 eq6.98 onwards
857 * return get_tau_vol(det_F) + get_tau_iso(b_bar);
862 * The fourth-order elasticity tensor in the spatial setting
863 * @f$\mathfrak{c}@f$ is calculated from the SEF @f$\Psi@f$ as @f$ J
864 * \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}@f$
865 * where @f$ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial
866 * \mathbf{C} \partial \mathbf{C}}@f$
869 * SymmetricTensor<4,dim,NumberType>
870 * get_Jc(const NumberType &det_F,
871 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
873 * return get_Jc_vol(det_F) + get_Jc_iso(b_bar);
879 * Define constitutive model parameters @f$\kappa@f$ (bulk modulus) and the
880 * neo-Hookean model parameter @f$c_1@f$:
883 * const double kappa;
888 * Value of the volumetric free energy
892 * get_Psi_vol(const NumberType &det_F) const
894 * return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
899 * Value of the isochoric free energy
903 * get_Psi_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
905 * return c_1 * (trace(b_bar) - dim);
910 * Derivative of the volumetric free energy with respect to
911 * @f$J@f$ return @f$\frac{\partial
912 * \Psi_{\text{vol}}(J)}{\partial J}@f$
916 * get_dPsi_vol_dJ(const NumberType &det_F) const
918 * return (kappa / 2.0) * (det_F - 1.0 / det_F);
923 * The following functions are used internally in determining the result
924 * of some of the public functions above. The first one determines the
925 * volumetric Kirchhoff stress @f$\boldsymbol{\tau}_{\textrm{vol}}@f$.
926 * Note the difference in its definition when compared to @ref step_44 "step-44".
929 * SymmetricTensor<2,dim,NumberType>
930 * get_tau_vol(const NumberType &det_F) const
932 * return NumberType(get_dPsi_vol_dJ(det_F) * det_F) * Physics::Elasticity::StandardTensors<dim>::I;
937 * Next, determine the isochoric Kirchhoff stress
938 * @f$\boldsymbol{\tau}_{\textrm{iso}} =
939 * \mathcal{P}:\overline{\boldsymbol{\tau}}@f$:
942 * SymmetricTensor<2,dim,NumberType>
943 * get_tau_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
945 * return Physics::Elasticity::StandardTensors<dim>::dev_P * get_tau_bar(b_bar);
950 * Then, determine the fictitious Kirchhoff stress
951 * @f$\overline{\boldsymbol{\tau}}@f$:
954 * SymmetricTensor<2,dim,NumberType>
955 * get_tau_bar(const SymmetricTensor<2,dim,NumberType> &b_bar) const
957 * return 2.0 * c_1 * b_bar;
962 * Second derivative of the volumetric free energy wrt @f$J@f$. We
963 * need the following computation explicitly in the tangent so we make it
964 * public. We calculate @f$\frac{\partial^2
965 * \Psi_{\textrm{vol}}(J)}{\partial J \partial
970 * get_d2Psi_vol_dJ2(const NumberType &det_F) const
972 * return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
977 * Calculate the volumetric part of the tangent @f$J
978 * \mathfrak{c}_\textrm{vol}@f$. Again, note the difference in its
979 * definition when compared to @ref step_44 "step-44". The extra terms result from two
980 * quantities in @f$\boldsymbol{\tau}_{\textrm{vol}}@f$ being dependent on
981 * @f$\boldsymbol{F}@f$.
984 * SymmetricTensor<4,dim,NumberType>
985 * get_Jc_vol(const NumberType &det_F) const
993 * * ( (get_dPsi_vol_dJ(det_F) + det_F * get_d2Psi_vol_dJ2(det_F))*Physics::Elasticity::StandardTensors<dim>::IxI
994 * - (2.0 * get_dPsi_vol_dJ(det_F))*Physics::Elasticity::StandardTensors<dim>::S );
999 * Calculate the isochoric part of the tangent @f$J
1000 * \mathfrak{c}_\textrm{iso}@f$:
1003 * SymmetricTensor<4,dim,NumberType>
1004 * get_Jc_iso(const SymmetricTensor<2,dim,NumberType> &b_bar) const
1006 * const SymmetricTensor<2, dim> tau_bar = get_tau_bar(b_bar);
1007 * const SymmetricTensor<2, dim> tau_iso = get_tau_iso(b_bar);
1008 * const SymmetricTensor<4, dim> tau_iso_x_I
1009 * = outer_product(tau_iso,
1010 * Physics::Elasticity::StandardTensors<dim>::I);
1011 * const SymmetricTensor<4, dim> I_x_tau_iso
1012 * = outer_product(Physics::Elasticity::StandardTensors<dim>::I,
1014 * const SymmetricTensor<4, dim> c_bar = get_c_bar();
1016 * return (2.0 / dim) * trace(tau_bar)
1017 * * Physics::Elasticity::StandardTensors<dim>::dev_P
1018 * - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso)
1019 * + Physics::Elasticity::StandardTensors<dim>::dev_P * c_bar
1020 * * Physics::Elasticity::StandardTensors<dim>::dev_P;
1025 * Calculate the fictitious elasticity tensor @f$\overline{\mathfrak{c}}@f$.
1026 * For the material model chosen this is simply zero:
1029 * SymmetricTensor<4,dim,double>
1032 * return SymmetricTensor<4, dim>();
1039 * <a name="cook_membrane.cc-Quadraturepointhistory"></a>
1040 * <h3>Quadrature point history</h3>
1044 * As seen in @ref step_18 "step-18", the <code> PointHistory </code> class offers a method
1045 * for storing data at the quadrature points. Here each quadrature point
1046 * holds a pointer to a material description. Thus, different material models
1047 * can be used in different regions of the domain. Among other data, we
1048 * choose to store the Kirchhoff stress @f$\boldsymbol{\tau}@f$ and the tangent
1049 * @f$J\mathfrak{c}@f$ for the quadrature points.
1052 * template <int dim,typename NumberType>
1053 * class PointHistory
1059 * virtual ~PointHistory()
1064 * The first function is used to create a material object and to
1065 * initialize all tensors correctly: The second one updates the stored
1066 * values and stresses based on the current deformation measure
1067 * @f$\textrm{Grad}\mathbf{u}_{\textrm{n}}@f$.
1070 * void setup_lqp (const Parameters::AllParameters ¶meters)
1072 * material.reset(new Material_Compressible_Neo_Hook_One_Field<dim,NumberType>(parameters.mu,
1078 * We offer an interface to retrieve certain data.
1079 * This is the strain energy:
1083 * get_Psi(const NumberType &det_F,
1084 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1086 * return material->get_Psi(det_F,b_bar);
1091 * Here are the kinetic variables. These are used in the material and
1092 * global tangent matrix and residual assembly operations:
1093 * First is the Kirchhoff stress:
1096 * SymmetricTensor<2,dim,NumberType>
1097 * get_tau(const NumberType &det_F,
1098 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1100 * return material->get_tau(det_F,b_bar);
1108 * SymmetricTensor<4,dim,NumberType>
1109 * get_Jc(const NumberType &det_F,
1110 * const SymmetricTensor<2,dim,NumberType> &b_bar) const
1112 * return material->get_Jc(det_F,b_bar);
1117 * In terms of member functions, this class stores for the quadrature
1118 * point it represents a copy of a material type in case different
1119 * materials are used in different regions of the domain, as well as the
1120 * inverse of the deformation gradient...
1124 * std::shared_ptr< Material_Compressible_Neo_Hook_One_Field<dim,NumberType> > material;
1131 * <a name="cook_membrane.cc-Quasistaticcompressiblefinitestrainsolid"></a>
1132 * <h3>Quasi-static compressible finite-strain solid</h3>
1136 * Forward declarations for classes that will
1137 * perform assembly of the linear system.
1140 * template <int dim,typename NumberType>
1141 * struct Assembler_Base;
1142 * template <int dim,typename NumberType>
1147 * The Solid class is the central class in that it represents the problem at
1148 * hand. It follows the usual scheme in that all it really has is a
1149 * constructor, destructor and a <code>run()</code> function that dispatches
1150 * all the work to private functions of this class:
1153 * template <int dim,typename NumberType>
1157 * Solid(const Parameters::AllParameters ¶meters);
1169 * We start the collection of member functions with one that builds the
1178 * Set up the finite element system to be solved:
1186 * Several functions to assemble the system and right hand side matrices
1187 * using multithreading. Each of them comes as a wrapper function, one
1188 * that is executed to do the work in the WorkStream model on one cell,
1189 * and one that copies the work done on this one cell into the global
1190 * object that represents it:
1194 * assemble_system(const BlockVector<double> &solution_delta);
1198 * We use a separate data structure to perform the assembly. It needs access
1199 * to some low-level data, so we simply befriend the class instead of
1200 * creating a complex interface to provide access as necessary.
1203 * friend struct Assembler_Base<dim,NumberType>;
1204 * friend struct Assembler<dim,NumberType>;
1208 * Apply Dirichlet boundary conditions on the displacement field
1212 * make_constraints(const int &it_nr);
1216 * Create and update the quadrature points. Here, no data needs to be
1217 * copied into a global object, so the copy_local_to_global function is
1226 * Solve for the displacement using a Newton-Raphson method. We break this
1227 * function into the nonlinear loop and the function that solves the
1228 * linearized Newton-Raphson step:
1232 * solve_nonlinear_timestep(BlockVector<double> &solution_delta);
1234 * std::pair<unsigned int, double>
1235 * solve_linear_system(BlockVector<double> &newton_update);
1239 * Solution retrieval as well as post-processing and writing data to file :
1242 * BlockVector<double>
1243 * get_total_solution(const BlockVector<double> &solution_delta) const;
1246 * output_results() const;
1250 * Finally, some member variables that describe the current state: A
1251 * collection of the parameters used to describe the problem setup...
1254 * const Parameters::AllParameters ¶meters;
1258 * ...the volume of the reference and current configurations...
1261 * double vol_reference;
1262 * double vol_current;
1266 * ...and description of the geometry on which the problem is solved:
1269 * Triangulation<dim> triangulation;
1273 * Also, keep track of the current time and the time spent evaluating
1278 * TimerOutput timer;
1282 * A storage object for quadrature point information. As opposed to
1297 *
const unsigned int degree;
1300 *
const unsigned int dofs_per_cell;
1314 *
static const unsigned int n_blocks = 1;
1315 *
static const unsigned int n_components = dim;
1333 *
const unsigned int n_q_points;
1413 * <a name=
"cook_membrane.cc-ImplementationofthecodeSolidcodeclass"></a>
1419 * <a name=
"cook_membrane.cc-Publicinterface"></a>
1427 *
template <
int dim,
typename NumberType>
1430 *
parameters(parameters),
1434 *
time(parameters.end_time, parameters.delta_t),
1438 *
degree(parameters.poly_degree),
1445 *
fe(
FE_Q<dim>(parameters.poly_degree), dim),
1447 *
dofs_per_cell (fe.dofs_per_cell),
1450 *
qf_cell(parameters.quad_order),
1451 *
qf_face(parameters.quad_order),
1452 *
n_q_points (
qf_cell.size()),
1463 *
template <
int dim,
typename NumberType>
1487 *
template <
int dim,
typename NumberType>
1507 *
while (time.current() <= time.end())
1544 * <a name=
"cook_membrane.cc-Privateinterface"></a>
1550 * <a name=
"cook_membrane.cc-Solidmake_grid"></a>
1570 *
template <
int dim>
1573 *
const double &
x =
pt_in[0];
1574 *
const double &
y =
pt_in[1];
1576 *
const double y_upper = 44.0 + (16.0/48.0)*
x;
1577 *
const double y_lower = 0.0 + (44.0/48.0)*
x;
1578 *
const double theta =
y/44.0;
1587 *
template <
int dim,
typename NumberType>
1595 *
std::vector< unsigned int >
repetitions(dim, parameters.elements_per_edge);
1627 *
for (; cell !=
endc; ++cell)
1628 *
for (
unsigned int face = 0;
1630 *
if (cell->face(face)->at_boundary() ==
true)
1633 *
cell->face(face)->set_boundary_id(1);
1635 *
cell->face(face)->set_boundary_id(11);
1637 *
cell->face(face)->set_boundary_id(2);
1651 *
std::cout <<
"Grid:\n\t Reference volume: " <<
vol_reference << std::endl;
1658 * <a name=
"cook_membrane.cc-Solidsystem_setup"></a>
1668 *
template <
int dim,
typename NumberType>
1671 *
timer.enter_subsection(
"Setup system");
1686 *
std::cout <<
"Triangulation:"
1687 *
<<
"\n\t Number of active cells: " <<
triangulation.n_active_cells()
1688 *
<<
"\n\t Number of degrees of freedom: " <<
dof_handler_ref.n_dofs()
1703 *
csp.collect_sizes();
1712 *
for (
unsigned int ii = 0;
ii < n_components; ++
ii)
1713 *
for (
unsigned int jj = 0;
jj < n_components; ++
jj)
1720 *
sparsity_pattern.copy_from(
csp);
1738 * ...and
finally set
up the quadrature
1744 *
timer.leave_subsection();
1751 * <a name=
"cook_membrane.cc-Solidsetup_qph"></a>
1762 *
template <
int dim,
typename NumberType>
1765 *
std::cout <<
" Setting up quadrature point data..." << std::endl;
1781 *
const std::vector<std::shared_ptr<PointHistory<dim,NumberType> > >
lqph =
1794 * <a name=
"cook_membrane.cc-Solidsolve_nonlinear_timestep"></a>
1804 *
template <
int dim,
typename NumberType>
1808 *
std::cout << std::endl <<
"Timestep " << time.get_timestep() <<
" @ "
1809 *
<< time.current() <<
"s" << std::endl;
1841 *
std::cout <<
" " << std::setw(2) <<
newton_iteration <<
" " << std::flush;
1870 *
std::cout <<
" CONVERGED! " << std::endl;
1876 *
const std::pair<unsigned int, double>
1896 *
std::cout <<
" | " << std::fixed << std::setprecision(3) << std::setw(7)
1901 *
<<
" " << std::endl;
1937 *
static const unsigned int l_width = 87;
1939 *
for (
unsigned int i = 0; i <
l_width; ++i)
1941 *
std::cout << std::endl;
1943 *
std::cout <<
" SOLVER STEP "
1944 *
<<
" | LIN_IT LIN_RES RES_NORM "
1945 *
<<
" RES_U NU_NORM "
1946 *
<<
" NU_U " << std::endl;
1948 *
for (
unsigned int i = 0; i <
l_width; ++i)
1950 *
std::cout << std::endl;
1955 *
template <
int dim,
typename NumberType>
1958 *
static const unsigned int l_width = 87;
1960 *
for (
unsigned int i = 0; i <
l_width; ++i)
1962 *
std::cout << std::endl;
1964 *
std::cout <<
"Relative errors:" << std::endl
1978 *
template <
int dim,
typename NumberType>
1981 *
static const unsigned int l_width = 87;
1983 *
for (
unsigned int i = 0; i <
l_width; ++i)
1985 *
std::cout << std::endl;
1994 *
Point<dim>(48.0*parameters.scale, 52.0*parameters.scale, 0.5*parameters.scale) :
2001 *
for (; cell !=
endc; ++cell)
2005 *
if (cell->point_inside(
soln_pt) ==
true)
2009 *
if (cell->vertex(v).distance(
soln_pt) < 1e-6)
2016 * that have support at the vertices
2019 * vertical_tip_displacement = solution_n(cell->vertex_dof_index(v,u_dof+1));
2023 * Sanity check using alternate method to extract the solution
2024 * at the given point. To do this, we must create an FEValues instance
2025 * to help us extract the solution value at the desired point
2028 * const MappingQ<dim> mapping (parameters.poly_degree);
2029 * const Point<dim> qp_unit = mapping.transform_real_to_unit_cell(cell,soln_pt);
2030 * const Quadrature<dim> soln_qrule (qp_unit);
2031 * AssertThrow(soln_qrule.size() == 1, ExcInternalError());
2032 * FEValues<dim> fe_values_soln (fe, soln_qrule, update_values);
2033 * fe_values_soln.reinit(cell);
2037 * Extract y-component of solution at given point
2040 * std::vector< Tensor<1,dim> > soln_values (soln_qrule.size());
2041 * fe_values_soln[u_fe].get_function_values(solution_n,
2043 * vertical_tip_displacement_check = soln_values[0][u_dof+1];
2048 * AssertThrow(vertical_tip_displacement > 0.0, ExcMessage("Found no cell with point inside!"));
2050 * std::cout << "Vertical tip displacement: " << vertical_tip_displacement
2051 * << "\t Check: " << vertical_tip_displacement_check
2059 * <a name="cook_membrane.cc-Solidget_error_residual"></a>
2060 * <h4>Solid::get_error_residual</h4>
2064 * Determine the true residual error for the problem. That is, determine the
2065 * error in the residual for the unconstrained degrees of freedom. Note that to
2066 * do so, we need to ignore constrained DOFs by setting the residual in these
2067 * vector components to zero.
2070 * template <int dim,typename NumberType>
2071 * void Solid<dim,NumberType>::get_error_residual(Errors &error_residual)
2073 * BlockVector<double> error_res(dofs_per_block);
2075 * for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
2076 * if (!constraints.is_constrained(i))
2077 * error_res(i) = system_rhs(i);
2079 * error_residual.norm = error_res.l2_norm();
2080 * error_residual.u = error_res.block(u_dof).l2_norm();
2087 * <a name="cook_membrane.cc-Solidget_error_udpate"></a>
2088 * <h4>Solid::get_error_udpate</h4>
2092 * Determine the true Newton update error for the problem
2095 * template <int dim,typename NumberType>
2096 * void Solid<dim,NumberType>::get_error_update(const BlockVector<double> &newton_update,
2097 * Errors &error_update)
2099 * BlockVector<double> error_ud(dofs_per_block);
2100 * for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
2101 * if (!constraints.is_constrained(i))
2102 * error_ud(i) = newton_update(i);
2104 * error_update.norm = error_ud.l2_norm();
2105 * error_update.u = error_ud.block(u_dof).l2_norm();
2113 * <a name="cook_membrane.cc-Solidget_total_solution"></a>
2114 * <h4>Solid::get_total_solution</h4>
2118 * This function provides the total solution, which is valid at any Newton step.
2119 * This is required as, to reduce computational error, the total solution is
2120 * only updated at the end of the timestep.
2123 * template <int dim,typename NumberType>
2124 * BlockVector<double>
2125 * Solid<dim,NumberType>::get_total_solution(const BlockVector<double> &solution_delta) const
2127 * BlockVector<double> solution_total(solution_n);
2128 * solution_total += solution_delta;
2129 * return solution_total;
2136 * <a name="cook_membrane.cc-Solidassemble_system"></a>
2137 * <h4>Solid::assemble_system</h4>
2143 * template <int dim,typename NumberType>
2144 * struct Assembler_Base
2146 * virtual ~Assembler_Base() {}
2150 * Here we deal with the tangent matrix assembly structures. The
2151 * PerTaskData object stores local contributions.
2154 * struct PerTaskData_ASM
2156 * const Solid<dim,NumberType> *solid;
2157 * FullMatrix<double> cell_matrix;
2158 * Vector<double> cell_rhs;
2159 * std::vector<types::global_dof_index> local_dof_indices;
2161 * PerTaskData_ASM(const Solid<dim,NumberType> *solid)
2164 * cell_matrix(solid->dofs_per_cell, solid->dofs_per_cell),
2165 * cell_rhs(solid->dofs_per_cell),
2166 * local_dof_indices(solid->dofs_per_cell)
2171 * cell_matrix = 0.0;
2178 * On the other hand, the ScratchData object stores the larger objects such as
2179 * the shape-function values array (<code>Nx</code>) and a shape function
2180 * gradient and symmetric gradient vector which we will use during the
2184 * struct ScratchData_ASM
2186 * const BlockVector<double> &solution_total;
2187 * std::vector<Tensor<2, dim,NumberType> > solution_grads_u_total;
2189 * FEValues<dim> fe_values_ref;
2190 * FEFaceValues<dim> fe_face_values_ref;
2192 * std::vector<std::vector<Tensor<2, dim,NumberType> > > grad_Nx;
2193 * std::vector<std::vector<SymmetricTensor<2,dim,NumberType> > > symm_grad_Nx;
2195 * ScratchData_ASM(const FiniteElement<dim> &fe_cell,
2196 * const QGauss<dim> &qf_cell,
2197 * const UpdateFlags uf_cell,
2198 * const QGauss<dim-1> & qf_face,
2199 * const UpdateFlags uf_face,
2200 * const BlockVector<double> &solution_total)
2202 * solution_total(solution_total),
2203 * solution_grads_u_total(qf_cell.size()),
2204 * fe_values_ref(fe_cell, qf_cell, uf_cell),
2205 * fe_face_values_ref(fe_cell, qf_face, uf_face),
2206 * grad_Nx(qf_cell.size(),
2207 * std::vector<Tensor<2,dim,NumberType> >(fe_cell.dofs_per_cell)),
2208 * symm_grad_Nx(qf_cell.size(),
2209 * std::vector<SymmetricTensor<2,dim,NumberType> >
2210 * (fe_cell.dofs_per_cell))
2213 * ScratchData_ASM(const ScratchData_ASM &rhs)
2215 * solution_total (rhs.solution_total),
2216 * solution_grads_u_total(rhs.solution_grads_u_total),
2217 * fe_values_ref(rhs.fe_values_ref.get_fe(),
2218 * rhs.fe_values_ref.get_quadrature(),
2219 * rhs.fe_values_ref.get_update_flags()),
2220 * fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
2221 * rhs.fe_face_values_ref.get_quadrature(),
2222 * rhs.fe_face_values_ref.get_update_flags()),
2223 * grad_Nx(rhs.grad_Nx),
2224 * symm_grad_Nx(rhs.symm_grad_Nx)
2229 * const unsigned int n_q_points = fe_values_ref.get_quadrature().size();
2230 * const unsigned int n_dofs_per_cell = fe_values_ref.dofs_per_cell;
2231 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2233 * Assert( grad_Nx[q_point].size() == n_dofs_per_cell,
2234 * ExcInternalError());
2235 * Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell,
2236 * ExcInternalError());
2238 * solution_grads_u_total[q_point] = Tensor<2,dim,NumberType>();
2239 * for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
2241 * grad_Nx[q_point][k] = Tensor<2,dim,NumberType>();
2242 * symm_grad_Nx[q_point][k] = SymmetricTensor<2,dim,NumberType>();
2251 * Of course, we still have to define how we assemble the tangent matrix
2252 * contribution for a single cell.
2256 * assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2257 * ScratchData_ASM &scratch,
2258 * PerTaskData_ASM &data)
2262 * Due to the C++ specialization rules, we need one more
2263 * level of indirection in order to define the assembly
2264 * routine for all different number. The next function call
2265 * is specialized for each NumberType, but to prevent having
2266 * to specialize the whole class along with it we have inlined
2267 * the definition of the other functions that are common to
2268 * all implementations.
2271 * assemble_system_tangent_residual_one_cell(cell, scratch, data);
2272 * assemble_neumann_contribution_one_cell(cell, scratch, data);
2277 * This function adds the local contribution to the system matrix.
2281 * copy_local_to_global_ASM(const PerTaskData_ASM &data)
2283 * const AffineConstraints<double> &constraints = data.solid->constraints;
2284 * BlockSparseMatrix<double> &tangent_matrix = const_cast<Solid<dim,NumberType> *>(data.solid)->tangent_matrix;
2285 * BlockVector<double> &system_rhs = const_cast<Solid<dim,NumberType> *>(data.solid)->system_rhs;
2287 * constraints.distribute_local_to_global(
2288 * data.cell_matrix, data.cell_rhs,
2289 * data.local_dof_indices,
2290 * tangent_matrix, system_rhs);
2297 * This function needs to exist in the base class for
2298 * Workstream to work with a reference to the base class.
2302 * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &/*cell*/,
2303 * ScratchData_ASM &/*scratch*/,
2304 * PerTaskData_ASM &/*data*/)
2306 * AssertThrow(false, ExcPureFunctionCalled());
2310 * assemble_neumann_contribution_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2311 * ScratchData_ASM &scratch,
2312 * PerTaskData_ASM &data)
2316 * Aliases for data referenced from the Solid class
2319 * const unsigned int &n_q_points_f = data.solid->n_q_points_f;
2320 * const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2321 * const Parameters::AllParameters ¶meters = data.solid->parameters;
2322 * const Time &time = data.solid->time;
2323 * const FESystem<dim> &fe = data.solid->fe;
2324 * const unsigned int &u_dof = data.solid->u_dof;
2328 * Next we assemble the Neumann contribution. We first check to see it the
2329 * cell face exists on a boundary on which a traction is applied and add
2330 * the contribution if this is the case.
2333 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2335 * if (cell->face(face)->at_boundary() == true
2336 * && cell->face(face)->boundary_id() == 11)
2338 * scratch.fe_face_values_ref.reinit(cell, face);
2340 * for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
2345 * We specify the traction in reference configuration.
2346 * For this problem, a defined total vertical force is applied
2347 * in the reference configuration.
2348 * The direction of the applied traction is assumed not to
2349 * evolve with the deformation of the domain.
2353 * Note that the contributions to the right hand side vector we
2354 * compute here only exist in the displacement components of the
2358 * const double time_ramp = (time.current() / time.end());
2359 * const double magnitude = (1.0/(16.0*parameters.scale*1.0*parameters.scale))*time_ramp; // (Total force) / (RHS surface area)
2360 * Tensor<1,dim> dir;
2362 * const Tensor<1, dim> traction = magnitude*dir;
2364 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2366 * const unsigned int i_group =
2367 * fe.system_to_base_index(i).first.first;
2369 * if (i_group == u_dof)
2371 * const unsigned int component_i =
2372 * fe.system_to_component_index(i).first;
2374 * scratch.fe_face_values_ref.shape_value(i,
2376 * const double JxW = scratch.fe_face_values_ref.JxW(
2379 * data.cell_rhs(i) += (Ni * traction[component_i])
2389 * template <int dim>
2390 * struct Assembler<dim,double> : Assembler_Base<dim,double>
2392 * typedef double NumberType;
2393 * using typename Assembler_Base<dim,NumberType>::ScratchData_ASM;
2394 * using typename Assembler_Base<dim,NumberType>::PerTaskData_ASM;
2396 * virtual ~Assembler() {}
2399 * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2400 * ScratchData_ASM &scratch,
2401 * PerTaskData_ASM &data) override
2405 * Aliases for data referenced from the Solid class
2408 * const unsigned int &n_q_points = data.solid->n_q_points;
2409 * const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2410 * const FESystem<dim> &fe = data.solid->fe;
2411 * const unsigned int &u_dof = data.solid->u_dof;
2412 * const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
2416 * scratch.fe_values_ref.reinit(cell);
2417 * cell->get_dof_indices(data.local_dof_indices);
2419 * const std::vector<std::shared_ptr<const PointHistory<dim,NumberType> > > lqph =
2420 * const_cast<const Solid<dim,NumberType> *>(data.solid)->quadrature_point_history.get_data(cell);
2421 * Assert(lqph.size() == n_q_points, ExcInternalError());
2425 * We first need to find the solution gradients at quadrature points
2426 * inside the current cell and then we update each local QP using the
2427 * displacement gradient:
2430 * scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
2431 * scratch.solution_grads_u_total);
2435 * Now we build the local cell stiffness matrix. Since the global and
2436 * local system matrices are symmetric, we can exploit this property by
2437 * building only the lower half of the local matrix and copying the values
2438 * to the upper half.
2442 * In doing so, we first extract some configuration dependent variables
2443 * from our QPH history objects for the current quadrature point.
2446 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2448 * const Tensor<2,dim,NumberType> &grad_u = scratch.solution_grads_u_total[q_point];
2449 * const Tensor<2,dim,NumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
2450 * const NumberType det_F = determinant(F);
2451 * const Tensor<2,dim,NumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
2452 * const SymmetricTensor<2,dim,NumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
2453 * const Tensor<2,dim,NumberType> F_inv = invert(F);
2454 * Assert(det_F > NumberType(0.0), ExcInternalError());
2456 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2458 * const unsigned int k_group = fe.system_to_base_index(k).first.first;
2460 * if (k_group == u_dof)
2462 * scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
2463 * scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
2466 * Assert(k_group <= u_dof, ExcInternalError());
2469 * const SymmetricTensor<2,dim,NumberType> tau = lqph[q_point]->get_tau(det_F,b_bar);
2470 * const SymmetricTensor<4,dim,NumberType> Jc = lqph[q_point]->get_Jc(det_F,b_bar);
2471 * const Tensor<2,dim,NumberType> tau_ns (tau);
2475 * Next we define some aliases to make the assembly process easier to
2479 * const std::vector<SymmetricTensor<2, dim> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
2480 * const std::vector<Tensor<2, dim> > &grad_Nx = scratch.grad_Nx[q_point];
2481 * const double JxW = scratch.fe_values_ref.JxW(q_point);
2483 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2485 * const unsigned int component_i = fe.system_to_component_index(i).first;
2486 * const unsigned int i_group = fe.system_to_base_index(i).first.first;
2488 * if (i_group == u_dof)
2489 * data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
2491 * Assert(i_group <= u_dof, ExcInternalError());
2493 * for (unsigned int j = 0; j <= i; ++j)
2495 * const unsigned int component_j = fe.system_to_component_index(j).first;
2496 * const unsigned int j_group = fe.system_to_base_index(j).first.first;
2500 * This is the @f$\mathsf{\mathbf{k}}_{\mathbf{u} \mathbf{u}}@f$
2501 * contribution. It comprises a material contribution, and a
2502 * geometrical stress contribution which is only added along
2503 * the local matrix diagonals:
2506 * if ((i_group == j_group) && (i_group == u_dof))
2508 * data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
2509 * * symm_grad_Nx[j] * JxW;
2510 * if (component_i == component_j) // geometrical stress contribution
2511 * data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau_ns
2512 * * grad_Nx[j][component_j] * JxW;
2515 * Assert((i_group <= u_dof) && (j_group <= u_dof),
2516 * ExcInternalError());
2524 * Finally, we need to copy the lower half of the local matrix into the
2528 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2529 * for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
2530 * data.cell_matrix(i, j) = data.cell_matrix(j, i);
2535 * #ifdef ENABLE_SACADO_FORMULATION
2538 * template <int dim>
2539 * struct Assembler<dim,Sacado::Fad::DFad<double> > : Assembler_Base<dim,Sacado::Fad::DFad<double> >
2541 * typedef Sacado::Fad::DFad<double> ADNumberType;
2542 * using typename Assembler_Base<dim,ADNumberType>::ScratchData_ASM;
2543 * using typename Assembler_Base<dim,ADNumberType>::PerTaskData_ASM;
2545 * virtual ~Assembler() {}
2548 * assemble_system_tangent_residual_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2549 * ScratchData_ASM &scratch,
2550 * PerTaskData_ASM &data) override
2554 * Aliases for data referenced from the Solid class
2557 * const unsigned int &n_q_points = data.solid->n_q_points;
2558 * const unsigned int &dofs_per_cell = data.solid->dofs_per_cell;
2559 * const FESystem<dim> &fe = data.solid->fe;
2560 * const unsigned int &u_dof = data.solid->u_dof;
2561 * const FEValuesExtractors::Vector &u_fe = data.solid->u_fe;
2565 * scratch.fe_values_ref.reinit(cell);
2566 * cell->get_dof_indices(data.local_dof_indices);
2568 * const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > > lqph =
2569 * const_cast<const Solid<dim,ADNumberType> *>(data.solid)->quadrature_point_history.get_data(cell);
2570 * Assert(lqph.size() == n_q_points, ExcInternalError());
2572 * const unsigned int n_independent_variables = data.local_dof_indices.size();
2573 * std::vector<double> local_dof_values(n_independent_variables);
2574 * cell->get_dof_values(scratch.solution_total,
2575 * local_dof_values.begin(),
2576 * local_dof_values.end());
2580 * We now retrieve a set of degree-of-freedom values that
2581 * have the operations that are performed with them tracked.
2584 * std::vector<ADNumberType> local_dof_values_ad (n_independent_variables);
2585 * for (unsigned int i=0; i<n_independent_variables; ++i)
2586 * local_dof_values_ad[i] = ADNumberType(n_independent_variables, i, local_dof_values[i]);
2590 * Compute all values, gradients etc. based on sensitive
2591 * AD degree-of-freedom values.
2594 * scratch.fe_values_ref[u_fe].get_function_gradients_from_local_dof_values(
2595 * local_dof_values_ad,
2596 * scratch.solution_grads_u_total);
2600 * Accumulate the residual value for each degree of freedom.
2601 * Note: Its important that the vectors is initialised (zero'd)
correctly.
2615 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
2617 *
const unsigned int k_group = fe.system_to_base_index(
k).first.first;
2636 *
const std::vector<SymmetricTensor<2, dim,ADNumberType> > &
symm_grad_Nx = scratch.symm_grad_Nx[
q_point];
2637 *
const double JxW = scratch.fe_values_ref.JxW(
q_point);
2639 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
2641 *
const unsigned int i_group = fe.system_to_base_index(i).first.first;
2650 *
for (
unsigned int I=0; I<n_independent_variables; ++I)
2654 *
for (
unsigned int J=0; J<n_independent_variables; ++J)
2669 *
template <
int dim>
2673 *
typedef Sacado::Rad::ADvar<ADDerivType>
ADNumberType;
2689 *
const unsigned int &n_q_points =
data.solid->n_q_points;
2694 *
scratch.fe_values_ref.reinit(cell);
2695 *
cell->get_dof_indices(
data.local_dof_indices);
2697 *
const std::vector<std::shared_ptr<const PointHistory<dim,ADNumberType> > >
lqph =
2698 *
data.solid->quadrature_point_history.get_data(cell);
2699 *
Assert(
lqph.size() == n_q_points, ExcInternalError());
2701 *
const unsigned int n_independent_variables =
data.local_dof_indices.size();
2702 *
std::vector<double> local_dof_values(n_independent_variables);
2703 *
cell->get_dof_values(scratch.solution_total,
2704 *
local_dof_values.begin(),
2705 *
local_dof_values.end());
2714 *
for (
unsigned int i=0; i<n_independent_variables; ++i)
2723 *
scratch.fe_values_ref[
u_fe].get_function_gradients_from_local_dof_values(
2725 *
scratch.solution_grads_u_total);
2735 * ADNumberType cell_energy_ad = ADNumberType(0.0);
2736 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2738 * const Tensor<2,dim,ADNumberType> &grad_u = scratch.solution_grads_u_total[q_point];
2739 * const Tensor<2,dim,ADNumberType> F = Physics::Elasticity::Kinematics::F(grad_u);
2740 * const ADNumberType det_F = determinant(F);
2741 * const Tensor<2,dim,ADNumberType> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
2742 * const SymmetricTensor<2,dim,ADNumberType> b_bar = Physics::Elasticity::Kinematics::b(F_bar);
2743 * Assert(det_F > ADNumberType(0.0), ExcInternalError());
2747 * Next we define some position-dependent aliases, again to
2748 * make the assembly process easier to follow.
2751 * const double JxW = scratch.fe_values_ref.JxW(q_point);
2753 * const ADNumberType Psi = lqph[q_point]->get_Psi(det_F,b_bar);
2757 * We extract the configuration-dependent material energy
2758 * from our QPH history objects for the current quadrature point
2759 * and integrate its contribution to increment the total
2763 * cell_energy_ad += Psi * JxW;
2768 * Compute derivatives of reverse-mode AD variables
2771 * ADNumberType::Gradcomp();
2773 * for (unsigned int I=0; I<n_independent_variables; ++I)
2777 * This computes the adjoint df/dX_{i} [reverse-mode]
2780 * const ADDerivType residual_I = local_dof_values_ad[I].adj();
2781 * data.cell_rhs(I) = -residual_I.val(); // RHS = - residual
2782 * for (unsigned int J=0; J<n_independent_variables; ++J)
2786 * Compute the gradients of the residual entry [forward-mode]
2789 * data.cell_matrix(I,J) = residual_I.dx(J); // linearisation_IJ
2802 * Since we use TBB for assembly, we simply setup a copy of the
2803 * data structures required for the process and pass them, along
2804 * with the memory addresses of the assembly functions to the
2805 * WorkStream object for processing. Note that we must ensure that
2806 * the matrix is reset before any assembly operations can occur.
2809 * template <int dim,typename NumberType>
2810 * void Solid<dim,NumberType>::assemble_system(const BlockVector<double> &solution_delta)
2812 * timer.enter_subsection("Assemble linear system");
2813 * std::cout << " ASM " << std::flush;
2815 * tangent_matrix = 0.0;
2818 * const UpdateFlags uf_cell(update_gradients |
2819 * update_JxW_values);
2820 * const UpdateFlags uf_face(update_values |
2821 * update_JxW_values);
2823 * const BlockVector<double> solution_total(get_total_solution(solution_delta));
2824 * typename Assembler_Base<dim,NumberType>::PerTaskData_ASM per_task_data(this);
2825 * typename Assembler_Base<dim,NumberType>::ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total);
2826 * Assembler<dim,NumberType> assembler;
2828 * WorkStream::run(dof_handler_ref.begin_active(),
2829 * dof_handler_ref.end(),
2830 * static_cast<Assembler_Base<dim,NumberType>&>(assembler),
2831 * &Assembler_Base<dim,NumberType>::assemble_system_one_cell,
2832 * &Assembler_Base<dim,NumberType>::copy_local_to_global_ASM,
2836 * timer.leave_subsection();
2843 * <a name="cook_membrane.cc-Solidmake_constraints"></a>
2844 * <h4>Solid::make_constraints</h4>
2845 * The constraints for this problem are simple to describe.
2846 * However, since we are dealing with an iterative Newton method,
2847 * it should be noted that any displacement constraints should only
2848 * be specified at the zeroth iteration and subsequently no
2849 * additional contributions are to be made since the constraints
2850 * are already exactly satisfied.
2853 * template <int dim,typename NumberType>
2854 * void Solid<dim,NumberType>::make_constraints(const int &it_nr)
2856 * std::cout << " CST " << std::flush;
2860 * Since the constraints are different at different Newton iterations, we
2861 * need to clear the constraints matrix and completely rebuild
2862 * it. However, after the first iteration, the constraints remain the same
2863 * and we can simply skip the rebuilding step if we do not clear it.
2868 * const bool apply_dirichlet_bc = (it_nr == 0);
2872 * The boundary conditions for the indentation problem are as follows: On
2873 * the -x face (ID = 1), we set up a zero-displacement condition, -y and +y traction
2897 *
constraints.
clear();
2910 *
fe.component_mask(
u_fe));
2932 *
if (constraints.has_inhomogeneities())
2938 *
constraints.
clear();
2943 *
constraints.close();
2949 * <a name=
"cook_membrane.cc-Solidsolve_linear_system"></a>
2955 *
template <
int dim,
typename NumberType>
2956 *
std::pair<unsigned int, double>
2962 *
unsigned int lin_it = 0;
2971 *
timer.enter_subsection(
"Linear solver");
2972 *
std::cout <<
" SLV " << std::flush;
2973 *
if (parameters.type_lin ==
"CG")
2975 *
const int solver_its =
static_cast<unsigned int>(
2977 *
* parameters.max_iterations_lin);
2978 *
const double tol_sol = parameters.tol_lin
2988 *
We've chosen by default a SSOR preconditioner as it appears to
2989 * provide the fastest solver convergence characteristics for this
2990 * problem on a single-thread machine. However, for multicore
2991 * computing, the Jacobi preconditioner which is multithreaded may
2992 * converge quicker for larger linear systems.
2995 * PreconditionSelector<SparseMatrix<double>, Vector<double> >
2996 * preconditioner (parameters.preconditioner_type,
2997 * parameters.preconditioner_relaxation);
2998 * preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
3000 * solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
3001 * newton_update.block(u_dof),
3002 * system_rhs.block(u_dof),
3005 * lin_it = solver_control.last_step();
3006 * lin_res = solver_control.last_value();
3008 * else if (parameters.type_lin == "Direct")
3012 * Otherwise if the problem is small
3013 * enough, a direct solver can be
3017 * SparseDirectUMFPACK A_direct;
3018 * A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
3019 * A_direct.vmult(newton_update.block(u_dof), system_rhs.block(u_dof));
3025 * Assert (false, ExcMessage("Linear solver type not implemented"));
3027 * timer.leave_subsection();
3032 * Now that we have the displacement update, distribute the constraints
3033 * back to the Newton update:
3036 * constraints.distribute(newton_update);
3038 * return std::make_pair(lin_it, lin_res);
3044 * <a name="cook_membrane.cc-Solidoutput_results"></a>
3045 * <h4>Solid::output_results</h4>
3046 * Here we present how the results are written to file to be viewed
3047 * using ParaView or Visit. The method is similar to that shown in the
3048 * tutorials so will not be discussed in detail.
3051 * template <int dim,typename NumberType>
3052 * void Solid<dim,NumberType>::output_results() const
3054 * DataOut<dim> data_out;
3055 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3056 * data_component_interpretation(dim,
3057 * DataComponentInterpretation::component_is_part_of_vector);
3059 * std::vector<std::string> solution_name(dim, "displacement");
3061 * data_out.attach_dof_handler(dof_handler_ref);
3062 * data_out.add_data_vector(solution_n,
3064 * DataOut<dim>::type_dof_data,
3065 * data_component_interpretation);
3069 * Since we are dealing with a large deformation problem, it would be nice
3070 * to display the result on a displaced grid! The MappingQEulerian class
3071 * linked with the DataOut class provides an interface through which this
3072 * can be achieved without physically moving the grid points in the
3073 * Triangulation object ourselves. We first need to copy the solution to
3074 * a temporary vector and then create the Eulerian mapping. We also
3075 * specify the polynomial degree to the DataOut object in order to produce
3076 * a more refined output data set when higher order polynomials are used.
3079 * Vector<double> soln(solution_n.size());
3080 * for (unsigned int i = 0; i < soln.size(); ++i)
3081 * soln(i) = solution_n(i);
3082 * MappingQEulerian<dim> q_mapping(degree, dof_handler_ref, soln);
3083 * data_out.build_patches(q_mapping, degree);
3085 * std::ostringstream filename;
3086 * filename << "solution-" << time.get_timestep() << ".vtk";
3088 * std::ofstream output(filename.str().c_str());
3089 * data_out.write_vtk(output);
3098 * <a name="cook_membrane.cc-Mainfunction"></a>
3099 * <h3>Main function</h3>
3100 * Lastly we provide the main driver function which appears
3101 * no different to the other tutorials.
3104 * int main (int argc, char *argv[])
3106 * using namespace dealii;
3107 * using namespace Cook_Membrane;
3109 * const unsigned int dim = 3;
3112 * deallog.depth_console(0);
3113 * Parameters::AllParameters parameters("parameters.prm");
3114 * if (parameters.automatic_differentiation_order == 0)
3116 * std::cout << "Assembly method: Residual and linearisation are computed manually." << std::endl;
3120 * Allow multi-threading
3123 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3124 * ::numbers::invalid_unsigned_int);
3126 * typedef double NumberType;
3127 * Solid<dim,NumberType> solid_3d(parameters);
3130 * #ifdef ENABLE_SACADO_FORMULATION
3131 * else if (parameters.automatic_differentiation_order == 1)
3133 * std::cout << "Assembly method: Residual computed manually; linearisation performed using AD." << std::endl;
3137 * Allow multi-threading
3140 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3141 * ::numbers::invalid_unsigned_int);
3143 * typedef Sacado::Fad::DFad<double> NumberType;
3144 * Solid<dim,NumberType> solid_3d(parameters);
3147 * else if (parameters.automatic_differentiation_order == 2)
3149 * std::cout << "Assembly method: Residual and linearisation computed using AD." << std::endl;
3153 * Sacado Rad-Fad is not thread-safe, so disable threading.
3154 * Parallelisation using MPI would be possible though.
3157 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
3160 * typedef Sacado::Rad::ADvar< Sacado::Fad::DFad<double> > NumberType;
3161 * Solid<dim,NumberType> solid_3d(parameters);
3167 * AssertThrow(false,
3168 * ExcMessage("The selected assembly method is not supported. "
3169 * "You need deal.II 9.0 and Trilinos with the Sacado package "
3170 * "to enable assembly using automatic differentiation."));
3173 * catch (std::exception &exc)
3175 * std::cerr << std::endl << std::endl
3176 * << "----------------------------------------------------"
3178 * std::cerr << "Exception on processing: " << std::endl << exc.what()
3179 * << std::endl << "Aborting!" << std::endl
3180 * << "----------------------------------------------------"
3187 * std::cerr << std::endl << std::endl
3188 * << "----------------------------------------------------"
3190 * std::cerr << "Unknown exception!" << std::endl << "Aborting!"
3192 * << "----------------------------------------------------"
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define DEAL_II_VERSION_MAJOR
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
std::vector< index_type > data
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F_iso(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)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
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
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)