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