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