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