deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-30.h
Go to the documentation of this file.
1 = 0) const override
602 *   {
603 *   (void)points;
604 *   AssertDimension(values.size(), points.size());
605 *  
606 *   std::fill(values.begin(), values.end(), 0.);
607 *   }
608 *   };
609 *  
610 *  
611 *   template <int dim>
612 *   class BoundaryValues : public Function<dim>
613 *   {
614 *   public:
615 *   virtual void value_list(const std::vector<Point<dim>> &points,
616 *   std::vector<double> &values,
617 *   const unsigned int /*component*/ = 0) const override
618 *   {
619 *   AssertDimension(values.size(), points.size());
620 *  
621 *   for (unsigned int i = 0; i < values.size(); ++i)
622 *   {
623 *   if (points[i][0] < 0.5)
624 *   values[i] = 1.;
625 *   else
626 *   values[i] = 0.;
627 *   }
628 *   }
629 *   };
630 *  
631 *  
632 *   template <int dim>
633 *   class Beta
634 *   {
635 *   public:
636 * @endcode
637 *
638 * The flow field is chosen to be a quarter circle with counterclockwise
639 * flow direction and with the origin as midpoint for the right half of the
640 * domain with positive @f$x@f$ values, whereas the flow simply goes to the left
641 * in the left part of the domain at a velocity that matches the one coming
642 * in from the right. In the circular part the magnitude of the flow
643 * velocity is proportional to the distance from the origin. This is a
644 * difference to @ref step_12 "step-12", where the magnitude was 1 everywhere. the new
645 * definition leads to a linear variation of @f$\beta@f$ along each given face
646 * of a cell. On the other hand, the solution @f$u(x,y)@f$ is exactly the same
647 * as before.
648 *
649 * @code
650 *   void value_list(const std::vector<Point<dim>> &points,
651 *   std::vector<Point<dim>> &values) const
652 *   {
653 *   AssertDimension(values.size(), points.size());
654 *  
655 *   for (unsigned int i = 0; i < points.size(); ++i)
656 *   {
657 *   if (points[i][0] > 0)
658 *   {
659 *   values[i][0] = -points[i][1];
660 *   values[i][1] = points[i][0];
661 *   }
662 *   else
663 *   {
664 *   values[i] = Point<dim>();
665 *   values[i][0] = -points[i][1];
666 *   }
667 *   }
668 *   }
669 *   };
670 *  
671 *  
672 *  
673 * @endcode
674 *
675 *
676 * <a name="step_30-ClassDGTransportEquation"></a>
677 * <h3>Class: DGTransportEquation</h3>
678 *
679
680 *
681 * This declaration of this class is utterly unaffected by our current
682 * changes.
683 *
684 * @code
685 *   template <int dim>
686 *   class DGTransportEquation
687 *   {
688 *   public:
689 *   DGTransportEquation();
690 *  
691 *   void assemble_cell_term(const FEValues<dim> &fe_v,
692 *   FullMatrix<double> &ui_vi_matrix,
693 *   Vector<double> &cell_vector) const;
694 *  
695 *   void assemble_boundary_term(const FEFaceValues<dim> &fe_v,
696 *   FullMatrix<double> &ui_vi_matrix,
697 *   Vector<double> &cell_vector) const;
698 *  
699 *   void assemble_face_term(const FEFaceValuesBase<dim> &fe_v,
700 *   const FEFaceValuesBase<dim> &fe_v_neighbor,
701 *   FullMatrix<double> &ui_vi_matrix,
702 *   FullMatrix<double> &ue_vi_matrix,
703 *   FullMatrix<double> &ui_ve_matrix,
704 *   FullMatrix<double> &ue_ve_matrix) const;
705 *  
706 *   private:
707 *   const Beta<dim> beta_function;
708 *   const RHS<dim> rhs_function;
709 *   const BoundaryValues<dim> boundary_function;
710 *   };
711 *  
712 *  
713 *  
714 * @endcode
715 *
716 * Likewise, the constructor of the class as well as the functions
717 * assembling the terms corresponding to cell interiors and boundary faces
718 * are unchanged from before. The function that assembles face terms between
719 * cells also did not change because all it does is operate on two objects
720 * of type FEFaceValuesBase (which is the base class of both FEFaceValues
721 * and FESubfaceValues). Where these objects come from, i.e. how they are
722 * initialized, is of no concern to this function: it simply assumes that
723 * the quadrature points on faces or subfaces represented by the two objects
724 * correspond to the same points in physical space.
725 *
726 * @code
727 *   template <int dim>
728 *   DGTransportEquation<dim>::DGTransportEquation()
729 *   : beta_function()
730 *   , rhs_function()
731 *   , boundary_function()
732 *   {}
733 *  
734 *  
735 *  
736 *   template <int dim>
737 *   void DGTransportEquation<dim>::assemble_cell_term(
738 *   const FEValues<dim> &fe_v,
739 *   FullMatrix<double> &ui_vi_matrix,
740 *   Vector<double> &cell_vector) const
741 *   {
742 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
743 *  
744 *   std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
745 *   std::vector<double> rhs(fe_v.n_quadrature_points);
746 *  
747 *   beta_function.value_list(fe_v.get_quadrature_points(), beta);
748 *   rhs_function.value_list(fe_v.get_quadrature_points(), rhs);
749 *  
750 *   for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
751 *   for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
752 *   {
753 *   for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
754 *   ui_vi_matrix(i, j) -= beta[point] * fe_v.shape_grad(i, point) *
755 *   fe_v.shape_value(j, point) * JxW[point];
756 *  
757 *   cell_vector(i) +=
758 *   rhs[point] * fe_v.shape_value(i, point) * JxW[point];
759 *   }
760 *   }
761 *  
762 *  
763 *   template <int dim>
764 *   void DGTransportEquation<dim>::assemble_boundary_term(
765 *   const FEFaceValues<dim> &fe_v,
766 *   FullMatrix<double> &ui_vi_matrix,
767 *   Vector<double> &cell_vector) const
768 *   {
769 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
770 *   const std::vector<Tensor<1, dim>> &normals = fe_v.get_normal_vectors();
771 *  
772 *   std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
773 *   std::vector<double> g(fe_v.n_quadrature_points);
774 *  
775 *   beta_function.value_list(fe_v.get_quadrature_points(), beta);
776 *   boundary_function.value_list(fe_v.get_quadrature_points(), g);
777 *  
778 *   for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
779 *   {
780 *   const double beta_n = beta[point] * normals[point];
781 *   if (beta_n > 0)
782 *   for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
783 *   for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
784 *   ui_vi_matrix(i, j) += beta_n * fe_v.shape_value(j, point) *
785 *   fe_v.shape_value(i, point) * JxW[point];
786 *   else
787 *   for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
788 *   cell_vector(i) -=
789 *   beta_n * g[point] * fe_v.shape_value(i, point) * JxW[point];
790 *   }
791 *   }
792 *  
793 *  
794 *   template <int dim>
795 *   void DGTransportEquation<dim>::assemble_face_term(
796 *   const FEFaceValuesBase<dim> &fe_v,
797 *   const FEFaceValuesBase<dim> &fe_v_neighbor,
798 *   FullMatrix<double> &ui_vi_matrix,
799 *   FullMatrix<double> &ue_vi_matrix,
800 *   FullMatrix<double> &ui_ve_matrix,
801 *   FullMatrix<double> &ue_ve_matrix) const
802 *   {
803 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
804 *   const std::vector<Tensor<1, dim>> &normals = fe_v.get_normal_vectors();
805 *  
806 *   std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
807 *  
808 *   beta_function.value_list(fe_v.get_quadrature_points(), beta);
809 *  
810 *   for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
811 *   {
812 *   const double beta_n = beta[point] * normals[point];
813 *   if (beta_n > 0)
814 *   {
815 *   for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
816 *   for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
817 *   ui_vi_matrix(i, j) += beta_n * fe_v.shape_value(j, point) *
818 *   fe_v.shape_value(i, point) * JxW[point];
819 *  
820 *   for (unsigned int k = 0; k < fe_v_neighbor.dofs_per_cell; ++k)
821 *   for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
822 *   ui_ve_matrix(k, j) -= beta_n * fe_v.shape_value(j, point) *
823 *   fe_v_neighbor.shape_value(k, point) *
824 *   JxW[point];
825 *   }
826 *   else
827 *   {
828 *   for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
829 *   for (unsigned int l = 0; l < fe_v_neighbor.dofs_per_cell; ++l)
830 *   ue_vi_matrix(i, l) += beta_n *
831 *   fe_v_neighbor.shape_value(l, point) *
832 *   fe_v.shape_value(i, point) * JxW[point];
833 *  
834 *   for (unsigned int k = 0; k < fe_v_neighbor.dofs_per_cell; ++k)
835 *   for (unsigned int l = 0; l < fe_v_neighbor.dofs_per_cell; ++l)
836 *   ue_ve_matrix(k, l) -=
837 *   beta_n * fe_v_neighbor.shape_value(l, point) *
838 *   fe_v_neighbor.shape_value(k, point) * JxW[point];
839 *   }
840 *   }
841 *   }
842 *  
843 *  
844 * @endcode
845 *
846 *
847 * <a name="step_30-ClassDGMethod"></a>
848 * <h3>Class: DGMethod</h3>
849 *
850
851 *
852 * This declaration is much like that of @ref step_12 "step-12". However, we introduce a
853 * new routine (set_anisotropic_flags) and modify another one (refine_grid).
854 *
855 * @code
856 *   template <int dim>
857 *   class DGMethod
858 *   {
859 *   public:
860 *   DGMethod(const bool anisotropic);
861 *  
862 *   void run();
863 *  
864 *   private:
865 *   void setup_system();
866 *   void assemble_system();
867 *   void solve(Vector<double> &solution);
868 *   void refine_grid();
869 *   void set_anisotropic_flags();
870 *   void output_results(const unsigned int cycle) const;
871 *  
873 *   const MappingQ1<dim> mapping;
874 * @endcode
875 *
876 * Again we want to use DG elements of degree 1 (but this is only
877 * specified in the constructor). If you want to use a DG method of a
878 * different degree replace 1 in the constructor by the new degree.
879 *
880 * @code
881 *   const unsigned int degree;
882 *   const FE_DGQ<dim> fe;
883 *   DoFHandler<dim> dof_handler;
884 *  
885 *   SparsityPattern sparsity_pattern;
886 *   SparseMatrix<double> system_matrix;
887 * @endcode
888 *
889 * This is new, the threshold value used in the evaluation of the
890 * anisotropic jump indicator explained in the introduction. Its value is
891 * set to 3.0 in the constructor, but it can easily be changed to a
892 * different value greater than 1.
893 *
894 * @code
895 *   const double anisotropic_threshold_ratio;
896 * @endcode
897 *
898 * This is a bool flag indicating whether anisotropic refinement shall be
899 * used or not. It is set by the constructor, which takes an argument of
900 * the same name.
901 *
902 * @code
903 *   const bool anisotropic;
904 *  
905 *   const QGauss<dim> quadrature;
906 *   const QGauss<dim - 1> face_quadrature;
907 *  
908 *   Vector<double> solution2;
909 *   Vector<double> right_hand_side;
910 *  
911 *   const DGTransportEquation<dim> dg;
912 *   };
913 *  
914 *  
915 *   template <int dim>
916 *   DGMethod<dim>::DGMethod(const bool anisotropic)
917 *   : mapping()
918 * @endcode
919 *
920 * Change here for DG methods of different degrees.
921 *
922 * @code
923 *   , degree(1)
924 *   , fe(degree)
925 *   , dof_handler(triangulation)
926 *   , anisotropic_threshold_ratio(3.)
927 *   , anisotropic(anisotropic)
928 * @endcode
929 *
930 * As beta is a linear function, we can choose the degree of the
931 * quadrature for which the resulting integration is correct. Thus, we
932 * choose to use <code>degree+1</code> Gauss points, which enables us to
933 * integrate exactly polynomials of degree <code>2*degree+1</code>, enough
934 * for all the integrals we will perform in this program.
935 *
936 * @code
937 *   , quadrature(degree + 1)
938 *   , face_quadrature(degree + 1)
939 *   , dg()
940 *   {}
941 *  
942 *  
943 *  
944 *   template <int dim>
945 *   void DGMethod<dim>::setup_system()
946 *   {
947 *   dof_handler.distribute_dofs(fe);
948 *   sparsity_pattern.reinit(dof_handler.n_dofs(),
949 *   dof_handler.n_dofs(),
952 *   1) *
953 *   fe.n_dofs_per_cell());
954 *  
955 *   DoFTools::make_flux_sparsity_pattern(dof_handler, sparsity_pattern);
956 *  
957 *   sparsity_pattern.compress();
958 *  
959 *   system_matrix.reinit(sparsity_pattern);
960 *  
961 *   solution2.reinit(dof_handler.n_dofs());
962 *   right_hand_side.reinit(dof_handler.n_dofs());
963 *   }
964 *  
965 *  
966 * @endcode
967 *
968 *
969 * <a name="step_30-Functionassemble_system"></a>
970 * <h4>Function: assemble_system</h4>
971 *
972
973 *
974 * We proceed with the <code>assemble_system</code> function that implements
975 * the DG discretization. This function does the same thing as the
976 * <code>assemble_system</code> function from @ref step_12 "step-12" (but without
977 * MeshWorker). The four cases considered for the neighbor-relations of a
978 * cell are the same as the isotropic case, namely a) cell is at the
979 * boundary, b) there are finer neighboring cells, c) the neighbor is
980 * neither coarser nor finer and d) the neighbor is coarser. However, the
981 * way in which we decide upon which case we have are modified in the way
982 * described in the introduction.
983 *
984 * @code
985 *   template <int dim>
986 *   void DGMethod<dim>::assemble_system()
987 *   {
988 *   const unsigned int dofs_per_cell = dof_handler.get_fe().n_dofs_per_cell();
989 *   std::vector<types::global_dof_index> dofs(dofs_per_cell);
990 *   std::vector<types::global_dof_index> dofs_neighbor(dofs_per_cell);
991 *  
992 *   const UpdateFlags update_flags = update_values | update_gradients |
995 *  
996 *   const UpdateFlags face_update_flags =
999 *  
1000 *   const UpdateFlags neighbor_face_update_flags = update_values;
1001 *  
1002 *   FEValues<dim> fe_v(mapping, fe, quadrature, update_flags);
1003 *   FEFaceValues<dim> fe_v_face(mapping,
1004 *   fe,
1005 *   face_quadrature,
1006 *   face_update_flags);
1007 *   FESubfaceValues<dim> fe_v_subface(mapping,
1008 *   fe,
1009 *   face_quadrature,
1010 *   face_update_flags);
1011 *   FEFaceValues<dim> fe_v_face_neighbor(mapping,
1012 *   fe,
1013 *   face_quadrature,
1014 *   neighbor_face_update_flags);
1015 *  
1016 *  
1017 *   FullMatrix<double> ui_vi_matrix(dofs_per_cell, dofs_per_cell);
1018 *   FullMatrix<double> ue_vi_matrix(dofs_per_cell, dofs_per_cell);
1019 *  
1020 *   FullMatrix<double> ui_ve_matrix(dofs_per_cell, dofs_per_cell);
1021 *   FullMatrix<double> ue_ve_matrix(dofs_per_cell, dofs_per_cell);
1022 *  
1023 *   Vector<double> cell_vector(dofs_per_cell);
1024 *  
1025 *   for (const auto &cell : dof_handler.active_cell_iterators())
1026 *   {
1027 *   ui_vi_matrix = 0;
1028 *   cell_vector = 0;
1029 *  
1030 *   fe_v.reinit(cell);
1031 *  
1032 *   dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector);
1033 *  
1034 *   cell->get_dof_indices(dofs);
1035 *  
1036 *   for (const auto face_no : cell->face_indices())
1037 *   {
1038 *   const auto face = cell->face(face_no);
1039 *  
1040 * @endcode
1041 *
1042 * Case (a): The face is at the boundary.
1043 *
1044 * @code
1045 *   if (face->at_boundary())
1046 *   {
1047 *   fe_v_face.reinit(cell, face_no);
1048 *  
1049 *   dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector);
1050 *   }
1051 *   else
1052 *   {
1053 *   Assert(cell->neighbor(face_no).state() == IteratorState::valid,
1054 *   ExcInternalError());
1055 *   const auto neighbor = cell->neighbor(face_no);
1056 *  
1057 * @endcode
1058 *
1059 * Case (b): This is an internal face and the neighbor
1060 * is refined (which we can test by asking whether the
1061 * face of the current cell has children). In this
1062 * case, we will need to integrate over the
1063 * "sub-faces", i.e., the children of the face of the
1064 * current cell.
1065 *
1066
1067 *
1068 * (There is a slightly confusing corner case: If we
1069 * are in 1d -- where admittedly the current program
1070 * and its demonstration of anisotropic refinement is
1071 * not particularly relevant -- then the faces between
1072 * cells are always the same: they are just
1073 * vertices. In other words, in 1d, we do not want to
1074 * treat faces between cells of different level
1075 * differently. The condition `face->has_children()`
1076 * we check here ensures this: in 1d, this function
1077 * always returns `false`, and consequently in 1d we
1078 * will not ever go into this `if` branch. But we will
1079 * have to come back to this corner case below in case
1080 * (c).)
1081 *
1082 * @code
1083 *   if (face->has_children())
1084 *   {
1085 * @endcode
1086 *
1087 * We need to know, which of the neighbors faces points in
1088 * the direction of our cell. Using the @p
1089 * neighbor_face_no function we get this information for
1090 * both coarser and non-coarser neighbors.
1091 *
1092 * @code
1093 *   const unsigned int neighbor2 =
1094 *   cell->neighbor_face_no(face_no);
1095 *  
1096 * @endcode
1097 *
1098 * Now we loop over all subfaces, i.e. the children and
1099 * possibly grandchildren of the current face.
1100 *
1101 * @code
1102 *   for (unsigned int subface_no = 0;
1103 *   subface_no < face->n_active_descendants();
1104 *   ++subface_no)
1105 *   {
1106 * @endcode
1107 *
1108 * To get the cell behind the current subface we can
1109 * use the @p neighbor_child_on_subface function. it
1110 * takes care of all the complicated situations of
1111 * anisotropic refinement and non-standard faces.
1112 *
1113 * @code
1114 *   const auto neighbor_child =
1115 *   cell->neighbor_child_on_subface(face_no, subface_no);
1116 *   Assert(!neighbor_child->has_children(),
1117 *   ExcInternalError());
1118 *  
1119 * @endcode
1120 *
1121 * The remaining part of this case is unchanged.
1122 *
1123 * @code
1124 *   ue_vi_matrix = 0;
1125 *   ui_ve_matrix = 0;
1126 *   ue_ve_matrix = 0;
1127 *  
1128 *   fe_v_subface.reinit(cell, face_no, subface_no);
1129 *   fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
1130 *  
1131 *   dg.assemble_face_term(fe_v_subface,
1132 *   fe_v_face_neighbor,
1133 *   ui_vi_matrix,
1134 *   ue_vi_matrix,
1135 *   ui_ve_matrix,
1136 *   ue_ve_matrix);
1137 *  
1138 *   neighbor_child->get_dof_indices(dofs_neighbor);
1139 *  
1140 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1141 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1142 *   {
1143 *   system_matrix.add(dofs[i],
1144 *   dofs_neighbor[j],
1145 *   ue_vi_matrix(i, j));
1146 *   system_matrix.add(dofs_neighbor[i],
1147 *   dofs[j],
1148 *   ui_ve_matrix(i, j));
1149 *   system_matrix.add(dofs_neighbor[i],
1150 *   dofs_neighbor[j],
1151 *   ue_ve_matrix(i, j));
1152 *   }
1153 *   }
1154 *   }
1155 *   else
1156 *   {
1157 * @endcode
1158 *
1159 * Case (c). We get here if this is an internal
1160 * face and if the neighbor is not further refined
1161 * (or, as mentioned above, we are in 1d in which
1162 * case we get here for every internal face). We
1163 * then need to decide whether we want to
1164 * integrate over the current face. If the
1165 * neighbor is in fact coarser, then we ignore the
1166 * face and instead handle it when we visit the
1167 * neighboring cell and look at the current face
1168 * (except in 1d, where as mentioned above this is
1169 * not happening):
1170 *
1171 * @code
1172 *   if (dim > 1 && cell->neighbor_is_coarser(face_no))
1173 *   continue;
1174 *  
1175 * @endcode
1176 *
1177 * On the other hand, if the neighbor is more
1178 * refined, then we have already handled the face
1179 * in case (b) above (except in 1d). So for 2d and
1180 * 3d, we just have to decide whether we want to
1181 * handle a face between cells at the same level
1182 * from the current side or from the neighboring
1183 * side. We do this by introducing a tie-breaker:
1184 * We'll just take the cell with the smaller index
1185 * (within the current refinement level). In 1d,
1186 * we take either the coarser cell, or if they are
1187 * on the same level, the one with the smaller
1188 * index within that level. This leads to a
1189 * complicated condition that, hopefully, makes
1190 * sense given the description above:
1191 *
1192 * @code
1193 *   if (((dim > 1) && (cell->index() < neighbor->index())) ||
1194 *   ((dim == 1) && ((cell->level() < neighbor->level()) ||
1195 *   ((cell->level() == neighbor->level()) &&
1196 *   (cell->index() < neighbor->index())))))
1197 *   {
1198 * @endcode
1199 *
1200 * Here we know, that the neighbor is not coarser so we
1201 * can use the usual @p neighbor_of_neighbor
1202 * function. However, we could also use the more
1203 * general @p neighbor_face_no function.
1204 *
1205 * @code
1206 *   const unsigned int neighbor2 =
1207 *   cell->neighbor_of_neighbor(face_no);
1208 *  
1209 *   ue_vi_matrix = 0;
1210 *   ui_ve_matrix = 0;
1211 *   ue_ve_matrix = 0;
1212 *  
1213 *   fe_v_face.reinit(cell, face_no);
1214 *   fe_v_face_neighbor.reinit(neighbor, neighbor2);
1215 *  
1216 *   dg.assemble_face_term(fe_v_face,
1217 *   fe_v_face_neighbor,
1218 *   ui_vi_matrix,
1219 *   ue_vi_matrix,
1220 *   ui_ve_matrix,
1221 *   ue_ve_matrix);
1222 *  
1223 *   neighbor->get_dof_indices(dofs_neighbor);
1224 *  
1225 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1226 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1227 *   {
1228 *   system_matrix.add(dofs[i],
1229 *   dofs_neighbor[j],
1230 *   ue_vi_matrix(i, j));
1231 *   system_matrix.add(dofs_neighbor[i],
1232 *   dofs[j],
1233 *   ui_ve_matrix(i, j));
1234 *   system_matrix.add(dofs_neighbor[i],
1235 *   dofs_neighbor[j],
1236 *   ue_ve_matrix(i, j));
1237 *   }
1238 *   }
1239 *  
1240 * @endcode
1241 *
1242 * We do not need to consider a case (d), as those
1243 * faces are treated 'from the other side within
1244 * case (b).
1245 *
1246 * @code
1247 *   }
1248 *   }
1249 *   }
1250 *  
1251 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1252 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1253 *   system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i, j));
1254 *  
1255 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1256 *   right_hand_side(dofs[i]) += cell_vector(i);
1257 *   }
1258 *   }
1259 *  
1260 *  
1261 * @endcode
1262 *
1263 *
1264 * <a name="step_30-Solver"></a>
1265 * <h3>Solver</h3>
1266 *
1267
1268 *
1269 * For this simple problem we use the simple Richardson iteration again. The
1270 * solver is completely unaffected by our anisotropic changes.
1271 *
1272 * @code
1273 *   template <int dim>
1274 *   void DGMethod<dim>::solve(Vector<double> &solution)
1275 *   {
1276 *   SolverControl solver_control(1000, 1e-12, false, false);
1277 *   SolverRichardson<Vector<double>> solver(solver_control);
1278 *  
1279 *   PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
1280 *  
1281 *   preconditioner.initialize(system_matrix, fe.n_dofs_per_cell());
1282 *  
1283 *   solver.solve(system_matrix, solution, right_hand_side, preconditioner);
1284 *   }
1285 *  
1286 *  
1287 * @endcode
1288 *
1289 *
1290 * <a name="step_30-Refinement"></a>
1291 * <h3>Refinement</h3>
1292 *
1293
1294 *
1295 * We refine the grid according to the same simple refinement criterion used
1296 * in @ref step_12 "step-12", namely an approximation to the gradient of the solution.
1297 *
1298 * @code
1299 *   template <int dim>
1300 *   void DGMethod<dim>::refine_grid()
1301 *   {
1302 *   Vector<float> gradient_indicator(triangulation.n_active_cells());
1303 *  
1304 * @endcode
1305 *
1306 * We approximate the gradient,
1307 *
1308 * @code
1310 *   dof_handler,
1311 *   solution2,
1312 *   gradient_indicator);
1313 *  
1314 * @endcode
1315 *
1316 * and scale it to obtain an error indicator.
1317 *
1318 * @code
1319 *   for (const auto &cell : triangulation.active_cell_iterators())
1320 *   gradient_indicator[cell->active_cell_index()] *=
1321 *   std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
1322 * @endcode
1323 *
1324 * Then we use this indicator to flag the 30 percent of the cells with
1325 * highest error indicator to be refined.
1326 *
1327 * @code
1329 *   gradient_indicator,
1330 *   0.3,
1331 *   0.1);
1332 * @endcode
1333 *
1334 * Now the refinement flags are set for those cells with a large error
1335 * indicator. If nothing is done to change this, those cells will be
1336 * refined isotropically. If the @p anisotropic flag given to this
1337 * function is set, we now call the set_anisotropic_flags() function,
1338 * which uses the jump indicator to reset some of the refinement flags to
1339 * anisotropic refinement.
1340 *
1341 * @code
1342 *   if (anisotropic)
1343 *   set_anisotropic_flags();
1344 * @endcode
1345 *
1346 * Now execute the refinement considering anisotropic as well as isotropic
1347 * refinement flags.
1348 *
1349 * @code
1350 *   triangulation.execute_coarsening_and_refinement();
1351 *   }
1352 *  
1353 * @endcode
1354 *
1355 * Once an error indicator has been evaluated and the cells with largest
1356 * error are flagged for refinement we want to loop over the flagged cells
1357 * again to decide whether they need isotropic refinement or whether
1358 * anisotropic refinement is more appropriate. This is the anisotropic jump
1359 * indicator explained in the introduction.
1360 *
1361 * @code
1362 *   template <int dim>
1363 *   void DGMethod<dim>::set_anisotropic_flags()
1364 *   {
1365 * @endcode
1366 *
1367 * We want to evaluate the jump over faces of the flagged cells, so we
1368 * need some objects to evaluate values of the solution on faces.
1369 *
1370 * @code
1371 *   UpdateFlags face_update_flags =
1373 *  
1374 *   FEFaceValues<dim> fe_v_face(mapping,
1375 *   fe,
1376 *   face_quadrature,
1377 *   face_update_flags);
1378 *   FESubfaceValues<dim> fe_v_subface(mapping,
1379 *   fe,
1380 *   face_quadrature,
1381 *   face_update_flags);
1382 *   FEFaceValues<dim> fe_v_face_neighbor(mapping,
1383 *   fe,
1384 *   face_quadrature,
1385 *   update_values);
1386 *  
1387 * @endcode
1388 *
1389 * Now we need to loop over all active cells.
1390 *
1391 * @code
1392 *   for (const auto &cell : dof_handler.active_cell_iterators())
1393 * @endcode
1394 *
1395 * We only need to consider cells which are flagged for refinement.
1396 *
1397 * @code
1398 *   if (cell->refine_flag_set())
1399 *   {
1400 *   std::array<double, dim> jump_in_coordinate_direction;
1401 *   std::array<double, dim> face_area_in_coordinate_direction;
1402 *  
1403 *   for (const auto face_no : cell->face_indices())
1404 *   {
1405 *   const auto face = cell->face(face_no);
1406 *  
1407 *   if (!face->at_boundary())
1408 *   {
1409 *   Assert(cell->neighbor(face_no).state() ==
1411 *   ExcInternalError());
1412 *   const auto neighbor = cell->neighbor(face_no);
1413 *  
1414 *   std::vector<double> u(fe_v_face.n_quadrature_points);
1415 *   std::vector<double> u_neighbor(fe_v_face.n_quadrature_points);
1416 *  
1417 * @endcode
1418 *
1419 * The four cases of different neighbor relations seen in
1420 * the assembly routines are repeated much in the same way
1421 * here.
1422 *
1423 * @code
1424 *   if (face->has_children())
1425 *   {
1426 * @endcode
1427 *
1428 * The neighbor is refined. First we store the
1429 * information, which of the neighbor's faces points in
1430 * the direction of our current cell. This property is
1431 * inherited to the children.
1432 *
1433 * @code
1434 *   unsigned int neighbor2 = cell->neighbor_face_no(face_no);
1435 * @endcode
1436 *
1437 * Now we loop over all subfaces,
1438 *
1439 * @code
1440 *   for (unsigned int subface_no = 0;
1441 *   subface_no < face->n_active_descendants();
1442 *   ++subface_no)
1443 *   {
1444 * @endcode
1445 *
1446 * get an iterator pointing to the cell behind the
1447 * present subface...
1448 *
1449 * @code
1450 *   const auto neighbor_child =
1451 *   cell->neighbor_child_on_subface(face_no,
1452 *   subface_no);
1453 *   Assert(!neighbor_child->has_children(),
1454 *   ExcInternalError());
1455 * @endcode
1456 *
1457 * ... and reinit the respective FEFaceValues and
1458 * FESubFaceValues objects.
1459 *
1460 * @code
1461 *   fe_v_subface.reinit(cell, face_no, subface_no);
1462 *   fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
1463 * @endcode
1464 *
1465 * We obtain the function values
1466 *
1467 * @code
1468 *   fe_v_subface.get_function_values(solution2, u);
1469 *   fe_v_face_neighbor.get_function_values(solution2,
1470 *   u_neighbor);
1471 * @endcode
1472 *
1473 * as well as the quadrature weights, multiplied by
1474 * the Jacobian determinant.
1475 *
1476 * @code
1477 *   const std::vector<double> &JxW =
1478 *   fe_v_subface.get_JxW_values();
1479 * @endcode
1480 *
1481 * Now we loop over all quadrature points
1482 *
1483 * @code
1484 *   for (unsigned int x = 0;
1485 *   x < fe_v_subface.n_quadrature_points;
1486 *   ++x)
1487 *   {
1488 * @endcode
1489 *
1490 * and integrate the absolute value of the jump
1491 * of the solution, i.e. the absolute value of
1492 * the difference between the function value
1493 * seen from the current cell and the
1494 * neighboring cell, respectively. We know, that
1495 * the first two faces are orthogonal to the
1496 * first coordinate direction on the unit cell,
1497 * the second two faces are orthogonal to the
1498 * second coordinate direction and so on, so we
1499 * accumulate these values into vectors with
1500 * <code>dim</code> components.
1501 *
1502 * @code
1503 *   jump_in_coordinate_direction[face_no / 2] +=
1504 *   std::abs(u[x] - u_neighbor[x]) * JxW[x];
1505 * @endcode
1506 *
1507 * We also sum up the scaled weights to obtain
1508 * the measure of the face.
1509 *
1510 * @code
1511 *   face_area_in_coordinate_direction[face_no / 2] +=
1512 *   JxW[x];
1513 *   }
1514 *   }
1515 *   }
1516 *   else
1517 *   {
1518 *   if (!cell->neighbor_is_coarser(face_no))
1519 *   {
1520 * @endcode
1521 *
1522 * Our current cell and the neighbor have the same
1523 * refinement along the face under
1524 * consideration. Apart from that, we do much the
1525 * same as with one of the subcells in the above
1526 * case.
1527 *
1528 * @code
1529 *   unsigned int neighbor2 =
1530 *   cell->neighbor_of_neighbor(face_no);
1531 *  
1532 *   fe_v_face.reinit(cell, face_no);
1533 *   fe_v_face_neighbor.reinit(neighbor, neighbor2);
1534 *  
1535 *   fe_v_face.get_function_values(solution2, u);
1536 *   fe_v_face_neighbor.get_function_values(solution2,
1537 *   u_neighbor);
1538 *  
1539 *   const std::vector<double> &JxW =
1540 *   fe_v_face.get_JxW_values();
1541 *  
1542 *   for (unsigned int x = 0;
1543 *   x < fe_v_face.n_quadrature_points;
1544 *   ++x)
1545 *   {
1546 *   jump_in_coordinate_direction[face_no / 2] +=
1547 *   std::abs(u[x] - u_neighbor[x]) * JxW[x];
1548 *   face_area_in_coordinate_direction[face_no / 2] +=
1549 *   JxW[x];
1550 *   }
1551 *   }
1552 *   else // i.e. neighbor is coarser than cell
1553 *   {
1554 * @endcode
1555 *
1556 * Now the neighbor is actually coarser. This case
1557 * is new, in that it did not occur in the assembly
1558 * routine. Here, we have to consider it, but this
1559 * is not overly complicated. We simply use the @p
1560 * neighbor_of_coarser_neighbor function, which
1561 * again takes care of anisotropic refinement and
1562 * non-standard face orientation by itself.
1563 *
1564 * @code
1565 *   std::pair<unsigned int, unsigned int>
1566 *   neighbor_face_subface =
1567 *   cell->neighbor_of_coarser_neighbor(face_no);
1568 *   Assert(neighbor_face_subface.first < cell->n_faces(),
1569 *   ExcInternalError());
1570 *   Assert(neighbor_face_subface.second <
1571 *   neighbor->face(neighbor_face_subface.first)
1572 *   ->n_active_descendants(),
1573 *   ExcInternalError());
1574 *   Assert(neighbor->neighbor_child_on_subface(
1575 *   neighbor_face_subface.first,
1576 *   neighbor_face_subface.second) == cell,
1577 *   ExcInternalError());
1578 *  
1579 *   fe_v_face.reinit(cell, face_no);
1580 *   fe_v_subface.reinit(neighbor,
1581 *   neighbor_face_subface.first,
1582 *   neighbor_face_subface.second);
1583 *  
1584 *   fe_v_face.get_function_values(solution2, u);
1585 *   fe_v_subface.get_function_values(solution2,
1586 *   u_neighbor);
1587 *  
1588 *   const std::vector<double> &JxW =
1589 *   fe_v_face.get_JxW_values();
1590 *  
1591 *   for (unsigned int x = 0;
1592 *   x < fe_v_face.n_quadrature_points;
1593 *   ++x)
1594 *   {
1595 *   jump_in_coordinate_direction[face_no / 2] +=
1596 *   std::abs(u[x] - u_neighbor[x]) * JxW[x];
1597 *   face_area_in_coordinate_direction[face_no / 2] +=
1598 *   JxW[x];
1599 *   }
1600 *   }
1601 *   }
1602 *   }
1603 *   }
1604 * @endcode
1605 *
1606 * Now we analyze the size of the mean jumps, which we get dividing
1607 * the jumps by the measure of the respective faces.
1608 *
1609 * @code
1610 *   std::array<double, dim> average_jumps;
1611 *   double sum_of_average_jumps = 0.;
1612 *   for (unsigned int i = 0; i < dim; ++i)
1613 *   {
1614 *   average_jumps[i] = jump_in_coordinate_direction[i] /
1615 *   face_area_in_coordinate_direction[i];
1616 *   sum_of_average_jumps += average_jumps[i];
1617 *   }
1618 *  
1619 * @endcode
1620 *
1621 * Now we loop over the <code>dim</code> coordinate directions of
1622 * the unit cell and compare the average jump over the faces
1623 * orthogonal to that direction with the average jumps over faces
1624 * orthogonal to the remaining direction(s). If the first is larger
1625 * than the latter by a given factor, we refine only along hat
1626 * axis. Otherwise we leave the refinement flag unchanged, resulting
1627 * in isotropic refinement.
1628 *
1629 * @code
1630 *   for (unsigned int i = 0; i < dim; ++i)
1631 *   if (average_jumps[i] > anisotropic_threshold_ratio *
1632 *   (sum_of_average_jumps - average_jumps[i]))
1633 *   cell->set_refine_flag(RefinementCase<dim>::cut_axis(i));
1634 *   }
1635 *   }
1636 *  
1637 * @endcode
1638 *
1639 *
1640 * <a name="step_30-TheRest"></a>
1641 * <h3>The Rest</h3>
1642 *
1643
1644 *
1645 * The remaining part of the program very much follows the scheme of
1646 * previous tutorial programs. We output the mesh in VTU format (just
1647 * as we did in @ref step_1 "step-1", for example), and the visualization output
1648 * in VTU format as we almost always do.
1649 *
1650 * @code
1651 *   template <int dim>
1652 *   void DGMethod<dim>::output_results(const unsigned int cycle) const
1653 *   {
1654 *   std::string refine_type;
1655 *   if (anisotropic)
1656 *   refine_type = ".aniso";
1657 *   else
1658 *   refine_type = ".iso";
1659 *  
1660 *   {
1661 *   const std::string filename =
1662 *   "grid-" + std::to_string(cycle) + refine_type + ".svg";
1663 *   std::cout << " Writing grid to <" << filename << ">..." << std::endl;
1664 *   std::ofstream svg_output(filename);
1665 *  
1666 *   GridOut grid_out;
1667 *   grid_out.write_svg(triangulation, svg_output);
1668 *   }
1669 *  
1670 *   {
1671 *   const std::string filename =
1672 *   "sol-" + std::to_string(cycle) + refine_type + ".vtu";
1673 *   std::cout << " Writing solution to <" << filename << ">..."
1674 *   << std::endl;
1675 *   std::ofstream gnuplot_output(filename);
1676 *  
1677 *   DataOut<dim> data_out;
1678 *   data_out.attach_dof_handler(dof_handler);
1679 *   data_out.add_data_vector(solution2, "u");
1680 *  
1681 *   data_out.build_patches(degree);
1682 *  
1683 *   data_out.write_vtu(gnuplot_output);
1684 *   }
1685 *   }
1686 *  
1687 *  
1688 *  
1689 *   template <int dim>
1690 *   void DGMethod<dim>::run()
1691 *   {
1692 *   for (unsigned int cycle = 0; cycle < 6; ++cycle)
1693 *   {
1694 *   std::cout << "Cycle " << cycle << ':' << std::endl;
1695 *  
1696 *   if (cycle == 0)
1697 *   {
1698 * @endcode
1699 *
1700 * Create the rectangular domain.
1701 *
1702 * @code
1703 *   Point<dim> p1, p2;
1704 *   p1[0] = -1.;
1705 *   for (unsigned int i = 0; i < dim; ++i)
1706 *   p2[i] = 1.;
1707 * @endcode
1708 *
1709 * Adjust the number of cells in different directions to obtain
1710 * completely isotropic cells for the original mesh.
1711 *
1712 * @code
1713 *   std::vector<unsigned int> repetitions(dim, 1);
1714 *   repetitions[0] = 2;
1715 *   GridGenerator::subdivided_hyper_rectangle(triangulation,
1716 *   repetitions,
1717 *   p1,
1718 *   p2);
1719 *  
1720 *   triangulation.refine_global(5 - dim);
1721 *   }
1722 *   else
1723 *   refine_grid();
1724 *  
1725 *  
1726 *   std::cout << " Number of active cells: "
1727 *   << triangulation.n_active_cells() << std::endl;
1728 *  
1729 *   setup_system();
1730 *  
1731 *   std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
1732 *   << std::endl;
1733 *  
1734 *   Timer assemble_timer;
1735 *   assemble_system();
1736 *   std::cout << " Time of assemble_system: " << assemble_timer.cpu_time()
1737 *   << std::endl;
1738 *   solve(solution2);
1739 *  
1740 *   output_results(cycle);
1741 *  
1742 *   std::cout << std::endl;
1743 *   }
1744 *   }
1745 *   } // namespace Step30
1746 *  
1747 *  
1748 *  
1749 *   int main()
1750 *   {
1751 *   try
1752 *   {
1753 *   using namespace Step30;
1754 *  
1755 * @endcode
1756 *
1757 * If you want to run the program in 3d, simply change the following
1758 * line to <code>const unsigned int dim = 3;</code>.
1759 *
1760 * @code
1761 *   const unsigned int dim = 2;
1762 *  
1763 *   {
1764 * @endcode
1765 *
1766 * First, we perform a run with isotropic refinement.
1767 *
1768 * @code
1769 *   std::cout << "Performing a " << dim
1770 *   << "D run with isotropic refinement..." << std::endl
1771 *   << "------------------------------------------------"
1772 *   << std::endl;
1773 *   DGMethod<dim> dgmethod_iso(false);
1774 *   dgmethod_iso.run();
1775 *   }
1776 *  
1777 *   {
1778 * @endcode
1779 *
1780 * Now we do a second run, this time with anisotropic refinement.
1781 *
1782 * @code
1783 *   std::cout << std::endl
1784 *   << "Performing a " << dim
1785 *   << "D run with anisotropic refinement..." << std::endl
1786 *   << "--------------------------------------------------"
1787 *   << std::endl;
1788 *   DGMethod<dim> dgmethod_aniso(true);
1789 *   dgmethod_aniso.run();
1790 *   }
1791 *   }
1792 *   catch (std::exception &exc)
1793 *   {
1794 *   std::cerr << std::endl
1795 *   << std::endl
1796 *   << "----------------------------------------------------"
1797 *   << std::endl;
1798 *   std::cerr << "Exception on processing: " << std::endl
1799 *   << exc.what() << std::endl
1800 *   << "Aborting!" << std::endl
1801 *   << "----------------------------------------------------"
1802 *   << std::endl;
1803 *   return 1;
1804 *   }
1805 *   catch (...)
1806 *   {
1807 *   std::cerr << std::endl
1808 *   << std::endl
1809 *   << "----------------------------------------------------"
1810 *   << std::endl;
1811 *   std::cerr << "Unknown exception!" << std::endl
1812 *   << "Aborting!" << std::endl
1813 *   << "----------------------------------------------------"
1814 *   << std::endl;
1815 *   return 1;
1816 *   };
1817 *  
1818 *   return 0;
1819 *   }
1820 * @endcode
1821<a name="step_30-Results"></a><h1>Results</h1>
1822
1823
1824
1825The output of this program consist of the console output, the SVG
1826files containing the grids, and the solutions given in VTU format.
1827@code
1828Performing a 2D run with isotropic refinement...
1829------------------------------------------------
1830Cycle 0:
1831 Number of active cells: 128
1832 Number of degrees of freedom: 512
1833 Time of assemble_system: 0.092049
1834 Writing grid to <grid-0.iso.svg>...
1835 Writing solution to <sol-0.iso.vtu>...
1836
1837Cycle 1:
1838 Number of active cells: 239
1839 Number of degrees of freedom: 956
1840 Time of assemble_system: 0.109519
1841 Writing grid to <grid-1.iso.svg>...
1842 Writing solution to <sol-1.iso.vtu>...
1843
1844Cycle 2:
1845 Number of active cells: 491
1846 Number of degrees of freedom: 1964
1847 Time of assemble_system: 0.08303
1848 Writing grid to <grid-2.iso.svg>...
1849 Writing solution to <sol-2.iso.vtu>...
1850
1851Cycle 3:
1852 Number of active cells: 1031
1853 Number of degrees of freedom: 4124
1854 Time of assemble_system: 0.278987
1855 Writing grid to <grid-3.iso.svg>...
1856 Writing solution to <sol-3.iso.vtu>...
1857
1858Cycle 4:
1859 Number of active cells: 2027
1860 Number of degrees of freedom: 8108
1861 Time of assemble_system: 0.305869
1862 Writing grid to <grid-4.iso.svg>...
1863 Writing solution to <sol-4.iso.vtu>...
1864
1865Cycle 5:
1866 Number of active cells: 4019
1867 Number of degrees of freedom: 16076
1868 Time of assemble_system: 0.47616
1869 Writing grid to <grid-5.iso.svg>...
1870 Writing solution to <sol-5.iso.vtu>...
1871
1872
1873Performing a 2D run with anisotropic refinement...
1874--------------------------------------------------
1875Cycle 0:
1876 Number of active cells: 128
1877 Number of degrees of freedom: 512
1878 Time of assemble_system: 0.052866
1879 Writing grid to <grid-0.aniso.svg>...
1880 Writing solution to <sol-0.aniso.vtu>...
1881
1882Cycle 1:
1883 Number of active cells: 171
1884 Number of degrees of freedom: 684
1885 Time of assemble_system: 0.050917
1886 Writing grid to <grid-1.aniso.svg>...
1887 Writing solution to <sol-1.aniso.vtu>...
1888
1889Cycle 2:
1890 Number of active cells: 255
1891 Number of degrees of freedom: 1020
1892 Time of assemble_system: 0.064132
1893 Writing grid to <grid-2.aniso.svg>...
1894 Writing solution to <sol-2.aniso.vtu>...
1895
1896Cycle 3:
1897 Number of active cells: 394
1898 Number of degrees of freedom: 1576
1899 Time of assemble_system: 0.119849
1900 Writing grid to <grid-3.aniso.svg>...
1901 Writing solution to <sol-3.aniso.vtu>...
1902
1903Cycle 4:
1904 Number of active cells: 648
1905 Number of degrees of freedom: 2592
1906 Time of assemble_system: 0.218244
1907 Writing grid to <grid-4.aniso.svg>...
1908 Writing solution to <sol-4.aniso.vtu>...
1909
1910Cycle 5:
1911 Number of active cells: 1030
1912 Number of degrees of freedom: 4120
1913 Time of assemble_system: 0.128121
1914 Writing grid to <grid-5.aniso.svg>...
1915 Writing solution to <sol-5.aniso.vtu>...
1916@endcode
1917
1918This text output shows the reduction in the number of cells which results from
1919the successive application of anisotropic refinement. After the last refinement
1920step the savings have accumulated so much that almost four times as many cells
1921and thus degrees of freedom are needed in the isotropic case. The time needed for assembly
1922scales with a similar factor.
1923
1924The first interesting part is of course to see how the meshes look like.
1925On the left are the isotropically refined ones, on the right the
1926anisotropic ones (colors indicate the refinement level of cells):
1927
1928<table width="80%" align="center">
1929 <tr>
1930 <td align="center">
1931 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-0.iso.9.2.png" alt="">
1932 </td>
1933 <td align="center">
1934 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-0.aniso.9.2.png" alt="">
1935 </td>
1936 </tr>
1937 <tr>
1938 <td align="center">
1939 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-1.iso.9.2.png" alt="">
1940 </td>
1941 <td align="center">
1942 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-1.aniso.9.2.png" alt="">
1943 </td>
1944 </tr>
1945 <tr>
1946 <td align="center">
1947 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-2.iso.9.2.png" alt="">
1948 </td>
1949 <td align="center">
1950 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-2.aniso.9.2.png" alt="">
1951 </td>
1952 </tr>
1953 <tr>
1954 <td align="center">
1955 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-3.iso.9.2.png" alt="">
1956 </td>
1957 <td align="center">
1958 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-3.aniso.9.2.png" alt="">
1959 </td>
1960 </tr>
1961 <tr>
1962 <td align="center">
1963 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-4.iso.9.2.png" alt="">
1964 </td>
1965 <td align="center">
1966 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-4.aniso.9.2.png" alt="">
1967 </td>
1968 </tr>
1969 <tr>
1970 <td align="center">
1971 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-5.iso.9.2.png" alt="">
1972 </td>
1973 <td align="center">
1974 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-5.aniso.9.2.png" alt="">
1975 </td>
1976 </tr>
1977</table>
1978
1979
1980The other interesting thing is, of course, to see the solution on these
1981two sequences of meshes. Here they are, on the refinement cycles 1 and 4,
1982clearly showing that the solution is indeed composed of <i>discontinuous</i> piecewise
1983polynomials:
1984
1985<table width="60%" align="center">
1986 <tr>
1987 <td align="center">
1988 <img src="https://www.dealii.org/images/steps/developer/step-30.sol-1.iso.9.2.png" alt="">
1989 </td>
1990 <td align="center">
1991 <img src="https://www.dealii.org/images/steps/developer/step-30.sol-1.aniso.9.2.png" alt="">
1992 </td>
1993 </tr>
1994 <tr>
1995 <td align="center">
1996 <img src="https://www.dealii.org/images/steps/developer/step-30.sol-4.iso.9.2.png" alt="">
1997 </td>
1998 <td align="center">
1999 <img src="https://www.dealii.org/images/steps/developer/step-30.sol-4.aniso.9.2.png" alt="">
2000 </td>
2001 </tr>
2002</table>
2003
2004We see, that the solution on the anisotropically refined mesh is very similar to
2005the solution obtained on the isotropically refined mesh. Thus the anisotropic
2006indicator seems to effectively select the appropriate cells for anisotropic
2007refinement.
2008
2009The pictures also explain why the mesh is refined as it is.
2010In the whole left part of the domain refinement is only performed along the
2011@f$y@f$-axis of cells. In the right part of the domain the refinement is dominated by
2012isotropic refinement, as the anisotropic feature of the solution - the jump from
2013one to zero - is not well aligned with the mesh where the advection direction
2014takes a turn. However, at the bottom and closest (to the observer) parts of the
2015quarter circle this jumps again becomes more and more aligned
2016with the mesh and the refinement algorithm reacts by creating anisotropic cells
2017of increasing aspect ratio.
2018
2019It might seem that the necessary alignment of anisotropic features and the
2020coarse mesh can decrease performance significantly for real world
2021problems. That is not wrong in general: If one were, for example, to apply
2022anisotropic refinement to problems in which shocks appear (e.g., the
2023equations solved in @ref step_69 "step-69"), then it many cases the shock is not aligned
2024with the mesh and anisotropic refinement will help little unless one also
2025introduces techniques to move the mesh in alignment with the shocks.
2026On the other hand, many steep features of solutions are due to boundary layers.
2027In those cases, the mesh is already aligned with the anisotropic features
2028because it is of course aligned with the boundary itself, and anisotropic
2029refinement will almost always increase the efficiency of computations on
2030adapted grids for these cases.
2031 *
2032 *
2033<a name="step_30-PlainProg"></a>
2034<h1> The plain program</h1>
2035@include "step-30.cc"
2036*/
const std::vector< double > & get_JxW_values() const
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
Definition point.h:111
void initialize(const MatrixType &A, const AdditionalData parameters)
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
Point< 3 > vertices[4]
unsigned int level
Definition grid_out.cc:4626
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
UpdateFlags
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void approximate_gradient(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InputVector &solution, Vector< float > &derivative_norm, const unsigned int component=0)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
double diameter(const Triangulation< dim, spacedim > &tria)
@ valid
Iterator points to a valid object.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
STL namespace.
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation