Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity.h
Go to the documentation of this file.
1 
228  *
229  * /*
230  * * Authors: Jean-Paul Pelteret, University of Erlangen-Nuremberg, 2017
231  * */
232  *
233  * #include <deal.II/base/mpi.h>
234  * #include <deal.II/base/function.h>
235  * #include <deal.II/base/parameter_handler.h>
236  * #include <deal.II/base/point.h>
237  * #include <deal.II/base/quadrature_lib.h>
238  * #include <deal.II/base/symmetric_tensor.h>
239  * #include <deal.II/base/tensor.h>
240  * #include <deal.II/base/timer.h>
241  * #include <deal.II/base/work_stream.h>
242  * #include <deal.II/dofs/dof_renumbering.h>
243  * #include <deal.II/dofs/dof_tools.h>
244  * #include <deal.II/base/quadrature_point_data.h>
245  * #include <deal.II/grid/filtered_iterator.h>
246  * #include <deal.II/grid/grid_generator.h>
247  * #include <deal.II/grid/grid_tools.h>
248  * #include <deal.II/grid/grid_in.h>
249  * #include <deal.II/grid/manifold_lib.h>
250  * #include <deal.II/grid/tria.h>
251  * #include <deal.II/grid/tria_boundary_lib.h>
252  * #include <deal.II/fe/fe_dgp_monomial.h>
253  * #include <deal.II/fe/fe_q.h>
254  * #include <deal.II/fe/fe_system.h>
255  * #include <deal.II/fe/fe_tools.h>
256  * #include <deal.II/fe/fe_values.h>
257  * #include <deal.II/fe/mapping_q_eulerian.h>
258  * #include <deal.II/lac/block_sparsity_pattern.h>
259  * #include <deal.II/lac/dynamic_sparsity_pattern.h>
260  * #include <deal.II/lac/full_matrix.h>
261  * #include <deal.II/lac/constraint_matrix.h>
262  * #include <deal.II/lac/solver_selector.h>
263  * #include <deal.II/lac/trilinos_block_sparse_matrix.h>
264  * #include <deal.II/lac/trilinos_precondition.h>
265  * #include <deal.II/lac/trilinos_sparsity_pattern.h>
266  * #include <deal.II/lac/trilinos_sparse_matrix.h>
267  * #include <deal.II/lac/trilinos_vector.h>
268  * #include <deal.II/lac/linear_operator.h>
269  * #include <deal.II/lac/packaged_operation.h>
270  * #include <deal.II/lac/trilinos_linear_operator.h>
271  * #include <deal.II/numerics/data_out.h>
272  * #include <deal.II/numerics/vector_tools.h>
273  * #include <deal.II/physics/transformations.h>
274  * #include <deal.II/physics/elasticity/kinematics.h>
275  * #include <deal.II/physics/elasticity/standard_tensors.h>
276  * #include <iostream>
277  * #include <fstream>
278  * #include <numeric>
279  *
280  * #include <deal.II/grid/grid_out.h>
281  *
282  * namespace ViscoElasStripHole
283  * {
284  * using namespace dealii;
285  * namespace LA = TrilinosWrappers;
286  * namespace Parameters
287  * {
288  * struct BoundaryConditions
289  * {
290  * BoundaryConditions();
291  *
292  * std::string driver;
293  * double stretch;
294  * double pressure;
295  * double load_time;
296  *
297  * const types::boundary_id boundary_id_minus_X;
298  * const types::boundary_id boundary_id_plus_X;
299  * const types::boundary_id boundary_id_minus_Y;
300  * const types::boundary_id boundary_id_plus_Y;
301  * const types::boundary_id boundary_id_minus_Z;
302  * const types::boundary_id boundary_id_plus_Z;
303  * const types::boundary_id boundary_id_hole;
304  * const types::manifold_id manifold_id_hole;
305  *
306  * static void
307  * declare_parameters(ParameterHandler &prm);
308  * void
309  * parse_parameters(ParameterHandler &prm);
310  * };
311  * BoundaryConditions::BoundaryConditions()
312  * :
313  * driver ("Neumann"),
314  * stretch (2.0),
315  * pressure(0.0),
316  * load_time(2.5),
317  * boundary_id_minus_X (1),
318  * boundary_id_plus_X (2),
319  * boundary_id_minus_Y (3),
320  * boundary_id_plus_Y (4),
321  * boundary_id_minus_Z (5),
322  * boundary_id_plus_Z (6),
323  * boundary_id_hole (10),
324  * manifold_id_hole (10)
325  * { }
326  * void BoundaryConditions::declare_parameters(ParameterHandler &prm)
327  * {
328  * prm.enter_subsection("Boundary conditions");
329  * {
330  * prm.declare_entry("Driver", "Dirichlet",
331  * Patterns::Selection("Dirichlet|Neumann"),
332  * "Driver boundary condition for the problem");
333  * prm.declare_entry("Final stretch", "2.0",
334  * Patterns::Double(1.0),
335  * "Positive stretch applied length-ways to the strip");
336  * prm.declare_entry("Applied pressure", "0.0",
337  * Patterns::Double(-1e3,1e3),
338  * "Hydrostatic pressure applied (in the referential configuration) to the interior surface of the hole");
339  * prm.declare_entry("Load time", "2.5",
340  * Patterns::Double(0.0),
341  * "Total time over which the stretch/pressure is ramped up");
342  * }
343  * prm.leave_subsection();
344  * }
345  * void BoundaryConditions::parse_parameters(ParameterHandler &prm)
346  * {
347  * prm.enter_subsection("Boundary conditions");
348  * {
349  * driver = prm.get("Driver");
350  * stretch = prm.get_double("Final stretch");
351  * pressure = prm.get_double("Applied pressure");
352  * load_time = prm.get_double("Load time");
353  * }
354  * prm.leave_subsection();
355  * }
356  * struct FESystem
357  * {
358  * unsigned int poly_degree;
359  * unsigned int quad_order;
360  * static void
361  * declare_parameters(ParameterHandler &prm);
362  * void
363  * parse_parameters(ParameterHandler &prm);
364  * };
365  * void FESystem::declare_parameters(ParameterHandler &prm)
366  * {
367  * prm.enter_subsection("Finite element system");
368  * {
369  * prm.declare_entry("Polynomial degree", "2",
370  * Patterns::Integer(0),
371  * "Displacement system polynomial order");
372  * prm.declare_entry("Quadrature order", "3",
373  * Patterns::Integer(0),
374  * "Gauss quadrature order");
375  * }
376  * prm.leave_subsection();
377  * }
378  * void FESystem::parse_parameters(ParameterHandler &prm)
379  * {
380  * prm.enter_subsection("Finite element system");
381  * {
382  * poly_degree = prm.get_integer("Polynomial degree");
383  * quad_order = prm.get_integer("Quadrature order");
384  * }
385  * prm.leave_subsection();
386  * }
387  * struct Geometry
388  * {
389  * double length;
390  * double width;
391  * double thickness;
392  * double hole_diameter;
393  * double hole_division_fraction;
394  * unsigned int n_repetitions_xy;
395  * unsigned int n_repetitions_z;
396  * unsigned int global_refinement;
397  * double scale;
398  * static void
399  * declare_parameters(ParameterHandler &prm);
400  * void
401  * parse_parameters(ParameterHandler &prm);
402  * };
403  * void Geometry::declare_parameters(ParameterHandler &prm)
404  * {
405  * prm.enter_subsection("Geometry");
406  * {
407  * prm.declare_entry("Length", "100.0",
408  * Patterns::Double(0.0),
409  * "Total sample length");
410  * prm.declare_entry("Width", "50.0",
411  * Patterns::Double(0.0),
412  * "Total sample width");
413  * prm.declare_entry("Thickness", "5.0",
414  * Patterns::Double(0.0),
415  * "Total sample thickness");
416  * prm.declare_entry("Hole diameter", "20.0",
417  * Patterns::Double(0.0),
418  * "Hole diameter");
419  * prm.declare_entry("Hole division fraction", "0.5",
420  * Patterns::Double(0.0,1.0),
421  * "A geometric factor affecting the discretisation near the hole");
422  * prm.declare_entry("Number of subdivisions in cross-section", "2",
423  * Patterns::Integer(1.0),
424  * "A factor defining the number of initial grid subdivisions in the cross-section");
425  * prm.declare_entry("Number of subdivisions thickness", "6",
426  * Patterns::Integer(1.0),
427  * "A factor defining the number of initial grid subdivisions through the thickness");
428  * prm.declare_entry("Global refinement", "2",
429  * Patterns::Integer(0),
430  * "Global refinement level");
431  * prm.declare_entry("Grid scale", "1e-3",
432  * Patterns::Double(0.0),
433  * "Global grid scaling factor");
434  * }
435  * prm.leave_subsection();
436  * }
437  * void Geometry::parse_parameters(ParameterHandler &prm)
438  * {
439  * prm.enter_subsection("Geometry");
440  * {
441  * length = prm.get_double("Length");
442  * width = prm.get_double("Width");
443  * thickness = prm.get_double("Thickness");
444  * hole_diameter = prm.get_double("Hole diameter");
445  * hole_division_fraction = prm.get_double("Hole division fraction");
446  * n_repetitions_xy = prm.get_integer("Number of subdivisions in cross-section");
447  * n_repetitions_z = prm.get_integer("Number of subdivisions thickness");
448  * global_refinement = prm.get_integer("Global refinement");
449  * scale = prm.get_double("Grid scale");
450  * }
451  * prm.leave_subsection();
452  * }
453  * struct Materials
454  * {
455  * double nu_e;
456  * double mu_e;
457  * double mu_v;
458  * double tau_v;
459  * static void
460  * declare_parameters(ParameterHandler &prm);
461  * void
462  * parse_parameters(ParameterHandler &prm);
463  * };
464  * void Materials::declare_parameters(ParameterHandler &prm)
465  * {
466  * prm.enter_subsection("Material properties");
467  * {
468  * prm.declare_entry("Poisson's ratio", "0.4999",
469  * Patterns::Double(-1.0,0.5),
470  * "Poisson's ratio");
471  * prm.declare_entry("Elastic shear modulus", "80.194e6",
472  * Patterns::Double(0.0),
473  * "Elastic shear modulus");
474  * prm.declare_entry("Viscous shear modulus", "80.194e6",
475  * Patterns::Double(0.0),
476  * "Viscous shear modulus");
477  * prm.declare_entry("Viscous relaxation time", "2.0",
478  * Patterns::Double(0.0),
479  * "Viscous relaxation time");
480  * }
481  * prm.leave_subsection();
482  * }
483  * void Materials::parse_parameters(ParameterHandler &prm)
484  * {
485  * prm.enter_subsection("Material properties");
486  * {
487  * nu_e = prm.get_double("Poisson's ratio");
488  * mu_e = prm.get_double("Elastic shear modulus");
489  * mu_v = prm.get_double("Viscous shear modulus");
490  * tau_v = prm.get_double("Viscous relaxation time");
491  * }
492  * prm.leave_subsection();
493  * }
494  * struct LinearSolver
495  * {
496  * std::string type_lin;
497  * double tol_lin;
498  * double max_iterations_lin;
499  * static void
500  * declare_parameters(ParameterHandler &prm);
501  * void
502  * parse_parameters(ParameterHandler &prm);
503  * };
504  * void LinearSolver::declare_parameters(ParameterHandler &prm)
505  * {
506  * prm.enter_subsection("Linear solver");
507  * {
508  * prm.declare_entry("Solver type", "cg",
510  * "Type of solver used to solve the linear system");
511  * prm.declare_entry("Residual", "1e-6",
512  * Patterns::Double(0.0),
513  * "Linear solver residual (scaled by residual norm)");
514  * prm.declare_entry("Max iteration multiplier", "1",
515  * Patterns::Double(1.0),
516  * "Linear solver iterations (multiples of the system matrix size)");
517  * }
518  * prm.leave_subsection();
519  * }
520  * void LinearSolver::parse_parameters(ParameterHandler &prm)
521  * {
522  * prm.enter_subsection("Linear solver");
523  * {
524  * type_lin = prm.get("Solver type");
525  * tol_lin = prm.get_double("Residual");
526  * max_iterations_lin = prm.get_double("Max iteration multiplier");
527  * }
528  * prm.leave_subsection();
529  * }
530  * struct NonlinearSolver
531  * {
532  * unsigned int max_iterations_NR;
533  * double tol_f;
534  * double tol_u;
535  * static void
536  * declare_parameters(ParameterHandler &prm);
537  * void
538  * parse_parameters(ParameterHandler &prm);
539  * };
540  * void NonlinearSolver::declare_parameters(ParameterHandler &prm)
541  * {
542  * prm.enter_subsection("Nonlinear solver");
543  * {
544  * prm.declare_entry("Max iterations Newton-Raphson", "10",
545  * Patterns::Integer(0),
546  * "Number of Newton-Raphson iterations allowed");
547  * prm.declare_entry("Tolerance displacement", "1.0e-6",
548  * Patterns::Double(0.0),
549  * "Displacement error tolerance");
550  * prm.declare_entry("Tolerance force", "1.0e-9",
551  * Patterns::Double(0.0),
552  * "Force residual tolerance");
553  * }
554  * prm.leave_subsection();
555  * }
556  * void NonlinearSolver::parse_parameters(ParameterHandler &prm)
557  * {
558  * prm.enter_subsection("Nonlinear solver");
559  * {
560  * max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
561  * tol_f = prm.get_double("Tolerance force");
562  * tol_u = prm.get_double("Tolerance displacement");
563  * }
564  * prm.leave_subsection();
565  * }
566  * struct Time
567  * {
568  * double delta_t;
569  * double end_time;
570  * static void
571  * declare_parameters(ParameterHandler &prm);
572  * void
573  * parse_parameters(ParameterHandler &prm);
574  * };
575  * void Time::declare_parameters(ParameterHandler &prm)
576  * {
577  * prm.enter_subsection("Time");
578  * {
579  * prm.declare_entry("End time", "1",
580  * Patterns::Double(),
581  * "End time");
582  * prm.declare_entry("Time step size", "0.1",
583  * Patterns::Double(),
584  * "Time step size");
585  * }
586  * prm.leave_subsection();
587  * }
588  * void Time::parse_parameters(ParameterHandler &prm)
589  * {
590  * prm.enter_subsection("Time");
591  * {
592  * end_time = prm.get_double("End time");
593  * delta_t = prm.get_double("Time step size");
594  * }
595  * prm.leave_subsection();
596  * }
597  * struct AllParameters
598  * : public BoundaryConditions,
599  * public FESystem,
600  * public Geometry,
601  * public Materials,
602  * public LinearSolver,
603  * public NonlinearSolver,
604  * public Time
605  * {
606  * AllParameters(const std::string &input_file);
607  * static void
608  * declare_parameters(ParameterHandler &prm);
609  * void
610  * parse_parameters(ParameterHandler &prm);
611  * };
612  * AllParameters::AllParameters(const std::string &input_file)
613  * {
614  * ParameterHandler prm;
615  * declare_parameters(prm);
616  * prm.parse_input(input_file);
617  * parse_parameters(prm);
618  * }
619  * void AllParameters::declare_parameters(ParameterHandler &prm)
620  * {
621  * BoundaryConditions::declare_parameters(prm);
622  * FESystem::declare_parameters(prm);
623  * Geometry::declare_parameters(prm);
624  * Materials::declare_parameters(prm);
625  * LinearSolver::declare_parameters(prm);
626  * NonlinearSolver::declare_parameters(prm);
627  * Time::declare_parameters(prm);
628  * }
629  * void AllParameters::parse_parameters(ParameterHandler &prm)
630  * {
631  * BoundaryConditions::parse_parameters(prm);
632  * FESystem::parse_parameters(prm);
633  * Geometry::parse_parameters(prm);
634  * Materials::parse_parameters(prm);
635  * LinearSolver::parse_parameters(prm);
636  * NonlinearSolver::parse_parameters(prm);
637  * Time::parse_parameters(prm);
638  * }
639  * }
640  * class Time
641  * {
642  * public:
643  * Time (const double time_end,
644  * const double delta_t)
645  * :
646  * timestep(0),
647  * time_current(0.0),
648  * time_end(time_end),
649  * delta_t(delta_t)
650  * {}
651  * virtual ~Time()
652  * {}
653  * double current() const
654  * {
655  * return time_current;
656  * }
657  * double end() const
658  * {
659  * return time_end;
660  * }
661  * double get_delta_t() const
662  * {
663  * return delta_t;
664  * }
665  * unsigned int get_timestep() const
666  * {
667  * return timestep;
668  * }
669  * void increment()
670  * {
671  * time_current += delta_t;
672  * ++timestep;
673  * }
674  * private:
675  * unsigned int timestep;
676  * double time_current;
677  * const double time_end;
678  * const double delta_t;
679  * };
680  * template <int dim>
681  * class Material_Compressible_Three_Field_Linear_Viscoelastic
682  * {
683  * public:
684  * Material_Compressible_Three_Field_Linear_Viscoelastic(const double mu_e,
685  * const double nu_e,
686  * const double mu_v,
687  * const double tau_v,
688  * const Time &time)
689  * :
690  * kappa((2.0 * mu_e * (1.0 + nu_e)) / (3.0 * (1.0 - 2.0 * nu_e))),
691  * mu_e(mu_e),
692  * mu_v(mu_v),
693  * tau_v(tau_v),
694  * time(time),
697  * {
698  * Assert(kappa > 0, ExcInternalError());
699  * }
700  * ~Material_Compressible_Three_Field_Linear_Viscoelastic()
701  * {}
702  *
704  * get_tau(const Tensor<2,dim> &F,
705  * const double &p_tilde) const
706  * {
707  * return get_tau_iso(F) + get_tau_vol(F,p_tilde);
708  * }
709  * SymmetricTensor<4,dim> get_Jc(const Tensor<2,dim> &F,
710  * const double &p_tilde) const
711  * {
712  * return get_Jc_iso(F) + get_Jc_vol(F,p_tilde);
713  * }
714  * double
715  * get_dPsi_vol_dJ(const double &J_tilde) const
716  * {
717  * return (kappa / 2.0) * (J_tilde - 1.0 / J_tilde);
718  * }
719  * double
720  * get_d2Psi_vol_dJ2(const double &J_tilde) const
721  * {
722  * return ( (kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde)));
723  * }
724  * void
725  * update_internal_equilibrium(const Tensor<2, dim> &F,
726  * const double &/*p_tilde*/,
727  * const double &/*J_tilde*/)
728  * {
729  * const double det_F = determinant(F);
730  * const SymmetricTensor<2,dim> C_bar = std::pow(det_F, -2.0 / dim) * Physics::Elasticity::Kinematics::C(F);
731  * @endcode
732  *
733  * Linder2011 eq 54
734  * Assumes first-oder backward Euler time discretisation
735  *
736  * @code
737  * Q_n_t = (1.0/(1.0 + time.get_delta_t()/tau_v))*(Q_t1 + (time.get_delta_t()/tau_v)*invert(C_bar));
738  * }
739  * void
740  * update_end_timestep()
741  * {
742  * Q_t1 = Q_n_t;
743  * }
744  *
745  * protected:
746  * const double kappa;
747  * const double mu_e;
748  * const double mu_v;
749  * const double tau_v;
750  * const Time &time;
751  * SymmetricTensor<2,dim> Q_n_t; // Value of internal variable at this Newton step and timestep
752  * SymmetricTensor<2,dim> Q_t1; // Value of internal variable at the previous timestep
753  *
755  * get_tau_vol(const Tensor<2,dim> &F,
756  * const double &p_tilde) const
757  * {
758  * const double det_F = determinant(F);
759  *
760  * return p_tilde * det_F * Physics::Elasticity::StandardTensors<dim>::I;
761  * }
763  * get_tau_iso(const Tensor<2,dim> &F) const
764  * {
765  * return Physics::Elasticity::StandardTensors<dim>::dev_P * get_tau_bar(F);
766  * }
768  * get_tau_bar(const Tensor<2,dim> &F) const
769  * {
770  * const double det_F = determinant(F);
771  * const Tensor<2,dim> F_bar = std::pow(det_F, -1.0 / dim) * F;
772  * const SymmetricTensor<2,dim> b_bar = std::pow(det_F, -2.0 / dim) * symmetrize(F * transpose(F));
773  * @endcode
774  *
775  * Elastic Neo-Hookean + Linder2011 eq 47
776  *
777  * @code
778  * return mu_e * b_bar
779  * + mu_v * symmetrize(F_bar*static_cast<Tensor<2,dim> >(Q_n_t)*transpose(F_bar));
780  * }
781  * SymmetricTensor<4, dim> get_Jc_vol(const Tensor<2,dim> &F,
782  * const double &p_tilde) const
783  * {
784  * const double det_F = determinant(F);
785  * return p_tilde * det_F
788  * }
789  * SymmetricTensor<4, dim> get_Jc_iso(const Tensor<2,dim> &F) const
790  * {
791  * const SymmetricTensor<2, dim> tau_bar = get_tau_bar(F);
792  * const SymmetricTensor<2, dim> tau_iso = get_tau_iso(F);
793  * const SymmetricTensor<4, dim> tau_iso_x_I
794  * = outer_product(tau_iso,
796  * const SymmetricTensor<4, dim> I_x_tau_iso
798  * tau_iso);
799  * const SymmetricTensor<4, dim> c_bar = get_c_bar(F);
800  * return (2.0 / dim) * trace(tau_bar)
802  * - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso)
805  * }
806  * SymmetricTensor<4, dim> get_c_bar(const Tensor<2,dim> &/*F*/) const
807  * {
808  * @endcode
809  *
810  * Elastic Neo-Hookean + Linder2011 eq 56
811  *
812  * @code
813  * return -2.0*mu_v*((time.get_delta_t()/tau_v)/(1.0 + time.get_delta_t()/tau_v))*Physics::Elasticity::StandardTensors<dim>::S;
814  * }
815  * };
816  * template <int dim>
817  * class PointHistory
818  * {
819  * public:
820  * PointHistory()
821  * {}
822  * virtual ~PointHistory()
823  * {}
824  * void
825  * setup_lqp (const Parameters::AllParameters &parameters,
826  * const Time &time)
827  * {
828  * material.reset(new Material_Compressible_Three_Field_Linear_Viscoelastic<dim>(
829  * parameters.mu_e, parameters.nu_e,
830  * parameters.mu_v, parameters.tau_v,
831  * time));
832  * }
833  *
835  * get_tau(const Tensor<2, dim> &F,
836  * const double &p_tilde) const
837  * {
838  * return material->get_tau(F, p_tilde);
839  * }
841  * get_Jc(const Tensor<2, dim> &F,
842  * const double &p_tilde) const
843  * {
844  * return material->get_Jc(F, p_tilde);
845  * }
846  * double
847  * get_dPsi_vol_dJ(const double &J_tilde) const
848  * {
849  * return material->get_dPsi_vol_dJ(J_tilde);
850  * }
851  * double
852  * get_d2Psi_vol_dJ2(const double &J_tilde) const
853  * {
854  * return material->get_d2Psi_vol_dJ2(J_tilde);
855  * }
856  * void
857  * update_internal_equilibrium(const Tensor<2, dim> &F,
858  * const double &p_tilde,
859  * const double &J_tilde)
860  * {
861  * material->update_internal_equilibrium(F,p_tilde,J_tilde);
862  * }
863  * void
864  * update_end_timestep()
865  * {
866  * material->update_end_timestep();
867  * }
868  * private:
869  * std::shared_ptr< Material_Compressible_Three_Field_Linear_Viscoelastic<dim> > material;
870  * };
871  * template <int dim>
872  * class Solid
873  * {
874  * public:
875  * Solid(const std::string &input_file);
876  * virtual
877  * ~Solid();
878  * void
879  * run();
880  * private:
881  * struct PerTaskData_ASM;
882  * struct ScratchData_ASM;
883  * void
884  * make_grid();
885  * void
886  * make_2d_quarter_plate_with_hole(Triangulation<2> &tria_2d,
887  * const double half_length,
888  * const double half_width,
889  * const double hole_radius,
890  * const unsigned int n_repetitions_xy = 1,
891  * const double hole_division_fraction = 0.25);
892  * void
893  * setup_system(LA::MPI::BlockVector &solution_delta);
894  * void
895  * determine_component_extractors();
896  * void
897  * assemble_system(const LA::MPI::BlockVector &solution_delta);
898  * void
899  * assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
900  * ScratchData_ASM &scratch,
901  * PerTaskData_ASM &data) const;
902  * void
903  * copy_local_to_global_system(const PerTaskData_ASM &data);
904  * void
905  * make_constraints(const int &it_nr);
906  * void
907  * setup_qph();
908  * void
909  * solve_nonlinear_timestep(LA::MPI::BlockVector &solution_delta);
910  * std::pair<unsigned int, double>
911  * solve_linear_system(LA::MPI::BlockVector &newton_update);
913  * get_solution_total(const LA::MPI::BlockVector &solution_delta) const;
914  * void
915  * update_end_timestep();
916  * void
917  * output_results(const unsigned int timestep,
918  * const double current_time) const;
919  * void
920  * compute_vertex_positions(std::vector<double> &real_time,
921  * std::vector<std::vector<Point<dim> > > &tracked_vertices,
922  * const LA::MPI::BlockVector &solution_total) const;
923  *
924  * @endcode
925  *
926  * Parallel communication
927  *
928  * @code
929  * MPI_Comm mpi_communicator;
930  * const unsigned int n_mpi_processes;
931  * const unsigned int this_mpi_process;
932  * mutable ConditionalOStream pcout;
933  *
934  * Parameters::AllParameters parameters;
936  * Time time;
937  * mutable TimerOutput timer;
939  * PointHistory<dim> > quadrature_point_history;
940  * const unsigned int degree;
941  * const FESystem<dim> fe;
942  * DoFHandler<dim> dof_handler;
943  * const unsigned int dofs_per_cell;
944  * const FEValuesExtractors::Vector u_fe;
945  * const FEValuesExtractors::Scalar p_fe;
946  * const FEValuesExtractors::Scalar J_fe;
947  * static const unsigned int n_blocks = 3;
948  * static const unsigned int n_components = dim + 2;
949  * static const unsigned int first_u_component = 0;
950  * static const unsigned int p_component = dim;
951  * static const unsigned int J_component = dim + 1;
952  * enum
953  * {
954  * u_block = 0,
955  * p_block = 1,
956  * J_block = 2
957  * };
958  * @endcode
959  *
960  * Block data
961  *
962  * @code
963  * std::vector<unsigned int> block_component;
964  *
965  * @endcode
966  *
967  * DoF index data
968  *
969  * @code
970  * std::vector<IndexSet> all_locally_owned_dofs;
971  * IndexSet locally_owned_dofs;
972  * IndexSet locally_relevant_dofs;
973  * std::vector<IndexSet> locally_owned_partitioning;
974  * std::vector<IndexSet> locally_relevant_partitioning;
975  * std::vector<types::global_dof_index> dofs_per_block;
976  * std::vector<types::global_dof_index> element_indices_u;
977  * std::vector<types::global_dof_index> element_indices_p;
978  * std::vector<types::global_dof_index> element_indices_J;
979  * const QGauss<dim> qf_cell;
980  * const QGauss<dim - 1> qf_face;
981  * const unsigned int n_q_points;
982  * const unsigned int n_q_points_f;
983  * ConstraintMatrix constraints;
984  * LA::BlockSparseMatrix tangent_matrix;
985  * LA::MPI::BlockVector system_rhs;
986  * LA::MPI::BlockVector solution_n;
987  * struct Errors
988  * {
989  * Errors()
990  * :
991  * norm(1.0), u(1.0), p(1.0), J(1.0)
992  * {}
993  * void reset()
994  * {
995  * norm = 1.0;
996  * u = 1.0;
997  * p = 1.0;
998  * J = 1.0;
999  * }
1000  * void normalise(const Errors &rhs)
1001  * {
1002  * if (rhs.norm != 0.0)
1003  * norm /= rhs.norm;
1004  * if (rhs.u != 0.0)
1005  * u /= rhs.u;
1006  * if (rhs.p != 0.0)
1007  * p /= rhs.p;
1008  * if (rhs.J != 0.0)
1009  * J /= rhs.J;
1010  * }
1011  * double norm, u, p, J;
1012  * };
1013  * Errors error_residual, error_residual_0, error_residual_norm, error_update,
1014  * error_update_0, error_update_norm;
1015  * void
1016  * get_error_residual(Errors &error_residual);
1017  * void
1018  * get_error_update(const LA::MPI::BlockVector &newton_update,
1019  * Errors &error_update);
1020  * std::pair<double, std::pair<double,double> >
1021  * get_error_dilation(const LA::MPI::BlockVector &solution_total) const;
1022  * void
1023  * print_conv_header();
1024  * void
1025  * print_conv_footer(const LA::MPI::BlockVector &solution_delta);
1026  * };
1027  * template <int dim>
1028  * Solid<dim>::Solid(const std::string &input_file)
1029  * :
1030  * mpi_communicator(MPI_COMM_WORLD),
1031  * n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
1032  * this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
1033  * pcout(std::cout, this_mpi_process == 0),
1034  * parameters(input_file),
1036  * time(parameters.end_time, parameters.delta_t),
1037  * timer(mpi_communicator,
1038  * pcout,
1041  * degree(parameters.poly_degree),
1042  * fe(FE_Q<dim>(parameters.poly_degree), dim, // displacement
1043  * FE_DGPMonomial<dim>(parameters.poly_degree - 1), 1, // pressure
1044  * FE_DGPMonomial<dim>(parameters.poly_degree - 1), 1), // dilatation
1045  * dof_handler(triangulation),
1046  * dofs_per_cell (fe.dofs_per_cell),
1047  * u_fe(first_u_component),
1048  * p_fe(p_component),
1049  * J_fe(J_component),
1050  * dofs_per_block(n_blocks),
1051  * qf_cell(parameters.quad_order),
1052  * qf_face(parameters.quad_order),
1053  * n_q_points (qf_cell.size()),
1054  * n_q_points_f (qf_face.size())
1055  * {
1056  * Assert(dim==2 || dim==3, ExcMessage("This problem only works in 2 or 3 space dimensions."));
1057  * determine_component_extractors();
1058  * }
1059  * template <int dim>
1060  * Solid<dim>::~Solid()
1061  * {
1062  * dof_handler.clear();
1063  * }
1064  * template <int dim>
1065  * void Solid<dim>::run()
1066  * {
1067  * LA::MPI::BlockVector solution_delta;
1068  *
1069  * make_grid();
1070  * setup_system(solution_delta);
1071  * {
1072  * ConstraintMatrix constraints;
1073  * constraints.close();
1075  * J_mask (J_component, n_components);
1076  * VectorTools::project (dof_handler,
1077  * constraints,
1078  * QGauss<dim>(degree+2),
1079  * J_mask,
1080  * solution_n);
1081  * }
1082  * output_results(time.get_timestep(), time.current());
1083  * time.increment();
1084  *
1085  * @endcode
1086  *
1087  * Some points for post-processing
1088  *
1089  * @code
1090  * std::vector<double> real_time;
1091  * real_time.push_back(0);
1092  * std::vector<std::vector<Point<dim> > > tracked_vertices (4);
1093  * {
1094  * Point<dim> p;
1095  * p[1] = parameters.length/2.0;
1096  * tracked_vertices[0].push_back(p*parameters.scale);
1097  * }
1098  * {
1099  * Point<dim> p;
1100  * p[1] = parameters.hole_diameter/2.0;
1101  * tracked_vertices[1].push_back(p*parameters.scale);
1102  * }
1103  * {
1104  * Point<dim> p;
1105  * p[0] = parameters.hole_diameter/2.0;
1106  * tracked_vertices[2].push_back(p*parameters.scale);
1107  * }
1108  * {
1109  * Point<dim> p;
1110  * p[0] = parameters.width/2.0;
1111  * tracked_vertices[3].push_back(p*parameters.scale);
1112  * }
1113  *
1114  * while (time.current() < time.end()+0.01*time.get_delta_t())
1115  * {
1116  * solve_nonlinear_timestep(solution_delta);
1117  * solution_n += solution_delta;
1118  * solution_delta = 0.0;
1119  * output_results(time.get_timestep(), time.current());
1120  * compute_vertex_positions(real_time,
1121  * tracked_vertices,
1122  * get_solution_total(solution_delta));
1123  * update_end_timestep();
1124  * time.increment();
1125  * }
1126  *
1127  * pcout << "\n\n*** Spatial position history for tracked vertices ***" << std::endl;
1128  * for (unsigned int t=0; t<real_time.size(); ++t)
1129  * {
1130  * if (t == 0)
1131  * {
1132  * pcout << "Time,";
1133  * for (unsigned int p=0; p<tracked_vertices.size(); ++p)
1134  * {
1135  * for (unsigned int d=0; d<dim; ++d)
1136  * {
1137  * pcout << "Point " << p << " [" << d << "]";
1138  * if (!(p == tracked_vertices.size()-1 && d == dim-1))
1139  * pcout << ",";
1140  * }
1141  * }
1142  * pcout << std::endl;
1143  * }
1144  *
1145  * pcout << std::setprecision(6);
1146  * pcout << real_time[t] << ",";
1147  * for (unsigned int p=0; p<tracked_vertices.size(); ++p)
1148  * {
1149  * Assert(tracked_vertices[p].size() == real_time.size(),
1150  * ExcMessage("Vertex not tracked at each timestep"));
1151  * for (unsigned int d=0; d<dim; ++d)
1152  * {
1153  * pcout << tracked_vertices[p][t][d];
1154  * if (!(p == tracked_vertices.size()-1 && d == dim-1))
1155  * pcout << ",";
1156  * }
1157  * }
1158  * pcout << std::endl;
1159  * }
1160  * }
1161  * template <int dim>
1162  * struct Solid<dim>::PerTaskData_ASM
1163  * {
1165  * Vector<double> cell_rhs;
1166  * std::vector<types::global_dof_index> local_dof_indices;
1167  * PerTaskData_ASM(const unsigned int dofs_per_cell)
1168  * :
1169  * cell_matrix(dofs_per_cell, dofs_per_cell),
1170  * cell_rhs(dofs_per_cell),
1171  * local_dof_indices(dofs_per_cell)
1172  * {}
1173  * void reset()
1174  * {
1175  * cell_matrix = 0.0;
1176  * cell_rhs = 0.0;
1177  * }
1178  * };
1179  * template <int dim>
1180  * struct Solid<dim>::ScratchData_ASM
1181  * {
1182  * const LA::MPI::BlockVector &solution_total;
1183  *
1184  * @endcode
1185  *
1186  * Integration helper
1187  *
1188  * @code
1189  * FEValues<dim> fe_values_ref;
1190  * FEFaceValues<dim> fe_face_values_ref;
1191  *
1192  * @endcode
1193  *
1194  * Quadrature point solution
1195  *
1196  * @code
1197  * std::vector<Tensor<2, dim> > solution_grads_u_total;
1198  * std::vector<double> solution_values_p_total;
1199  * std::vector<double> solution_values_J_total;
1200  *
1201  * @endcode
1202  *
1203  * Shape function values and gradients
1204  *
1205  * @code
1206  * std::vector<std::vector<double> > Nx;
1207  * std::vector<std::vector<Tensor<2, dim> > > grad_Nx;
1208  * std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
1209  *
1210  * ScratchData_ASM(const FiniteElement<dim> &fe_cell,
1211  * const QGauss<dim> &qf_cell, const UpdateFlags uf_cell,
1212  * const QGauss<dim - 1> & qf_face, const UpdateFlags uf_face,
1213  * const LA::MPI::BlockVector &solution_total)
1214  * :
1215  * solution_total (solution_total),
1216  * fe_values_ref(fe_cell, qf_cell, uf_cell),
1217  * fe_face_values_ref(fe_cell, qf_face, uf_face),
1218  * solution_grads_u_total(qf_cell.size()),
1219  * solution_values_p_total(qf_cell.size()),
1220  * solution_values_J_total(qf_cell.size()),
1221  * Nx(qf_cell.size(),
1222  * std::vector<double>(fe_cell.dofs_per_cell)),
1223  * grad_Nx(qf_cell.size(),
1224  * std::vector<Tensor<2, dim> >(fe_cell.dofs_per_cell)),
1225  * symm_grad_Nx(qf_cell.size(),
1226  * std::vector<SymmetricTensor<2, dim> >
1227  * (fe_cell.dofs_per_cell))
1228  * {}
1229  * ScratchData_ASM(const ScratchData_ASM &rhs)
1230  * :
1231  * solution_total (rhs.solution_total),
1232  * fe_values_ref(rhs.fe_values_ref.get_fe(),
1233  * rhs.fe_values_ref.get_quadrature(),
1234  * rhs.fe_values_ref.get_update_flags()),
1235  * fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
1236  * rhs.fe_face_values_ref.get_quadrature(),
1237  * rhs.fe_face_values_ref.get_update_flags()),
1238  * solution_grads_u_total(rhs.solution_grads_u_total),
1239  * solution_values_p_total(rhs.solution_values_p_total),
1240  * solution_values_J_total(rhs.solution_values_J_total),
1241  * Nx(rhs.Nx),
1242  * grad_Nx(rhs.grad_Nx),
1243  * symm_grad_Nx(rhs.symm_grad_Nx)
1244  * {}
1245  * void reset()
1246  * {
1247  * const unsigned int n_q_points = solution_grads_u_total.size();
1248  * const unsigned int n_dofs_per_cell = Nx[0].size();
1249  *
1250  * Assert(solution_grads_u_total.size() == n_q_points,
1251  * ExcInternalError());
1252  * Assert(solution_values_p_total.size() == n_q_points,
1253  * ExcInternalError());
1254  * Assert(solution_values_J_total.size() == n_q_points,
1255  * ExcInternalError());
1256  * Assert(Nx.size() == n_q_points,
1257  * ExcInternalError());
1258  * Assert(grad_Nx.size() == n_q_points,
1259  * ExcInternalError());
1260  * Assert(symm_grad_Nx.size() == n_q_points,
1261  * ExcInternalError());
1262  *
1263  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1264  * {
1265  * Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
1266  * Assert( grad_Nx[q_point].size() == n_dofs_per_cell,
1267  * ExcInternalError());
1268  * Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell,
1269  * ExcInternalError());
1270  *
1271  * solution_grads_u_total[q_point] = 0.0;
1272  * solution_values_p_total[q_point] = 0.0;
1273  * solution_values_J_total[q_point] = 0.0;
1274  * for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
1275  * {
1276  * Nx[q_point][k] = 0.0;
1277  * grad_Nx[q_point][k] = 0.0;
1278  * symm_grad_Nx[q_point][k] = 0.0;
1279  * }
1280  * }
1281  * }
1282  * };
1283  * template <>
1284  * void Solid<2>::make_grid()
1285  * {
1286  * const int dim = 2;
1287  * const double tol = 1e-12;
1288  * make_2d_quarter_plate_with_hole(triangulation,
1289  * parameters.length/2.0,
1290  * parameters.width/2.0,
1291  * parameters.hole_diameter/2.0,
1292  * parameters.n_repetitions_xy,
1293  * parameters.hole_division_fraction);
1294  *
1295  * @endcode
1296  *
1297  * Clear boundary ID's
1298  *
1299  * @code
1300  * for (typename Triangulation<dim>::active_cell_iterator
1301  * cell = triangulation.begin_active();
1302  * cell != triangulation.end(); ++cell)
1303  * {
1304  * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1305  * if (cell->face(face)->at_boundary())
1306  * {
1307  * cell->face(face)->set_all_boundary_ids(0);
1308  * }
1309  * }
1310  *
1311  * @endcode
1312  *
1313  * Set boundary IDs and and manifolds
1314  *
1315  * @code
1316  * const Point<dim> centre (0,0);
1317  * for (typename Triangulation<dim>::active_cell_iterator
1318  * cell = triangulation.begin_active();
1319  * cell != triangulation.end(); ++cell)
1320  * {
1321  * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1322  * if (cell->face(face)->at_boundary())
1323  * {
1324  * @endcode
1325  *
1326  * Set boundary IDs
1327  *
1328  * @code
1329  * if (std::abs(cell->face(face)->center()[0] - 0.0) < tol)
1330  * {
1331  * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X);
1332  * }
1333  * else if (std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol)
1334  * {
1335  * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X);
1336  * }
1337  * else if (std::abs(cell->face(face)->center()[1] - 0.0) < tol)
1338  * {
1339  * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y);
1340  * }
1341  * else if (std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol)
1342  * {
1343  * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y);
1344  * }
1345  * else
1346  * {
1347  * for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1348  * if (std::abs(cell->vertex(vertex).distance(centre) - parameters.hole_diameter/2.0) < tol)
1349  * {
1350  * cell->face(face)->set_boundary_id(parameters.boundary_id_hole);
1351  * break;
1352  * }
1353  * }
1354  *
1355  * @endcode
1356  *
1357  * Set manifold IDs
1358  *
1359  * @code
1360  * for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1361  * if (std::abs(cell->vertex(vertex).distance(centre) - parameters.hole_diameter/2.0) < tol)
1362  * {
1363  * cell->face(face)->set_manifold_id(parameters.manifold_id_hole);
1364  * break;
1365  * }
1366  * }
1367  * }
1368  * static SphericalManifold<dim> spherical_manifold (centre);
1369  * triangulation.set_manifold(parameters.manifold_id_hole,spherical_manifold);
1370  * triangulation.refine_global(parameters.global_refinement);
1371  * GridTools::scale(parameters.scale,triangulation);
1372  * }
1373  * template <>
1374  * void Solid<3>::make_grid()
1375  * {
1376  * const int dim = 3;
1377  * const double tol = 1e-12;
1378  * Triangulation<2> tria_2d;
1379  * make_2d_quarter_plate_with_hole(tria_2d,
1380  * parameters.length/2.0,
1381  * parameters.width/2.0,
1382  * parameters.hole_diameter/2.0,
1383  * parameters.n_repetitions_xy,
1384  * parameters.hole_division_fraction);
1385  * GridGenerator::extrude_triangulation(tria_2d,
1386  * parameters.n_repetitions_z+1,
1387  * parameters.thickness/2.0,
1388  * triangulation);
1389  *
1390  * @endcode
1391  *
1392  * Clear boundary ID's
1393  *
1394  * @code
1396  * cell = triangulation.begin_active();
1397  * cell != triangulation.end(); ++cell)
1398  * {
1399  * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1400  * if (cell->face(face)->at_boundary())
1401  * {
1402  * cell->face(face)->set_all_boundary_ids(0);
1403  * }
1404  * }
1405  *
1406  * @endcode
1407  *
1408  * Set boundary IDs and and manifolds
1409  *
1410  * @code
1411  * const Point<dim> direction (0,0,1);
1412  * const Point<dim> centre (0,0,0);
1414  * cell = triangulation.begin_active();
1415  * cell != triangulation.end(); ++cell)
1416  * {
1417  * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1418  * if (cell->face(face)->at_boundary())
1419  * {
1420  * @endcode
1421  *
1422  * Set boundary IDs
1423  *
1424  * @code
1425  * if (std::abs(cell->face(face)->center()[0] - 0.0) < tol)
1426  * {
1427  * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X);
1428  * }
1429  * else if (std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol)
1430  * {
1431  * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X);
1432  * }
1433  * else if (std::abs(cell->face(face)->center()[1] - 0.0) < tol)
1434  * {
1435  * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y);
1436  * }
1437  * else if (std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol)
1438  * {
1439  * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y);
1440  * }
1441  * else if (std::abs(cell->face(face)->center()[2] - 0.0) < tol)
1442  * {
1443  * cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Z);
1444  * }
1445  * else if (std::abs(cell->face(face)->center()[2] - parameters.thickness/2.0) < tol)
1446  * {
1447  * cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Z);
1448  * }
1449  * else
1450  * {
1451  * for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1452  * {
1453  * @endcode
1454  *
1455  * Project the cell vertex to the XY plane and
1456  * test the distance from the cylinder axis
1457  *
1458  * @code
1459  * Point<dim> vertex_proj = cell->vertex(vertex);
1460  * vertex_proj[2] = 0.0;
1461  * if (std::abs(vertex_proj.distance(centre) - parameters.hole_diameter/2.0) < tol)
1462  * {
1463  * cell->face(face)->set_boundary_id(parameters.boundary_id_hole);
1464  * break;
1465  * }
1466  * }
1467  * }
1468  *
1469  * @endcode
1470  *
1471  * Set manifold IDs
1472  *
1473  * @code
1474  * for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
1475  * {
1476  * @endcode
1477  *
1478  * Project the cell vertex to the XY plane and
1479  * test the distance from the cylinder axis
1480  *
1481  * @code
1482  * Point<dim> vertex_proj = cell->vertex(vertex);
1483  * vertex_proj[2] = 0.0;
1484  * if (std::abs(vertex_proj.distance(centre) - parameters.hole_diameter/2.0) < 1e-12)
1485  * {
1486  * @endcode
1487  *
1488  * Set manifold ID on face and edges
1489  *
1490  * @code
1491  * cell->face(face)->set_all_manifold_ids(parameters.manifold_id_hole);
1492  * break;
1493  * }
1494  * }
1495  * }
1496  * }
1497  * static CylindricalManifold<dim> cylindrical_manifold (direction,centre);
1498  * triangulation.set_manifold(parameters.manifold_id_hole,cylindrical_manifold);
1499  * triangulation.refine_global(parameters.global_refinement);
1500  * GridTools::scale(parameters.scale,triangulation);
1501  * }
1502  * template <int dim>
1503  * void Solid<dim>::make_2d_quarter_plate_with_hole(Triangulation<2> &tria_2d,
1504  * const double half_length,
1505  * const double half_width,
1506  * const double hole_radius,
1507  * const unsigned int n_repetitions_xy,
1508  * const double hole_division_fraction)
1509  * {
1510  * const double length = 2.0*half_length;
1511  * const double width = 2.0*half_width;
1512  * const double hole_diameter = 2.0*hole_radius;
1513  *
1514  * const double internal_width = hole_diameter + hole_division_fraction*(width - hole_diameter);
1515  * Triangulation<2> tria_quarter_plate_hole;
1516  * {
1517  * Triangulation<2> tria_plate_hole;
1519  * hole_diameter/2.0,
1520  * internal_width/2.0);
1521  *
1522  * std::set<typename Triangulation<2>::active_cell_iterator > cells_to_remove;
1524  * cell = tria_plate_hole.begin_active();
1525  * cell != tria_plate_hole.end(); ++cell)
1526  * {
1527  * @endcode
1528  *
1529  * Remove all cells that are not in the first quadrant
1530  *
1531  * @code
1532  * if (cell->center()[0] < 0.0 || cell->center()[1] < 0.0)
1533  * cells_to_remove.insert(cell);
1534  * }
1535  * Assert(cells_to_remove.size() > 0, ExcInternalError());
1536  * Assert(cells_to_remove.size() != tria_plate_hole.n_active_cells(), ExcInternalError());
1537  * GridGenerator::create_triangulation_with_removed_cells(tria_plate_hole,cells_to_remove,tria_quarter_plate_hole);
1538  * }
1539  *
1540  * Triangulation<2> tria_cut_plate;
1541  * {
1542  * Triangulation<2> tria_plate;
1543  * @endcode
1544  *
1545  * Subdivide the plate so that we're left one
1546  * cell to remove (we'll replace this with the
1547  * plate with the hole) and then make the
1548  * rest of the subdivisions so that we're left
1549  * with cells with a decent aspect ratio
1550  *
1551  * @code
1552  * std::vector<std::vector<double> > step_sizes;
1553  * {
1554  * std::vector<double> subdivision_width;
1555  * subdivision_width.push_back(internal_width/2.0);
1556  * const double width_remaining = (width - internal_width)/2.0;
1557  * const unsigned int n_subs = std::max(1.0,std::ceil(width_remaining/(internal_width/2.0)));
1558  * Assert(n_subs>0, ExcInternalError());
1559  * for (unsigned int s=0; s<n_subs; ++s)
1560  * subdivision_width.push_back(width_remaining/n_subs);
1561  * step_sizes.push_back(subdivision_width);
1562  *
1563  * const double sum_half_width = std::accumulate(subdivision_width.begin(), subdivision_width.end(), 0.0);
1564  * Assert(std::abs(sum_half_width-width/2.0) < 1e-12, ExcInternalError());
1565  * }
1566  * {
1567  * std::vector<double> subdivision_length;
1568  * subdivision_length.push_back(internal_width/2.0);
1569  * const double length_remaining = (length - internal_width)/2.0;
1570  * const unsigned int n_subs = std::max(1.0,std::ceil(length_remaining/(internal_width/2.0)));
1571  * Assert(n_subs>0, ExcInternalError());
1572  * for (unsigned int s=0; s<n_subs; ++s)
1573  * subdivision_length.push_back(length_remaining/n_subs);
1574  * step_sizes.push_back(subdivision_length);
1575  *
1576  * const double sum_half_length = std::accumulate(subdivision_length.begin(), subdivision_length.end(), 0.0);
1577  * Assert(std::abs(sum_half_length-length/2.0) < 1e-12, ExcInternalError());
1578  * }
1579  *
1580  * GridGenerator::subdivided_hyper_rectangle(tria_plate,
1581  * step_sizes,
1582  * Point<2>(0.0, 0.0),
1583  * Point<2>(width/2.0, length/2.0));
1584  *
1585  * std::set<typename Triangulation<2>::active_cell_iterator > cells_to_remove;
1586  * for (typename Triangulation<2>::active_cell_iterator
1587  * cell = tria_plate.begin_active();
1588  * cell != tria_plate.end(); ++cell)
1589  * {
1590  * @endcode
1591  *
1592  * Remove all cells that are in the first quadrant
1593  *
1594  * @code
1595  * if (cell->center()[0] < internal_width/2.0 && cell->center()[1] < internal_width/2.0)
1596  * cells_to_remove.insert(cell);
1597  * }
1598  * Assert(cells_to_remove.size() > 0, ExcInternalError());
1599  * Assert(cells_to_remove.size() != tria_plate.n_active_cells(), ExcInternalError());
1600  * GridGenerator::create_triangulation_with_removed_cells(tria_plate,cells_to_remove,tria_cut_plate);
1601  * }
1602  *
1603  * Triangulation<2> tria_2d_not_flat;
1604  * GridGenerator::merge_triangulations(tria_quarter_plate_hole,
1605  * tria_cut_plate,
1606  * tria_2d_not_flat);
1607  *
1608  * @endcode
1609  *
1610  * Attach a manifold to the curved boundary and refine
1611  * Note: We can only guarentee that the vertices sit on
1612  * the curve, so we must test with their position instead
1613  * of the cell centre.
1614  *
1615  * @code
1616  * const Point<2> centre_2d (0,0);
1617  * for (typename Triangulation<2>::active_cell_iterator
1618  * cell = tria_2d_not_flat.begin_active();
1619  * cell != tria_2d_not_flat.end(); ++cell)
1620  * {
1621  * for (unsigned int face=0; face<GeometryInfo<2>::faces_per_cell; ++face)
1622  * if (cell->face(face)->at_boundary())
1623  * for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
1624  * if (std::abs(cell->vertex(vertex).distance(centre_2d) - hole_diameter/2.0) < 1e-12)
1625  * {
1626  * cell->face(face)->set_manifold_id(10);
1627  * break;
1628  * }
1629  * }
1630  * SphericalManifold<2> spherical_manifold_2d (centre_2d);
1631  * tria_2d_not_flat.set_manifold(10,spherical_manifold_2d);
1632  * tria_2d_not_flat.refine_global(std::max (1U, n_repetitions_xy));
1633  * tria_2d_not_flat.reset_manifold(10); // Clear manifold
1634  *
1635  * GridGenerator::flatten_triangulation(tria_2d_not_flat,tria_2d);
1636  * }
1637  * template <int dim>
1638  * void Solid<dim>::setup_system(LA::MPI::BlockVector &solution_delta)
1639  * {
1640  * timer.enter_subsection("Setup system");
1641  * pcout << "Setting up linear system..." << std::endl;
1642  *
1643  * @endcode
1644  *
1645  * Partition triangulation
1646  *
1647  * @code
1648  * GridTools::partition_triangulation (n_mpi_processes,
1649  * triangulation);
1650  *
1651  * block_component = std::vector<unsigned int> (n_components, u_block); // Displacement
1652  * block_component[p_component] = p_block; // Pressure
1653  * block_component[J_component] = J_block; // Dilatation
1654  * dof_handler.distribute_dofs(fe);
1655  * DoFRenumbering::Cuthill_McKee(dof_handler);
1656  * DoFRenumbering::component_wise(dof_handler, block_component);
1657  *
1658  * @endcode
1659  *
1660  * Count DoFs in each block
1661  *
1662  * @code
1663  * dofs_per_block.clear();
1664  * dofs_per_block.resize(n_blocks);
1665  * DoFTools::count_dofs_per_block(dof_handler, dofs_per_block,
1666  * block_component);
1667  *
1668  * all_locally_owned_dofs = DoFTools::locally_owned_dofs_per_subdomain (dof_handler);
1669  * std::vector<IndexSet> all_locally_relevant_dofs
1670  * = DoFTools::locally_relevant_dofs_per_subdomain (dof_handler);
1671  *
1672  * locally_owned_dofs.clear();
1673  * locally_owned_partitioning.clear();
1674  * Assert(all_locally_owned_dofs.size() > this_mpi_process, ExcInternalError());
1675  * locally_owned_dofs = all_locally_owned_dofs[this_mpi_process];
1676  *
1677  * locally_relevant_dofs.clear();
1678  * locally_relevant_partitioning.clear();
1679  * Assert(all_locally_relevant_dofs.size() > this_mpi_process, ExcInternalError());
1680  * locally_relevant_dofs = all_locally_relevant_dofs[this_mpi_process];
1681  *
1682  * locally_owned_partitioning.reserve(n_blocks);
1683  * locally_relevant_partitioning.reserve(n_blocks);
1684  * for (unsigned int b=0; b<n_blocks; ++b)
1685  * {
1686  * const types::global_dof_index idx_begin
1687  * = std::accumulate(dofs_per_block.begin(),
1688  * std::next(dofs_per_block.begin(),b), 0);
1689  * const types::global_dof_index idx_end
1690  * = std::accumulate(dofs_per_block.begin(),
1691  * std::next(dofs_per_block.begin(),b+1), 0);
1692  * locally_owned_partitioning.push_back(locally_owned_dofs.get_view(idx_begin, idx_end));
1693  * locally_relevant_partitioning.push_back(locally_relevant_dofs.get_view(idx_begin, idx_end));
1694  * }
1695  *
1696  * pcout
1697  * << " Number of active cells: " << triangulation.n_active_cells()
1698  * << " (by partition:";
1699  * for (unsigned int p=0; p<n_mpi_processes; ++p)
1700  * pcout
1701  * << (p==0 ? ' ' : '+')
1702  * << (GridTools::count_cells_with_subdomain_association (triangulation,p));
1703  * pcout << ")" << std::endl;
1704  *
1705  * pcout
1706  * << " Number of degrees of freedom: " << dof_handler.n_dofs()
1707  * << " (by partition:";
1708  * for (unsigned int p=0; p<n_mpi_processes; ++p)
1709  * pcout
1710  * << (p==0 ? ' ' : '+')
1711  * << (DoFTools::count_dofs_with_subdomain_association (dof_handler,p));
1712  * pcout << ")" << std::endl;
1713  * pcout
1714  * << " Number of degrees of freedom per block: "
1715  * << "[n_u, n_p, n_J] = ["
1716  * << dofs_per_block[u_block] << ", "
1717  * << dofs_per_block[p_block] << ", "
1718  * << dofs_per_block[J_block] << "]"
1719  * << std::endl;
1720  *
1721  *
1722  * Table<2, DoFTools::Coupling> coupling(n_components, n_components);
1723  * for (unsigned int ii = 0; ii < n_components; ++ii)
1724  * for (unsigned int jj = 0; jj < n_components; ++jj)
1725  * if (((ii < p_component) && (jj == J_component))
1726  * || ((ii == J_component) && (jj < p_component))
1727  * || ((ii == p_component) && (jj == p_component)))
1728  * coupling[ii][jj] = DoFTools::none;
1729  * else
1730  * coupling[ii][jj] = DoFTools::always;
1731  *
1732  * TrilinosWrappers::BlockSparsityPattern bsp (locally_owned_partitioning,
1733  * locally_owned_partitioning,
1734  * locally_relevant_partitioning,
1735  * mpi_communicator);
1736  * DoFTools::make_sparsity_pattern (dof_handler, bsp,
1737  * constraints, false,
1738  * this_mpi_process);
1739  * bsp.compress();
1740  * tangent_matrix.reinit (bsp);
1741  *
1742  * @endcode
1743  *
1744  * We then set up storage vectors
1745  *
1746  * @code
1747  * system_rhs.reinit(locally_owned_partitioning,
1748  * mpi_communicator);
1749  * solution_n.reinit(locally_owned_partitioning,
1750  * mpi_communicator);
1751  * solution_delta.reinit(locally_owned_partitioning,
1752  * mpi_communicator);
1753  * setup_qph();
1754  * timer.leave_subsection();
1755  * }
1756  * template <int dim>
1757  * void
1758  * Solid<dim>::determine_component_extractors()
1759  * {
1760  * element_indices_u.clear();
1761  * element_indices_p.clear();
1762  * element_indices_J.clear();
1763  * for (unsigned int k = 0; k < fe.dofs_per_cell; ++k)
1764  * {
1765  * const unsigned int k_group = fe.system_to_base_index(k).first.first;
1766  * if (k_group == u_block)
1767  * element_indices_u.push_back(k);
1768  * else if (k_group == p_block)
1769  * element_indices_p.push_back(k);
1770  * else if (k_group == J_block)
1771  * element_indices_J.push_back(k);
1772  * else
1773  * {
1774  * Assert(k_group <= J_block, ExcInternalError());
1775  * }
1776  * }
1777  * }
1778  * template <int dim>
1779  * void Solid<dim>::setup_qph()
1780  * {
1781  * pcout << "Setting up quadrature point data..." << std::endl;
1782  * quadrature_point_history.initialize(triangulation.begin_active(),
1783  * triangulation.end(),
1784  * n_q_points);
1785  * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
1786  * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1787  * dof_handler.begin_active()),
1788  * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1789  * dof_handler.end());
1790  * for (; cell!=endc; ++cell)
1791  * {
1792  * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
1793  * const std::vector<std::shared_ptr<PointHistory<dim> > > lqph =
1794  * quadrature_point_history.get_data(cell);
1795  * Assert(lqph.size() == n_q_points, ExcInternalError());
1796  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1797  * lqph[q_point]->setup_lqp(parameters, time);
1798  * }
1799  * }
1800  * template <int dim>
1801  * void
1802  * Solid<dim>::solve_nonlinear_timestep(LA::MPI::BlockVector &solution_delta)
1803  * {
1804  * pcout << std::endl
1805  * << "Timestep " << time.get_timestep() << " @ "
1806  * << time.current() << "s of "
1807  * << time.end() << "s" << std::endl;
1808  * LA::MPI::BlockVector newton_update(locally_owned_partitioning,
1809  * mpi_communicator);
1810  * error_residual.reset();
1811  * error_residual_0.reset();
1812  * error_residual_norm.reset();
1813  * error_update.reset();
1814  * error_update_0.reset();
1815  * error_update_norm.reset();
1816  * print_conv_header();
1817  * unsigned int newton_iteration = 0;
1818  * for (; newton_iteration < parameters.max_iterations_NR;
1819  * ++newton_iteration)
1820  * {
1821  * pcout << " " << std::setw(2) << newton_iteration << " " << std::flush;
1822  * make_constraints(newton_iteration);
1823  * assemble_system(solution_delta);
1824  * get_error_residual(error_residual);
1825  * if (newton_iteration == 0)
1826  * error_residual_0 = error_residual;
1827  * error_residual_norm = error_residual;
1828  * error_residual_norm.normalise(error_residual_0);
1829  * if (newton_iteration > 0 &&
1830  * (error_update_norm.u <= parameters.tol_u &&
1831  * error_residual_norm.u <= parameters.tol_f) )
1832  * {
1833  * pcout << " CONVERGED! " << std::endl;
1834  * print_conv_footer(solution_delta);
1835  * break;
1836  * }
1837  * const std::pair<unsigned int, double>
1838  * lin_solver_output = solve_linear_system(newton_update);
1839  * get_error_update(newton_update, error_update);
1840  * if (newton_iteration == 0)
1841  * error_update_0 = error_update;
1842  * error_update_norm = error_update;
1843  * error_update_norm.normalise(error_update_0);
1844  * solution_delta += newton_update;
1845  * newton_update = 0.0;
1846  * pcout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
1847  * << std::scientific << lin_solver_output.first << " "
1848  * << lin_solver_output.second << " " << error_residual_norm.norm
1849  * << " " << error_residual_norm.u << " "
1850  * << error_residual_norm.p << " " << error_residual_norm.J
1851  * << " " << error_update_norm.norm << " " << error_update_norm.u
1852  * << " " << error_update_norm.p << " " << error_update_norm.J
1853  * << " " << std::endl;
1854  * }
1855  * AssertThrow (newton_iteration <= parameters.max_iterations_NR,
1856  * ExcMessage("No convergence in nonlinear solver!"));
1857  * }
1858  * template <int dim>
1859  * void Solid<dim>::print_conv_header()
1860  * {
1861  * pcout << std::string(132,'_') << std::endl;
1862  * pcout << " SOLVER STEP "
1863  * << " | LIN_IT LIN_RES RES_NORM "
1864  * << " RES_U RES_P RES_J NU_NORM "
1865  * << " NU_U NU_P NU_J " << std::endl;
1866  * pcout << std::string(132,'_') << std::endl;
1867  * }
1868  * template <int dim>
1869  * void Solid<dim>::print_conv_footer(const LA::MPI::BlockVector &solution_delta)
1870  * {
1871  * pcout << std::string(132,'_') << std::endl;
1872  * const std::pair<double,std::pair<double,double> > error_dil = get_error_dilation(get_solution_total(solution_delta));
1873  * pcout << "Relative errors:" << std::endl
1874  * << "Displacement:\t" << error_update.u / error_update_0.u << std::endl
1875  * << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
1876  * << "Dilatation:\t" << error_dil.first << std::endl
1877  * << "v / V_0:\t" << error_dil.second.second << " / " << error_dil.second.first
1878  * << " = " << (error_dil.second.second/error_dil.second.first) << std::endl;
1879  * }
1880  * template <int dim>
1881  * std::pair<double,std::pair<double,double> >
1882  * Solid<dim>::get_error_dilation(const LA::MPI::BlockVector &solution_total) const
1883  * {
1884  * double vol_reference = 0.0;
1885  * double vol_current = 0.0;
1886  * double dil_L2_error = 0.0;
1887  * FEValues<dim> fe_values_ref(fe, qf_cell,
1888  * update_values | update_gradients | update_JxW_values);
1889  * std::vector<Tensor<2, dim> > solution_grads_u_total (qf_cell.size());
1890  * std::vector<double> solution_values_J_total (qf_cell.size());
1891  * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
1892  * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1893  * dof_handler.begin_active()),
1894  * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
1895  * dof_handler.end());
1896  * for (; cell != endc; ++cell)
1897  * {
1898  * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
1899  * fe_values_ref.reinit(cell);
1900  * fe_values_ref[u_fe].get_function_gradients(solution_total,
1901  * solution_grads_u_total);
1902  * fe_values_ref[J_fe].get_function_values(solution_total,
1903  * solution_values_J_total);
1904  * const std::vector<std::shared_ptr<const PointHistory<dim> > > lqph =
1905  * quadrature_point_history.get_data(cell);
1906  * Assert(lqph.size() == n_q_points, ExcInternalError());
1907  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1908  * {
1909  * const double det_F_qp = determinant(Physics::Elasticity::Kinematics::F(solution_grads_u_total[q_point]));
1910  * const double J_tilde_qp = solution_values_J_total[q_point];
1911  * const double the_error_qp_squared = std::pow((det_F_qp - J_tilde_qp),
1912  * 2);
1913  * const double JxW = fe_values_ref.JxW(q_point);
1914  * dil_L2_error += the_error_qp_squared * JxW;
1915  * vol_reference += JxW;
1916  * vol_current += det_F_qp * JxW;
1917  * }
1918  * }
1919  * Assert(vol_current > 0.0, ExcInternalError());
1920  * @endcode
1921  *
1922  * Sum across all porcessors
1923  *
1924  * @code
1925  * dil_L2_error = Utilities::MPI::sum(dil_L2_error,mpi_communicator);
1926  * vol_reference = Utilities::MPI::sum(vol_reference,mpi_communicator);
1927  * vol_current = Utilities::MPI::sum(vol_current,mpi_communicator);
1928  *
1929  * return std::make_pair(std::sqrt(dil_L2_error),
1930  * std::make_pair(vol_reference,vol_current));
1931  * }
1932  * template <int dim>
1933  * void Solid<dim>::get_error_residual(Errors &error_residual)
1934  * {
1935  * @endcode
1936  *
1937  * Construct a residual vector that has the values for all of its
1938  * constrained DoFs set to zero.
1939  *
1940  * @code
1941  * LA::MPI::BlockVector error_res (system_rhs);
1942  * constraints.set_zero(error_res);
1943  * error_residual.norm = error_res.l2_norm();
1944  * error_residual.u = error_res.block(u_block).l2_norm();
1945  * error_residual.p = error_res.block(p_block).l2_norm();
1946  * error_residual.J = error_res.block(J_block).l2_norm();
1947  * }
1948  * template <int dim>
1949  * void Solid<dim>::get_error_update(const LA::MPI::BlockVector &newton_update,
1950  * Errors &error_update)
1951  * {
1952  * @endcode
1953  *
1954  * Construct a update vector that has the values for all of its
1955  * constrained DoFs set to zero.
1956  *
1957  * @code
1958  * LA::MPI::BlockVector error_ud (newton_update);
1959  * constraints.set_zero(error_ud);
1960  * error_update.norm = error_ud.l2_norm();
1961  * error_update.u = error_ud.block(u_block).l2_norm();
1962  * error_update.p = error_ud.block(p_block).l2_norm();
1963  * error_update.J = error_ud.block(J_block).l2_norm();
1964  * }
1965  * template <int dim>
1966  * LA::MPI::BlockVector
1967  * Solid<dim>::get_solution_total(const LA::MPI::BlockVector &solution_delta) const
1968  * {
1969  * @endcode
1970  *
1971  * Cell interpolation -> Ghosted vector
1972  *
1973  * @code
1974  * LA::MPI::BlockVector solution_total (locally_owned_partitioning,
1975  * locally_relevant_partitioning,
1976  * mpi_communicator,
1977  * /*vector_writable = */ false);
1978  * LA::MPI::BlockVector tmp (solution_total);
1979  * solution_total = solution_n;
1980  * tmp = solution_delta;
1981  * solution_total += tmp;
1982  * return solution_total;
1983  * }
1984  * template <int dim>
1985  * void Solid<dim>::assemble_system(const LA::MPI::BlockVector &solution_delta)
1986  * {
1987  * timer.enter_subsection("Assemble system");
1988  * pcout << " ASM_SYS " << std::flush;
1989  * tangent_matrix = 0.0;
1990  * system_rhs = 0.0;
1991  * const LA::MPI::BlockVector solution_total(get_solution_total(solution_delta));
1992  * const UpdateFlags uf_cell(update_values |
1993  * update_gradients |
1994  * update_JxW_values);
1995  * const UpdateFlags uf_face(update_values |
1996  * update_normal_vectors |
1997  * update_JxW_values);
1998  * PerTaskData_ASM per_task_data(dofs_per_cell);
1999  * ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total);
2000  *
2001  * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
2002  * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2003  * dof_handler.begin_active()),
2004  * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2005  * dof_handler.end());
2006  * for (; cell != endc; ++cell)
2007  * {
2008  * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2009  * assemble_system_one_cell(cell, scratch_data, per_task_data);
2010  * copy_local_to_global_system(per_task_data);
2011  * }
2012  * tangent_matrix.compress(VectorOperation::add);
2013  * system_rhs.compress(VectorOperation::add);
2014  * timer.leave_subsection();
2015  * }
2016  * template <int dim>
2017  * void Solid<dim>::copy_local_to_global_system(const PerTaskData_ASM &data)
2018  * {
2019  * constraints.distribute_local_to_global(data.cell_matrix, data.cell_rhs,
2020  * data.local_dof_indices,
2021  * tangent_matrix, system_rhs);
2022  * }
2023  * template <int dim>
2024  * void
2025  * Solid<dim>::assemble_system_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
2026  * ScratchData_ASM &scratch,
2027  * PerTaskData_ASM &data) const
2028  * {
2029  * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2030  *
2031  * data.reset();
2032  * scratch.reset();
2033  * scratch.fe_values_ref.reinit(cell);
2034  * cell->get_dof_indices(data.local_dof_indices);
2035  * const std::vector<std::shared_ptr<const PointHistory<dim> > > lqph =
2036  * quadrature_point_history.get_data(cell);
2037  * Assert(lqph.size() == n_q_points, ExcInternalError());
2038  *
2039  * @endcode
2040  *
2041  * Update quadrature point solution
2042  *
2043  * @code
2044  * scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
2045  * scratch.solution_grads_u_total);
2046  * scratch.fe_values_ref[p_fe].get_function_values(scratch.solution_total,
2047  * scratch.solution_values_p_total);
2048  * scratch.fe_values_ref[J_fe].get_function_values(scratch.solution_total,
2049  * scratch.solution_values_J_total);
2050  *
2051  * @endcode
2052  *
2053  * Update shape functions and their gradients (push-forward)
2054  *
2055  * @code
2056  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2057  * {
2058  * const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]);
2059  * const Tensor<2, dim> F_inv = invert(F);
2060  *
2061  * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2062  * {
2063  * const unsigned int k_group = fe.system_to_base_index(k).first.first;
2064  * if (k_group == u_block)
2065  * {
2066  * scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point)
2067  * * F_inv;
2068  * scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
2069  * }
2070  * else if (k_group == p_block)
2071  * scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k,
2072  * q_point);
2073  * else if (k_group == J_block)
2074  * scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k,
2075  * q_point);
2076  * else
2077  * Assert(k_group <= J_block, ExcInternalError());
2078  * }
2079  * }
2080  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2081  * {
2082  * const SymmetricTensor<2, dim> &I = Physics::Elasticity::StandardTensors<dim>::I;
2083  * const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]);
2084  * const double det_F = determinant(F);
2085  * const double &p_tilde = scratch.solution_values_p_total[q_point];
2086  * const double &J_tilde = scratch.solution_values_J_total[q_point];
2087  * Assert(det_F > 0, ExcInternalError());
2088  *
2089  * {
2090  * PointHistory<dim> *lqph_q_point_nc = const_cast<PointHistory<dim>*>(lqph[q_point].get());
2091  * lqph_q_point_nc->update_internal_equilibrium(F,p_tilde,J_tilde);
2092  * }
2093  *
2094  * const SymmetricTensor<2, dim> tau = lqph[q_point]->get_tau(F,p_tilde);
2095  * const Tensor<2, dim> tau_ns (tau);
2096  * const SymmetricTensor<4, dim> Jc = lqph[q_point]->get_Jc(F,p_tilde);
2097  * const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ(J_tilde);
2098  * const double d2Psi_vol_dJ2 = lqph[q_point]->get_d2Psi_vol_dJ2(J_tilde);
2099  *
2100  * const std::vector<double> &Nx = scratch.Nx[q_point];
2101  * const std::vector<Tensor<2, dim> > &grad_Nx = scratch.grad_Nx[q_point];
2102  * const std::vector<SymmetricTensor<2, dim> > &symm_grad_Nx = scratch.symm_grad_Nx[q_point];
2103  * const double JxW = scratch.fe_values_ref.JxW(q_point);
2104  *
2105  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2106  * {
2107  * const unsigned int component_i = fe.system_to_component_index(i).first;
2108  * const unsigned int i_group = fe.system_to_base_index(i).first.first;
2109  * if (i_group == u_block)
2110  * data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
2111  * else if (i_group == p_block)
2112  * data.cell_rhs(i) -= Nx[i] * (det_F - J_tilde) * JxW;
2113  * else if (i_group == J_block)
2114  * data.cell_rhs(i) -= Nx[i] * (dPsi_vol_dJ - p_tilde) * JxW;
2115  * else
2116  * Assert(i_group <= J_block, ExcInternalError());
2117  *
2118  * for (unsigned int j = 0; j <= i; ++j)
2119  * {
2120  * const unsigned int component_j = fe.system_to_component_index(j).first;
2121  * const unsigned int j_group = fe.system_to_base_index(j).first.first;
2122  * if ((i_group == u_block) && (j_group == u_block))
2123  * {
2124  * data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
2125  * * symm_grad_Nx[j] * JxW;
2126  * if (component_i == component_j) // geometrical stress contribution
2127  * data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau_ns
2128  * * grad_Nx[j][component_j] * JxW;
2129  * }
2130  * else if ((i_group == u_block) && (j_group == p_block))
2131  * {
2132  * data.cell_matrix(i, j) += (symm_grad_Nx[i] * I)
2133  * * Nx[j] * det_F
2134  * * JxW;
2135  * }
2136  * else if ((i_group == p_block) && (j_group == u_block))
2137  * {
2138  * data.cell_matrix(i, j) += Nx[i] * det_F
2139  * * (symm_grad_Nx[j] * I)
2140  * * JxW;
2141  * }
2142  * else if ((i_group == p_block) && (j_group == J_block))
2143  * data.cell_matrix(i, j) -= Nx[i] * Nx[j] * JxW;
2144  * else if ((i_group == J_block) && (j_group == p_block))
2145  * data.cell_matrix(i, j) -= Nx[i] * Nx[j] * JxW;
2146  * else if ((i_group == J_block) && (j_group == J_block))
2147  * data.cell_matrix(i, j) += Nx[i] * d2Psi_vol_dJ2 * Nx[j] * JxW;
2148  * else
2149  * Assert((i_group <= J_block) && (j_group <= J_block),
2150  * ExcInternalError());
2151  * }
2152  * }
2153  * }
2154  *
2155  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2156  * for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
2157  * data.cell_matrix(i, j) = data.cell_matrix(j, i);
2158  *
2159  * if (parameters.driver == "Neumann")
2160  * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2161  * ++face)
2162  * if (cell->face(face)->at_boundary() == true
2163  * && cell->face(face)->boundary_id() == parameters.boundary_id_plus_Y)
2164  * {
2165  * scratch.fe_face_values_ref.reinit(cell, face);
2166  * for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
2167  * ++f_q_point)
2168  * {
2169  * const Tensor<1, dim> &N =
2170  * scratch.fe_face_values_ref.normal_vector(f_q_point);
2171  * static const double pressure_nom = parameters.pressure
2172  * / (parameters.scale * parameters.scale);
2173  * const double time_ramp = (time.current() < parameters.load_time ?
2174  * time.current() / parameters.load_time : 1.0);
2175  * const double pressure = -pressure_nom * time_ramp;
2176  * const Tensor<1, dim> traction = pressure * N;
2177  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2178  * {
2179  * const unsigned int i_group =
2180  * fe.system_to_base_index(i).first.first;
2181  * if (i_group == u_block)
2182  * {
2183  * const unsigned int component_i =
2184  * fe.system_to_component_index(i).first;
2185  * const double Ni =
2186  * scratch.fe_face_values_ref.shape_value(i,
2187  * f_q_point);
2188  * const double JxW = scratch.fe_face_values_ref.JxW(
2189  * f_q_point);
2190  * data.cell_rhs(i) += (Ni * traction[component_i])
2191  * * JxW;
2192  * }
2193  * }
2194  * }
2195  * }
2196  * }
2197  * template <int dim>
2198  * void Solid<dim>::make_constraints(const int &it_nr)
2199  * {
2200  * pcout << " CST " << std::flush;
2201  * if (it_nr > 1)
2202  * return;
2203  * constraints.clear();
2204  * const bool apply_dirichlet_bc = (it_nr == 0);
2205  * const FEValuesExtractors::Scalar x_displacement(0);
2206  * const FEValuesExtractors::Scalar y_displacement(1);
2207  * {
2208  * const int boundary_id = parameters.boundary_id_minus_X;
2209  * if (apply_dirichlet_bc == true)
2210  * VectorTools::interpolate_boundary_values(dof_handler,
2211  * boundary_id,
2212  * ZeroFunction<dim>(n_components),
2213  * constraints,
2214  * fe.component_mask(x_displacement));
2215  * else
2216  * VectorTools::interpolate_boundary_values(dof_handler,
2217  * boundary_id,
2218  * ZeroFunction<dim>(n_components),
2219  * constraints,
2220  * fe.component_mask(x_displacement));
2221  * }
2222  * {
2223  * const int boundary_id = parameters.boundary_id_minus_Y;
2224  * if (apply_dirichlet_bc == true)
2225  * VectorTools::interpolate_boundary_values(dof_handler,
2226  * boundary_id,
2227  * ZeroFunction<dim>(n_components),
2228  * constraints,
2229  * fe.component_mask(y_displacement));
2230  * else
2231  * VectorTools::interpolate_boundary_values(dof_handler,
2232  * boundary_id,
2233  * ZeroFunction<dim>(n_components),
2234  * constraints,
2235  * fe.component_mask(y_displacement));
2236  * }
2237  * if (dim==3)
2238  * {
2239  * const FEValuesExtractors::Scalar z_displacement(2);
2240  * {
2241  * const int boundary_id = parameters.boundary_id_minus_Z;
2242  * if (apply_dirichlet_bc == true)
2243  * VectorTools::interpolate_boundary_values(dof_handler,
2244  * boundary_id,
2245  * ZeroFunction<dim>(n_components),
2246  * constraints,
2247  * fe.component_mask(z_displacement));
2248  * else
2249  * VectorTools::interpolate_boundary_values(dof_handler,
2250  * boundary_id,
2251  * ZeroFunction<dim>(n_components),
2252  * constraints,
2253  * fe.component_mask(z_displacement));
2254  * }
2255  * {
2256  * const int boundary_id = parameters.boundary_id_plus_Z;
2257  * if (apply_dirichlet_bc == true)
2258  * VectorTools::interpolate_boundary_values(dof_handler,
2259  * boundary_id,
2260  * ZeroFunction<dim>(n_components),
2261  * constraints,
2262  * fe.component_mask(z_displacement));
2263  * else
2264  * VectorTools::interpolate_boundary_values(dof_handler,
2265  * boundary_id,
2266  * ZeroFunction<dim>(n_components),
2267  * constraints,
2268  * fe.component_mask(z_displacement));
2269  * }
2270  * }
2271  * if (parameters.driver == "Dirichlet")
2272  * {
2273  * const int boundary_id = parameters.boundary_id_plus_Y;
2274  * if (apply_dirichlet_bc == true)
2275  * {
2276  *
2277  * if (time.current() < parameters.load_time+0.01*time.get_delta_t())
2278  * {
2279  * const double delta_length = parameters.length*(parameters.stretch - 1.0)*parameters.scale;
2280  * const unsigned int n_stretch_steps = parameters.load_time/time.get_delta_t();
2281  * const double delta_u_y = delta_length/2.0/n_stretch_steps;
2282  * VectorTools::interpolate_boundary_values(dof_handler,
2283  * boundary_id,
2284  * ConstantFunction<dim>(delta_u_y,n_components),
2285  * constraints,
2286  * fe.component_mask(y_displacement));
2287  * }
2288  * else
2289  * VectorTools::interpolate_boundary_values(dof_handler,
2290  * boundary_id,
2291  * ZeroFunction<dim>(n_components),
2292  * constraints,
2293  * fe.component_mask(y_displacement));
2294  * }
2295  * else
2296  * VectorTools::interpolate_boundary_values(dof_handler,
2297  * boundary_id,
2298  * ZeroFunction<dim>(n_components),
2299  * constraints,
2300  * fe.component_mask(y_displacement));
2301  * }
2302  * constraints.close();
2303  * }
2304  * template <int dim>
2305  * std::pair<unsigned int, double>
2306  * Solid<dim>::solve_linear_system(LA::MPI::BlockVector &newton_update)
2307  * {
2308  * unsigned int lin_it = 0;
2309  * double lin_res = 0.0;
2310  *
2311  * timer.enter_subsection("Linear solver");
2312  * pcout << " SLV " << std::flush;
2313  *
2314  * const LA::MPI::Vector &f_u = system_rhs.block(u_block);
2315  * const LA::MPI::Vector &f_p = system_rhs.block(p_block);
2316  * const LA::MPI::Vector &f_J = system_rhs.block(J_block);
2317  * LA::MPI::Vector &d_u = newton_update.block(u_block);
2318  * LA::MPI::Vector &d_p = newton_update.block(p_block);
2319  * LA::MPI::Vector &d_J = newton_update.block(J_block);
2320  * const auto K_uu = linear_operator<LA::MPI::Vector>(tangent_matrix.block(u_block, u_block));
2321  * const auto K_up = linear_operator<LA::MPI::Vector>(tangent_matrix.block(u_block, p_block));
2322  * const auto K_pu = linear_operator<LA::MPI::Vector>(tangent_matrix.block(p_block, u_block));
2323  * const auto K_Jp = linear_operator<LA::MPI::Vector>(tangent_matrix.block(J_block, p_block));
2324  * const auto K_JJ = linear_operator<LA::MPI::Vector>(tangent_matrix.block(J_block, J_block));
2325  *
2326  * LA::PreconditionJacobi preconditioner_K_Jp_inv;
2327  * preconditioner_K_Jp_inv.initialize(
2328  * tangent_matrix.block(J_block, p_block),
2329  * LA::PreconditionJacobi::AdditionalData());
2330  * ReductionControl solver_control_K_Jp_inv (
2331  * tangent_matrix.block(J_block, p_block).m() * parameters.max_iterations_lin,
2332  * 1.0e-30, 1e-6);
2333  * ::SolverCG<LA::MPI::Vector> solver_K_Jp_inv (solver_control_K_Jp_inv);
2334  *
2335  * const auto K_Jp_inv = inverse_operator(K_Jp,
2336  * solver_K_Jp_inv,
2337  * preconditioner_K_Jp_inv);
2338  * const auto K_pJ_inv = transpose_operator(K_Jp_inv);
2339  * const auto K_pp_bar = K_Jp_inv * K_JJ * K_pJ_inv;
2340  * const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu;
2341  * const auto K_uu_con = K_uu + K_uu_bar_bar;
2342  *
2343  * LA::PreconditionAMG preconditioner_K_con_inv;
2344  * preconditioner_K_con_inv.initialize(
2345  * tangent_matrix.block(u_block, u_block),
2346  * LA::PreconditionAMG::AdditionalData(
2347  * true /*elliptic*/,
2348  * (parameters.poly_degree > 1 /*higher_order_elements*/)) );
2349  * ReductionControl solver_control_K_con_inv (
2350  * tangent_matrix.block(u_block, u_block).m() * parameters.max_iterations_lin,
2351  * 1.0e-30, parameters.tol_lin);
2352  * ::SolverSelector<LA::MPI::Vector> solver_K_con_inv;
2353  * solver_K_con_inv.select(parameters.type_lin);
2354  * solver_K_con_inv.set_control(solver_control_K_con_inv);
2355  * const auto K_uu_con_inv = inverse_operator(K_uu_con,
2356  * solver_K_con_inv,
2357  * preconditioner_K_con_inv);
2358  *
2359  * d_u = K_uu_con_inv*(f_u - K_up*(K_Jp_inv*f_J - K_pp_bar*f_p));
2360  * lin_it = solver_control_K_con_inv.last_step();
2361  * lin_res = solver_control_K_con_inv.last_value();
2362  * timer.leave_subsection();
2363  *
2364  * timer.enter_subsection("Linear solver postprocessing");
2365  * d_J = K_pJ_inv*(f_p - K_pu*d_u);
2366  * d_p = K_Jp_inv*(f_J - K_JJ*d_J);
2367  * timer.leave_subsection();
2368  *
2369  * constraints.distribute(newton_update);
2370  * return std::make_pair(lin_it, lin_res);
2371  * }
2372  * template <int dim>
2373  * void
2374  * Solid<dim>::update_end_timestep ()
2375  * {
2376  * FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
2377  * cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2378  * dof_handler.begin_active()),
2379  * endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),
2380  * dof_handler.end());
2381  * for (; cell != endc; ++cell)
2382  * {
2383  * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2384  * const std::vector<std::shared_ptr<PointHistory<dim> > > lqph =
2385  * quadrature_point_history.get_data(cell);
2386  * Assert(lqph.size() == n_q_points, ExcInternalError());
2387  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2388  * lqph[q_point]->update_end_timestep();
2389  * }
2390  * }
2391  *
2392  * template<int dim, class DH=DoFHandler<dim> >
2393  * class FilteredDataOut : public DataOut<dim, DH>
2394  * {
2395  * public:
2396  * FilteredDataOut (const unsigned int subdomain_id)
2397  * :
2398  * subdomain_id (subdomain_id)
2399  * {}
2400  *
2401  * virtual ~FilteredDataOut() {}
2402  *
2403  * virtual typename DataOut<dim, DH>::cell_iterator
2404  * first_cell ()
2405  * {
2406  * typename DataOut<dim, DH>::active_cell_iterator
2407  * cell = this->dofs->begin_active();
2408  * while ((cell != this->dofs->end()) &&
2409  * (cell->subdomain_id() != subdomain_id))
2410  * ++cell;
2411  * return cell;
2412  * }
2413  *
2414  * virtual typename DataOut<dim, DH>::cell_iterator
2415  * next_cell (const typename DataOut<dim, DH>::cell_iterator &old_cell)
2416  * {
2417  * if (old_cell != this->dofs->end())
2418  * {
2419  * const IteratorFilters::SubdomainEqualTo predicate(subdomain_id);
2420  * return
2421  * ++(FilteredIterator
2422  * <typename DataOut<dim, DH>::active_cell_iterator>
2423  * (predicate,old_cell));
2424  * }
2425  * else
2426  * return old_cell;
2427  * }
2428  *
2429  * private:
2430  * const unsigned int subdomain_id;
2431  * };
2432  *
2433  * template <int dim>
2434  * void Solid<dim>::output_results(const unsigned int timestep,
2435  * const double current_time) const
2436  * {
2437  * @endcode
2438  *
2439  * Output -> Ghosted vector
2440  *
2441  * @code
2442  * LA::MPI::BlockVector solution_total (locally_owned_partitioning,
2443  * locally_relevant_partitioning,
2444  * mpi_communicator,
2445  * /*vector_writable = */ false);
2446  * LA::MPI::BlockVector residual (locally_owned_partitioning,
2447  * locally_relevant_partitioning,
2448  * mpi_communicator,
2449  * /*vector_writable = */ false);
2450  * solution_total = solution_n;
2451  * residual = system_rhs;
2452  * residual *= -1.0;
2453  *
2454  * @endcode
2455  *
2456  * --- Additional data ---
2457  *
2458  * @code
2459  * Vector<double> material_id;
2460  * Vector<double> polynomial_order;
2461  * material_id.reinit(triangulation.n_active_cells());
2462  * polynomial_order.reinit(triangulation.n_active_cells());
2463  * std::vector<types::subdomain_id> partition_int (triangulation.n_active_cells());
2464  *
2465  * FilteredDataOut<dim> data_out(this_mpi_process);
2466  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2467  * data_component_interpretation(dim,
2468  * DataComponentInterpretation::component_is_part_of_vector);
2469  * data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
2470  * data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
2471  *
2472  * GridTools::get_subdomain_association (triangulation, partition_int);
2473  *
2474  * @endcode
2475  *
2476  * Can't use filtered iterators here because the cell
2477  * count "c" is incorrect for the parallel case
2478  *
2479  * @code
2480  * unsigned int c = 0;
2482  * cell = dof_handler.begin_active(),
2483  * endc = dof_handler.end();
2484  * for (; cell!=endc; ++cell, ++c)
2485  * {
2486  * if (cell->subdomain_id() != this_mpi_process) continue;
2487  *
2488  * material_id(c) = static_cast<int>(cell->material_id());
2489  * }
2490  *
2491  * std::vector<std::string> solution_name(n_components, "solution_");
2492  * std::vector<std::string> residual_name(n_components, "residual_");
2493  * for (unsigned int c=0; c<n_components; ++c)
2494  * {
2495  * if (block_component[c] == u_block)
2496  * {
2497  * solution_name[c] += "u";
2498  * residual_name[c] += "u";
2499  * }
2500  * else if (block_component[c] == p_block)
2501  * {
2502  * solution_name[c] += "p";
2503  * residual_name[c] += "p";
2504  * }
2505  * else if (block_component[c] == J_block)
2506  * {
2507  * solution_name[c] += "J";
2508  * residual_name[c] += "J";
2509  * }
2510  * else
2511  * {
2512  * Assert(c <= J_block, ExcInternalError());
2513  * }
2514  * }
2515  *
2516  * data_out.attach_dof_handler(dof_handler);
2517  * data_out.add_data_vector(solution_total,
2518  * solution_name,
2520  * data_component_interpretation);
2521  * data_out.add_data_vector(residual,
2522  * residual_name,
2524  * data_component_interpretation);
2525  * const Vector<double> partitioning(partition_int.begin(),
2526  * partition_int.end());
2527  * data_out.add_data_vector (material_id, "material_id");
2528  * data_out.add_data_vector (partitioning, "partitioning");
2529  * data_out.build_patches(degree);
2530  *
2531  * struct Filename
2532  * {
2533  * static std::string get_filename_vtu (unsigned int process,
2534  * unsigned int timestep,
2535  * const unsigned int n_digits = 4)
2536  * {
2537  * std::ostringstream filename_vtu;
2538  * filename_vtu
2539  * << "solution-"
2540  * << (std::to_string(dim) + "d")
2541  * << "."
2542  * << Utilities::int_to_string (process, n_digits)
2543  * << "."
2544  * << Utilities::int_to_string(timestep, n_digits)
2545  * << ".vtu";
2546  * return filename_vtu.str();
2547  * }
2548  *
2549  * static std::string get_filename_pvtu (unsigned int timestep,
2550  * const unsigned int n_digits = 4)
2551  * {
2552  * std::ostringstream filename_vtu;
2553  * filename_vtu
2554  * << "solution-"
2555  * << (std::to_string(dim) + "d")
2556  * << "."
2557  * << Utilities::int_to_string(timestep, n_digits)
2558  * << ".pvtu";
2559  * return filename_vtu.str();
2560  * }
2561  *
2562  * static std::string get_filename_pvd (void)
2563  * {
2564  * std::ostringstream filename_vtu;
2565  * filename_vtu
2566  * << "solution-"
2567  * << (std::to_string(dim) + "d")
2568  * << ".pvd";
2569  * return filename_vtu.str();
2570  * }
2571  * };
2572  *
2573  * @endcode
2574  *
2575  * Write out main data file
2576  *
2577  * @code
2578  * const std::string filename_vtu = Filename::get_filename_vtu(this_mpi_process, timestep);
2579  * std::ofstream output(filename_vtu.c_str());
2580  * data_out.write_vtu(output);
2581  *
2582  * @endcode
2583  *
2584  * Collection of files written in parallel
2585  * This next set of steps should only be performed
2586  * by master process
2587  *
2588  * @code
2589  * if (this_mpi_process == 0)
2590  * {
2591  * @endcode
2592  *
2593  * List of all files written out at this timestep by all processors
2594  *
2595  * @code
2596  * std::vector<std::string> parallel_filenames_vtu;
2597  * for (unsigned int p=0; p < n_mpi_processes; ++p)
2598  * {
2599  * parallel_filenames_vtu.push_back(Filename::get_filename_vtu(p, timestep));
2600  * }
2601  *
2602  * const std::string filename_pvtu (Filename::get_filename_pvtu(timestep));
2603  * std::ofstream pvtu_master(filename_pvtu.c_str());
2604  * data_out.write_pvtu_record(pvtu_master,
2605  * parallel_filenames_vtu);
2606  *
2607  * @endcode
2608  *
2609  * Time dependent data master file
2610  *
2611  * @code
2612  * static std::vector<std::pair<double,std::string> > time_and_name_history;
2613  * time_and_name_history.push_back (std::make_pair (current_time,
2614  * filename_pvtu));
2615  * const std::string filename_pvd (Filename::get_filename_pvd());
2616  * std::ofstream pvd_output (filename_pvd.c_str());
2617  * DataOutBase::write_pvd_record (pvd_output, time_and_name_history);
2618  * }
2619  * }
2620  * template <int dim>
2621  * void Solid<dim>::compute_vertex_positions(std::vector<double> &real_time,
2622  * std::vector<std::vector<Point<dim> > > &tracked_vertices,
2623  * const LA::MPI::BlockVector &solution_total) const
2624  * {
2625  * real_time.push_back(time.current());
2626  *
2627  * std::vector<bool> vertex_found (tracked_vertices.size(), false);
2628  * std::vector<Tensor<1,dim> > vertex_update (tracked_vertices.size());
2629  *
2632  * dof_handler.begin_active()),
2634  * dof_handler.end());
2635  * for (; cell != endc; ++cell)
2636  * {
2637  * Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
2638  * for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
2639  * {
2640  * for (unsigned int p=0; p<tracked_vertices.size(); ++p)
2641  * {
2642  * if (vertex_found[p] == true) continue;
2643  *
2644  * const Point<dim> pt_ref = tracked_vertices[p][0];
2645  * if (cell->vertex(v).distance(pt_ref) < 1e-6*parameters.scale)
2646  * {
2647  * for (unsigned int d=0; d<dim; ++d)
2648  * vertex_update[p][d] = solution_total(cell->vertex_dof_index(v,u_block+d));
2649  *
2650  * vertex_found[p] = true;
2651  * }
2652  * }
2653  * }
2654  * }
2655  *
2656  * for (unsigned int p=0; p<tracked_vertices.size(); ++p)
2657  * {
2658  * const int found_on_n_processes = Utilities::MPI::sum(int(vertex_found[p]), mpi_communicator);
2659  * Assert(found_on_n_processes>0, ExcMessage("Vertex not found on any processor"));
2660  * Tensor<1,dim> update;
2661  * for (unsigned int d=0; d<dim; ++d)
2662  * update[d] = Utilities::MPI::sum(vertex_update[p][d], mpi_communicator);
2663  * update /= found_on_n_processes;
2664  * tracked_vertices[p].push_back(tracked_vertices[p][0] + update);
2665  * }
2666  *
2667  * }
2668  * }
2669  * int main (int argc, char *argv[])
2670  * {
2671  * using namespace dealii;
2672  * using namespace ViscoElasStripHole;
2673  *
2674  * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
2675  *
2676  * try
2677  * {
2678  * const unsigned int dim = 2; // Works in both 2d and 3d
2679  * Solid<dim> solid("parameters.prm");
2680  * solid.run();
2681  * }
2682  * catch (std::exception &exc)
2683  * {
2684  * if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
2685  * {
2686  * std::cerr << std::endl << std::endl
2687  * << "----------------------------------------------------"
2688  * << std::endl;
2689  * std::cerr << "Exception on processing: " << std::endl << exc.what()
2690  * << std::endl << "Aborting!" << std::endl
2691  * << "----------------------------------------------------"
2692  * << std::endl;
2693  * return 1;
2694  * }
2695  * }
2696  * catch (...)
2697  * {
2698  * if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
2699  * {
2700  * std::cerr << std::endl << std::endl
2701  * << "----------------------------------------------------"
2702  * << std::endl;
2703  * std::cerr << "Unknown exception!" << std::endl << "Aborting!"
2704  * << std::endl
2705  * << "----------------------------------------------------"
2706  * << std::endl;
2707  * return 1;
2708  * }
2709  * }
2710  * return 0;
2711  * }
2712  * @endcode
2713 
2714 
2715 */
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
IndexSet
Definition: index_set.h:74
GridGenerator::hyper_cube_with_cylindrical_hole
void hyper_cube_with_cylindrical_hole(Triangulation< dim > &triangulation, const double inner_radius=.25, const double outer_radius=.5, const double L=.5, const unsigned int repetitions=1, const bool colorize=false)
FE_Q
Definition: fe_q.h:554
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
SolverSelector
Definition: solver_selector.h:91
SymmetricTensor< 2, dim >
dealii
Definition: namespace_dealii.h:25
Triangulation
Definition: tria.h:1109
determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2645
Physics::Elasticity::StandardTensors
Definition: standard_tensors.h:47
FEValuesExtractors::Scalar
Definition: fe_values_extractors.h:95
FE_DGPMonomial
Definition: fe_dgp_monomial.h:286
TimerOutput::summary
@ summary
Definition: timer.h:605
TrilinosWrappers
Definition: types.h:161
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2360
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
MPI_Comm
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorTools::project
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >(0)), const bool project_to_boundary_first=false)
TimerOutput::wall_times
@ wall_times
Definition: timer.h:649
LocalIntegrators::Advection::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double >> &velocity, const double factor=1.)
Definition: advection.h:80
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
DoFHandler
Definition: dof_handler.h:205
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
FEValues< dim >
CylindricalManifold
Definition: manifold_lib.h:387
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
TimerOutput
Definition: timer.h:546
LinearOperator::linear_operator
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
Definition: linear_operator.h:1373
GridGenerator::cylinder
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
GridGenerator::create_triangulation_with_removed_cells
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
LinearAlgebraDealII::BlockVector
BlockVector< double > BlockVector
Definition: generic_linear_algebra.h:48
DataOutBase::write_pvd_record
void write_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string >> &times_and_names)
Definition: data_out_base.cc:5586
FEValuesExtractors::Vector
Definition: fe_values_extractors.h:150
types::material_id
unsigned int material_id
Definition: types.h:152
FilteredIterator
Definition: filtered_iterator.h:529
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
Triangulation::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12013
Tensor< 2, dim >
ComponentSelectFunction
Definition: function.h:560
outer_product
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Definition: symmetric_tensor.h:3520
LinearAlgebraDealII::BlockSparseMatrix
BlockSparseMatrix< double > BlockSparseMatrix
Definition: generic_linear_algebra.h:58
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
trace
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
FiniteElement< dim >
TriaActiveIterator
Definition: tria_iterator.h:759
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:66
QGauss
Definition: quadrature_lib.h:40
DoFHandler::end
cell_iterator end() const
Definition: dof_handler.cc:951
unsigned int
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AffineConstraints< double >
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
IteratorFilters::SubdomainEqualTo
Definition: filtered_iterator.h:153
DoFHandler::clear
virtual void clear()
Definition: dof_handler.cc:1352
std::pow
inline ::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &x, const Number p)
Definition: vectorization.h:5428
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
ConditionalOStream
Definition: conditional_ostream.h:81
CellDataStorage
Definition: quadrature_point_data.h:64
Point< dim >
Triangulation::n_active_cells
unsigned int n_active_cells() const
Definition: tria.cc:12675
ParameterHandler
Definition: parameter_handler.h:845
Quadrature::size
unsigned int size() const
FullMatrix< double >
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
FEFaceValues< dim >
Patterns::Selection
Definition: patterns.h:383
Quadrature
Definition: quadrature.h:85
first
Point< 2 > first
Definition: grid_out.cc:4352
DataOut
Definition: data_out.h:148
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:935
Vector< double >
AffineConstraints::close
void close()
FESystem
Definition: fe.h:44
ParameterHandler::parse_input
virtual void parse_input(std::istream &input, const std::string &filename="input file", const std::string &last_line="", const bool skip_undefined=false)
Definition: parameter_handler.cc:399
parallel
Definition: distributed.h:416
Patterns::Double
Definition: patterns.h:293
Triangulation::end
cell_iterator end() const
Definition: tria.cc:12079
invert
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:3467
Utilities::MPI::MPI_InitFinalize
Definition: mpi.h:828
symmetrize
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3547
Patterns::Integer
Definition: patterns.h:190