Reference documentation for deal.II version 9.6.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
Swift-Hohenberg-Solver.h
Go to the documentation of this file.
1
378 *  
379 *  
380 *   #include <deal.II/base/utilities.h>
381 *   #include <deal.II/base/quadrature_lib.h>
382 *   #include <deal.II/base/function.h>
383 *   #include <deal.II/base/logstream.h>
384 *   #include <deal.II/lac/vector.h>
385 *   #include <deal.II/lac/full_matrix.h>
386 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
387 *   #include <deal.II/lac/sparse_matrix.h>
388 *   #include <deal.II/lac/solver_cg.h>
389 *   #include <deal.II/lac/precondition.h>
390 *   #include <deal.II/lac/affine_constraints.h>
391 *   #include <deal.II/grid/tria.h>
392 *   #include <deal.II/grid/grid_generator.h>
393 *   #include <deal.II/grid/grid_refinement.h>
394 *   #include <deal.II/grid/grid_out.h>
395 *   #include <deal.II/dofs/dof_handler.h>
396 *   #include <deal.II/dofs/dof_tools.h>
397 *   #include <deal.II/fe/fe_q.h>
398 *   #include <deal.II/fe/fe_system.h>
399 *   #include <deal.II/fe/fe_values.h>
400 *   #include <deal.II/numerics/data_out.h>
401 *   #include <deal.II/numerics/vector_tools.h>
402 *   #include <deal.II/numerics/error_estimator.h>
403 *   #include <deal.II/numerics/solution_transfer.h>
404 *   #include <deal.II/numerics/matrix_tools.h>
405 *   #include <deal.II/lac/sparse_direct.h>
406 *   #include <deal.II/base/timer.h>
407 *  
408 *   #include <deal.II/grid/manifold_lib.h>
409 *   #include <deal.II/grid/grid_tools.h>
410 *  
411 *   #include <boost/math/special_functions/ellint_1.hpp>
412 *  
413 *   #include <fstream>
414 *   #include <iostream>
415 *   #include <random>
416 *  
417 *   namespace SwiftHohenbergSolver
418 *   {
419 *   using namespace dealii;
420 *  
421 *  
422 *  
423 *  
428 *   enum MeshType {HYPERCUBE, CYLINDER, SPHERE, TORUS, SINUSOID};
429 *  
430 *  
431 *  
437 *   enum InitialConditionType {HOTSPOT, PSUEDORANDOM, RANDOM};
438 *  
439 *  
440 *  
441 *  
442 *  
449 *   template<int spacedim>
450 *   Point<spacedim> transform_function(const Point<spacedim>&p)
451 *   {
452 * @endcode
453 *
454 * Currently this only works for a 3-dimensional embedding space
455 * because we are explicitly referencing the x, y, and z coordinates
456 *
457 * @code
458 *   Assert(spacedim == 3, ExcNotImplemented());
459 *  
460 * @endcode
461 *
462 * Returns a point where the x-coordinate is unchanged but the y and z coordinates are adjusted
463 * by a cos wave of period 20, amplitude .5, and vertical shift 1
464 *
465 * @code
466 *   return Point<spacedim>(p(0), p(1)*(1 + .5*std::cos((3.14159/10)*p(0))), p(2)*(1 + .5*std::cos((3.14159/10)*p(0))));
467 *   }
468 *  
469 *  
470 *  
471 *  
472 *  
480 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
481 *   class SHEquation
482 *   {
483 *   public:
484 * @endcode
485 *
486 * / @brief Default constructor, initializes all variables and objects with default values
487 *
488 * @code
489 *   SHEquation();
490 *  
491 *  
492 *  
501 *   SHEquation(const unsigned int degree
502 *   , double time_step_denominator
503 *   , unsigned int ref_num
504 *   , double r_constant = 0.5
505 *   , double g1_constant = 0.5
506 *   , std::string output_file_name = "solution-"
507 *   , double end_time = 0.5);
508 *   void run();
509 *  
510 *   private:
511 *   void setup_system();
512 *   void solve_time_step();
513 *   void output_results() const;
514 *  
517 *   void make_grid();
518 *  
519 *  
522 *   void make_cylinder();
523 * @endcode
524 *
525 * / @brief Uses the same process as creating a cylinder, but then also warps the boundary of the cylinder by the function (1 + 0.5*cos(pi*x/10))
526 *
527 * @code
528 *   void make_sinusoid();
529 * @endcode
530 *
531 * / @brief Generates a spherical mesh of radius 6*pi using GridGenerator and refines it refinement_number times.
532 *
533 * @code
534 *   void make_sphere();
535 * @endcode
536 *
537 * / @brief Generates a torus mesh with inner radius 4 and outer radius 9 using GridGenerator and refines it refinement_number times.
538 *
539 * @code
540 *   void make_torus();
541 * @endcode
542 *
543 * / @brief Generates a hypercube mesh with sidelenth 12*pi using GridGenerator and refines it refinement_number times.
544 *
545 * @code
546 *   void make_hypercube();
547 *  
548 *  
549 * @endcode
550 *
551 * / @brief The degree of finite element to be used, default 1
552 *
553 * @code
554 *   const unsigned int degree;
555 *  
556 * @endcode
557 *
558 * / @brief Object holding the mesh
559 *
560 * @code
562 *  
566 * @endcode
567 *
568 * / @brief Object which understands which finite elements are at each node
569 *
570 * @code
571 *   DoFHandler<dim, spacedim> dof_handler;
572 *  
573 * @endcode
574 *
575 * / @brief Describes the sparsity of the system matrix, allows for more efficient storage
576 *
577 * @code
578 *   SparsityPattern sparsity_pattern;
579 *  
580 * @endcode
581 *
582 * / @brief Object holding the system matrix, stored as a sparse matrix
583 *
584 * @code
585 *   SparseMatrix<double> system_matrix;
586 *  
587 *  
590 *   Vector<double> solution;
591 * @endcode
592 *
593 * / @brief Stores the solution from the previous timestep. Used to compute non-linear terms
594 *
595 * @code
596 *   Vector<double> old_solution;
597 *  
600 *   Vector<double> system_rhs;
601 *  
602 * @endcode
603 *
604 * / @brief Stores the current time, in the units of the problem
605 *
606 * @code
607 *   double time;
608 * @endcode
609 *
610 * / @brief The amount time is increased each iteration/ the denominator of the discretized time derivative
611 *
612 * @code
613 *   double time_step;
614 * @endcode
615 *
616 * / @brief Counts the number of iterations that have elapsed
617 *
618 * @code
619 *   unsigned int timestep_number;
620 * @endcode
621 *
622 * / @brief Used to compute the time_step: time_step = 1/timestep_denominator
623 *
624 * @code
625 *   unsigned int timestep_denominator;
626 * @endcode
627 *
628 * / @brief Determines how much to globally refine each mesh
629 *
630 * @code
631 *   unsigned int refinement_number;
632 *  
633 * @endcode
634 *
635 * / @brief Coefficient of the linear term in the SH equation. This is often taken to be constant and g_1 allowed to vary
636 *
637 * @code
638 *   const double r;
639 * @endcode
640 *
641 * / @brief Coefficient of the quadratic term in the SH equation. Determines whether hexagonal lattices can form
642 *
643 * @code
644 *   const double g1;
645 * @endcode
646 *
647 * / @brief A control parameter for the cubic term. Can be useful for testing, in this code we let k=1 in all cases
648 *
649 * @code
650 *   const double k;
651 *  
652 * @endcode
653 *
654 * / @brief Name used to create output file. Should not include extension
655 *
656 * @code
657 *   const std::string output_file_name;
658 *  
659 * @endcode
660 *
661 * / @brief Determines when the solver terminates, endtime of ~100 are useful to see equilibrium results
662 *
663 * @code
664 *   const double end_time;
665 *   };
666 *  
667 *  
668 *  
673 *   template <int spacedim>
674 *   class BoundaryValues : public Function<spacedim>
675 *   {
676 *   public:
677 *   BoundaryValues()
678 *   : Function<spacedim>(2)
679 *   {}
680 *  
681 *   virtual double value(const Point<spacedim> & p,
682 *   const unsigned int component = 0) const override;
683 *   };
684 *  
685 *  
686 *  
687 *  
694 *   template <int spacedim>
695 *   double BoundaryValues<spacedim>::value(const Point<spacedim> & p,
696 *   const unsigned int component) const
697 *   {
698 *   (void)component;
699 *   AssertIndexRange(component, 2);
700 *  
701 *   return 0.;
702 *   }
703 *  
704 *  
717 *   template<int spacedim, MeshType MESH, InitialConditionType ICTYPE>
718 *   class InitialCondition : public Function<spacedim>
719 *   {
720 *   private:
721 * @endcode
722 *
723 * / @brief The value of the parameter r, used to determine a bound for the magnitude of the initial conditions
724 *
725 * @code
726 *   const double r;
727 * @endcode
728 *
729 * / @brief A center point, used to determine the location of the hot spot for the HotSpot initial condition
730 *
731 * @code
733 * @endcode
734 *
735 * / @brief Radius of the hot spot
736 *
737 * @code
738 *   double radius;
739 * @endcode
740 *
741 * / @brief Stores the randomly generated coefficients for planar sine waves along the x-axis, used for psuedorandom initial conditions
742 *
743 * @code
744 *   double x_sin_coefficients[10];
745 * @endcode
746 *
747 * / @brief Stores the randomly generated coefficients for planar sine waves along the y-axis, used for psuedorandom initial conditions
748 *
749 * @code
750 *   double y_sin_coefficients[10];
751 *  
752 *   public:
753 *  
756 *   InitialCondition()
757 *   : Function<spacedim>(2),
758 *   r(0.5),
759 *   radius(.5)
760 *   {
761 *   for(int i = 0; i < 10; ++i){
762 *   x_sin_coefficients[i] = 2*std::sqrt(r)*(std::rand()%1001)/1000 - std::sqrt(r);
763 *   y_sin_coefficients[i] = 2*std::sqrt(r)*(std::rand()%1001)/1000 - std::sqrt(r);
764 *   }
765 *   }
766 *  
767 *  
772 *   InitialCondition(const double r,
773 *   const double radius)
774 *   : Function<spacedim>(2),
775 *   r(r),
776 *   radius(radius)
777 *   {
778 *   for(int i = 0; i < 10; ++i){
779 *   x_sin_coefficients[i] = 2*std::sqrt(r)*(std::rand()%1001)/1000 - std::sqrt(r);
780 *   y_sin_coefficients[i] = 2*std::sqrt(r)*(std::rand()%1001)/1000 - std::sqrt(r);
781 *   }
782 *   }
783 *  
784 *  
804 *   virtual double value(const Point<spacedim> &p, const unsigned int component) const override;
805 *   };
806 *  
807 *  
812 *   template <>
813 *   double InitialCondition<2, HYPERCUBE, HOTSPOT>::value(
814 *   const Point<2> &p,
815 *   const unsigned int component) const
816 *   {
817 *   if(component == 0){
818 *   if(p.square() <= radius){
819 *   return std::sqrt(r);
820 *   }
821 *   else{
822 *   return -std::sqrt(r);
823 *   }
824 *   }
825 *   else{
826 *   return 1e18;
827 *   }
828 *   }
829 *  
830 *  
835 *   template <>
836 *   double InitialCondition<3, CYLINDER, HOTSPOT>::value(
837 *   const Point<3> &p,
838 *   const unsigned int component) const
839 *   {
840 *   if(component == 0){
841 *   const Point<3> center(0, 0, 6);
842 *   const Point<3> compare(p - center);
843 *   if(compare.square() <= radius){
844 *   return std::sqrt(r);
845 *   }
846 *   else{
847 *   return -std::sqrt(r);
848 *   }
849 *   }
850 *   else{
851 *   return 1e18;
852 *   }
853 *   }
854 *  
855 *  
860 *   template <>
861 *   double InitialCondition<3, SPHERE, HOTSPOT>::value(
862 *   const Point<3> &p,
863 *   const unsigned int component) const
864 *   {
865 *   if(component == 0){
866 *   const Point<3> center(18.41988074, 0, 0);
867 *   const Point<3> compare(p - center);
868 *   if(compare.square() <= radius){
869 *   return std::sqrt(r);
870 *   }
871 *   else{
872 *   return -std::sqrt(r);
873 *   }
874 *   }
875 *   else{
876 *   return 1e18;
877 *   }
878 *   }
879 *  
880 *  
885 *   template <>
886 *   double InitialCondition<3, TORUS, HOTSPOT>::value(
887 *   const Point<3> &p,
888 *   const unsigned int component) const
889 *   {
890 *   if(component == 0){
891 *   const Point<3> center(13., 0, 0);
892 *   const Point<3> compare(p - center);
893 *   if(compare.square() <= radius){
894 *   return std::sqrt(r);
895 *   }
896 *   else{
897 *   return -std::sqrt(r);
898 *   }
899 *   }
900 *   else{
901 *   return 1e18;
902 *   }
903 *   }
904 *  
905 *  
910 *   template <>
911 *   double InitialCondition<3, SINUSOID, HOTSPOT>::value(
912 *   const Point<3> &p,
913 *   const unsigned int component) const
914 *   {
915 *   if(component == 0){
916 *   const Point<3> center(0, 0, 9.);
917 *   const Point<3> compare(p - center);
918 *   if(compare.square() <= radius){
919 *   return std::sqrt(r);
920 *   }
921 *   else{
922 *   return -std::sqrt(r);
923 *   }
924 *   }
925 *   else{
926 *   return 1e18;
927 *   }
928 *   }
929 *  
930 *  
935 *   template <>
936 *   double InitialCondition<2, HYPERCUBE, PSUEDORANDOM>::value(
937 *   const Point<2> &p,
938 *   const unsigned int component) const
939 *   {
940 *   if(component == 0){
941 *   double x_val = 0;
942 *   double y_val = 0;
943 *   for(int i=0; i < 10; ++i){
944 *   x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
945 *   y_val += y_sin_coefficients[i]*std::sin(2*3.141592653*p(1)/((i+1)*1.178097245));
946 *   }
947 *  
948 *   return x_val*y_val;
949 *   }
950 *   else{
951 *   return 1e18;
952 *   }
953 *   }
954 *  
955 *  
960 *   template <>
961 *   double InitialCondition<3, CYLINDER, PSUEDORANDOM>::value(
962 *   const Point<3> &p,
963 *   const unsigned int component) const
964 *   {
965 *   if(component == 0){
966 *   double x_val = 0;
967 *   double w_val = 0;
968 *   double width = ((std::atan2(p(1),p(2)) - 3.1415926)/3.1415926)*18.84955592;
969 *   for(int i=0; i < 10; ++i){
970 *   x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
971 *   w_val += y_sin_coefficients[i]*std::sin(2*3.141592653*width/((i+1)*1.178097245));
972 *   }
973 *  
974 *   return x_val*w_val;
975 *   }
976 *   else{
977 *   return 1e18;
978 *   }
979 *   }
980 *  
981 *  
986 *   template <>
987 *   double InitialCondition<3, SPHERE, PSUEDORANDOM>::value(
988 *   const Point<3> &p,
989 *   const unsigned int component) const
990 *   {
991 *   if(component == 0){
992 *   double x_val = 0;
993 *   double y_val = 0;
994 *   for(int i=0; i < 10; ++i){
995 *   x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
996 *   y_val += y_sin_coefficients[i]*std::sin(2*3.141592653*p(1)/((i+1)*1.178097245));
997 *   }
998 *  
999 *   return x_val*y_val;
1000 *   }
1001 *   else{
1002 *   return 1e18;
1003 *   }
1004 *   }
1005 *  
1006 *  
1011 *   template <>
1012 *   double InitialCondition<3, TORUS, PSUEDORANDOM>::value(
1013 *   const Point<3> &p,
1014 *   const unsigned int component) const
1015 *   {
1016 *   if(component == 0){
1017 *   double x_val = 0;
1018 *   double z_val = 0;
1019 *   for(int i=0; i < 10; ++i){
1020 *   x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
1021 *   z_val += y_sin_coefficients[i]*std::sin(2*3.141592653*p(2)/((i+1)*1.178097245));
1022 *   }
1023 *  
1024 *   return x_val*z_val;
1025 *   }
1026 *   else{
1027 *   return 1e18;
1028 *   }
1029 *   }
1030 *  
1031 *  
1036 *   template <>
1037 *   double InitialCondition<3, SINUSOID, PSUEDORANDOM>::value(
1038 *   const Point<3> &p,
1039 *   const unsigned int component) const
1040 *   {
1041 *   if(component == 0){
1042 *   double x_val = 0;
1043 *   double w_val = 0;
1044 *   double width = ((std::atan2(p(1),p(2)) - 3.1415926)/3.1415926)*18.84955592;
1045 *   for(int i=0; i < 10; ++i){
1046 *   x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
1047 *   w_val += y_sin_coefficients[i]*std::sin(2*3.141592653*width/((i+1)*1.178097245));
1048 *   }
1049 *  
1050 *   return x_val*w_val;
1051 *   }
1052 *   else{
1053 *   return 1e18;
1054 *   }
1055 *   }
1056 *  
1057 *  
1062 *   template <>
1063 *   double InitialCondition<2, HYPERCUBE, RANDOM>::value(
1064 *   const Point<2> &/*p*/,
1065 *   const unsigned int component) const
1066 *   {
1067 *   if(component == 0){
1068 *   return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
1069 *   }
1070 *   else{
1071 *   return 1e18;
1072 *   }
1073 *   }
1074 *  
1075 *  
1080 *   template <>
1081 *   double InitialCondition<3, CYLINDER, RANDOM>::value(
1082 *   const Point<3> &/*p*/,
1083 *   const unsigned int component) const
1084 *   {
1085 *   if(component == 0){
1086 *   return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
1087 *   }
1088 *   else{
1089 *   return 1e18;
1090 *   }
1091 *   }
1092 *  
1093 *  
1098 *   template <>
1099 *   double InitialCondition<3, SPHERE, RANDOM>::value(
1100 *   const Point<3> &/*p*/,
1101 *   const unsigned int component) const
1102 *   {
1103 *   if(component == 0){
1104 *   return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
1105 *   }
1106 *   else{
1107 *   return 1e18;
1108 *   }
1109 *   }
1110 *  
1111 *  
1116 *   template <>
1117 *   double InitialCondition<3, TORUS, RANDOM>::value(
1118 *   const Point<3> &/*p*/,
1119 *   const unsigned int component) const
1120 *   {
1121 *   if(component == 0){
1122 *   return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
1123 *   }
1124 *   else{
1125 *   return 1e18;
1126 *   }
1127 *   }
1128 *  
1129 *  
1134 *   template <>
1135 *   double InitialCondition<3, SINUSOID, RANDOM>::value(
1136 *   const Point<3> &/*p*/,
1137 *   const unsigned int component) const
1138 *   {
1139 *   if(component == 0){
1140 *   return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
1141 *   }
1142 *   else{
1143 *   return 1e18;
1144 *   }
1145 *   }
1146 *  
1147 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1148 *   SHEquation<dim, spacedim, MESH, ICTYPE>::SHEquation()
1149 *   : degree(1)
1150 *   , fe(FE_Q<dim, spacedim>(degree), 2)
1151 *   , dof_handler(triangulation)
1152 *   , time_step(1. / 1500)
1153 *   , timestep_denominator(1500)
1154 *   , refinement_number(4)
1155 *   , r(0.5)
1156 *   , g1(0.5)
1157 *   , k(1.)
1158 *   , output_file_name("solution-")
1159 *   , end_time(0.5)
1160 *   {}
1161 *  
1162 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1163 *   SHEquation<dim, spacedim, MESH, ICTYPE>::SHEquation(const unsigned int degree,
1164 *   double time_step_denominator,
1165 *   unsigned int ref_num,
1166 *   double r_constant,
1167 *   double g1_constant,
1168 *   std::string output_file_name,
1169 *   double end_time)
1170 *   : degree(degree)
1171 *   , fe(FE_Q<dim, spacedim>(degree), 2)
1172 *   , dof_handler(triangulation)
1173 *   , time_step(1. / time_step_denominator)
1174 *   , timestep_denominator(time_step_denominator)
1175 *   , refinement_number(ref_num)
1176 *   , r(r_constant)
1177 *   , g1(g1_constant)
1178 *   , k(1.)
1179 *   , output_file_name(output_file_name)
1180 *   , end_time(end_time)
1181 *   {}
1182 *  
1183 *  
1190 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1191 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::setup_system()
1192 *   {
1193 *   dof_handler.distribute_dofs(fe);
1194 *  
1195 * @endcode
1196 *
1197 * Counts the DoF's for outputting to consolse
1198 *
1199 * @code
1200 *   const std::vector<types::global_dof_index> dofs_per_component =
1201 *   DoFTools::count_dofs_per_fe_component(dof_handler);
1202 *   const unsigned int n_u = dofs_per_component[0],
1203 *   n_v = dofs_per_component[1];
1204 *  
1205 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
1206 *   << std::endl
1207 *   << "Total number of cells: " << triangulation.n_cells()
1208 *   << std::endl
1209 *   << "Number of degrees of freedom: " << dof_handler.n_dofs()
1210 *   << " (" << n_u << '+' << n_v << ')' << std::endl;
1211 *  
1212 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
1213 *  
1214 *   DoFTools::make_sparsity_pattern(dof_handler,
1215 *   dsp);
1216 *   sparsity_pattern.copy_from(dsp);
1217 *  
1218 *   system_matrix.reinit(sparsity_pattern);
1219 *  
1220 *   solution.reinit(dof_handler.n_dofs());
1221 *   old_solution.reinit(dof_handler.n_dofs());
1222 *   system_rhs.reinit(dof_handler.n_dofs());
1223 *   }
1224 *  
1225 *  
1226 *   /** @brief Uses a direct solver to invert the system matrix, then multiplies the RHS vector by the inverted matrix to get the solution.
1227 *   * Also includes a timer feature, which is currently commented out, but can be helpful to compute how long a run will take
1228 *   * @tparam dim The dimension of the manifold
1229 *   * @tparam spacedim The dimension of the ambient space
1230 *   * @tparam MESH The type of mesh being used, doesn't change how this function works
1231 *   * @tparam ICTYPE The type of initial condition used, doesn't change how this function works
1232 *   */
1233 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1234 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::solve_time_step()
1235 *   {
1236 * @endcode
1237 *
1238 * std::cout << "Solving linear system" << std::endl;
1239 * Timer timer;
1240 *
1241
1242 *
1243 *
1244 * @code
1245 *   SparseDirectUMFPACK direct_solver;
1246 *  
1247 *   direct_solver.initialize(system_matrix);
1248 *  
1249 *   direct_solver.vmult(solution, system_rhs);
1250 *  
1251 * @endcode
1252 *
1253 * timer.stop();
1254 * std::cout << "done (" << timer.cpu_time() << " s)" << std::endl;
1255 *
1256 * @code
1257 *   }
1258 *  
1259 *  
1260 *  
1261 *   /** @brief Converts the solution vector into a .vtu file and labels the outputs as u and v
1262 *   * @tparam dim The dimension of the manifold
1263 *   * @tparam spacedim The dimension of the ambient space
1264 *   * @tparam MESH The type of mesh being used, doesn't change how this function works
1265 *   * @tparam ICTYPE The type of initial condition used, doesn't change how this function works
1266 *   */
1267 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1268 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::output_results() const
1269 *   {
1270 *   std::vector<std::string> solution_names(1, "u");
1271 *   solution_names.emplace_back("v");
1272 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
1273 *   interpretation(1,
1274 *   DataComponentInterpretation::component_is_scalar);
1275 *   interpretation.push_back(DataComponentInterpretation::component_is_scalar);
1276 *  
1277 *   DataOut<dim, spacedim> data_out;
1278 *   data_out.add_data_vector(dof_handler,
1279 *   solution,
1280 *   solution_names,
1281 *   interpretation /*,
1282 *   DataOut<dim, spacedim>::type_dof_data*/);
1283 *  
1284 *   data_out.build_patches(degree + 1);
1285 *  
1286 * @endcode
1287 *
1288 * Takes the output_file_name string and appends timestep_number with up to three leading 0's
1289 *
1290 * @code
1291 *   const std::string filename =
1292 *   output_file_name + Utilities::int_to_string(timestep_number, 3) + ".vtu";
1293 *  
1294 *   std::ofstream output(filename);
1295 *   data_out.write_vtu(output);
1296 *   }
1297 *  
1298 * @endcode
1299 *
1300 * Below are all the different template cases for the make_grid() function
1301 *
1302 * @code
1303 *   template <>
1304 *   void SHEquation<2, 2, HYPERCUBE, HOTSPOT>::make_grid()
1305 *   {
1306 *   make_hypercube();
1307 *   }
1308 *  
1309 *   template <>
1310 *   void SHEquation<2, 3, CYLINDER, HOTSPOT>::make_grid()
1311 *   {
1312 *   make_cylinder();
1313 *   }
1314 *  
1315 *   template <>
1316 *   void SHEquation<2, 3, SPHERE, HOTSPOT>::make_grid()
1317 *   {
1318 *   make_sphere();
1319 *   }
1320 *  
1321 *   template <>
1322 *   void SHEquation<2, 3, TORUS, HOTSPOT>::make_grid()
1323 *   {
1324 *   make_torus();
1325 *   }
1326 *  
1327 *   template <>
1328 *   void SHEquation<2, 3, SINUSOID, HOTSPOT>::make_grid()
1329 *   {
1330 *   make_sinusoid();
1331 *   }
1332 *  
1333 *   template <>
1334 *   void SHEquation<2, 2, HYPERCUBE, PSUEDORANDOM>::make_grid()
1335 *   {
1336 *   make_hypercube();
1337 *   }
1338 *  
1339 *   template <>
1340 *   void SHEquation<2, 3, CYLINDER, PSUEDORANDOM>::make_grid()
1341 *   {
1342 *   make_cylinder();
1343 *   }
1344 *  
1345 *   template <>
1346 *   void SHEquation<2, 3, SPHERE, PSUEDORANDOM>::make_grid()
1347 *   {
1348 *   make_sphere();
1349 *   }
1350 *  
1351 *   template <>
1352 *   void SHEquation<2, 3, TORUS, PSUEDORANDOM>::make_grid()
1353 *   {
1354 *   make_torus();
1355 *   }
1356 *  
1357 *   template <>
1358 *   void SHEquation<2, 3, SINUSOID, PSUEDORANDOM>::make_grid()
1359 *   {
1360 *   make_sinusoid();
1361 *   }
1362 *  
1363 *   template <>
1364 *   void SHEquation<2, 2, HYPERCUBE, RANDOM>::make_grid()
1365 *   {
1366 *   make_hypercube();
1367 *   }
1368 *  
1369 *   template <>
1370 *   void SHEquation<2, 3, CYLINDER, RANDOM>::make_grid()
1371 *   {
1372 *   make_cylinder();
1373 *   }
1374 *  
1375 *   template <>
1376 *   void SHEquation<2, 3, SPHERE, RANDOM>::make_grid()
1377 *   {
1378 *   make_sphere();
1379 *   }
1380 *  
1381 *   template <>
1382 *   void SHEquation<2, 3, TORUS, RANDOM>::make_grid()
1383 *   {
1384 *   make_torus();
1385 *   }
1386 *  
1387 *   template <>
1388 *   void SHEquation<2, 3, SINUSOID, RANDOM>::make_grid()
1389 *   {
1390 *   make_sinusoid();
1391 *   }
1392 *  
1393 *  
1394 *  
1401 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1402 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::run()
1403 *   {
1404 *   make_grid();
1405 *  
1406 *   setup_system();
1407 *  
1408 * @endcode
1409 *
1410 * Counts total time elapsed
1411 *
1412 * @code
1413 *   time = 0.0;
1414 * @endcode
1415 *
1416 * Counts number of iterations
1417 *
1418 * @code
1419 *   timestep_number = 0;
1420 *  
1421 * @endcode
1422 *
1423 * Sets the random seed so runs are repeatable, remove for varying random initial conditions
1424 *
1425 * @code
1426 *   std::srand(314);
1427 *  
1428 *   InitialCondition<spacedim, MESH, ICTYPE> initial_conditions(r, 0.5);
1429 *  
1430 * @endcode
1431 *
1432 * Applies the initial conditions to the old_solution
1433 *
1434 * @code
1435 *   VectorTools::interpolate(dof_handler,
1436 *   initial_conditions,
1437 *   old_solution);
1438 *   solution = old_solution;
1439 *  
1440 * @endcode
1441 *
1442 * Outputs initial solution
1443 *
1444 * @code
1445 *   output_results();
1446 *  
1447 * @endcode
1448 *
1449 * Sets up the quadrature formula and FEValues object
1450 *
1451 * @code
1452 *   const QGauss<dim> quadrature_formula(degree + 2);
1453 *  
1454 *   FEValues<dim, spacedim> fe_values(fe, quadrature_formula,
1457 *  
1458 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1459 *  
1460 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1461 *   Vector<double> cell_rhs(dofs_per_cell);
1462 *  
1463 * @endcode
1464 *
1465 * The vector which stores the global indices that each local index connects to
1466 *
1467 * @code
1468 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1469 *  
1470 * @endcode
1471 *
1472 * Extracts the finite elements associated to u and v
1473 *
1474 * @code
1475 *   const FEValuesExtractors::Scalar u(0);
1476 *   const FEValuesExtractors::Scalar v(1);
1477 *  
1478 * @endcode
1479 *
1480 * Loops over the cells to create the system matrix. We do this only once because the timestep is constant
1481 *
1482 * @code
1483 *   for(const auto &cell : dof_handler.active_cell_iterators()){
1484 *   cell_matrix = 0;
1485 *   cell_rhs = 0;
1486 *  
1487 *   fe_values.reinit(cell);
1488 *  
1489 *   cell->get_dof_indices(local_dof_indices);
1490 *  
1491 *   for(const unsigned int q_index : fe_values.quadrature_point_indices()){
1492 *  
1493 *   for(const unsigned int i : fe_values.dof_indices()){
1494 * @endcode
1495 *
1496 * These are the ith finite elements associated to u and v
1497 *
1498 * @code
1499 *   const double phi_i_u = fe_values[u].value(i, q_index);
1500 *   const Tensor<1, spacedim> grad_phi_i_u = fe_values[u].gradient(i, q_index);
1501 *   const double phi_i_v = fe_values[v].value(i, q_index);
1502 *   const Tensor<1, spacedim> grad_phi_i_v = fe_values[v].gradient(i, q_index);
1503 *  
1504 *   for(const unsigned int j : fe_values.dof_indices())
1505 *   {
1506 * @endcode
1507 *
1508 * These are the jth finite elements associated to u and v
1509 *
1510 * @code
1511 *   const double phi_j_u = fe_values[u].value(j, q_index);
1512 *   const Tensor<1, spacedim> grad_phi_j_u = fe_values[u].gradient(j, q_index);
1513 *   const double phi_j_v = fe_values[v].value(j, q_index);
1514 *   const Tensor<1, spacedim> grad_phi_j_v = fe_values[v].gradient(j, q_index);
1515 *  
1516 * @endcode
1517 *
1518 * This formula comes from expanding the PDE system
1519 *
1520 * @code
1521 *   cell_matrix(i, j) += (phi_i_u*phi_j_u - time_step*r*phi_i_u*phi_j_u
1522 *   + time_step*phi_i_u*phi_j_v - time_step*grad_phi_i_u*grad_phi_j_v
1523 *   + phi_i_v*phi_j_u - grad_phi_i_v*grad_phi_j_u
1524 *   - phi_i_v*phi_j_v)*fe_values.JxW(q_index);
1525 *   }
1526 *   }
1527 *   }
1528 *  
1529 * @endcode
1530 *
1531 * Loops over the dof indices to fill the entries of the system_matrix with the local data
1532 *
1533 * @code
1534 *   for(unsigned int i : fe_values.dof_indices()){
1535 *   for(unsigned int j : fe_values.dof_indices()){
1536 *   system_matrix.add(local_dof_indices[i],
1537 *   local_dof_indices[j],
1538 *   cell_matrix(i, j));
1539 *   }
1540 *   }
1541 *   }
1542 *  
1543 * @endcode
1544 *
1545 * Loops over time, incrementing by timestep, to create the RHS, solve the linear system, then output the result
1546 *
1547 * @code
1548 *   while (time <= end_time)
1549 *   {
1550 * @endcode
1551 *
1552 * Increments time and timestep_number
1553 *
1554 * @code
1555 *   time += time_step;
1556 *   ++timestep_number;
1557 *  
1558 * @endcode
1559 *
1560 * Outputs to console the number of iterations and current time. Currently outputs once every "second"
1561 *
1562 * @code
1563 *   if(timestep_number%timestep_denominator == 0){
1564 *   std::cout << "Time step " << timestep_number << " at t=" << time
1565 *   << std::endl;
1566 *   }
1567 *  
1568 * @endcode
1569 *
1570 * Resets the system_rhs vector. THIS IS VERY IMPORTANT TO ENSURE THE SYSTEM IS SOLVED CORRECTLY AT EACH TIMESTEP
1571 *
1572 * @code
1573 *   system_rhs = 0;
1574 *  
1575 * @endcode
1576 *
1577 * Loops over cells, then quadrature points, then dof indices to construct the RHS
1578 *
1579 * @code
1580 *   for(const auto &cell : dof_handler.active_cell_iterators()){
1581 * @endcode
1582 *
1583 * Resets the cell_rhs. THIS IS ALSO VERY IMPORTANT TO ENSURE THE SYSTEM IS SOLVED CORRECTLY
1584 *
1585 * @code
1586 *   cell_rhs = 0;
1587 *  
1588 * @endcode
1589 *
1590 * Resets the FEValues object to only the current cell
1591 *
1592 * @code
1593 *   fe_values.reinit(cell);
1594 *  
1595 *   cell->get_dof_indices(local_dof_indices);
1596 *  
1597 * @endcode
1598 *
1599 * Loop over the quadrature points
1600 *
1601 * @code
1602 *   for(const unsigned int q_index : fe_values.quadrature_point_indices()){
1603 * @endcode
1604 *
1605 * Stores the value of the previous solution at the quadrature point
1606 *
1607 * @code
1608 *   double Un1 = 0;
1609 *  
1610 * @endcode
1611 *
1612 * Loops over the dof indices to get the value of Un1
1613 *
1614 * @code
1615 *   for(const unsigned int i : fe_values.dof_indices()){
1616 *   Un1 += old_solution(local_dof_indices[i])*fe_values[u].value(i, q_index);
1617 *   }
1618 *  
1619 * @endcode
1620 *
1621 * Loops over the dof indices, using Un1 to construct the RHS for the current timestep. Un1 is used to account for the nonlinear terms in the SH equation
1622 *
1623 * @code
1624 *   for(const unsigned int i : fe_values.dof_indices()){
1625 *   cell_rhs(i) += (Un1 + time_step*g1*std::pow(Un1, 2) - time_step*k*std::pow(Un1, 3))
1626 *   *fe_values[u].value(i, q_index)*fe_values.JxW(q_index);
1627 *   }
1628 *   }
1629 *  
1630 * @endcode
1631 *
1632 * Loops over the dof indices to store the local data in the global RHS vector
1633 *
1634 * @code
1635 *   for(unsigned int i : fe_values.dof_indices()){
1636 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
1637 *   }
1638 *  
1639 *  
1640 *   }
1641 * @endcode
1642 *
1643 * This is where Dirichlet conditions are applied, or Neumann conditions if the code is commented out
1644 *
1645 * @code
1646 *   /* {
1647 *   BoundaryValues<spacedim> boundary_values_function;
1648 *   boundary_values_function.set_time(time);
1649 *  
1650 *   std::map<types::global_dof_index, double> boundary_values;
1651 *   VectorTools::interpolate_boundary_values(dof_handler,
1652 *   0,
1653 *   boundary_values_function,
1654 *   boundary_values);
1655 *  
1656 *   MatrixTools::apply_boundary_values(boundary_values,
1657 *   system_matrix,
1658 *   solution,
1659 *   system_rhs);
1660 *   } */
1661 *  
1662 *   solve_time_step();
1663 *  
1664 * @endcode
1665 *
1666 * Outputs the solution at regular intervals, currently once every "second" The SH equation evolves slowly in time, so this saves disk space
1667 *
1668 * @code
1669 *   if(timestep_number%timestep_denominator == 0){
1670 *   output_results();
1671 *   }
1672 *  
1673 *   old_solution = solution;
1674 *   }
1675 *   }
1676 *  
1677 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1678 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::make_cylinder()
1679 *   {
1680 * @endcode
1681 *
1682 * Creates a volumetric cylinder
1683 *
1684 * @code
1686 *   GridGenerator::cylinder(cylinder, 6, 18.84955592);
1687 *  
1688 * @endcode
1689 *
1690 * Extracts the boundary mesh with ID 0, which happens to be the tube part of the cylinder
1691 *
1692 * @code
1694 *  
1695 * @endcode
1696 *
1697 * The manifold information is lost upon boundary extraction. This sets the mesh boundary type to be a cylinder again
1698 *
1699 * @code
1700 *   const CylindricalManifold<dim, spacedim> boundary;
1702 *   triangulation.set_manifold(0, boundary);
1703 *  
1704 *   triangulation.refine_global(refinement_number);
1705 *   }
1706 *  
1707 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1708 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::make_sinusoid()
1709 *   {
1710 * @endcode
1711 *
1712 * Same process as above
1713 *
1714 * @code
1716 *   GridGenerator::cylinder(cylinder, 6, 18.84955592);
1717 *  
1719 *  
1720 *   const CylindricalManifold<dim, spacedim> boundary;
1722 *   triangulation.set_manifold(0, boundary);
1723 *  
1724 *   triangulation.refine_global(refinement_number);
1725 *  
1726 * @endcode
1727 *
1728 * We warp the mesh after refinement to avoid a jagged mesh. We can't tell the code that the boundary should be a perfect sine wave, so we only warp after the
1729 * mesh is fine enough to resolve this
1730 *
1731 * @code
1732 *   GridTools::transform(transform_function<spacedim>, triangulation);
1733 *   }
1734 *  
1735 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1736 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::make_sphere()
1737 *   {
1738 *   GridGenerator::hyper_sphere(triangulation, Point<3>(0, 0, 0), 18.41988074);
1739 *   triangulation.refine_global(refinement_number);
1740 *   }
1741 *  
1742 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1743 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::make_torus()
1744 *   {
1745 *   GridGenerator::torus(triangulation, 9., 4.);
1746 *   triangulation.refine_global(refinement_number);
1747 *   }
1748 *   template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1749 *   void SHEquation<dim, spacedim, MESH, ICTYPE>::make_hypercube()
1750 *   {
1751 *   GridGenerator::hyper_cube(triangulation, -18.84955592, 18.84955592);
1752 *   triangulation.refine_global(refinement_number);
1753 *   }
1754 *   } // namespace SwiftHohenbergSolver
1755 *  
1756 *  
1757 *  
1758 *   int main()
1759 *   {
1760 *   using namespace SwiftHohenbergSolver;
1761 *  
1762 * @endcode
1763 *
1764 * An array of mesh types. We iterate over this to allow for longer runs without having to stop the code
1765 *
1766 * @code
1767 *   MeshType mesh_types[5] = {HYPERCUBE, CYLINDER, SPHERE, TORUS, SINUSOID};
1768 * @endcode
1769 *
1770 * An array of initial condition types. We iterate this as well, for the same reason
1771 *
1772 * @code
1773 *   InitialConditionType ic_types[3] = {HOTSPOT, PSUEDORANDOM, RANDOM};
1774 *  
1775 * @endcode
1776 *
1777 * Controls how long the code runs
1778 *
1779 * @code
1780 *   const double end_time = 100.;
1781 *  
1782 * @endcode
1783 *
1784 * The number of times we refine the hypercube mesh
1785 *
1786 * @code
1787 *   const unsigned int ref_num = 6;
1788 *  
1789 * @endcode
1790 *
1791 * The timestep will be 1/timestep_denominator
1792 *
1793 * @code
1794 *   const unsigned int timestep_denominator = 25;
1795 *  
1796 * @endcode
1797 *
1798 * Loops over mesh types, then initial condition types, then loops over values of g_1
1799 *
1800 * @code
1801 *   for(const auto MESH : mesh_types){
1802 *   for(const auto ICTYPE: ic_types){
1803 *   for(int i = 0; i < 8; ++i){
1804 * @endcode
1805 *
1806 * The value of g_1 passed to the solver object
1807 *
1808 * @code
1809 *   const double g_constant = 0.2*i;
1810 *  
1811 * @endcode
1812 *
1813 * Used to distinguish the start of each run
1814 *
1815 * @code
1816 *   std::cout<< std::endl << std::endl;
1817 *  
1818 *   try{
1819 * @endcode
1820 *
1821 * Switch statement that determines what template parameters are used by the solver object. Template parameters must be known at compile time, so we cannot
1822 * pass this as a variable unfortunately. In each case, we create a filename string (named appropriately for the particular case), output to the console what
1823 * we are running, create the solver object, and call run(). Note that for the cylinder, sphere, and sinusoid we decrease the refinement number by 1. This keeps
1824 * the number of dofs used in these cases comparable to the number of dofs on the 2D hypercube (otherwise the number of dofs is much larger). For the torus, we
1825 * decrease the refinement number by 2.
1826 *
1827 * @code
1828 *   switch (MESH)
1829 *   {
1830 *   case HYPERCUBE:
1831 *   switch (ICTYPE){
1832 *   case HOTSPOT:
1833 *   {
1834 *   std::string filename = "HYPERCUBE-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1835 *   std::cout << "Running: " << filename << std::endl << std::endl;
1836 *  
1837 *   SHEquation<2, 2, HYPERCUBE, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1838 *   ref_num, 0.3, g_constant,
1839 *   filename, end_time);
1840 *   heat_equation_solver.run();
1841 *   }
1842 *   break;
1843 *  
1844 *   case PSUEDORANDOM:
1845 *   {
1846 *   std::string filename = "HYPERCUBE-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1847 *   std::cout << "Running: " << filename << std::endl << std::endl;
1848 *  
1849 *   SHEquation<2, 2, HYPERCUBE, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
1850 *   ref_num, 0.3, g_constant,
1851 *   filename, end_time);
1852 *   heat_equation_solver.run();
1853 *   }
1854 *   break;
1855 *  
1856 *   case RANDOM:
1857 *   {
1858 *   std::string filename = "HYPERCUBE-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1859 *   std::cout << "Running: " << filename << std::endl << std::endl;
1860 *  
1861 *   SHEquation<2, 2, HYPERCUBE, RANDOM> heat_equation_solver(1, timestep_denominator,
1862 *   ref_num, 0.3, g_constant,
1863 *   filename, end_time);
1864 *   heat_equation_solver.run();
1865 *   }
1866 *   break;
1867 *   }
1868 *   break;
1869 *   case CYLINDER:
1870 *   switch (ICTYPE){
1871 *   case HOTSPOT:
1872 *   {
1873 *   std::string filename = "CYLINDER-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1874 *   std::cout << "Running: " << filename << std::endl << std::endl;
1875 *  
1876 *   SHEquation<2, 3, CYLINDER, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1877 *   ref_num-1, 0.3, g_constant,
1878 *   filename, end_time);
1879 *   heat_equation_solver.run();
1880 *   }
1881 *   break;
1882 *  
1883 *   case PSUEDORANDOM:
1884 *   {
1885 *   std::string filename = "CYLINDER-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1886 *   std::cout << "Running: " << filename << std::endl << std::endl;
1887 *  
1888 *   SHEquation<2, 3, CYLINDER, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
1889 *   ref_num-1, 0.3, g_constant,
1890 *   filename, end_time);
1891 *   heat_equation_solver.run();
1892 *   }
1893 *   break;
1894 *  
1895 *   case RANDOM:
1896 *   {
1897 *   std::string filename = "CYLINDER-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1898 *   std::cout << "Running: " << filename << std::endl << std::endl;
1899 *  
1900 *   SHEquation<2, 3, CYLINDER, RANDOM> heat_equation_solver(1, timestep_denominator,
1901 *   ref_num-1, 0.3, g_constant,
1902 *   filename, end_time);
1903 *   heat_equation_solver.run();
1904 *   }
1905 *   break;
1906 *   }
1907 *   break;
1908 *   case SPHERE:
1909 *   switch (ICTYPE){
1910 *   case HOTSPOT:
1911 *   {
1912 *   std::string filename = "SPHERE-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1913 *   std::cout << "Running: " << filename << std::endl << std::endl;
1914 *  
1915 *   SHEquation<2, 3, SPHERE, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1916 *   ref_num-1, 0.3, g_constant,
1917 *   filename, end_time);
1918 *   heat_equation_solver.run();
1919 *   }
1920 *   break;
1921 *  
1922 *   case PSUEDORANDOM:
1923 *   {
1924 *   std::string filename = "SPHERE-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1925 *   std::cout << "Running: " << filename << std::endl << std::endl;
1926 *  
1927 *   SHEquation<2, 3, SPHERE, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
1928 *   ref_num-1, 0.3, g_constant,
1929 *   filename, end_time);
1930 *   heat_equation_solver.run();
1931 *   }
1932 *   break;
1933 *  
1934 *   case RANDOM:
1935 *   {
1936 *   std::string filename = "SPHERE-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1937 *   std::cout << "Running: " << filename << std::endl << std::endl;
1938 *  
1939 *   SHEquation<2, 3, SPHERE, RANDOM> heat_equation_solver(1, timestep_denominator,
1940 *   ref_num-1, 0.3, g_constant,
1941 *   filename, end_time);
1942 *   heat_equation_solver.run();
1943 *   }
1944 *   break;
1945 *   }
1946 *   break;
1947 *   case TORUS:
1948 *   switch (ICTYPE){
1949 *   case HOTSPOT:
1950 *   {
1951 *   std::string filename = "TORUS-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1952 *   std::cout << "Running: " << filename << std::endl << std::endl;
1953 *  
1954 *   SHEquation<2, 3, TORUS, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1955 *   ref_num-2, 0.3, g_constant,
1956 *   filename, end_time);
1957 *   heat_equation_solver.run();
1958 *   }
1959 *   break;
1960 *  
1961 *   case PSUEDORANDOM:
1962 *   {
1963 *   std::string filename = "TORUS-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1964 *   std::cout << "Running: " << filename << std::endl << std::endl;
1965 *  
1966 *   SHEquation<2, 3, TORUS, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
1967 *   ref_num-2, 0.3, g_constant,
1968 *   filename, end_time);
1969 *   heat_equation_solver.run();
1970 *   }
1971 *   break;
1972 *  
1973 *   case RANDOM:
1974 *   {
1975 *   std::string filename = "TORUS-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1976 *   std::cout << "Running: " << filename << std::endl << std::endl;
1977 *  
1978 *   SHEquation<2, 3, TORUS, RANDOM> heat_equation_solver(1, timestep_denominator,
1979 *   ref_num-2, 0.3, g_constant,
1980 *   filename, end_time);
1981 *   heat_equation_solver.run();
1982 *   }
1983 *   break;
1984 *   }
1985 *   break;
1986 *   case SINUSOID:
1987 *   switch (ICTYPE){
1988 *   case HOTSPOT:
1989 *   {
1990 *   std::string filename = "SINUSOID-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1991 *   std::cout << "Running: " << filename << std::endl << std::endl;
1992 *  
1993 *   SHEquation<2, 3, SINUSOID, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1994 *   ref_num-1, 0.3, g_constant,
1995 *   filename, end_time);
1996 *   heat_equation_solver.run();
1997 *   }
1998 *   break;
1999 *  
2000 *   case PSUEDORANDOM:
2001 *   {
2002 *   std::string filename = "SINUSOID-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
2003 *   std::cout << "Running: " << filename << std::endl << std::endl;
2004 *  
2005 *   SHEquation<2, 3, SINUSOID, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
2006 *   ref_num-1, 0.3, g_constant,
2007 *   filename, end_time);
2008 *   heat_equation_solver.run();
2009 *   }
2010 *   break;
2011 *  
2012 *   case RANDOM:
2013 *   {
2014 *   std::string filename = "SINUSOID-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
2015 *   std::cout << "Running: " << filename << std::endl << std::endl;
2016 *  
2017 *   SHEquation<2, 3, SINUSOID, RANDOM> heat_equation_solver(1, timestep_denominator,
2018 *   ref_num-1, 0.3, g_constant,
2019 *   filename, end_time);
2020 *   heat_equation_solver.run();
2021 *   }
2022 *   break;
2023 *   }
2024 *   break;
2025 *   default:
2026 *   break;
2027 *   }
2028 *   }
2029 *   catch (std::exception &exc)
2030 *   {
2031 *   std::cout << "An error occurred" << std::endl;
2032 *   std::cerr << std::endl
2033 *   << std::endl
2034 *   << "----------------------------------------------------"
2035 *   << std::endl;
2036 *   std::cerr << "Exception on processing: " << std::endl
2037 *   << exc.what() << std::endl
2038 *   << "Aborting!" << std::endl
2039 *   << "----------------------------------------------------"
2040 *   << std::endl;
2041 *  
2042 *   return 1;
2043 *   }
2044 *   catch (...)
2045 *   {
2046 *   std::cout << "Error occurred, made it past first catch" << std::endl;
2047 *   std::cerr << std::endl
2048 *   << std::endl
2049 *   << "----------------------------------------------------"
2050 *   << std::endl;
2051 *   std::cerr << "Unknown exception!" << std::endl
2052 *   << "Aborting!" << std::endl
2053 *   << "----------------------------------------------------"
2054 *   << std::endl;
2055 *   return 1;
2056 *   }
2057 *   }
2058 *   }
2059 *   }
2060 *   return 0;
2061 *   }
2062 * @endcode
2063
2064
2065*/
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
Definition fe_q.h:554
Function(const unsigned int n_components=1, const time_type initial_time=0.0)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:111
void refine_global(const unsigned int times=1)
Point< 3 > center
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
void set_all_manifold_ids(const types::manifold_id number)
const Event initial
Definition event.cc:64
void random(DoFHandler< dim, spacedim > &dof_handler)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
return_type extract_boundary_mesh(const MeshType< dim, spacedim > &volume_mesh, MeshType< dim - 1, spacedim > &surface_mesh, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
void torus(Triangulation< dim, spacedim > &tria, const double centerline_radius, const double inner_radius, const unsigned int n_cells_toroidal=6, const double phi=2.0 *numbers::PI)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
@ matrix
Contents is actually a matrix.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={})
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)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation