Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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/affine_constraints.h>
294 *   #include <deal.II/lac/vector.h>
295 *   #include <deal.II/lac/full_matrix.h>
296 *   #include <deal.II/lac/sparse_matrix.h>
297 *   #include <deal.II/lac/solver_cg.h>
298 *   #include <deal.II/lac/precondition.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/dofs/dof_handler.h>
308 *   #include <deal.II/dofs/dof_accessor.h>
309 *   #include <deal.II/dofs/dof_tools.h>
310 *   #include <deal.II/fe/fe_values.h>
311 *   #include <deal.II/fe/fe_system.h>
312 *   #include <deal.II/fe/fe_q.h>
313 *   #include <deal.II/numerics/vector_tools.h>
314 *   #include <deal.II/numerics/matrix_tools.h>
315 *   #include <deal.II/numerics/data_out.h>
316 *   #include <deal.II/numerics/error_estimator.h>
317 *   #include <deal.II/physics/transformations.h>
318 *  
319 *   #include <fstream>
320 *   #include <iostream>
321 *   #include <vector>
322 *  
323 *   namespace LMM
324 *   {
325 *   using namespace dealii;
326 *  
327 * @endcode
328 *
329 *
330 * <a name="Runtimeparameters"></a>
331 * <h3>Run-time parameters</h3>
332 *
333
334 *
335 * There are several parameters that can be set in the code so we set up a
336 * ParameterHandler object to read in the choices at run-time.
337 *
338 * @code
339 *   namespace Parameters
340 *   {
341 * @endcode
342 *
343 *
344 * <a name="FiniteElementsystem"></a>
345 * <h4>Finite Element system</h4>
346 *
347
348 *
349 * Here we specify the polynomial order used to approximate the solution.
350 * The quadrature order should be adjusted accordingly.
351 *
352 * @code
353 *   struct FESystem
354 *   {
355 *   unsigned int poly_degree;
356 *   unsigned int quad_order;
357 *  
358 *   static void
359 *   declare_parameters(ParameterHandler &prm);
360 *  
361 *   void
362 *   parse_parameters(ParameterHandler &prm);
363 *   };
364 *  
365 *   void FESystem::declare_parameters(ParameterHandler &prm)
366 *   {
367 *   prm.enter_subsection("Finite element system");
368 *   {
369 *   prm.declare_entry("Polynomial degree", "1",
370 *   Patterns::Integer(0),
371 *   "Displacement system polynomial order");
372 *  
373 *   prm.declare_entry("Quadrature order", "2",
374 *   Patterns::Integer(0),
375 *   "Gauss quadrature order");
376 *   }
377 *   prm.leave_subsection();
378 *   }
379 *  
380 *   void FESystem::parse_parameters(ParameterHandler &prm)
381 *   {
382 *   prm.enter_subsection("Finite element system");
383 *   {
384 *   poly_degree = prm.get_integer("Polynomial degree");
385 *   quad_order = prm.get_integer("Quadrature order");
386 *   }
387 *   prm.leave_subsection();
388 *   }
389 *  
390 * @endcode
391 *
392 *
393 * <a name="Problem"></a>
394 * <h4>Problem</h4>
395 *
396
397 *
398 * Choose which problem is going to be solved
399 *
400 * @code
401 *   struct Problem
402 *   {
403 *   std::string problem;
404 *  
405 *   static void
406 *   declare_parameters(ParameterHandler &prm);
407 *  
408 *   void
409 *   parse_parameters(ParameterHandler &prm);
410 *   };
411 *  
412 *   void Problem::declare_parameters(ParameterHandler &prm)
413 *   {
414 *   prm.enter_subsection("Problem");
415 *   {
416 *   prm.declare_entry("Problem", "IsotonicContraction",
417 *   Patterns::Selection("IsotonicContraction|BicepsBrachii"),
418 *   "The problem that is to be solved");
419 *   }
420 *   prm.leave_subsection();
421 *   }
422 *  
423 *   void Problem::parse_parameters(ParameterHandler &prm)
424 *   {
425 *   prm.enter_subsection("Problem");
426 *   {
427 *   problem = prm.get("Problem");
428 *   }
429 *   prm.leave_subsection();
430 *   }
431 *  
432 * @endcode
433 *
434 *
435 * <a name="IsotonicContractionGeometry"></a>
436 * <h4>IsotonicContractionGeometry</h4>
437 *
438
439 *
440 * Make adjustments to the geometry and discretisation of the
441 * isotonic contraction model from Martins2006.
442 *
443
444 *
445 *
446 * @code
447 *   struct IsotonicContraction
448 *   {
449 *   const double half_length_x = 10e-3/2.0;
450 *   const double half_length_y = 10e-3/2.0;
451 *   const double half_length_z = 1e-3/2.0;
452 *   const types::boundary_id bid_CC_dirichlet_symm_X = 1;
453 *   const types::boundary_id bid_CC_dirichlet_symm_Z = 2;
454 *   const types::boundary_id bid_CC_neumann = 10;
455 *  
456 *   static void
457 *   declare_parameters(ParameterHandler &prm);
458 *  
459 *   void
460 *   parse_parameters(ParameterHandler &prm);
461 *   };
462 *  
463 *   void IsotonicContraction::declare_parameters(ParameterHandler &/*prm*/)
464 *   {
465 *  
466 *   }
467 *  
468 *   void IsotonicContraction::parse_parameters(ParameterHandler &/*prm*/)
469 *   {
470 *  
471 *   }
472 *  
473 * @endcode
474 *
475 *
476 * <a name="BicepsBrachiiGeometry"></a>
477 * <h4>BicepsBrachiiGeometry</h4>
478 *
479
480 *
481 * Make adjustments to the geometry and discretisation of the
482 * biceps model.
483 *
484
485 *
486 *
487 * @code
488 *   struct BicepsBrachii
489 *   {
490 *   double axial_length;
491 *   double radius_insertion_origin;
492 *   double radius_midpoint;
493 *   double scale;
494 *   unsigned int elements_along_axis;
495 *   unsigned int n_refinements_radial;
496 *   bool include_gravity;
497 *   double axial_force;
498 *  
499 *   const types::boundary_id bid_BB_dirichlet_X = 1;
500 *   const types::boundary_id bid_BB_neumann = 10;
501 *   const types::boundary_id mid_BB_radial = 100;
502 *  
503 *   static void
504 *   declare_parameters(ParameterHandler &prm);
505 *  
506 *   void
507 *   parse_parameters(ParameterHandler &prm);
508 *   };
509 *  
510 *   void BicepsBrachii::declare_parameters(ParameterHandler &prm)
511 *   {
512 *   prm.enter_subsection("Biceps Brachii geometry");
513 *   {
514 *   prm.declare_entry("Axial length", "250",
515 *   Patterns::Double(0),
516 *   "Axial length of the muscle");
517 *  
518 *   prm.declare_entry("Radius insertion and origin", "5",
519 *   Patterns::Double(0),
520 *   "Insertion and origin radius");
521 *  
522 *   prm.declare_entry("Radius midpoint", "7.5",
523 *   Patterns::Double(0),
524 *   "Radius at the midpoint of the muscle");
525 *  
526 *   prm.declare_entry("Grid scale", "1e-3",
527 *   Patterns::Double(0.0),
528 *   "Global grid scaling factor");
529 *  
530 *   prm.declare_entry("Elements along axis", "32",
531 *   Patterns::Integer(2),
532 *   "Number of elements along the muscle axis");
533 *  
534 *   prm.declare_entry("Radial refinements", "4",
535 *   Patterns::Integer(0),
536 *   "Control the discretisation in the radial direction");
537 *  
538 *   prm.declare_entry("Gravity", "false",
539 *   Patterns::Bool(),
540 *   "Include the effects of gravity (in the y-direction; "
541 *   " perpendicular to the muscle axis)");
542 *  
543 *   prm.declare_entry("Axial force", "1",
544 *   Patterns::Double(),
545 *   "Applied distributed axial force (in Newtons)");
546 *   }
547 *   prm.leave_subsection();
548 *   }
549 *  
550 *   void BicepsBrachii::parse_parameters(ParameterHandler &prm)
551 *   {
552 *   prm.enter_subsection("Biceps Brachii geometry");
553 *   {
554 *   axial_length = prm.get_double("Axial length");
555 *   radius_insertion_origin = prm.get_double("Radius insertion and origin");
556 *   radius_midpoint = prm.get_double("Radius midpoint");
557 *   scale = prm.get_double("Grid scale");
558 *   elements_along_axis = prm.get_integer("Elements along axis");
559 *   n_refinements_radial = prm.get_integer("Radial refinements");
560 *   include_gravity = prm.get_bool("Gravity");
561 *   axial_force = prm.get_double("Axial force");
562 *   }
563 *   prm.leave_subsection();
564 *  
565 *   AssertThrow(radius_midpoint >= radius_insertion_origin,
566 *   ExcMessage("Unrealistic geometry"));
567 *   }
568 *  
569 * @endcode
570 *
571 *
572 * <a name="Neurologicalsignal"></a>
573 * <h4>Neurological signal</h4>
574 *
575
576 *
577 *
578 * @code
579 *   struct NeurologicalSignal
580 *   {
581 *   double neural_signal_start_time;
582 *   double neural_signal_end_time;
583 *  
584 *   static void
585 *   declare_parameters(ParameterHandler &prm);
586 *  
587 *   void
588 *   parse_parameters(ParameterHandler &prm);
589 *   };
590 *  
591 *   void NeurologicalSignal::declare_parameters(ParameterHandler &prm)
592 *   {
593 *   prm.enter_subsection("Neurological signal");
594 *   {
595 *   prm.declare_entry("Start time", "1.0",
596 *   Patterns::Double(0),
597 *   "Time at which to start muscle activation");
598 *  
599 *   prm.declare_entry("End time", "2.0",
600 *   Patterns::Double(0),
601 *   "Time at which to remove muscle activation signal");
602 *   }
603 *   prm.leave_subsection();
604 *   }
605 *  
606 *   void NeurologicalSignal::parse_parameters(ParameterHandler &prm)
607 *   {
608 *   prm.enter_subsection("Neurological signal");
609 *   {
610 *   neural_signal_start_time = prm.get_double("Start time");
611 *   neural_signal_end_time = prm.get_double("End time");
612 *   }
613 *   prm.leave_subsection();
614 *  
615 *   Assert(neural_signal_start_time < neural_signal_end_time,
616 *   ExcMessage("Invalid neural signal times."));
617 *   }
618 *  
619 * @endcode
620 *
621 *
622 * <a name="Time"></a>
623 * <h4>Time</h4>
624 *
625
626 *
627 * Set the timestep size @f$ \varDelta t @f$ and the simulation end-time.
628 *
629 * @code
630 *   struct Time
631 *   {
632 *   double delta_t;
633 *   double end_time;
634 *   double end_ramp_time;
635 *  
636 *   static void
637 *   declare_parameters(ParameterHandler &prm);
638 *  
639 *   void
640 *   parse_parameters(ParameterHandler &prm);
641 *   };
642 *  
643 *   void Time::declare_parameters(ParameterHandler &prm)
644 *   {
645 *   prm.enter_subsection("Time");
646 *   {
647 *   prm.declare_entry("End time", "3",
648 *   Patterns::Double(0),
649 *   "End time");
650 *  
651 *   prm.declare_entry("End ramp time", "1",
652 *   Patterns::Double(0),
653 *   "Force ramp end time");
654 *  
655 *   prm.declare_entry("Time step size", "0.1",
656 *   Patterns::Double(0),
657 *   "Time step size");
658 *   }
659 *   prm.leave_subsection();
660 *   }
661 *  
662 *   void Time::parse_parameters(ParameterHandler &prm)
663 *   {
664 *   prm.enter_subsection("Time");
665 *   {
666 *   end_time = prm.get_double("End time");
667 *   end_ramp_time = prm.get_double("End ramp time");
668 *   delta_t = prm.get_double("Time step size");
669 *   }
670 *   prm.leave_subsection();
671 *   }
672 *  
673 * @endcode
674 *
675 *
676 * <a name="Allparameters"></a>
677 * <h4>All parameters</h4>
678 *
679
680 *
681 * Finally we consolidate all of the above structures into a single container
682 * that holds all of our run-time selections.
683 *
684 * @code
685 *   struct AllParameters : public FESystem,
686 *   public Problem,
687 *   public IsotonicContraction,
688 *   public BicepsBrachii,
689 *   public NeurologicalSignal,
690 *   public Time
691 *   {
692 *   AllParameters(const std::string &input_file);
693 *  
694 *   static void
695 *   declare_parameters(ParameterHandler &prm);
696 *  
697 *   void
698 *   parse_parameters(ParameterHandler &prm);
699 *   };
700 *  
701 *   AllParameters::AllParameters(const std::string &input_file)
702 *   {
703 *   ParameterHandler prm;
704 *   declare_parameters(prm);
705 *   prm.parse_input(input_file);
706 *   parse_parameters(prm);
707 *   }
708 *  
709 *   void AllParameters::declare_parameters(ParameterHandler &prm)
710 *   {
711 *   FESystem::declare_parameters(prm);
712 *   Problem::declare_parameters(prm);
713 *   IsotonicContraction::declare_parameters(prm);
714 *   BicepsBrachii::declare_parameters(prm);
715 *   NeurologicalSignal::declare_parameters(prm);
716 *   Time::declare_parameters(prm);
717 *   }
718 *  
719 *   void AllParameters::parse_parameters(ParameterHandler &prm)
720 *   {
721 *   FESystem::parse_parameters(prm);
722 *   Problem::parse_parameters(prm);
723 *   IsotonicContraction::parse_parameters(prm);
724 *   BicepsBrachii::parse_parameters(prm);
725 *   NeurologicalSignal::parse_parameters(prm);
726 *   Time::parse_parameters(prm);
727 *  
728 * @endcode
729 *
730 * Override time setting for test defined
731 * in the literature
732 *
733 * @code
734 *   if (problem == "IsotonicContraction")
735 *   {
736 *   end_time = 3.0;
737 *   end_ramp_time = 1.0;
738 *   delta_t = 0.1;
739 *  
740 *   neural_signal_start_time = 1.0;
741 *   neural_signal_end_time = 2.0;
742 *   }
743 *   }
744 *   }
745 *  
746 * @endcode
747 *
748 *
749 * <a name="Bodyforcevalues"></a>
750 * <h3>Body force values</h3>
751 *
752
753 *
754 *
755 * @code
756 *   template <int dim>
757 *   class BodyForce : public Function<dim>
758 *   {
759 *   public:
760 *   BodyForce (const double rho,
761 *   const Tensor<1,dim> direction);
762 *   virtual ~BodyForce () {}
763 *  
764 *   virtual void vector_value (const Point<dim> &p,
765 *   Vector<double> &values) const override;
766 *  
767 *   virtual void vector_value_list (const std::vector<Point<dim> > &points,
768 *   std::vector<Vector<double> > &value_list) const override;
769 *  
770 *   const double rho;
771 *   const double g;
772 *   const Tensor<1,dim> M;
773 *   };
774 *  
775 *  
776 *   template <int dim>
777 *   BodyForce<dim>::BodyForce (const double rho,
778 *   const Tensor<1,dim> direction)
779 *   :
780 *   Function<dim> (dim),
781 *   rho (rho),
782 *   g (9.81),
783 *   M (direction)
784 *   {
785 *   Assert(M.norm() == 1.0, ExcMessage("Direction vector is not a unit vector"));
786 *   }
787 *  
788 *  
789 *   template <int dim>
790 *   inline
791 *   void BodyForce<dim>::vector_value (const Point<dim> &/*p*/,
792 *   Vector<double> &values) const
793 *   {
794 *   Assert (values.size() == dim,
795 *   ExcDimensionMismatch (values.size(), dim));
796 *   Assert (dim >= 2, ExcNotImplemented());
797 *   for (unsigned int d=0; d<dim; ++d)
798 *   {
799 *   values(d) = rho*g*M[d];
800 *   }
801 *   }
802 *  
803 *  
804 *   template <int dim>
805 *   void BodyForce<dim>::vector_value_list (const std::vector<Point<dim> > &points,
806 *   std::vector<Vector<double> > &value_list) const
807 *   {
808 *   Assert (value_list.size() == points.size(),
809 *   ExcDimensionMismatch (value_list.size(), points.size()));
810 *  
811 *   const unsigned int n_points = points.size();
812 *  
813 *   for (unsigned int p=0; p<n_points; ++p)
814 *   BodyForce<dim>::vector_value (points[p],
815 *   value_list[p]);
816 *   }
817 *  
818 *   template <int dim>
819 *   class Traction : public Function<dim>
820 *   {
821 *   public:
822 *   Traction (const double force,
823 *   const double area);
824 *   virtual ~Traction () {}
825 *  
826 *   virtual void vector_value (const Point<dim> &p,
827 *   Vector<double> &values) const override;
828 *  
829 *   virtual void vector_value_list (const std::vector<Point<dim> > &points,
830 *   std::vector<Vector<double> > &value_list) const override;
831 *  
832 *   const double t;
833 *   };
834 *  
835 *  
836 *   template <int dim>
837 *   Traction<dim>::Traction (const double force,
838 *   const double area)
839 *   :
840 *   Function<dim> (dim),
841 *   t (force/area)
842 *   {}
843 *  
844 *  
845 *   template <int dim>
846 *   inline
847 *   void Traction<dim>::vector_value (const Point<dim> &/*p*/,
848 *   Vector<double> &values) const
849 *   {
850 *   Assert (values.size() == dim,
851 *   ExcDimensionMismatch (values.size(), dim));
852 *   Assert (dim == 3, ExcNotImplemented());
853 *  
854 * @endcode
855 *
856 * Assume uniform distributed load
857 *
858 * @code
859 *   values(0) = t;
860 *   values(1) = 0.0;
861 *   values(2) = 0.0;
862 *   }
863 *  
864 *  
865 *   template <int dim>
866 *   void Traction<dim>::vector_value_list (const std::vector<Point<dim> > &points,
867 *   std::vector<Vector<double> > &value_list) const
868 *   {
869 *   Assert (value_list.size() == points.size(),
870 *   ExcDimensionMismatch (value_list.size(), points.size()));
871 *  
872 *   const unsigned int n_points = points.size();
873 *  
874 *   for (unsigned int p=0; p<n_points; ++p)
875 *   Traction<dim>::vector_value (points[p],
876 *   value_list[p]);
877 *   }
878 *  
879 * @endcode
880 *
881 *
882 * <a name="Utilityfunctions"></a>
883 * <h3>Utility functions</h3>
884 *
885
886 *
887 *
888 * @code
889 *   template <int dim>
890 *   inline
891 *   Tensor<2,dim> get_deformation_gradient (std::vector<Tensor<1,dim> > &grad)
892 *   {
893 *   Assert (grad.size() == dim, ExcInternalError());
894 *  
895 *   Tensor<2,dim> F (unit_symmetric_tensor<dim>());
896 *   for (unsigned int i=0; i<dim; ++i)
897 *   for (unsigned int j=0; j<dim; ++j)
898 *   F[i][j] += grad[i][j];
899 *   return F;
900 *   }
901 *  
902 *   template <int dim>
903 *   inline
904 *   SymmetricTensor<2,dim> get_small_strain (std::vector<Tensor<1,dim> > &grad)
905 *   {
906 *   Assert (grad.size() == dim, ExcInternalError());
907 *  
908 *   SymmetricTensor<2,dim> strain;
909 *   for (unsigned int i=0; i<dim; ++i)
910 *   strain[i][i] = grad[i][i];
911 *  
912 *   for (unsigned int i=0; i<dim; ++i)
913 *   for (unsigned int j=i+1; j<dim; ++j)
914 *   strain[i][j] = (grad[i][j] + grad[j][i]) / 2;
915 *   return strain;
916 *   }
917 *  
918 * @endcode
919 *
920 *
921 * <a name="Propertiesformusclematrix"></a>
922 * <h3>Properties for muscle matrix</h3>
923 *
924
925 *
926 *
927 * @code
928 *   struct MuscleMatrix
929 *   {
930 *   static const double E; // Young's modulus
931 *   static const double nu; // Poisson ratio
932 *  
933 *   static const double mu; // Shear modulus
934 *   static const double lambda; // Lame parameter
935 *   };
936 *  
937 *   const double MuscleMatrix::E = 26e3;
938 *   const double MuscleMatrix::nu = 0.45;
939 *   const double MuscleMatrix::mu = MuscleMatrix::E/(2.0*(1.0 + MuscleMatrix::nu));
940 *   const double MuscleMatrix::lambda = 2.0*MuscleMatrix::mu *MuscleMatrix::nu/(1.0 - 2.0*MuscleMatrix::nu);
941 *  
942 * @endcode
943 *
944 *
945 * <a name="Localdataformusclefibres"></a>
946 * <h3>Local data for muscle fibres</h3>
947 *
948
949 *
950 *
951 * @code
952 *   #define convert_gf_to_N 1.0/101.97
953 *   #define convert_gf_per_cm2_to_N_per_m2 convert_gf_to_N*1e2*1e2
954 *   #define T0 6280.0*convert_gf_per_cm2_to_N_per_m2
955 *  
956 * @endcode
957 *
958 * A struct that governs the functioning of a single muscle fibre
959 *
960 * @code
961 *   template <int dim>
962 *   struct MuscleFibre
963 *   {
964 *   MuscleFibre (void)
965 *   : alpha (0.0),
966 *   alpha_t1 (0.0),
967 *   epsilon_f (0.0),
968 *   epsilon_c (0.0),
969 *   epsilon_c_t1 (0.0),
970 *   epsilon_c_dot (0.0)
971 *   {
972 *  
973 *   }
974 *  
975 *   MuscleFibre(const Tensor<1,dim> &direction)
976 *   : M (direction),
977 *   alpha (0.0),
978 *   alpha_t1 (0.0),
979 *   epsilon_f (0.0),
980 *   epsilon_c (0.0),
981 *   epsilon_c_t1 (0.0),
982 *   epsilon_c_dot (0.0)
983 *   {
984 *   Assert(M.norm() == 1.0,
985 *   ExcMessage("Fibre direction is not a unit vector"));
986 *   }
987 *  
988 *   void update_alpha (const double u,
989 *   const double dt);
990 *  
991 *   void update_state(const SymmetricTensor<2,dim> &strain_tensor,
992 *   const double dt);
993 *  
994 *   const Tensor<1,dim> &get_M () const
995 *   {
996 *   return M;
997 *   }
998 *   double get_m_p () const;
999 *   double get_m_s () const;
1000 *   double get_beta (const double dt) const;
1001 *   double get_gamma (const double dt) const;
1002 *  
1003 * @endcode
1004 *
1005 * Postprocessing
1006 *
1007 * @code
1008 *   const double &get_alpha() const
1009 *   {
1010 *   return alpha;
1011 *   }
1012 *   const double &get_epsilon_f() const
1013 *   {
1014 *   return epsilon_f;
1015 *   }
1016 *   const double &get_epsilon_c() const
1017 *   {
1018 *   return epsilon_c;
1019 *   }
1020 *   const double &get_epsilon_c_dot() const
1021 *   {
1022 *   return epsilon_c_dot;
1023 *   }
1024 *  
1025 *   private:
1026 *   Tensor<1,dim> M; // Direction
1027 *  
1028 *   double alpha; // Activation level at current timestep
1029 *   double alpha_t1; // Activation level at previous timestep
1030 *  
1031 *   double epsilon_f; // Fibre strain at current timestep
1032 *   double epsilon_c; // Contractile strain at current timestep
1033 *   double epsilon_c_t1; // Contractile strain at previous timestep
1034 *   double epsilon_c_dot; // Contractile velocity at previous timestep
1035 *  
1036 *   double get_f_c_L () const;
1037 *   double get_m_c_V () const;
1038 *   double get_c_c_V () const;
1039 *   };
1040 *  
1041 *   template <int dim>
1042 *   void MuscleFibre<dim>::update_alpha (const double u,
1043 *   const double dt)
1044 *   {
1045 *   static const double tau_r = 0.15; // s
1046 *   static const double tau_f = 0.15; // s
1047 *   static const double alpha_min = 0;
1048 *  
1049 *   if (u == 1.0)
1050 *   alpha = (alpha_t1*tau_r*tau_f + dt*tau_f) / (tau_r*tau_f + dt*tau_f);
1051 *   else if (u == 0)
1052 *   alpha = (alpha_t1*tau_r*tau_f + dt*alpha_min*tau_r) / (tau_r*tau_f + dt*tau_r);
1053 *   else
1054 *   {
1055 *   const double b = 1.0/tau_r - 1.0/tau_f;
1056 *   const double c = 1.0/tau_f;
1057 *   const double d = alpha_min/tau_f;
1058 *   const double f1 = 1.0/tau_r - alpha_min/tau_f;
1059 *   const double p = b*u + c;
1060 *   const double q = f1*u + d;
1061 *  
1062 *   alpha = (q*dt + alpha_t1)/(1.0 + p*dt);
1063 *   }
1064 *   }
1065 *  
1066 *  
1067 *   template <int dim>
1068 *   double MuscleFibre<dim>::get_m_p () const
1069 *   {
1070 *   static const double A = 8.568e-4*convert_gf_per_cm2_to_N_per_m2;
1071 *   static const double a = 12.43;
1072 *   if (epsilon_f >= 0.0)
1073 *   {
1074 * @endcode
1075 *
1076 * 100 times more compliant than Martins2006
1077 *
1078 * @code
1079 *   static const double m_p = 2.0*A*a/1e2;
1080 *   return m_p;
1081 *   }
1082 *   else
1083 *   return 0.0;
1084 *   }
1085 *  
1086 *   template <int dim>
1087 *   double MuscleFibre<dim>::get_m_s (void) const
1088 *   {
1089 *   const double epsilon_s = epsilon_f - epsilon_c; // Small strain assumption
1090 *   if (epsilon_s >= -1e-6) // Tolerant check
1091 *   return 10.0;
1092 *   else
1093 *   return 0.0;
1094 *   }
1095 *  
1096 *   template <int dim>
1097 *   double MuscleFibre<dim>::get_f_c_L (void) const
1098 *   {
1099 *   if (epsilon_c <= 0.5 && epsilon_c >= -0.5)
1100 *   return 1.0;
1101 *   else
1102 *   return 0.0;
1103 *   }
1104 *  
1105 *   template <int dim>
1106 *   double MuscleFibre<dim>::get_m_c_V (void) const
1107 *   {
1108 *   if (epsilon_c_dot < -5.0)
1109 *   return 0.0;
1110 *   else if (epsilon_c_dot <= 3.0)
1111 *   return 1.0/5.0;
1112 *   else
1113 *   return 0.0;
1114 *   }
1115 *  
1116 *   template <int dim>
1117 *   double MuscleFibre<dim>::get_c_c_V (void) const
1118 *   {
1119 *   if (epsilon_c_dot < -5.0)
1120 *   return 0.0;
1121 *   else if (epsilon_c_dot <= 3.0)
1122 *   return 1.0;
1123 *   else
1124 *   return 1.6;
1125 *   }
1126 *  
1127 *   template <int dim>
1128 *   double MuscleFibre<dim>::get_beta(const double dt) const
1129 *   {
1130 *   return get_f_c_L()*get_m_c_V()*alpha/dt + get_m_s();
1131 *   }
1132 *  
1133 *   template <int dim>
1134 *   double MuscleFibre<dim>::get_gamma(const double dt) const
1135 *   {
1136 *   return get_f_c_L()*alpha*(get_m_c_V()*epsilon_c_t1/dt - get_c_c_V());
1137 *   }
1138 *  
1139 *   template <int dim>
1140 *   void MuscleFibre<dim>::update_state(const SymmetricTensor<2,dim> &strain_tensor,
1141 *   const double dt)
1142 *   {
1143 * @endcode
1144 *
1145 * Values from previous state
1146 * These were the values that were used in the assembly,
1147 * so we must use them in the update step to be consistant.
1148 * Need to compute these before we overwrite epsilon_c_t1
1149 *
1150 * @code
1151 *   const double m_s = get_m_s();
1152 *   const double beta = get_beta(dt);
1153 *   const double gamma = get_gamma(dt);
1154 *  
1155 * @endcode
1156 *
1157 * Update current state
1158 *
1159 * @code
1160 *   alpha_t1 = alpha;
1161 *   epsilon_f = M*static_cast< Tensor<2,dim> >(strain_tensor)*M;
1162 *   epsilon_c_t1 = epsilon_c;
1163 *   epsilon_c = (m_s*epsilon_f + gamma)/beta;
1164 *   epsilon_c_dot = (epsilon_c - epsilon_c_t1)/dt;
1165 *   }
1166 *  
1167 *  
1168 * @endcode
1169 *
1170 *
1171 * <a name="ThecodeLinearMuscleModelProblemcodeclasstemplate"></a>
1172 * <h3>The <code>LinearMuscleModelProblem</code> class template</h3>
1173 *
1174
1175 *
1176 *
1177 * @code
1178 *   template <int dim>
1179 *   class LinearMuscleModelProblem
1180 *   {
1181 *   public:
1182 *   LinearMuscleModelProblem (const std::string &input_file);
1183 *   ~LinearMuscleModelProblem ();
1184 *   void run ();
1185 *  
1186 *   private:
1187 *   void make_grid ();
1188 *   void setup_muscle_fibres ();
1189 *   double get_neural_signal (const double time);
1190 *   void update_fibre_activation (const double time);
1191 *   void update_fibre_state ();
1192 *   void setup_system ();
1193 *   void assemble_system (const double time);
1194 *   void apply_boundary_conditions ();
1195 *   void solve ();
1196 *   void output_results (const unsigned int timestep,
1197 *   const double time) const;
1198 *  
1199 *   Parameters::AllParameters parameters;
1200 *  
1201 *   Triangulation<dim> triangulation;
1202 *   DoFHandler<dim> dof_handler;
1203 *  
1204 *   FESystem<dim> fe;
1205 *   QGauss<dim> qf_cell;
1206 *   QGauss<dim-1> qf_face;
1207 *  
1208 *   AffineConstraints<double> hanging_node_constraints;
1209 *  
1210 *   SparsityPattern sparsity_pattern;
1211 *   SparseMatrix<double> system_matrix;
1212 *  
1213 *   Vector<double> solution;
1214 *   Vector<double> system_rhs;
1215 *  
1216 * @endcode
1217 *
1218 * Time
1219 *
1220 * @code
1221 *   const double t_end;
1222 *   const double dt;
1223 *   const double t_ramp_end; // Force ramp end time
1224 *  
1225 * @endcode
1226 *
1227 * Loading
1228 *
1229 * @code
1230 *   const BodyForce<dim> body_force;
1231 *   const Traction<dim> traction;
1232 *  
1233 * @endcode
1234 *
1235 * Local data
1236 *
1237 * @code
1238 *   std::vector< std::vector<MuscleFibre<dim> > > fibre_data;
1239 *  
1240 * @endcode
1241 *
1242 * Constitutive functions for assembly
1243 *
1244 * @code
1245 *   SymmetricTensor<4,dim> get_stiffness_tensor (const unsigned int cell,
1246 *   const unsigned int q_point_cell) const;
1247 *   SymmetricTensor<2,dim> get_rhs_tensor (const unsigned int cell,
1248 *   const unsigned int q_point_cell) const;
1249 *   };
1250 *  
1251 * @endcode
1252 *
1253 *
1254 * <a name="LinearMuscleModelProblemLinearMuscleModelProblem"></a>
1255 * <h4>LinearMuscleModelProblem::LinearMuscleModelProblem</h4>
1256 *
1257
1258 *
1259 *
1260 * @code
1261 *   template <int dim>
1262 *   LinearMuscleModelProblem<dim>::LinearMuscleModelProblem (const std::string &input_file)
1263 *   :
1264 *   parameters(input_file),
1265 *   dof_handler (triangulation),
1266 *   fe (FE_Q<dim>(parameters.poly_degree), dim),
1267 *   qf_cell (parameters.quad_order),
1268 *   qf_face (parameters.quad_order),
1269 *   t_end (parameters.end_time),
1270 *   dt (parameters.delta_t),
1271 *   t_ramp_end(parameters.end_ramp_time),
1272 *   body_force ((parameters.problem == "BicepsBrachii" &&parameters.include_gravity == true) ?
1273 *   BodyForce<dim>(0.375*1000.0, Tensor<1,dim>({0,-1,0})) : // (reduced) Density and direction
1274 *   BodyForce<dim>(0.0, Tensor<1,dim>({0,0,1})) ),
1275 *   traction (parameters.problem == "BicepsBrachii" ?
1276 *   Traction<dim>(parameters.axial_force, // Force, area
1277 *   M_PI*std::pow(parameters.radius_insertion_origin *parameters.scale,2.0) ) :
1278 *   Traction<dim>(4.9*convert_gf_to_N, // Force; Conversion of gf to N,
1279 *   (2.0*parameters.half_length_y)*(2.0*parameters.half_length_z)) ) // Area
1280 *   {
1281 *   Assert(dim==3, ExcNotImplemented());
1282 *   }
1283 *  
1284 *  
1285 * @endcode
1286 *
1287 *
1288 * <a name="LinearMuscleModelProblemLinearMuscleModelProblem"></a>
1289 * <h4>LinearMuscleModelProblem::~LinearMuscleModelProblem</h4>
1290 *
1291
1292 *
1293 *
1294 * @code
1295 *   template <int dim>
1296 *   LinearMuscleModelProblem<dim>::~LinearMuscleModelProblem ()
1297 *   {
1298 *   dof_handler.clear ();
1299 *   }
1300 *  
1301 *  
1302 * @endcode
1303 *
1304 *
1305 * <a name="LinearMuscleModelProblemmake_grid"></a>
1306 * <h4>LinearMuscleModelProblem::make_grid</h4>
1307 *
1308
1309 *
1310 *
1311 * @code
1312 *   template<int dim>
1313 *   struct BicepsGeometry
1314 *   {
1315 *   BicepsGeometry(const double axial_length,
1316 *   const double radius_ins_orig,
1317 *   const double radius_mid)
1318 *   :
1319 *   ax_lgth (axial_length),
1320 *   r_ins_orig (radius_ins_orig),
1321 *   r_mid (radius_mid)
1322 *   {}
1323 *  
1324 * @endcode
1325 *
1326 * The radial profile of the muscle
1327 * This provides the new coordinates for points @p pt
1328 * on a cylinder of radius r_ins_orig and length
1329 * ax_lgth to be moved to in order to create the
1330 * physiologically representative geometry of
1331 * the muscle
1332 *
1333 * @code
1334 *   Point<dim> profile (const Point<dim> &pt_0) const
1335 *   {
1336 *   Assert(pt_0[0] > -1e-6,
1337 *   ExcMessage("All points must have x-coordinate > 0"));
1338 *  
1339 *   const double r_scale = get_radial_scaling_factor(pt_0[0]);
1340 *   return pt_0 + Point<dim>(0.0, r_scale*pt_0[1], r_scale*pt_0[2]);
1341 *   }
1342 *  
1343 *   Point<dim> operator() (const Point<dim> &pt) const
1344 *   {
1345 *   return profile(pt);
1346 *   }
1347 *  
1348 * @endcode
1349 *
1350 * Provides the muscle direction at the point @p pt
1351 * in the real geometry (one that has undergone the
1352 * transformation given by the profile() function)
1353 * and subequent grid rescaling.
1354 * The directions are given by the gradient of the
1355 * transformation function (i.e. the fibres are
1356 * orientated by the curvature of the muscle).
1357 *
1358
1359 *
1360 * So, being lazy, we transform the current point back
1361 * to the original point on the completely unscaled
1362 * cylindrical grid. We then evaluate the transformation
1363 * at two points (axially displaced) very close to the
1364 * point of interest. The normalised vector joining the
1365 * transformed counterparts of the perturbed points is
1366 * the gradient of the transformation function and,
1367 * thus, defines the fibre direction.
1368 *
1369 * @code
1370 *   Tensor<1,dim> direction (const Point<dim> &pt_scaled,
1371 *   const double &grid_scale) const
1372 *   {
1373 *   const Point<dim> pt = (1.0/grid_scale)*pt_scaled;
1374 *   const Point<dim> pt_0 = inv_profile(pt);
1375 *  
1376 *   static const double eps = 1e-6;
1377 *   const Point<dim> pt_0_eps_p = pt_0 + Point<dim>(+eps,0,0);
1378 *   const Point<dim> pt_0_eps_m = pt_0 + Point<dim>(-eps,0,0);
1379 *   const Point<dim> pt_eps_p = profile(pt_0_eps_p);
1380 *   const Point<dim> pt_eps_m = profile(pt_0_eps_m);
1381 *  
1382 *   static const double tol = 1e-9;
1383 *   (void)tol;
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 *  
1440 *   GridGenerator::hyper_rectangle (triangulation, p1, p2);
1441 *  
1442 *   typename Triangulation<dim>::active_cell_iterator cell =
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);
1469 *   for (typename Triangulation<2>::active_cell_iterator
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(Tensor<1,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 *   Functions::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 *   Functions::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 *   Functions::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 *   Functions::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 *   Functions::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,
2355 *   update_values | update_gradients | update_quadrature_points);
2356 *   unsigned int cell_no = 0;
2357 *   unsigned int fibre_no = 0;
2358 *   for (typename DoFHandler<dim>::active_cell_iterator
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;
2490 *   types::global_dof_index dof_of_interest = numbers::invalid_dof_index;
2491 *  
2492 *   for (typename DoFHandler<dim>::active_cell_iterator
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>
2555 * <h4>LinearMuscleModelProblem::run</h4>
2556 *
2557
2558 *
2559 *
2560 * @code
2561 *   template <int dim>
2562 *   void LinearMuscleModelProblem<dim>::run ()
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*/
FE_Q
Definition fe_q.h:551
Function::vector_value_list
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< RangeNumberType > > &values) const
Function::vector_value
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
Point
Definition point.h:112
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
level
unsigned int level
Definition grid_out.cc:4618
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1703
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition dof_handler.h:441
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:76
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:82
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition fe_update_flags.h:128
DerivativeApproximation::internal::approximate
void approximate(SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Definition derivative_approximation.cc:921
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)
GridGenerator::hyper_rectangle
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
GridGenerator::cylinder
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
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={})
GridGenerator::flatten_triangulation
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:2202
GridTools::transform
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:58
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
Physics::Elasticity::Kinematics::E
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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 Tensor< 1, 3, Number > &axis, const Number &angle)
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition trilinos_sparse_matrix.cc:65
WorkStream::internal::tbb_no_coloring::run
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition work_stream.h:472
internal::QGaussLobatto::gamma
long double gamma(const unsigned int n)
Definition quadrature_lib.cc:104
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:13826
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition p4est_wrappers.cc:339
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:71
numbers::invalid_dof_index
const types::global_dof_index invalid_dof_index
Definition types.h:233
parallel::transform
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Function &function, const unsigned int grainsize)
Definition parallel.h:148
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6024
std::cos
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:5858
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:5946
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition vectorization.h:5962
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6008
triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition p4est_wrappers.cc:69
tria
const ::Triangulation< dim, spacedim > & tria
Definition tria_description.cc:704
Generated by doxygen 1.9.7