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