Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Linear_Elastic_Active_Skeletal_Muscle_Model.h
Go to the documentation of this file.
1 
269  *
270  * /*
271  * * Authors: Jean-Paul Pelteret, Tim Hamann,
272  * * University of Erlangen-Nuremberg, 2017
273  * *
274  * * The support of this work by the European Research Council (ERC) through
275  * * the Advanced Grant 289049 MOCOPOLY is gratefully acknowledged by the
276  * * first author.
277  * */
278  *
279  * @endcode
280  *
281  *
282  * <a name="Includefiles"></a>
283  * <h3>Include files</h3>
284  *
285 
286  *
287  *
288  * @code
289  * #include <deal.II/base/quadrature_lib.h>
290  * #include <deal.II/base/function.h>
291  * #include <deal.II/base/logstream.h>
292  * #include <deal.II/base/parameter_handler.h>
293  * #include <deal.II/lac/vector.h>
294  * #include <deal.II/lac/full_matrix.h>
295  * #include <deal.II/lac/sparse_matrix.h>
296  * #include <deal.II/lac/solver_cg.h>
297  * #include <deal.II/lac/precondition.h>
298  * #include <deal.II/lac/constraint_matrix.h>
299  * #include <deal.II/grid/grid_out.h>
300  * #include <deal.II/grid/manifold_lib.h>
301  * #include <deal.II/grid/tria.h>
302  * #include <deal.II/grid/grid_generator.h>
303  * #include <deal.II/grid/grid_refinement.h>
304  * #include <deal.II/grid/grid_tools.h>
305  * #include <deal.II/grid/tria_accessor.h>
306  * #include <deal.II/grid/tria_iterator.h>
307  * #include <deal.II/grid/tria_boundary_lib.h>
308  * #include <deal.II/dofs/dof_handler.h>
309  * #include <deal.II/dofs/dof_accessor.h>
310  * #include <deal.II/dofs/dof_tools.h>
311  * #include <deal.II/fe/fe_values.h>
312  * #include <deal.II/fe/fe_system.h>
313  * #include <deal.II/fe/fe_q.h>
314  * #include <deal.II/numerics/vector_tools.h>
315  * #include <deal.II/numerics/matrix_tools.h>
316  * #include <deal.II/numerics/data_out.h>
317  * #include <deal.II/numerics/error_estimator.h>
318  * #include <deal.II/physics/transformations.h>
319  *
320  * #include <fstream>
321  * #include <iostream>
322  * #include <vector>
323  *
324  * namespace LMM
325  * {
326  * using namespace dealii;
327  *
328  * @endcode
329  *
330  *
331  * <a name="Runtimeparameters"></a>
332  * <h3>Run-time parameters</h3>
333  *
334 
335  *
336  * There are several parameters that can be set in the code so we set up a
337  * ParameterHandler object to read in the choices at run-time.
338  *
339  * @code
340  * namespace Parameters
341  * {
342  * @endcode
343  *
344  *
345  * <a name="FiniteElementsystem"></a>
346  * <h4>Finite Element system</h4>
347  *
348 
349  *
350  * Here we specify the polynomial order used to approximate the solution.
351  * The quadrature order should be adjusted accordingly.
352  *
353  * @code
354  * struct FESystem
355  * {
356  * unsigned int poly_degree;
357  * unsigned int quad_order;
358  *
359  * static void
360  * declare_parameters(ParameterHandler &prm);
361  *
362  * void
363  * parse_parameters(ParameterHandler &prm);
364  * };
365  *
366  * void FESystem::declare_parameters(ParameterHandler &prm)
367  * {
368  * prm.enter_subsection("Finite element system");
369  * {
370  * prm.declare_entry("Polynomial degree", "1",
371  * Patterns::Integer(0),
372  * "Displacement system polynomial order");
373  *
374  * prm.declare_entry("Quadrature order", "2",
375  * Patterns::Integer(0),
376  * "Gauss quadrature order");
377  * }
378  * prm.leave_subsection();
379  * }
380  *
381  * void FESystem::parse_parameters(ParameterHandler &prm)
382  * {
383  * prm.enter_subsection("Finite element system");
384  * {
385  * poly_degree = prm.get_integer("Polynomial degree");
386  * quad_order = prm.get_integer("Quadrature order");
387  * }
388  * prm.leave_subsection();
389  * }
390  *
391  * @endcode
392  *
393  *
394  * <a name="Problem"></a>
395  * <h4>Problem</h4>
396  *
397 
398  *
399  * Choose which problem is going to be solved
400  *
401  * @code
402  * struct Problem
403  * {
404  * std::string problem;
405  *
406  * static void
407  * declare_parameters(ParameterHandler &prm);
408  *
409  * void
410  * parse_parameters(ParameterHandler &prm);
411  * };
412  *
413  * void Problem::declare_parameters(ParameterHandler &prm)
414  * {
415  * prm.enter_subsection("Problem");
416  * {
417  * prm.declare_entry("Problem", "IsotonicContraction",
418  * Patterns::Selection("IsotonicContraction|BicepsBrachii"),
419  * "The problem that is to be solved");
420  * }
421  * prm.leave_subsection();
422  * }
423  *
424  * void Problem::parse_parameters(ParameterHandler &prm)
425  * {
426  * prm.enter_subsection("Problem");
427  * {
428  * problem = prm.get("Problem");
429  * }
430  * prm.leave_subsection();
431  * }
432  *
433  * @endcode
434  *
435  *
436  * <a name="IsotonicContractionGeometry"></a>
437  * <h4>IsotonicContractionGeometry</h4>
438  *
439 
440  *
441  * Make adjustments to the geometry and discretisation of the
442  * isotonic contraction model from Martins2006.
443  *
444 
445  *
446  *
447  * @code
448  * struct IsotonicContraction
449  * {
450  * const double half_length_x = 10e-3/2.0;
451  * const double half_length_y = 10e-3/2.0;
452  * const double half_length_z = 1e-3/2.0;
453  * const types::boundary_id bid_CC_dirichlet_symm_X = 1;
454  * const types::boundary_id bid_CC_dirichlet_symm_Z = 2;
455  * const types::boundary_id bid_CC_neumann = 10;
456  *
457  * static void
458  * declare_parameters(ParameterHandler &prm);
459  *
460  * void
461  * parse_parameters(ParameterHandler &prm);
462  * };
463  *
464  * void IsotonicContraction::declare_parameters(ParameterHandler &/*prm*/)
465  * {
466  *
467  * }
468  *
469  * void IsotonicContraction::parse_parameters(ParameterHandler &/*prm*/)
470  * {
471  *
472  * }
473  *
474  * @endcode
475  *
476  *
477  * <a name="BicepsBrachiiGeometry"></a>
478  * <h4>BicepsBrachiiGeometry</h4>
479  *
480 
481  *
482  * Make adjustments to the geometry and discretisation of the
483  * biceps model.
484  *
485 
486  *
487  *
488  * @code
489  * struct BicepsBrachii
490  * {
491  * double axial_length;
492  * double radius_insertion_origin;
493  * double radius_midpoint;
494  * double scale;
495  * unsigned int elements_along_axis;
496  * unsigned int n_refinements_radial;
497  * bool include_gravity;
498  * double axial_force;
499  *
500  * const types::boundary_id bid_BB_dirichlet_X = 1;
501  * const types::boundary_id bid_BB_neumann = 10;
502  * const types::boundary_id mid_BB_radial = 100;
503  *
504  * static void
505  * declare_parameters(ParameterHandler &prm);
506  *
507  * void
508  * parse_parameters(ParameterHandler &prm);
509  * };
510  *
511  * void BicepsBrachii::declare_parameters(ParameterHandler &prm)
512  * {
513  * prm.enter_subsection("Biceps Brachii geometry");
514  * {
515  * prm.declare_entry("Axial length", "250",
516  * Patterns::Double(0),
517  * "Axial length of the muscle");
518  *
519  * prm.declare_entry("Radius insertion and origin", "5",
520  * Patterns::Double(0),
521  * "Insertion and origin radius");
522  *
523  * prm.declare_entry("Radius midpoint", "7.5",
524  * Patterns::Double(0),
525  * "Radius at the midpoint of the muscle");
526  *
527  * prm.declare_entry("Grid scale", "1e-3",
528  * Patterns::Double(0.0),
529  * "Global grid scaling factor");
530  *
531  * prm.declare_entry("Elements along axis", "32",
532  * Patterns::Integer(2),
533  * "Number of elements along the muscle axis");
534  *
535  * prm.declare_entry("Radial refinements", "4",
536  * Patterns::Integer(0),
537  * "Control the discretisation in the radial direction");
538  *
539  * prm.declare_entry("Gravity", "false",
540  * Patterns::Bool(),
541  * "Include the effects of gravity (in the y-direction; "
542  * " perpendicular to the muscle axis)");
543  *
544  * prm.declare_entry("Axial force", "1",
545  * Patterns::Double(),
546  * "Applied distributed axial force (in Newtons)");
547  * }
548  * prm.leave_subsection();
549  * }
550  *
551  * void BicepsBrachii::parse_parameters(ParameterHandler &prm)
552  * {
553  * prm.enter_subsection("Biceps Brachii geometry");
554  * {
555  * axial_length = prm.get_double("Axial length");
556  * radius_insertion_origin = prm.get_double("Radius insertion and origin");
557  * radius_midpoint = prm.get_double("Radius midpoint");
558  * scale = prm.get_double("Grid scale");
559  * elements_along_axis = prm.get_integer("Elements along axis");
560  * n_refinements_radial = prm.get_integer("Radial refinements");
561  * include_gravity = prm.get_bool("Gravity");
562  * axial_force = prm.get_double("Axial force");
563  * }
564  * prm.leave_subsection();
565  *
566  * AssertThrow(radius_midpoint >= radius_insertion_origin,
567  * ExcMessage("Unrealistic geometry"));
568  * }
569  *
570  * @endcode
571  *
572  *
573  * <a name="Neurologicalsignal"></a>
574  * <h4>Neurological signal</h4>
575  *
576 
577  *
578  *
579  * @code
580  * struct NeurologicalSignal
581  * {
582  * double neural_signal_start_time;
583  * double neural_signal_end_time;
584  *
585  * static void
586  * declare_parameters(ParameterHandler &prm);
587  *
588  * void
589  * parse_parameters(ParameterHandler &prm);
590  * };
591  *
592  * void NeurologicalSignal::declare_parameters(ParameterHandler &prm)
593  * {
594  * prm.enter_subsection("Neurological signal");
595  * {
596  * prm.declare_entry("Start time", "1.0",
597  * Patterns::Double(0),
598  * "Time at which to start muscle activation");
599  *
600  * prm.declare_entry("End time", "2.0",
601  * Patterns::Double(0),
602  * "Time at which to remove muscle activation signal");
603  * }
604  * prm.leave_subsection();
605  * }
606  *
607  * void NeurologicalSignal::parse_parameters(ParameterHandler &prm)
608  * {
609  * prm.enter_subsection("Neurological signal");
610  * {
611  * neural_signal_start_time = prm.get_double("Start time");
612  * neural_signal_end_time = prm.get_double("End time");
613  * }
614  * prm.leave_subsection();
615  *
616  * Assert(neural_signal_start_time < neural_signal_end_time,
617  * ExcMessage("Invalid neural signal times."));
618  * }
619  *
620  * @endcode
621  *
622  *
623  * <a name="Time"></a>
624  * <h4>Time</h4>
625  *
626 
627  *
628  * Set the timestep size @f$ \varDelta t @f$ and the simulation end-time.
629  *
630  * @code
631  * struct Time
632  * {
633  * double delta_t;
634  * double end_time;
635  * double end_ramp_time;
636  *
637  * static void
638  * declare_parameters(ParameterHandler &prm);
639  *
640  * void
641  * parse_parameters(ParameterHandler &prm);
642  * };
643  *
644  * void Time::declare_parameters(ParameterHandler &prm)
645  * {
646  * prm.enter_subsection("Time");
647  * {
648  * prm.declare_entry("End time", "3",
649  * Patterns::Double(0),
650  * "End time");
651  *
652  * prm.declare_entry("End ramp time", "1",
653  * Patterns::Double(0),
654  * "Force ramp end time");
655  *
656  * prm.declare_entry("Time step size", "0.1",
657  * Patterns::Double(0),
658  * "Time step size");
659  * }
660  * prm.leave_subsection();
661  * }
662  *
663  * void Time::parse_parameters(ParameterHandler &prm)
664  * {
665  * prm.enter_subsection("Time");
666  * {
667  * end_time = prm.get_double("End time");
668  * end_ramp_time = prm.get_double("End ramp time");
669  * delta_t = prm.get_double("Time step size");
670  * }
671  * prm.leave_subsection();
672  * }
673  *
674  * @endcode
675  *
676  *
677  * <a name="Allparameters"></a>
678  * <h4>All parameters</h4>
679  *
680 
681  *
682  * Finally we consolidate all of the above structures into a single container
683  * that holds all of our run-time selections.
684  *
685  * @code
686  * struct AllParameters : public FESystem,
687  * public Problem,
688  * public IsotonicContraction,
689  * public BicepsBrachii,
690  * public NeurologicalSignal,
691  * public Time
692  * {
693  * AllParameters(const std::string &input_file);
694  *
695  * static void
696  * declare_parameters(ParameterHandler &prm);
697  *
698  * void
699  * parse_parameters(ParameterHandler &prm);
700  * };
701  *
702  * AllParameters::AllParameters(const std::string &input_file)
703  * {
704  * ParameterHandler prm;
705  * declare_parameters(prm);
706  * prm.parse_input(input_file);
707  * parse_parameters(prm);
708  * }
709  *
710  * void AllParameters::declare_parameters(ParameterHandler &prm)
711  * {
712  * FESystem::declare_parameters(prm);
713  * Problem::declare_parameters(prm);
714  * IsotonicContraction::declare_parameters(prm);
715  * BicepsBrachii::declare_parameters(prm);
716  * NeurologicalSignal::declare_parameters(prm);
717  * Time::declare_parameters(prm);
718  * }
719  *
720  * void AllParameters::parse_parameters(ParameterHandler &prm)
721  * {
722  * FESystem::parse_parameters(prm);
723  * Problem::parse_parameters(prm);
724  * IsotonicContraction::parse_parameters(prm);
725  * BicepsBrachii::parse_parameters(prm);
726  * NeurologicalSignal::parse_parameters(prm);
727  * Time::parse_parameters(prm);
728  *
729  * @endcode
730  *
731  * Override time setting for test defined
732  * in the literature
733  *
734  * @code
735  * if (problem == "IsotonicContraction")
736  * {
737  * end_time = 3.0;
738  * end_ramp_time = 1.0;
739  * delta_t = 0.1;
740  *
741  * neural_signal_start_time = 1.0;
742  * neural_signal_end_time = 2.0;
743  * }
744  * }
745  * }
746  *
747  * @endcode
748  *
749  *
750  * <a name="Bodyforcevalues"></a>
751  * <h3>Body force values</h3>
752  *
753 
754  *
755  *
756  * @code
757  * template <int dim>
758  * class BodyForce : public Function<dim>
759  * {
760  * public:
761  * BodyForce (const double rho,
762  * const Tensor<1,dim> direction);
763  * virtual ~BodyForce () {}
764  *
765  * virtual void vector_value (const Point<dim> &p,
766  * Vector<double> &values) const;
767  *
768  * virtual void vector_value_list (const std::vector<Point<dim> > &points,
769  * std::vector<Vector<double> > &value_list) const;
770  *
771  * const double rho;
772  * const double g;
773  * const Tensor<1,dim> M;
774  * };
775  *
776  *
777  * template <int dim>
778  * BodyForce<dim>::BodyForce (const double rho,
779  * const Tensor<1,dim> direction)
780  * :
781  * Function<dim> (dim),
782  * rho (rho),
783  * g (9.81),
784  * M (direction)
785  * {
786  * Assert(M.norm() == 1.0, ExcMessage("Direction vector is not a unit vector"));
787  * }
788  *
789  *
790  * template <int dim>
791  * inline
792  * void BodyForce<dim>::vector_value (const Point<dim> &/*p*/,
793  * Vector<double> &values) const
794  * {
795  * Assert (values.size() == dim,
796  * ExcDimensionMismatch (values.size(), dim));
797  * Assert (dim >= 2, ExcNotImplemented());
798  * for (unsigned int d=0; d<dim; ++d)
799  * {
800  * values(d) = rho*g*M[d];
801  * }
802  * }
803  *
804  *
805  * template <int dim>
806  * void BodyForce<dim>::vector_value_list (const std::vector<Point<dim> > &points,
807  * std::vector<Vector<double> > &value_list) const
808  * {
809  * Assert (value_list.size() == points.size(),
810  * ExcDimensionMismatch (value_list.size(), points.size()));
811  *
812  * const unsigned int n_points = points.size();
813  *
814  * for (unsigned int p=0; p<n_points; ++p)
815  * BodyForce<dim>::vector_value (points[p],
816  * value_list[p]);
817  * }
818  *
819  * template <int dim>
820  * class Traction : public Function<dim>
821  * {
822  * public:
823  * Traction (const double force,
824  * const double area);
825  * virtual ~Traction () {}
826  *
827  * virtual void vector_value (const Point<dim> &p,
828  * Vector<double> &values) const;
829  *
830  * virtual void vector_value_list (const std::vector<Point<dim> > &points,
831  * std::vector<Vector<double> > &value_list) const;
832  *
833  * const double t;
834  * };
835  *
836  *
837  * template <int dim>
838  * Traction<dim>::Traction (const double force,
839  * const double area)
840  * :
841  * Function<dim> (dim),
842  * t (force/area)
843  * {}
844  *
845  *
846  * template <int dim>
847  * inline
848  * void Traction<dim>::vector_value (const Point<dim> &/*p*/,
849  * Vector<double> &values) const
850  * {
851  * Assert (values.size() == dim,
852  * ExcDimensionMismatch (values.size(), dim));
853  * Assert (dim == 3, ExcNotImplemented());
854  *
855  * @endcode
856  *
857  * Assume uniform distributed load
858  *
859  * @code
860  * values(0) = t;
861  * values(1) = 0.0;
862  * values(2) = 0.0;
863  * }
864  *
865  *
866  * template <int dim>
867  * void Traction<dim>::vector_value_list (const std::vector<Point<dim> > &points,
868  * std::vector<Vector<double> > &value_list) const
869  * {
870  * Assert (value_list.size() == points.size(),
871  * ExcDimensionMismatch (value_list.size(), points.size()));
872  *
873  * const unsigned int n_points = points.size();
874  *
875  * for (unsigned int p=0; p<n_points; ++p)
876  * Traction<dim>::vector_value (points[p],
877  * value_list[p]);
878  * }
879  *
880  * @endcode
881  *
882  *
883  * <a name="Utilityfunctions"></a>
884  * <h3>Utility functions</h3>
885  *
886 
887  *
888  *
889  * @code
890  * template <int dim>
891  * inline
892  * Tensor<2,dim> get_deformation_gradient (std::vector<Tensor<1,dim> > &grad)
893  * {
894  * Assert (grad.size() == dim, ExcInternalError());
895  *
896  * Tensor<2,dim> F (unit_symmetric_tensor<dim>());
897  * for (unsigned int i=0; i<dim; ++i)
898  * for (unsigned int j=0; j<dim; ++j)
899  * F[i][j] += grad[i][j];
900  * return F;
901  * }
902  *
903  * template <int dim>
904  * inline
905  * SymmetricTensor<2,dim> get_small_strain (std::vector<Tensor<1,dim> > &grad)
906  * {
907  * Assert (grad.size() == dim, ExcInternalError());
908  *
909  * SymmetricTensor<2,dim> strain;
910  * for (unsigned int i=0; i<dim; ++i)
911  * strain[i][i] = grad[i][i];
912  *
913  * for (unsigned int i=0; i<dim; ++i)
914  * for (unsigned int j=i+1; j<dim; ++j)
915  * strain[i][j] = (grad[i][j] + grad[j][i]) / 2;
916  * return strain;
917  * }
918  *
919  * @endcode
920  *
921  *
922  * <a name="Propertiesformusclematrix"></a>
923  * <h3>Properties for muscle matrix</h3>
924  *
925 
926  *
927  *
928  * @code
929  * struct MuscleMatrix
930  * {
931  * static const double E; // Young's modulus
932  * static const double nu; // Poisson ratio
933  *
934  * static const double mu; // Shear modulus
935  * static const double lambda; // Lame parameter
936  * };
937  *
938  * const double MuscleMatrix::E = 26e3;
939  * const double MuscleMatrix::nu = 0.45;
940  * const double MuscleMatrix::mu = MuscleMatrix::E/(2.0*(1.0 + MuscleMatrix::nu));
941  * const double MuscleMatrix::lambda = 2.0*MuscleMatrix::mu *MuscleMatrix::nu/(1.0 - 2.0*MuscleMatrix::nu);
942  *
943  * @endcode
944  *
945  *
946  * <a name="Localdataformusclefibres"></a>
947  * <h3>Local data for muscle fibres</h3>
948  *
949 
950  *
951  *
952  * @code
953  * #define convert_gf_to_N 1.0/101.97
954  * #define convert_gf_per_cm2_to_N_per_m2 convert_gf_to_N*1e2*1e2
955  * #define T0 6280.0*convert_gf_per_cm2_to_N_per_m2
956  *
957  * @endcode
958  *
959  * A struct that governs the functioning of a single muscle fibre
960  *
961  * @code
962  * template <int dim>
963  * struct MuscleFibre
964  * {
965  * MuscleFibre (void)
966  * : alpha (0.0),
967  * alpha_t1 (0.0),
968  * epsilon_f (0.0),
969  * epsilon_c (0.0),
970  * epsilon_c_t1 (0.0),
971  * epsilon_c_dot (0.0)
972  * {
973  *
974  * }
975  *
976  * MuscleFibre(const Tensor<1,dim> &direction)
977  * : M (direction),
978  * alpha (0.0),
979  * alpha_t1 (0.0),
980  * epsilon_f (0.0),
981  * epsilon_c (0.0),
982  * epsilon_c_t1 (0.0),
983  * epsilon_c_dot (0.0)
984  * {
985  * Assert(M.norm() == 1.0,
986  * ExcMessage("Fibre direction is not a unit vector"));
987  * }
988  *
989  * void update_alpha (const double u,
990  * const double dt);
991  *
992  * void update_state(const SymmetricTensor<2,dim> &strain_tensor,
993  * const double dt);
994  *
995  * const Tensor<1,dim> &get_M () const
996  * {
997  * return M;
998  * }
999  * double get_m_p () const;
1000  * double get_m_s () const;
1001  * double get_beta (const double dt) const;
1002  * double get_gamma (const double dt) const;
1003  *
1004  * @endcode
1005  *
1006  * Postprocessing
1007  *
1008  * @code
1009  * const double &get_alpha() const
1010  * {
1011  * return alpha;
1012  * }
1013  * const double &get_epsilon_f() const
1014  * {
1015  * return epsilon_f;
1016  * }
1017  * const double &get_epsilon_c() const
1018  * {
1019  * return epsilon_c;
1020  * }
1021  * const double &get_epsilon_c_dot() const
1022  * {
1023  * return epsilon_c_dot;
1024  * }
1025  *
1026  * private:
1027  * Tensor<1,dim> M; // Direction
1028  *
1029  * double alpha; // Activation level at current timestep
1030  * double alpha_t1; // Activation level at previous timestep
1031  *
1032  * double epsilon_f; // Fibre strain at current timestep
1033  * double epsilon_c; // Contractile strain at current timestep
1034  * double epsilon_c_t1; // Contractile strain at previous timestep
1035  * double epsilon_c_dot; // Contractile velocity at previous timestep
1036  *
1037  * double get_f_c_L () const;
1038  * double get_m_c_V () const;
1039  * double get_c_c_V () const;
1040  * };
1041  *
1042  * template <int dim>
1043  * void MuscleFibre<dim>::update_alpha (const double u,
1044  * const double dt)
1045  * {
1046  * static const double tau_r = 0.15; // s
1047  * static const double tau_f = 0.15; // s
1048  * static const double alpha_min = 0;
1049  *
1050  * if (u == 1.0)
1051  * alpha = (alpha_t1*tau_r*tau_f + dt*tau_f) / (tau_r*tau_f + dt*tau_f);
1052  * else if (u == 0)
1053  * alpha = (alpha_t1*tau_r*tau_f + dt*alpha_min*tau_r) / (tau_r*tau_f + dt*tau_r);
1054  * else
1055  * {
1056  * const double b = 1.0/tau_r - 1.0/tau_f;
1057  * const double c = 1.0/tau_f;
1058  * const double d = alpha_min/tau_f;
1059  * const double f1 = 1.0/tau_r - alpha_min/tau_f;
1060  * const double p = b*u + c;
1061  * const double q = f1*u + d;
1062  *
1063  * alpha = (q*dt + alpha_t1)/(1.0 + p*dt);
1064  * }
1065  * }
1066  *
1067  *
1068  * template <int dim>
1069  * double MuscleFibre<dim>::get_m_p () const
1070  * {
1071  * static const double A = 8.568e-4*convert_gf_per_cm2_to_N_per_m2;
1072  * static const double a = 12.43;
1073  * if (epsilon_f >= 0.0)
1074  * {
1075  * @endcode
1076  *
1077  * 100 times more compliant than Martins2006
1078  *
1079  * @code
1080  * static const double m_p = 2.0*A*a/1e2;
1081  * return m_p;
1082  * }
1083  * else
1084  * return 0.0;
1085  * }
1086  *
1087  * template <int dim>
1088  * double MuscleFibre<dim>::get_m_s (void) const
1089  * {
1090  * const double epsilon_s = epsilon_f - epsilon_c; // Small strain assumption
1091  * if (epsilon_s >= -1e-6) // Tolerant check
1092  * return 10.0;
1093  * else
1094  * return 0.0;
1095  * }
1096  *
1097  * template <int dim>
1098  * double MuscleFibre<dim>::get_f_c_L (void) const
1099  * {
1100  * if (epsilon_c <= 0.5 && epsilon_c >= -0.5)
1101  * return 1.0;
1102  * else
1103  * return 0.0;
1104  * }
1105  *
1106  * template <int dim>
1107  * double MuscleFibre<dim>::get_m_c_V (void) const
1108  * {
1109  * if (epsilon_c_dot < -5.0)
1110  * return 0.0;
1111  * else if (epsilon_c_dot <= 3.0)
1112  * return 1.0/5.0;
1113  * else
1114  * return 0.0;
1115  * }
1116  *
1117  * template <int dim>
1118  * double MuscleFibre<dim>::get_c_c_V (void) const
1119  * {
1120  * if (epsilon_c_dot < -5.0)
1121  * return 0.0;
1122  * else if (epsilon_c_dot <= 3.0)
1123  * return 1.0;
1124  * else
1125  * return 1.6;
1126  * }
1127  *
1128  * template <int dim>
1129  * double MuscleFibre<dim>::get_beta(const double dt) const
1130  * {
1131  * return get_f_c_L()*get_m_c_V()*alpha/dt + get_m_s();
1132  * }
1133  *
1134  * template <int dim>
1135  * double MuscleFibre<dim>::get_gamma(const double dt) const
1136  * {
1137  * return get_f_c_L()*alpha*(get_m_c_V()*epsilon_c_t1/dt - get_c_c_V());
1138  * }
1139  *
1140  * template <int dim>
1141  * void MuscleFibre<dim>::update_state(const SymmetricTensor<2,dim> &strain_tensor,
1142  * const double dt)
1143  * {
1144  * @endcode
1145  *
1146  * Values from previous state
1147  * These were the values that were used in the assembly,
1148  * so we must use them in the update step to be consistant.
1149  * Need to compute these before we overwrite epsilon_c_t1
1150  *
1151  * @code
1152  * const double m_s = get_m_s();
1153  * const double beta = get_beta(dt);
1154  * const double gamma = get_gamma(dt);
1155  *
1156  * @endcode
1157  *
1158  * Update current state
1159  *
1160  * @code
1161  * alpha_t1 = alpha;
1162  * epsilon_f = M*static_cast< Tensor<2,dim> >(strain_tensor)*M;
1163  * epsilon_c_t1 = epsilon_c;
1164  * epsilon_c = (m_s*epsilon_f + gamma)/beta;
1165  * epsilon_c_dot = (epsilon_c - epsilon_c_t1)/dt;
1166  * }
1167  *
1168  *
1169  * @endcode
1170  *
1171  *
1172  * <a name="ThecodeLinearMuscleModelProblemcodeclasstemplate"></a>
1173  * <h3>The <code>LinearMuscleModelProblem</code> class template</h3>
1174  *
1175 
1176  *
1177  *
1178  * @code
1179  * template <int dim>
1180  * class LinearMuscleModelProblem
1181  * {
1182  * public:
1183  * LinearMuscleModelProblem (const std::string &input_file);
1184  * ~LinearMuscleModelProblem ();
1185  * void run ();
1186  *
1187  * private:
1188  * void make_grid ();
1189  * void setup_muscle_fibres ();
1190  * double get_neural_signal (const double time);
1191  * void update_fibre_activation (const double time);
1192  * void update_fibre_state ();
1193  * void setup_system ();
1194  * void assemble_system (const double time);
1195  * void apply_boundary_conditions ();
1196  * void solve ();
1197  * void output_results (const unsigned int timestep,
1198  * const double time) const;
1199  *
1200  * Parameters::AllParameters parameters;
1201  *
1203  * DoFHandler<dim> dof_handler;
1204  *
1205  * FESystem<dim> fe;
1206  * QGauss<dim> qf_cell;
1207  * QGauss<dim-1> qf_face;
1208  *
1209  * ConstraintMatrix hanging_node_constraints;
1210  *
1211  * SparsityPattern sparsity_pattern;
1212  * SparseMatrix<double> system_matrix;
1213  *
1214  * Vector<double> solution;
1215  * Vector<double> system_rhs;
1216  *
1217  * @endcode
1218  *
1219  * Time
1220  *
1221  * @code
1222  * const double t_end;
1223  * const double dt;
1224  * const double t_ramp_end; // Force ramp end time
1225  *
1226  * @endcode
1227  *
1228  * Loading
1229  *
1230  * @code
1231  * const BodyForce<dim> body_force;
1232  * const Traction<dim> traction;
1233  *
1234  * @endcode
1235  *
1236  * Local data
1237  *
1238  * @code
1239  * std::vector< std::vector<MuscleFibre<dim> > > fibre_data;
1240  *
1241  * @endcode
1242  *
1243  * Constitutive functions for assembly
1244  *
1245  * @code
1246  * SymmetricTensor<4,dim> get_stiffness_tensor (const unsigned int cell,
1247  * const unsigned int q_point_cell) const;
1248  * SymmetricTensor<2,dim> get_rhs_tensor (const unsigned int cell,
1249  * const unsigned int q_point_cell) const;
1250  * };
1251  *
1252  * @endcode
1253  *
1254  *
1255  * <a name="LinearMuscleModelProblemLinearMuscleModelProblem"></a>
1256  * <h4>LinearMuscleModelProblem::LinearMuscleModelProblem</h4>
1257  *
1258 
1259  *
1260  *
1261  * @code
1262  * template <int dim>
1263  * LinearMuscleModelProblem<dim>::LinearMuscleModelProblem (const std::string &input_file)
1264  * :
1265  * parameters(input_file),
1266  * dof_handler (triangulation),
1267  * fe (FE_Q<dim>(parameters.poly_degree), dim),
1268  * qf_cell (parameters.quad_order),
1269  * qf_face (parameters.quad_order),
1270  * t_end (parameters.end_time),
1271  * dt (parameters.delta_t),
1272  * t_ramp_end(parameters.end_ramp_time),
1273  * body_force ((parameters.problem == "BicepsBrachii" &&parameters.include_gravity == true) ?
1274  * BodyForce<dim>(0.375*1000.0, Tensor<1,dim>({0,-1,0})) : // (reduced) Density and direction
1275  * BodyForce<dim>(0.0, Tensor<1,dim>({0,0,1})) ),
1276  * traction (parameters.problem == "BicepsBrachii" ?
1277  * Traction<dim>(parameters.axial_force, // Force, area
1278  * M_PI*std::pow(parameters.radius_insertion_origin *parameters.scale,2.0) ) :
1279  * Traction<dim>(4.9*convert_gf_to_N, // Force; Conversion of gf to N,
1280  * (2.0*parameters.half_length_y)*(2.0*parameters.half_length_z)) ) // Area
1281  * {
1282  * Assert(dim==3, ExcNotImplemented());
1283  * }
1284  *
1285  *
1286  * @endcode
1287  *
1288  *
1289  * <a name="LinearMuscleModelProblemLinearMuscleModelProblem"></a>
1290  * <h4>LinearMuscleModelProblem::~LinearMuscleModelProblem</h4>
1291  *
1292 
1293  *
1294  *
1295  * @code
1296  * template <int dim>
1297  * LinearMuscleModelProblem<dim>::~LinearMuscleModelProblem ()
1298  * {
1299  * dof_handler.clear ();
1300  * }
1301  *
1302  *
1303  * @endcode
1304  *
1305  *
1306  * <a name="LinearMuscleModelProblemmake_grid"></a>
1307  * <h4>LinearMuscleModelProblem::make_grid</h4>
1308  *
1309 
1310  *
1311  *
1312  * @code
1313  * template<int dim>
1314  * struct BicepsGeometry
1315  * {
1316  * BicepsGeometry(const double axial_length,
1317  * const double radius_ins_orig,
1318  * const double radius_mid)
1319  * :
1320  * ax_lgth (axial_length),
1321  * r_ins_orig (radius_ins_orig),
1322  * r_mid (radius_mid)
1323  * {}
1324  *
1325  * @endcode
1326  *
1327  * The radial profile of the muscle
1328  * This provides the new coordinates for points @p pt
1329  * on a cylinder of radius r_ins_orig and length
1330  * ax_lgth to be moved to in order to create the
1331  * physiologically representative geometry of
1332  * the muscle
1333  *
1334  * @code
1335  * Point<dim> profile (const Point<dim> &pt_0) const
1336  * {
1337  * Assert(pt_0[0] > -1e-6,
1338  * ExcMessage("All points must have x-coordinate > 0"));
1339  *
1340  * const double r_scale = get_radial_scaling_factor(pt_0[0]);
1341  * return pt_0 + Point<dim>(0.0, r_scale*pt_0[1], r_scale*pt_0[2]);
1342  * }
1343  *
1344  * Point<dim> operator() (const Point<dim> &pt) const
1345  * {
1346  * return profile(pt);
1347  * }
1348  *
1349  * @endcode
1350  *
1351  * Provides the muscle direction at the point @p pt
1352  * in the real geometry (one that has undergone the
1353  * transformation given by the profile() function)
1354  * and subequent grid rescaling.
1355  * The directions are given by the gradient of the
1356  * transformation function (i.e. the fibres are
1357  * orientated by the curvature of the muscle).
1358  *
1359 
1360  *
1361  * So, being lazy, we transform the current point back
1362  * to the original point on the completely unscaled
1363  * cylindrical grid. We then evaluate the transformation
1364  * at two points (axially displaced) very close to the
1365  * point of interest. The normalised vector joining the
1366  * transformed counterparts of the perturbed points is
1367  * the gradient of the transformation function and,
1368  * thus, defines the fibre direction.
1369  *
1370  * @code
1371  * Tensor<1,dim> direction (const Point<dim> &pt_scaled,
1372  * const double &grid_scale) const
1373  * {
1374  * const Point<dim> pt = (1.0/grid_scale)*pt_scaled;
1375  * const Point<dim> pt_0 = inv_profile(pt);
1376  *
1377  * static const double eps = 1e-6;
1378  * const Point<dim> pt_0_eps_p = pt_0 + Point<dim>(+eps,0,0);
1379  * const Point<dim> pt_0_eps_m = pt_0 + Point<dim>(-eps,0,0);
1380  * const Point<dim> pt_eps_p = profile(pt_0_eps_p);
1381  * const Point<dim> pt_eps_m = profile(pt_0_eps_m);
1382  *
1383  * static const double tol = 1e-9;
1384  * Assert(profile(pt_0).distance(pt) < tol, ExcInternalError());
1385  * Assert(inv_profile(pt_eps_p).distance(pt_0_eps_p) < tol, ExcInternalError());
1386  * Assert(inv_profile(pt_eps_m).distance(pt_0_eps_m) < tol, ExcInternalError());
1387  *
1388  * Tensor<1,dim> dir = pt_eps_p-pt_eps_m;
1389  * dir /= dir.norm();
1390  * return dir;
1391  * }
1392  *
1393  * private:
1394  * const double ax_lgth;
1395  * const double r_ins_orig;
1396  * const double r_mid;
1397  *
1398  * double get_radial_scaling_factor (const double &x) const
1399  * {
1400  * @endcode
1401  *
1402  * Expect all grid points with X>=0, but we provide a
1403  * tolerant location for points "on" the Cartesian plane X=0
1404  *
1405  * @code
1406  * const double lgth_frac = std::max(x/ax_lgth,0.0);
1407  * const double amplitude = 0.25*(r_mid - r_ins_orig);
1408  * const double phase_shift = M_PI;
1409  * const double y_shift = 1.0;
1410  * const double wave_func = y_shift + std::cos(phase_shift + 2.0*M_PI*lgth_frac);
1411  * Assert(wave_func >= 0.0, ExcInternalError());
1412  * return std::sqrt(amplitude*wave_func);
1413  * }
1414  *
1415  * Point<dim> inv_profile (const Point<dim> &pt) const
1416  * {
1417  * Assert(pt[0] > -1e-6,
1418  * ExcMessage("All points must have x-coordinate > 0"));
1419  *
1420  * const double r_scale = get_radial_scaling_factor(pt[0]);
1421  * const double trans_inv_scale = 1.0/(1.0+r_scale);
1422  * return Point<dim>(pt[0], trans_inv_scale*pt[1], trans_inv_scale*pt[2]);
1423  * }
1424  * };
1425  *
1426  * template <int dim>
1427  * void LinearMuscleModelProblem<dim>::make_grid ()
1428  * {
1429  * Assert (dim == 3, ExcNotImplemented());
1430  *
1431  * if (parameters.problem == "IsotonicContraction")
1432  * {
1433  * const Point<dim> p1(-parameters.half_length_x,
1434  * -parameters.half_length_y,
1435  * -parameters.half_length_z);
1436  * const Point<dim> p2( parameters.half_length_x,
1437  * parameters.half_length_y,
1438  * parameters.half_length_z);
1439  *
1441  *
1443  * triangulation.begin_active(), endc = triangulation.end();
1444  * for (; cell != endc; ++cell)
1445  * {
1446  * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
1447  * {
1448  * if (cell->face(face)->at_boundary() == true)
1449  * {
1450  * if (cell->face(face)->center()[0] == -parameters.half_length_x) // -X oriented face
1451  * cell->face(face)->set_boundary_id(parameters.bid_CC_dirichlet_symm_X); // Dirichlet
1452  * else if (cell->face(face)->center()[0] == parameters.half_length_x) // +X oriented face
1453  * cell->face(face)->set_boundary_id(parameters.bid_CC_neumann); // Neumann
1454  * else if (std::abs(cell->face(face)->center()[2]) == parameters.half_length_z) // -Z/+Z oriented face
1455  * cell->face(face)->set_boundary_id(parameters.bid_CC_dirichlet_symm_Z); // Dirichlet
1456  * }
1457  * }
1458  * }
1459  *
1460  * triangulation.refine_global (1);
1461  * }
1462  * else if (parameters.problem == "BicepsBrachii")
1463  * {
1464  * SphericalManifold<2> manifold_cap;
1465  * Triangulation<2> tria_cap;
1466  * GridGenerator::hyper_ball(tria_cap,
1467  * Point<2>(),
1468  * parameters.radius_insertion_origin);
1470  * cell = tria_cap.begin_active();
1471  * cell != tria_cap.end(); ++cell)
1472  * {
1473  * for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell; ++face)
1474  * {
1475  * if (cell->face(face)->at_boundary() == true)
1476  * cell->face(face)->set_all_manifold_ids(0);
1477  * }
1478  * }
1479  * tria_cap.set_manifold (0, manifold_cap);
1480  * tria_cap.refine_global(parameters.n_refinements_radial);
1481  *
1482  * Triangulation<2> tria_cap_flat;
1483  * GridGenerator::flatten_triangulation(tria_cap, tria_cap_flat);
1484  *
1485  * GridGenerator::extrude_triangulation(tria_cap_flat,
1486  * parameters.elements_along_axis,
1487  * parameters.axial_length,
1488  * triangulation);
1489  *
1490  * struct GridRotate
1491  * {
1492  * Point<dim> operator() (const Point<dim> &in) const
1493  * {
1494  * static const Tensor<2,dim> rot_mat = Physics::Transformations::Rotations::rotation_matrix_3d(Point<dim>(0,1,0), M_PI/2.0);
1495  * return Point<dim>(rot_mat*in);
1496  * }
1497  * };
1498  *
1499  * @endcode
1500  *
1501  * Rotate grid so that the length is axially
1502  * coincident and aligned with the X-axis
1503  *
1504  * @code
1505  * GridTools::transform (GridRotate(), triangulation);
1506  *
1507  * @endcode
1508  *
1509  * Deform the grid into something that vaguely
1510  * resemble's a Biceps Brachii
1511  *
1512  * @code
1513  * GridTools::transform (BicepsGeometry<dim>(parameters.axial_length,
1514  * parameters.radius_insertion_origin,
1515  * parameters.radius_midpoint), triangulation);
1516  *
1517  * @endcode
1518  *
1519  * Set boundary IDs
1520  *
1521  * @code
1522  * typename Triangulation<dim>::active_cell_iterator cell =
1523  * triangulation.begin_active(), endc = triangulation.end();
1524  * for (; cell != endc; ++cell)
1525  * {
1526  * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
1527  * {
1528  * if (cell->face(face)->at_boundary() == true)
1529  * {
1530  * static const double tol =1e-6;
1531  * if (std::abs(cell->face(face)->center()[0]) < tol) // -X oriented face
1532  * cell->face(face)->set_boundary_id(parameters.bid_BB_dirichlet_X); // Dirichlet
1533  * else if (std::abs(cell->face(face)->center()[0] - parameters.axial_length) < tol) // +X oriented face
1534  * cell->face(face)->set_boundary_id(parameters.bid_BB_neumann); // Neumann
1535  * }
1536  * }
1537  * }
1538  *
1539  * @endcode
1540  *
1541  * Finally resize the grid
1542  *
1543  * @code
1544  * GridTools::scale (parameters.scale, triangulation);
1545  * }
1546  * else
1547  * AssertThrow(false, ExcNotImplemented());
1548  * }
1549  *
1550  * @endcode
1551  *
1552  *
1553  * <a name="LinearMuscleModelProblemsetup_muscle_fibres"></a>
1554  * <h4>LinearMuscleModelProblem::setup_muscle_fibres</h4>
1555  *
1556 
1557  *
1558  *
1559  * @code
1560  * template <int dim>
1561  * void LinearMuscleModelProblem<dim>::setup_muscle_fibres ()
1562  * {
1563  * fibre_data.clear();
1564  * const unsigned int n_cells = triangulation.n_active_cells();
1565  * fibre_data.resize(n_cells);
1566  * const unsigned int n_q_points_cell = qf_cell.size();
1567  *
1568  * if (parameters.problem == "IsotonicContraction")
1569  * {
1570  * MuscleFibre<dim> fibre_template (Tensor<1,dim>({1,0,0}));
1571  *
1572  * for (unsigned int cell_no=0; cell_no<triangulation.n_active_cells(); ++cell_no)
1573  * {
1574  * fibre_data[cell_no].resize(n_q_points_cell);
1575  * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1576  * {
1577  * fibre_data[cell_no][q_point_cell] = fibre_template;
1578  * }
1579  * }
1580  * }
1581  * else if (parameters.problem == "BicepsBrachii")
1582  * {
1583  * FEValues<dim> fe_values (fe, qf_cell, update_quadrature_points);
1584  * BicepsGeometry<dim> bicep_geom (parameters.axial_length,
1585  * parameters.radius_insertion_origin,
1586  * parameters.radius_midpoint);
1587  *
1588  * unsigned int cell_no = 0;
1589  * for (typename Triangulation<dim>::active_cell_iterator
1590  * cell = triangulation.begin_active();
1591  * cell != triangulation.end();
1592  * ++cell, ++cell_no)
1593  * {
1594  * Assert(cell_no<fibre_data.size(), ExcMessage("Trying to access fibre data not stored for this cell index"));
1595  * fe_values.reinit(cell);
1596  *
1597  * fibre_data[cell_no].resize(n_q_points_cell);
1598  * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1599  * {
1600  * const Point<dim> pt = fe_values.get_quadrature_points()[q_point_cell];
1601  * fibre_data[cell_no][q_point_cell] = MuscleFibre<dim>(bicep_geom.direction(pt,parameters.scale));
1602  * }
1603  * }
1604  * }
1605  * else
1606  * AssertThrow(false, ExcNotImplemented());
1607  * }
1608  *
1609  * @endcode
1610  *
1611  *
1612  * <a name="LinearMuscleModelProblemupdate_fibre_state"></a>
1613  * <h4>LinearMuscleModelProblem::update_fibre_state</h4>
1614  *
1615 
1616  *
1617  *
1618  * @code
1619  * template <int dim>
1620  * double LinearMuscleModelProblem<dim>::get_neural_signal (const double time)
1621  * {
1622  * @endcode
1623  *
1624  * Note: 40 times less force generated than Martins2006
1625  * This is necessary due to the (compliant) linear tissue model
1626  *
1627  * @code
1628  * return (time > parameters.neural_signal_start_time && time < parameters.neural_signal_end_time ?
1629  * 1.0/40.0 :
1630  * 0.0);
1631  * }
1632  *
1633  * template <int dim>
1634  * void LinearMuscleModelProblem<dim>::update_fibre_activation (const double time)
1635  * {
1636  * const double u = get_neural_signal(time);
1637  *
1638  * const unsigned int n_q_points_cell = qf_cell.size();
1639  * for (unsigned int cell=0; cell<triangulation.n_active_cells(); ++cell)
1640  * {
1641  * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1642  * {
1643  * MuscleFibre<dim> &fibre = fibre_data[cell][q_point_cell];
1644  * fibre.update_alpha(u,dt);
1645  * }
1646  * }
1647  * }
1648  *
1649  * template <int dim>
1650  * void LinearMuscleModelProblem<dim>::update_fibre_state ()
1651  * {
1652  * const unsigned int n_q_points_cell = qf_cell.size();
1653  *
1654  * FEValues<dim> fe_values (fe, qf_cell, update_gradients);
1655  *
1656  * @endcode
1657  *
1658  * Displacement gradient
1659  *
1660  * @code
1661  * std::vector< std::vector< Tensor<1,dim> > > u_grads (n_q_points_cell,
1662  * std::vector<Tensor<1,dim> >(dim));
1663  *
1664  * unsigned int cell_no = 0;
1665  * for (typename DoFHandler<dim>::active_cell_iterator
1666  * cell = dof_handler.begin_active();
1667  * cell!=dof_handler.end(); ++cell, ++cell_no)
1668  * {
1669  * Assert(cell_no<fibre_data.size(), ExcMessage("Trying to access fibre data not stored for this cell index"));
1670  * fe_values.reinit(cell);
1671  * fe_values.get_function_gradients (solution, u_grads);
1672  *
1673  * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1674  * {
1675  * Assert(q_point_cell<fibre_data[cell_no].size(), ExcMessage("Trying to access fibre data not stored for this cell and qp index"));
1676  *
1677  * const SymmetricTensor<2,dim> strain_tensor = get_small_strain (u_grads[q_point_cell]);
1678  * MuscleFibre<dim> &fibre = fibre_data[cell_no][q_point_cell];
1679  * fibre.update_state(strain_tensor, dt);
1680  * }
1681  * }
1682  * }
1683  *
1684  * @endcode
1685  *
1686  *
1687  * <a name="LinearMuscleModelProblemsetup_system"></a>
1688  * <h4>LinearMuscleModelProblem::setup_system</h4>
1689  *
1690 
1691  *
1692  *
1693  * @code
1694  * template <int dim>
1695  * void LinearMuscleModelProblem<dim>::setup_system ()
1696  * {
1697  * dof_handler.distribute_dofs (fe);
1698  * hanging_node_constraints.clear ();
1699  * DoFTools::make_hanging_node_constraints (dof_handler,
1700  * hanging_node_constraints);
1701  * hanging_node_constraints.close ();
1702  * sparsity_pattern.reinit (dof_handler.n_dofs(),
1703  * dof_handler.n_dofs(),
1704  * dof_handler.max_couplings_between_dofs());
1705  * DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
1706  *
1707  * hanging_node_constraints.condense (sparsity_pattern);
1708  *
1709  * sparsity_pattern.compress();
1710  *
1711  * system_matrix.reinit (sparsity_pattern);
1712  *
1713  * solution.reinit (dof_handler.n_dofs());
1714  * system_rhs.reinit (dof_handler.n_dofs());
1715  *
1716  * std::cout << " Number of active cells: "
1717  * << triangulation.n_active_cells()
1718  * << std::endl;
1719  *
1720  * std::cout << " Number of degrees of freedom: "
1721  * << dof_handler.n_dofs()
1722  * << std::endl;
1723  * }
1724  *
1725  * @endcode
1726  *
1727  *
1728  * <a name="LinearMuscleModelProblemassemble_system"></a>
1729  * <h4>LinearMuscleModelProblem::assemble_system</h4>
1730  *
1731 
1732  *
1733  *
1734  * @code
1735  * template <int dim>
1736  * SymmetricTensor<4,dim>
1737  * LinearMuscleModelProblem<dim>::get_stiffness_tensor (const unsigned int cell,
1738  * const unsigned int q_point_cell) const
1739  * {
1740  * static const SymmetricTensor<2,dim> I = unit_symmetric_tensor<dim>();
1741  *
1742  * Assert(cell<fibre_data.size(), ExcMessage("Trying to access fibre data not stored for this cell index"));
1743  * Assert(q_point_cell<fibre_data[cell].size(), ExcMessage("Trying to access fibre data not stored for this cell and qp index"));
1744  * const MuscleFibre<dim> &fibre = fibre_data[cell][q_point_cell];
1745  *
1746  * @endcode
1747  *
1748  * Matrix
1749  *
1750  * @code
1751  * const double lambda = MuscleMatrix::lambda;
1752  * const double mu = MuscleMatrix::mu;
1753  * @endcode
1754  *
1755  * Fibre
1756  *
1757  * @code
1758  * const double m_p = fibre.get_m_p();
1759  * const double m_s = fibre.get_m_s();
1760  * const double beta = fibre.get_beta(dt);
1761  * AssertThrow(beta != 0.0, ExcInternalError());
1762  * const double Cf = T0*(m_p + m_s*(1.0 - m_s/beta));
1763  * const Tensor<1,dim> &M = fibre.get_M();
1764  *
1765  * SymmetricTensor<4,dim> C;
1766  * for (unsigned int i=0; i < dim; ++i)
1767  * for (unsigned int j=i; j < dim; ++j)
1768  * for (unsigned int k=0; k < dim; ++k)
1769  * for (unsigned int l=k; l < dim; ++l)
1770  * {
1771  * @endcode
1772  *
1773  * Matrix contribution
1774  *
1775  * @code
1776  * C[i][j][k][l] = lambda * I[i][j]*I[k][l]
1777  * + mu * (I[i][k]*I[j][l] + I[i][l]*I[j][k]);
1778  *
1779  * @endcode
1780  *
1781  * Fibre contribution (Passive + active branches)
1782  *
1783  * @code
1784  * C[i][j][k][l] += Cf * M[i]*M[j]*M[k]*M[l];
1785  * }
1786  *
1787  * return C;
1788  * }
1789  *
1790  * template <int dim>
1791  * SymmetricTensor<2,dim>
1792  * LinearMuscleModelProblem<dim>::get_rhs_tensor (const unsigned int cell,
1793  * const unsigned int q_point_cell) const
1794  * {
1795  * Assert(cell<fibre_data.size(), ExcMessage("Trying to access fibre data not stored for this cell index"));
1796  * Assert(q_point_cell<fibre_data[cell].size(), ExcMessage("Trying to access fibre data not stored for this cell and qp index"));
1797  * const MuscleFibre<dim> &fibre = fibre_data[cell][q_point_cell];
1798  *
1799  * const double m_s = fibre.get_m_s();
1800  * const double beta = fibre.get_beta(dt);
1801  * const double gamma = fibre.get_gamma(dt);
1802  * AssertThrow(beta != 0.0, ExcInternalError());
1803  * const double Sf = T0*(m_s*gamma/beta);
1804  * const Tensor<1,dim> &M = fibre.get_M();
1805  *
1806  * SymmetricTensor<2,dim> S;
1807  * for (unsigned int i=0; i < dim; ++i)
1808  * for (unsigned int j=i; j < dim; ++j)
1809  * {
1810  * @endcode
1811  *
1812  * Fibre contribution (Active branch)
1813  *
1814  * @code
1815  * S[i][j] = Sf * M[i]*M[j];
1816  * }
1817  *
1818  * return S;
1819  * }
1820  *
1821  * @endcode
1822  *
1823  *
1824  * <a name="LinearMuscleModelProblemassemble_system"></a>
1825  * <h4>LinearMuscleModelProblem::assemble_system</h4>
1826  *
1827 
1828  *
1829  *
1830  * @code
1831  * template <int dim>
1832  * void LinearMuscleModelProblem<dim>::assemble_system (const double time)
1833  * {
1834  * @endcode
1835  *
1836  * Reset system
1837  *
1838  * @code
1839  * system_matrix = 0;
1840  * system_rhs = 0;
1841  *
1842  * FEValues<dim> fe_values (fe, qf_cell,
1843  * update_values | update_gradients |
1844  * update_quadrature_points | update_JxW_values);
1845  * FEFaceValues<dim> fe_face_values (fe, qf_face,
1846  * update_values |
1847  * update_quadrature_points | update_JxW_values);
1848  *
1849  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1850  * const unsigned int n_q_points_cell = qf_cell.size();
1851  * const unsigned int n_q_points_face = qf_face.size();
1852  *
1853  * FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
1854  * Vector<double> cell_rhs (dofs_per_cell);
1855  *
1856  * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1857  *
1858  * @endcode
1859  *
1860  * Loading
1861  *
1862  * @code
1863  * std::vector<Vector<double> > body_force_values (n_q_points_cell,
1864  * Vector<double>(dim));
1865  * std::vector<Vector<double> > traction_values (n_q_points_face,
1866  * Vector<double>(dim));
1867  *
1868  * unsigned int cell_no = 0;
1869  * for (typename DoFHandler<dim>::active_cell_iterator
1870  * cell = dof_handler.begin_active();
1871  * cell!=dof_handler.end(); ++cell, ++cell_no)
1872  * {
1873  * cell_matrix = 0;
1874  * cell_rhs = 0;
1875  *
1876  * fe_values.reinit (cell);
1877  * body_force.vector_value_list (fe_values.get_quadrature_points(),
1878  * body_force_values);
1879  *
1880  * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1881  * {
1882  * const SymmetricTensor<4,dim> C = get_stiffness_tensor (cell_no, q_point_cell);
1883  * const SymmetricTensor<2,dim> R = get_rhs_tensor(cell_no, q_point_cell);
1884  *
1885  * for (unsigned int I=0; I<dofs_per_cell; ++I)
1886  * {
1887  * const unsigned int
1888  * component_I = fe.system_to_component_index(I).first;
1889  *
1890  * for (unsigned int J=0; J<dofs_per_cell; ++J)
1891  * {
1892  * const unsigned int
1893  * component_J = fe.system_to_component_index(J).first;
1894  *
1895  * for (unsigned int k=0; k < dim; ++k)
1896  * for (unsigned int l=0; l < dim; ++l)
1897  * cell_matrix(I,J)
1898  * += (fe_values.shape_grad(I,q_point_cell)[k] *
1899  * C[component_I][k][component_J][l] *
1900  * fe_values.shape_grad(J,q_point_cell)[l]) *
1901  * fe_values.JxW(q_point_cell);
1902  * }
1903  * }
1904  *
1905  * for (unsigned int I=0; I<dofs_per_cell; ++I)
1906  * {
1907  * const unsigned int
1908  * component_I = fe.system_to_component_index(I).first;
1909  *
1910  * cell_rhs(I)
1911  * += fe_values.shape_value(I,q_point_cell) *
1912  * body_force_values[q_point_cell](component_I) *
1913  * fe_values.JxW(q_point_cell);
1914  *
1915  * for (unsigned int k=0; k < dim; ++k)
1916  * cell_rhs(I)
1917  * += (fe_values.shape_grad(I,q_point_cell)[k] *
1918  * R[component_I][k]) *
1919  * fe_values.JxW(q_point_cell);
1920  * }
1921  * }
1922  *
1923  * for (unsigned int face = 0; face <GeometryInfo<dim>::faces_per_cell; ++face)
1924  * {
1925  * if (cell->face(face)->at_boundary() == true &&
1926  * ((parameters.problem == "IsotonicContraction" &&
1927  * cell->face(face)->boundary_id() == parameters.bid_CC_neumann) ||
1928  * (parameters.problem == "BicepsBrachii" &&
1929  * cell->face(face)->boundary_id() == parameters.bid_BB_neumann)) )
1930  * {
1931  * fe_face_values.reinit(cell, face);
1932  * traction.vector_value_list (fe_face_values.get_quadrature_points(),
1933  * traction_values);
1934  *
1935  * @endcode
1936  *
1937  * Scale applied traction according to time
1938  *
1939  * @code
1940  * const double ramp = (time <= t_ramp_end ? time/t_ramp_end : 1.0);
1941  * Assert(ramp >= 0.0 && ramp <= 1.0, ExcMessage("Invalid force ramp"));
1942  * for (unsigned int q_point_face = 0; q_point_face < n_q_points_face; ++q_point_face)
1943  * traction_values[q_point_face] *= ramp;
1944  *
1945  * for (unsigned int q_point_face = 0; q_point_face < n_q_points_face; ++q_point_face)
1946  * {
1947  * for (unsigned int I=0; I<dofs_per_cell; ++I)
1948  * {
1949  * const unsigned int
1950  * component_I = fe.system_to_component_index(I).first;
1951  *
1952  * cell_rhs(I)
1953  * += fe_face_values.shape_value(I,q_point_face)*
1954  * traction_values[q_point_face][component_I]*
1955  * fe_face_values.JxW(q_point_face);
1956  * }
1957  * }
1958  * }
1959  * }
1960  *
1961  * cell->get_dof_indices (local_dof_indices);
1962  * for (unsigned int i=0; i<dofs_per_cell; ++i)
1963  * {
1964  * for (unsigned int j=0; j<dofs_per_cell; ++j)
1965  * system_matrix.add (local_dof_indices[i],
1966  * local_dof_indices[j],
1967  * cell_matrix(i,j));
1968  *
1969  * system_rhs(local_dof_indices[i]) += cell_rhs(i);
1970  * }
1971  * }
1972  *
1973  * hanging_node_constraints.condense (system_matrix);
1974  * hanging_node_constraints.condense (system_rhs);
1975  * }
1976  *
1977  * template <int dim>
1978  * void LinearMuscleModelProblem<dim>::apply_boundary_conditions ()
1979  * {
1980  * std::map<types::global_dof_index,double> boundary_values;
1981  *
1982  * if (parameters.problem == "IsotonicContraction")
1983  * {
1984  * @endcode
1985  *
1986  * Symmetry condition on -X faces
1987  *
1988  * @code
1989  * {
1990  * ComponentMask component_mask_x (dim, false);
1991  * component_mask_x.set(0, true);
1992  * VectorTools::interpolate_boundary_values (dof_handler,
1993  * parameters.bid_CC_dirichlet_symm_X,
1994  * ZeroFunction<dim>(dim),
1995  * boundary_values,
1996  * component_mask_x);
1997  * }
1998  * @endcode
1999  *
2000  * Symmetry condition on -Z/+Z faces
2001  *
2002  * @code
2003  * {
2004  * ComponentMask component_mask_z (dim, false);
2005  * component_mask_z.set(2, true);
2006  * VectorTools::interpolate_boundary_values (dof_handler,
2007  * parameters.bid_CC_dirichlet_symm_Z,
2008  * ZeroFunction<dim>(dim),
2009  * boundary_values,
2010  * component_mask_z);
2011  * }
2012  * @endcode
2013  *
2014  * Fixed point on -X face
2015  *
2016  * @code
2017  * {
2018  * const Point<dim> fixed_point (-parameters.half_length_x,0.0,0.0);
2019  * std::vector<types::global_dof_index> fixed_dof_indices;
2020  * bool found_point_of_interest = false;
2021  *
2022  * for (typename DoFHandler<dim>::active_cell_iterator
2023  * cell = dof_handler.begin_active(),
2024  * endc = dof_handler.end(); cell != endc; ++cell)
2025  * {
2026  * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
2027  * {
2028  * @endcode
2029  *
2030  * We know that the fixed point is on the -X Dirichlet boundary
2031  *
2032  * @code
2033  * if (cell->face(face)->at_boundary() == true &&
2034  * cell->face(face)->boundary_id() == parameters.bid_CC_dirichlet_symm_X)
2035  * {
2036  * for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo<dim>::vertices_per_face; ++face_vertex_index)
2037  * {
2038  * if (cell->face(face)->vertex(face_vertex_index).distance(fixed_point) < 1e-6)
2039  * {
2040  * found_point_of_interest = true;
2041  * for (unsigned int index_component = 0; index_component < dim; ++index_component)
2042  * fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index,
2043  * index_component));
2044  * }
2045  *
2046  * if (found_point_of_interest == true) break;
2047  * }
2048  * }
2049  * if (found_point_of_interest == true) break;
2050  * }
2051  * if (found_point_of_interest == true) break;
2052  * }
2053  *
2054  * Assert(found_point_of_interest == true, ExcMessage("Didn't find point of interest"));
2055  * AssertThrow(fixed_dof_indices.size() == dim, ExcMessage("Didn't find the correct number of DoFs to fix"));
2056  *
2057  * for (unsigned int i=0; i < fixed_dof_indices.size(); ++i)
2058  * boundary_values[fixed_dof_indices[i]] = 0.0;
2059  * }
2060  * }
2061  * else if (parameters.problem == "BicepsBrachii")
2062  * {
2063  * if (parameters.include_gravity == false)
2064  * {
2065  * @endcode
2066  *
2067  * Symmetry condition on -X surface
2068  *
2069  * @code
2070  * {
2071  * ComponentMask component_mask_x (dim, false);
2072  * component_mask_x.set(0, true);
2073  * VectorTools::interpolate_boundary_values (dof_handler,
2074  * parameters.bid_BB_dirichlet_X,
2075  * ZeroFunction<dim>(dim),
2076  * boundary_values,
2077  * component_mask_x);
2078  * }
2079  *
2080  * @endcode
2081  *
2082  * Fixed central point on -X surface
2083  *
2084  * @code
2085  * {
2086  * const Point<dim> fixed_point (0.0,0.0,0.0);
2087  * std::vector<types::global_dof_index> fixed_dof_indices;
2088  * bool found_point_of_interest = false;
2089  *
2090  * for (typename DoFHandler<dim>::active_cell_iterator
2091  * cell = dof_handler.begin_active(),
2092  * endc = dof_handler.end(); cell != endc; ++cell)
2093  * {
2094  * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
2095  * {
2096  * @endcode
2097  *
2098  * We know that the fixed point is on the -X Dirichlet boundary
2099  *
2100  * @code
2101  * if (cell->face(face)->at_boundary() == true &&
2102  * cell->face(face)->boundary_id() == parameters.bid_BB_dirichlet_X)
2103  * {
2104  * for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo<dim>::vertices_per_face; ++face_vertex_index)
2105  * {
2106  * if (cell->face(face)->vertex(face_vertex_index).distance(fixed_point) < 1e-6)
2107  * {
2108  * found_point_of_interest = true;
2109  * for (unsigned int index_component = 0; index_component < dim; ++index_component)
2110  * fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index,
2111  * index_component));
2112  * }
2113  *
2114  * if (found_point_of_interest == true) break;
2115  * }
2116  * }
2117  * if (found_point_of_interest == true) break;
2118  * }
2119  * if (found_point_of_interest == true) break;
2120  * }
2121  *
2122  * Assert(found_point_of_interest == true, ExcMessage("Didn't find point of interest"));
2123  * AssertThrow(fixed_dof_indices.size() == dim, ExcMessage("Didn't find the correct number of DoFs to fix"));
2124  *
2125  * for (unsigned int i=0; i < fixed_dof_indices.size(); ++i)
2126  * boundary_values[fixed_dof_indices[i]] = 0.0;
2127  * }
2128  * }
2129  * else
2130  * {
2131  * @endcode
2132  *
2133  * When we apply gravity, some additional constraints
2134  * are required to support the load of the muscle, as
2135  * the material response is more compliant than would
2136  * be the case in reality.
2137  *
2138 
2139  *
2140  * Symmetry condition on -X surface
2141  *
2142  * @code
2143  * {
2144  * ComponentMask component_mask_x (dim, true);
2145  * VectorTools::interpolate_boundary_values (dof_handler,
2146  * parameters.bid_BB_dirichlet_X,
2147  * ZeroFunction<dim>(dim),
2148  * boundary_values,
2149  * component_mask_x);
2150  * }
2151  * @endcode
2152  *
2153  * Symmetry condition on -X surface
2154  *
2155  * @code
2156  * {
2157  * ComponentMask component_mask_x (dim, false);
2158  * component_mask_x.set(1, true);
2159  * component_mask_x.set(2, true);
2160  * VectorTools::interpolate_boundary_values (dof_handler,
2161  * parameters.bid_BB_neumann,
2162  * ZeroFunction<dim>(dim),
2163  * boundary_values,
2164  * component_mask_x);
2165  * }
2166  * }
2167  *
2168  * @endcode
2169  *
2170  * Roller condition at central point on +X face
2171  *
2172  * @code
2173  * {
2174  * const Point<dim> roller_point (parameters.axial_length*parameters.scale,0.0,0.0);
2175  * std::vector<types::global_dof_index> fixed_dof_indices;
2176  * bool found_point_of_interest = false;
2177  *
2178  * for (typename DoFHandler<dim>::active_cell_iterator
2179  * cell = dof_handler.begin_active(),
2180  * endc = dof_handler.end(); cell != endc; ++cell)
2181  * {
2182  * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
2183  * {
2184  * @endcode
2185  *
2186  * We know that the fixed point is on the +X Neumann boundary
2187  *
2188  * @code
2189  * if (cell->face(face)->at_boundary() == true &&
2190  * cell->face(face)->boundary_id() == parameters.bid_BB_neumann)
2191  * {
2192  * for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo<dim>::vertices_per_face; ++face_vertex_index)
2193  * {
2194  * if (cell->face(face)->vertex(face_vertex_index).distance(roller_point) < 1e-6)
2195  * {
2196  * found_point_of_interest = true;
2197  * for (unsigned int index_component = 1; index_component < dim; ++index_component)
2198  * fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index,
2199  * index_component));
2200  * }
2201  *
2202  * if (found_point_of_interest == true) break;
2203  * }
2204  * }
2205  * if (found_point_of_interest == true) break;
2206  * }
2207  * if (found_point_of_interest == true) break;
2208  * }
2209  *
2210  * Assert(found_point_of_interest == true, ExcMessage("Didn't find point of interest"));
2211  * AssertThrow(fixed_dof_indices.size() == dim-1, ExcMessage("Didn't find the correct number of DoFs to fix"));
2212  *
2213  * for (unsigned int i=0; i < fixed_dof_indices.size(); ++i)
2214  * boundary_values[fixed_dof_indices[i]] = 0.0;
2215  * }
2216  * }
2217  * else
2218  * AssertThrow(false, ExcNotImplemented());
2219  *
2220  * MatrixTools::apply_boundary_values (boundary_values,
2221  * system_matrix,
2222  * solution,
2223  * system_rhs);
2224  * }
2225  *
2226  *
2227  * @endcode
2228  *
2229  *
2230  * <a name="LinearMuscleModelProblemsolve"></a>
2231  * <h4>LinearMuscleModelProblem::solve</h4>
2232  *
2233 
2234  *
2235  *
2236  * @code
2237  * template <int dim>
2238  * void LinearMuscleModelProblem<dim>::solve ()
2239  * {
2240  * SolverControl solver_control (system_matrix.m(), 1e-12);
2241  * SolverCG<> cg (solver_control);
2242  *
2243  * PreconditionSSOR<> preconditioner;
2244  * preconditioner.initialize(system_matrix, 1.2);
2245  *
2246  * cg.solve (system_matrix, solution, system_rhs,
2247  * preconditioner);
2248  *
2249  * hanging_node_constraints.distribute (solution);
2250  * }
2251  *
2252  *
2253  * @endcode
2254  *
2255  *
2256  * <a name="LinearMuscleModelProblemoutput_results"></a>
2257  * <h4>LinearMuscleModelProblem::output_results</h4>
2258  *
2259 
2260  *
2261  *
2262 
2263  *
2264  *
2265  * @code
2266  * template <int dim>
2267  * void LinearMuscleModelProblem<dim>::output_results (const unsigned int timestep,
2268  * const double time) const
2269  * {
2270  * @endcode
2271  *
2272  * Visual output: FEM results
2273  *
2274  * @code
2275  * {
2276  * std::string filename = "solution-";
2277  * filename += Utilities::int_to_string(timestep,4);
2278  * filename += ".vtk";
2279  * std::ofstream output (filename.c_str());
2280  *
2281  * DataOut<dim> data_out;
2282  * data_out.attach_dof_handler (dof_handler);
2283  *
2284  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2285  * data_component_interpretation(dim,
2286  * DataComponentInterpretation::component_is_part_of_vector);
2287  * std::vector<std::string> solution_name(dim, "displacement");
2288  *
2289  * data_out.add_data_vector (solution, solution_name,
2290  * DataOut<dim>::type_dof_data,
2291  * data_component_interpretation);
2292  * data_out.build_patches ();
2293  * data_out.write_vtk (output);
2294  * }
2295  *
2296  * @endcode
2297  *
2298  * Visual output: FEM data
2299  *
2300  * @code
2301  * {
2302  * std::string filename = "fibres-";
2303  * filename += Utilities::int_to_string(timestep,4);
2304  * filename += ".vtk";
2305  * std::ofstream output (filename.c_str());
2306  *
2307  * output
2308  * << "# vtk DataFile Version 3.0" << std::endl
2309  * << "# " << std::endl
2310  * << "ASCII"<< std::endl
2311  * << "DATASET POLYDATA"<< std::endl << std::endl;
2312  *
2313  * @endcode
2314  *
2315  * Extract fibre data from quadrature points
2316  *
2317  * @code
2318  * const unsigned int n_cells = triangulation.n_active_cells();
2319  * const unsigned int n_q_points_cell = qf_cell.size();
2320  *
2321  * @endcode
2322  *
2323  * Data that we'll be outputting
2324  *
2325  * @code
2326  * std::vector<std::string> results_fibre_names;
2327  * results_fibre_names.push_back("alpha");
2328  * results_fibre_names.push_back("epsilon_f");
2329  * results_fibre_names.push_back("epsilon_c");
2330  * results_fibre_names.push_back("epsilon_c_dot");
2331  *
2332  * const unsigned int n_results = results_fibre_names.size();
2333  * const unsigned int n_data_points = n_cells*n_q_points_cell;
2334  * std::vector< Point<dim> > output_points(n_data_points);
2335  * std::vector< Tensor<1,dim> > output_displacements(n_data_points);
2336  * std::vector< Tensor<1,dim> > output_directions(n_data_points);
2337  * std::vector< std::vector<double> > output_values(n_results, std::vector<double>(n_data_points));
2338  *
2339  * @endcode
2340  *
2341  * Displacement
2342  *
2343  * @code
2344  * std::vector< Vector<double> > u_values (n_q_points_cell,
2345  * Vector<double>(dim));
2346  * @endcode
2347  *
2348  * Displacement gradient
2349  *
2350  * @code
2351  * std::vector< std::vector< Tensor<1,dim> > > u_grads (n_q_points_cell,
2352  * std::vector<Tensor<1,dim> >(dim));
2353  *
2354  * FEValues<dim> fe_values (fe, qf_cell,
2356  * unsigned int cell_no = 0;
2357  * unsigned int fibre_no = 0;
2359  * cell = dof_handler.begin_active();
2360  * cell != dof_handler.end();
2361  * ++cell, ++cell_no)
2362  * {
2363  * fe_values.reinit (cell);
2364  * fe_values.get_function_values (solution, u_values);
2365  * fe_values.get_function_gradients (solution, u_grads);
2366  *
2367  * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell, ++fibre_no)
2368  * {
2369  * const MuscleFibre<dim> &fibre = fibre_data[cell_no][q_point_cell];
2370  * output_points[fibre_no] = fe_values.get_quadrature_points()[q_point_cell]; // Position
2371  * for (unsigned int d=0; d<dim; ++d)
2372  * output_displacements[fibre_no][d] = u_values[q_point_cell][d]; // Displacement
2373  * @endcode
2374  *
2375  * Direction (spatial configuration)
2376  *
2377  * @code
2378  * output_directions[fibre_no] = get_deformation_gradient(u_grads[q_point_cell])*fibre.get_M();
2379  * output_directions[fibre_no] /= output_directions[fibre_no].norm();
2380  *
2381  * @endcode
2382  *
2383  * Fibre values
2384  *
2385  * @code
2386  * output_values[0][fibre_no] = fibre.get_alpha();
2387  * output_values[1][fibre_no] = fibre.get_epsilon_f();
2388  * output_values[2][fibre_no] = fibre.get_epsilon_c();
2389  * output_values[3][fibre_no] = fibre.get_epsilon_c_dot();
2390  * }
2391  * }
2392  *
2393  * @endcode
2394  *
2395  * FIBRE POSITION
2396  *
2397  * @code
2398  * output
2399  * << "POINTS "
2400  * << n_data_points
2401  * << " float" << std::endl;
2402  * for (unsigned int i=0; i < n_data_points; ++i)
2403  * {
2404  * for (unsigned int j=0; j < dim; ++j)
2405  * {
2406  * output << (output_points)[i][j] << "\t";
2407  * }
2408  * output << std::endl;
2409  * }
2410  *
2411  * @endcode
2412  *
2413  * HEADER FOR POINT DATA
2414  *
2415  * @code
2416  * output << "\nPOINT_DATA "
2417  * << n_data_points
2418  * << std::endl << std::endl;
2419  *
2420  * @endcode
2421  *
2422  * FIBRE DISPLACEMENTS
2423  *
2424  * @code
2425  * output
2426  * << "VECTORS displacement float"
2427  * << std::endl;
2428  * for (unsigned int i = 0; i < n_data_points; ++i)
2429  * {
2430  * for (unsigned int j=0; j < dim; ++j)
2431  * {
2432  * output << (output_displacements)[i][j] << "\t";
2433  * }
2434  * output << std::endl;
2435  * }
2436  * output << std::endl;
2437  *
2438  * @endcode
2439  *
2440  * FIBRE DIRECTIONS
2441  *
2442  * @code
2443  * output
2444  * << "VECTORS direction float"
2445  * << std::endl;
2446  * for (unsigned int i = 0; i < n_data_points; ++i)
2447  * {
2448  * for (unsigned int j=0; j < dim; ++j)
2449  * {
2450  * output << (output_directions)[i][j] << "\t";
2451  * }
2452  * output << std::endl;
2453  * }
2454  * output << std::endl;
2455  *
2456  * @endcode
2457  *
2458  * POINT DATA
2459  *
2460  * @code
2461  * for (unsigned int v=0; v < n_results; ++v)
2462  * {
2463  * output
2464  * << "SCALARS "
2465  * << results_fibre_names[v]
2466  * << " float 1" << std::endl
2467  * << "LOOKUP_TABLE default "
2468  * << std::endl;
2469  * for (unsigned int i=0; i<n_data_points; ++i)
2470  * {
2471  * output << (output_values)[v][i] << " ";
2472  * }
2473  * output << std::endl;
2474  * }
2475  * }
2476  *
2477  * @endcode
2478  *
2479  * Output X-displacement at measured point
2480  *
2481  * @code
2482  * {
2483  * const Point<dim> meas_pt (parameters.problem == "IsotonicContraction" ?
2484  * Point<dim>(parameters.half_length_x, 0.0, 0.0) :
2485  * Point<dim>(parameters.axial_length*parameters.scale, 0.0, 0.0) );
2486  *
2487  *
2488  * const unsigned int index_of_interest = 0;
2489  * bool found_point_of_interest = false;
2491  *
2493  * cell = dof_handler.begin_active(),
2494  * endc = dof_handler.end(); cell != endc; ++cell)
2495  * {
2496  * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
2497  * {
2498  * @endcode
2499  *
2500  * We know that the measurement point is on the Neumann boundary
2501  *
2502  * @code
2503  * if (cell->face(face)->at_boundary() == true &&
2504  * ((parameters.problem == "IsotonicContraction" &&
2505  * cell->face(face)->boundary_id() == parameters.bid_CC_neumann) ||
2506  * (parameters.problem == "BicepsBrachii" &&
2507  * cell->face(face)->boundary_id() == parameters.bid_BB_neumann)) )
2508  * {
2509  * for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo<dim>::vertices_per_face; ++face_vertex_index)
2510  * {
2511  * if (cell->face(face)->vertex(face_vertex_index).distance(meas_pt) < 1e-6)
2512  * {
2513  * found_point_of_interest = true;
2514  * dof_of_interest = cell->face(face)->vertex_dof_index(face_vertex_index,
2515  * index_of_interest);
2516  * }
2517  *
2518  * if (found_point_of_interest == true) break;
2519  * }
2520  * }
2521  * if (found_point_of_interest == true) break;
2522  * }
2523  * if (found_point_of_interest == true) break;
2524  * }
2525  *
2526  * Assert(found_point_of_interest == true, ExcMessage("Didn't find point of interest"));
2527  * Assert(dof_of_interest != numbers::invalid_dof_index, ExcMessage("Didn't find DoF of interest"));
2528  * Assert(dof_of_interest < dof_handler.n_dofs(), ExcMessage("DoF index out of range"));
2529  *
2530  * const std::string filename = "displacement_POI.csv";
2531  * std::ofstream output;
2532  * if (timestep == 0)
2533  * {
2534  * output.open(filename.c_str(), std::ofstream::out);
2535  * output
2536  * << "Time [s]" << "," << "X-displacement [mm]" << std::endl;
2537  * }
2538  * else
2539  * output.open(filename.c_str(), std::ios_base::app);
2540  *
2541  * output
2542  * << time
2543  * << ","
2544  * << solution[dof_of_interest]*1e3
2545  * << std::endl;
2546  * }
2547  * }
2548  *
2549  *
2550  *
2551  * @endcode
2552  *
2553  *
2554  * <a name="LinearMuscleModelProblemrun"></a>
2556  *
2557 
2558  *
2559  *
2560  * @code
2561  * template <int dim>
2563  * {
2564  * make_grid();
2565  * setup_system ();
2566  * setup_muscle_fibres ();
2567  *
2568  * @endcode
2569  *
2570  * const bool do_grid_refinement = false;
2571  *
2572  * @code
2573  * double time = 0.0;
2574  * for (unsigned int timestep=0; time<=t_end; ++timestep, time+=dt)
2575  * {
2576  * std::cout
2577  * << "Timestep " << timestep
2578  * << " @ time " << time
2579  * << std::endl;
2580  *
2581  * @endcode
2582  *
2583  * First we update the fibre activation level
2584  * based on the current time
2585  *
2586  * @code
2587  * update_fibre_activation(time);
2588  *
2589  * @endcode
2590  *
2591  * Next we assemble the system and enforce boundary
2592  * conditions.
2593  * Here we assume that the system and fibres have
2594  * a fixed state, and we will assemble based on how
2595  * epsilon_c will update given the current state of
2596  * the body.
2597  *
2598  * @code
2599  * assemble_system (time);
2600  * apply_boundary_conditions ();
2601  *
2602  * @endcode
2603  *
2604  * Then we solve the linear system
2605  *
2606  * @code
2607  * solve ();
2608  *
2609  * @endcode
2610  *
2611  * Now we update the fibre state based on the new
2612  * displacement solution and the constitutive
2613  * parameters assumed to govern the stiffness of
2614  * the fibres at the previous state. i.e. We
2615  * follow through with assumed update conditions
2616  * used in the assembly phase.
2617  *
2618  * @code
2619  * update_fibre_state();
2620  *
2621  * @endcode
2622  *
2623  * Output some values to file
2624  *
2625  * @code
2626  * output_results (timestep, time);
2627  * }
2628  * }
2629  * }
2630  *
2631  * @endcode
2632  *
2633  *
2634  * <a name="Thecodemaincodefunction"></a>
2635  * <h3>The <code>main</code> function</h3>
2636  *
2637 
2638  *
2639  *
2640  * @code
2641  * int main ()
2642  * {
2643  * try
2644  * {
2645  * ::deallog.depth_console (0);
2646  * const unsigned int dim = 3;
2647  *
2648  * LMM::LinearMuscleModelProblem<dim> lmm_problem ("parameters.prm");
2649  * lmm_problem.run ();
2650  * }
2651  * catch (std::exception &exc)
2652  * {
2653  * std::cerr << std::endl << std::endl
2654  * << "----------------------------------------------------"
2655  * << std::endl;
2656  * std::cerr << "Exception on processing: " << std::endl
2657  * << exc.what() << std::endl
2658  * << "Aborting!" << std::endl
2659  * << "----------------------------------------------------"
2660  * << std::endl;
2661  *
2662  * return 1;
2663  * }
2664  * catch (...)
2665  * {
2666  * std::cerr << std::endl << std::endl
2667  * << "----------------------------------------------------"
2668  * << std::endl;
2669  * std::cerr << "Unknown exception!" << std::endl
2670  * << "Aborting!" << std::endl
2671  * << "----------------------------------------------------"
2672  * << std::endl;
2673  * return 1;
2674  * }
2675  *
2676  * return 0;
2677  * }
2678  * @endcode
2679 
2680 
2681 */
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Physics::Transformations::Rotations::rotation_matrix_3d
Tensor< 2, 3, Number > rotation_matrix_3d(const Point< 3, Number > &axis, const Number &angle)
internal::QGaussLobatto::gamma
long double gamma(const unsigned int n)
Definition: quadrature_lib.cc:96
numbers::invalid_dof_index
const types::global_dof_index invalid_dof_index
Definition: types.h:206
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
FE_Q
Definition: fe_q.h:554
SymmetricTensor< 2, dim >
dealii
Definition: namespace_dealii.h:25
GridGenerator::flatten_triangulation
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
Differentiation::SD::OptimizerType::lambda
@ lambda
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
mu
Definition: function_parser.h:32
Triangulation< dim >
parallel::transform
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:211
Triangulation::refine_global
void refine_global(const unsigned int times=1)
Definition: tria.cc:10851
SparseMatrix< double >
Patterns::Bool
Definition: patterns.h:984
SphericalManifold
Definition: manifold_lib.h:231
GridGenerator::hyper_rectangle
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler< dim >
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
std::cos
inline ::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5324
deallog
LogStream deallog
Definition: logstream.cc:37
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
FEValues< dim >
Triangulation::set_manifold
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
Definition: tria.cc:10245
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
level
unsigned int level
Definition: grid_out.cc:4355
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
GridGenerator::cylinder
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
Triangulation::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12013
Tensor< 1, dim >
Triangulation::set_all_manifold_ids
void set_all_manifold_ids(const types::manifold_id number)
Definition: tria.cc:10277
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
numbers::E
static constexpr double E
Definition: numbers.h:212
SparsityPattern
Definition: sparsity_pattern.h:865
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
TriaActiveIterator
Definition: tria_iterator.h:759
QGauss
Definition: quadrature_lib.h:40
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
LogStream::depth_console
unsigned int depth_console(const unsigned int n)
Definition: logstream.cc:349
unsigned int
GridGenerator::hyper_ball
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AffineConstraints< double >
DataOutBase::eps
@ eps
Definition: data_out_base.h:1582
std::sqrt
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
GridGenerator::extrude_triangulation
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
GridTools::transform
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
DerivativeApproximation::internal::approximate
void approximate(SynchronousIterators< std::tuple< TriaActiveIterator< ::DoFCellAccessor< DoFHandlerType< dim, spacedim >, false >>, Vector< float >::iterator >> const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Definition: derivative_approximation.cc:924
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
ParameterHandler
Definition: parameter_handler.h:845
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
Patterns::Selection
Definition: patterns.h:383
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12618
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:935
Vector< double >
FESystem
Definition: fe.h:44
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
ParameterHandler::parse_input
virtual void parse_input(std::istream &input, const std::string &filename="input file", const std::string &last_line="", const bool skip_undefined=false)
Definition: parameter_handler.cc:399
Patterns::Double
Definition: patterns.h:293
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
Triangulation::end
cell_iterator end() const
Definition: tria.cc:12079
Patterns::Integer
Definition: patterns.h:190