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