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