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