606 *
Assert(values.size() == points.size(),
609 * std::fill(values.begin(), values.end(), 0.);
615 *
class BoundaryValues :
public Function<dim>
619 * std::vector<double> & values,
620 *
const unsigned int = 0) const override
622 *
Assert(values.size() == points.size(),
625 *
for (
unsigned int i = 0; i < values.size(); ++i)
627 *
if (points[i](0) < 0.5)
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
654 *
void value_list(
const std::vector<
Point<dim>> &points,
657 *
Assert(values.size() == points.size(),
660 *
for (
unsigned int i = 0; i < points.size(); ++i)
662 *
if (points[i](0) > 0)
664 * values[i](0) = -points[i](1);
665 * values[i](1) = points[i](0);
670 * values[i](0) = -points[i](1);
681 * <a name=
"ClassDGTransportEquation"></a>
682 * <h3>Class: DGTransportEquation</h3>
686 * This declaration of
this class is utterly unaffected by our current
691 *
class DGTransportEquation
694 * DGTransportEquation();
712 *
const Beta<dim> beta_function;
713 *
const RHS<dim> rhs_function;
714 *
const BoundaryValues<dim> boundary_function;
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
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.
733 * DGTransportEquation<dim>::DGTransportEquation()
736 * , boundary_function()
742 *
void DGTransportEquation<dim>::assemble_cell_term(
749 * std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
750 * std::vector<double> rhs(fe_v.n_quadrature_points);
752 * beta_function.value_list(fe_v.get_quadrature_points(), beta);
753 * rhs_function.value_list(fe_v.get_quadrature_points(), rhs);
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)
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) *
769 *
void DGTransportEquation<dim>::assemble_boundary_term(
774 *
const std::vector<double> & JxW = fe_v.get_JxW_values();
775 *
const std::vector<Tensor<1, dim>> &normals = fe_v.get_normal_vectors();
777 * std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
778 * std::vector<double> g(fe_v.n_quadrature_points);
780 * beta_function.value_list(fe_v.get_quadrature_points(), beta);
781 * boundary_function.value_list(fe_v.get_quadrature_points(), g);
783 *
for (
unsigned int point = 0;
point < fe_v.n_quadrature_points; ++
point)
785 *
const double beta_n = beta[
point] * normals[
point];
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) *
792 *
for (
unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
800 *
void DGTransportEquation<dim>::assemble_face_term(
808 *
const std::vector<double> & JxW = fe_v.get_JxW_values();
809 *
const std::vector<Tensor<1, dim>> &normals = fe_v.get_normal_vectors();
811 * std::vector<Point<dim>> beta(fe_v.n_quadrature_points);
813 * beta_function.value_list(fe_v.get_quadrature_points(), beta);
815 *
for (
unsigned int point = 0;
point < fe_v.n_quadrature_points; ++
point)
817 *
const double beta_n = beta[
point] * normals[
point];
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) *
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) *
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) *
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];
852 * <a name=
"ClassDGMethod"></a>
853 * <h3>Class: DGMethod</h3>
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).
865 * DGMethod(
const bool anisotropic);
870 *
void setup_system();
871 *
void assemble_system();
873 *
void refine_grid();
874 *
void set_anisotropic_flags();
875 *
void output_results(
const unsigned int cycle)
const;
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.
886 *
const unsigned int degree;
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.
900 *
const double anisotropic_threshold_ratio;
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
908 *
const bool anisotropic;
916 *
const DGTransportEquation<dim> dg;
921 * DGMethod<dim>::DGMethod(
const bool anisotropic)
926 * Change here
for DG methods of different degrees.
932 * , anisotropic_threshold_ratio(3.)
933 * , anisotropic(anisotropic)
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.
944 * quadrature(degree + 1)
945 * , face_quadrature(degree + 1)
952 *
void DGMethod<dim>::setup_system()
954 * dof_handler.distribute_dofs(fe);
955 * sparsity_pattern.
reinit(dof_handler.n_dofs(),
956 * dof_handler.n_dofs(),
964 * sparsity_pattern.compress();
966 * system_matrix.reinit(sparsity_pattern);
968 * solution2.reinit(dof_handler.n_dofs());
969 * right_hand_side.reinit(dof_handler.n_dofs());
976 * <a name=
"Functionassemble_system"></a>
977 * <h4>
Function: assemble_system</h4>
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.
993 *
void DGMethod<dim>::assemble_system()
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);
1009 *
FEValues<dim> fe_v(mapping, fe, quadrature, update_flags);
1013 * face_update_flags);
1017 * face_update_flags);
1021 * neighbor_face_update_flags);
1032 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1037 * fe_v.reinit(cell);
1039 * dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector);
1041 * cell->get_dof_indices(dofs);
1045 *
const auto face = cell->face(face_no);
1049 * Case (a): The face is at the boundary.
1052 * if (face->at_boundary())
1054 * fe_v_face.reinit(cell, face_no);
1056 * dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector);
1062 *
const auto neighbor = cell->neighbor(face_no);
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
1075 * (There is a slightly confusing corner case: If we
1076 * are in 1
d -- 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 1
d, 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 1
d, this function
1084 *
always returns `false`, and consequently in 1
d 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
1090 * if (face->has_children())
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.
1100 *
const unsigned int neighbor2 =
1101 * cell->neighbor_face_no(face_no);
1105 * Now we
loop over all subfaces, i.e. the children and
1106 * possibly grandchildren of the current face.
1109 *
for (
unsigned int subface_no = 0;
1110 * subface_no < face->number_of_children();
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.
1121 *
const auto neighbor_child =
1122 * cell->neighbor_child_on_subface(face_no, subface_no);
1123 *
Assert(!neighbor_child->has_children(),
1128 * The remaining part of
this case is unchanged.
1135 * fe_v_subface.reinit(cell, face_no, subface_no);
1136 * fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
1138 * dg.assemble_face_term(fe_v_subface,
1139 * fe_v_face_neighbor,
1145 * neighbor_child->get_dof_indices(dofs_neighbor);
1147 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1148 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1150 * system_matrix.add(dofs[i],
1152 * ue_vi_matrix(i, j));
1153 * system_matrix.add(dofs_neighbor[i],
1155 * ui_ve_matrix(i, j));
1156 * system_matrix.add(dofs_neighbor[i],
1158 * ue_ve_matrix(i, j));
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 1
d 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 1
d, where as mentioned above
this is
1179 *
if (dim > 1 && cell->neighbor_is_coarser(face_no))
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 1
d). So
for 2
d and
1187 * 3
d, 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:
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())))))
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.
1213 * const unsigned int neighbor2 =
1214 * cell->neighbor_of_neighbor(face_no);
1220 * fe_v_face.reinit(cell, face_no);
1221 * fe_v_face_neighbor.reinit(neighbor, neighbor2);
1223 * dg.assemble_face_term(fe_v_face,
1224 * fe_v_face_neighbor,
1230 * neighbor->get_dof_indices(dofs_neighbor);
1232 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1233 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1235 * system_matrix.add(dofs[i],
1237 * ue_vi_matrix(i, j));
1238 * system_matrix.add(dofs_neighbor[i],
1240 * ui_ve_matrix(i, j));
1241 * system_matrix.add(dofs_neighbor[i],
1243 * ue_ve_matrix(i, j));
1249 * We do not need to consider a case (d), as those
1250 * faces are treated 'from the other side within
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));
1262 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1263 * right_hand_side(dofs[i]) += cell_vector(i);
1271 * <a name=
"Solver"></a>
1276 * For
this simple problem we use the simple Richardson iteration again. The
1277 * solver is completely unaffected by our anisotropic changes.
1280 *
template <
int dim>
1288 * preconditioner.
initialize(system_matrix, fe.dofs_per_cell);
1290 * solver.solve(system_matrix, solution, right_hand_side, preconditioner);
1297 * <a name=
"Refinement"></a>
1298 * <h3>Refinement</h3>
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.
1306 *
template <
int dim>
1307 *
void DGMethod<dim>::refine_grid()
1319 * gradient_indicator);
1323 * and
scale it to obtain an error indicator.
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);
1331 * Then we use
this indicator to flag the 30 percent of the cells with
1332 * highest error indicator to be refined.
1336 * gradient_indicator,
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.
1350 * set_anisotropic_flags();
1353 * Now execute the refinement considering anisotropic as well as isotropic
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.
1369 * template <
int dim>
1370 *
void DGMethod<dim>::set_anisotropic_flags()
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.
1384 * face_update_flags);
1388 * face_update_flags);
1396 * Now we need to
loop over all active cells.
1399 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1402 * We only need to consider cells which are flagged
for refinement.
1405 *
if (cell->refine_flag_set())
1412 *
const auto face = cell->face(face_no);
1414 *
if (!face->at_boundary())
1416 *
Assert(cell->neighbor(face_no).state() ==
1419 *
const auto neighbor = cell->neighbor(face_no);
1421 * std::vector<double> u(fe_v_face.n_quadrature_points);
1422 * std::vector<double> u_neighbor(fe_v_face.n_quadrature_points);
1426 * The four cases of different neighbor relations seen in
1427 * the assembly routines are repeated much in the same way
1431 *
if (face->has_children())
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.
1441 * unsigned int neighbor2 = cell->neighbor_face_no(face_no);
1444 * Now we loop over all subfaces,
1447 * for (unsigned int subface_no = 0;
1448 * subface_no < face->number_of_children();
1453 * get an iterator pointing to the cell behind the
1454 * present subface...
1457 * const auto neighbor_child =
1458 * cell->neighbor_child_on_subface(face_no,
1460 * Assert(!neighbor_child->has_children(),
1461 * ExcInternalError());
1464 * ... and reinit the respective FEFaceValues and
1465 * FESubFaceValues objects.
1468 * fe_v_subface.reinit(cell, face_no, subface_no);
1469 * fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
1472 * We obtain the function values
1475 * fe_v_subface.get_function_values(solution2, u);
1476 * fe_v_face_neighbor.get_function_values(solution2,
1480 * as well as the quadrature weights, multiplied by
1481 * the Jacobian determinant.
1484 * const std::vector<double> &JxW =
1485 * fe_v_subface.get_JxW_values();
1488 * Now we loop over all quadrature points
1491 * for (unsigned int x = 0;
1492 * x < fe_v_subface.n_quadrature_points;
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.
1510 * jump[face_no / 2] +=
1511 * std::abs(u[x] - u_neighbor[x]) * JxW[x];
1514 * We also sum up the scaled weights to obtain
1515 * the measure of the face.
1518 * area[face_no / 2] += JxW[x];
1524 * if (!cell->neighbor_is_coarser(face_no))
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
1535 * unsigned int neighbor2 =
1536 * cell->neighbor_of_neighbor(face_no);
1538 * fe_v_face.reinit(cell, face_no);
1539 * fe_v_face_neighbor.reinit(neighbor, neighbor2);
1541 * fe_v_face.get_function_values(solution2, u);
1542 * fe_v_face_neighbor.get_function_values(solution2,
1545 * const std::vector<double> &JxW =
1546 * fe_v_face.get_JxW_values();
1548 * for (unsigned int x = 0;
1549 * x < fe_v_face.n_quadrature_points;
1552 * jump[face_no / 2] +=
1553 * std::abs(u[x] - u_neighbor[x]) * JxW[x];
1554 * area[face_no / 2] += JxW[x];
1557 * else // i.e. neighbor is coarser than cell
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.
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());
1585 * fe_v_face.reinit(cell, face_no);
1586 * fe_v_subface.reinit(neighbor,
1587 * neighbor_face_subface.first,
1588 * neighbor_face_subface.second);
1590 * fe_v_face.get_function_values(solution2, u);
1591 * fe_v_subface.get_function_values(solution2,
1594 * const std::vector<double> &JxW =
1595 * fe_v_face.get_JxW_values();
1597 * for (unsigned int x = 0;
1598 * x < fe_v_face.n_quadrature_points;
1601 * jump[face_no / 2] +=
1602 * std::abs(u[x] - u_neighbor[x]) * JxW[x];
1603 * area[face_no / 2] += JxW[x];
1611 * Now we analyze the size of the mean jumps, which we get dividing
1612 * the jumps by the measure of the respective faces.
1615 * std::array<double, dim> average_jumps;
1616 * double sum_of_average_jumps = 0.;
1617 * for (unsigned int i = 0; i < dim; ++i)
1619 * average_jumps[i] = jump(i) / area(i);
1620 * sum_of_average_jumps += average_jumps[i];
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.
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));
1644 * <a name="TheRest"></a>
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.
1655 * template <int dim>
1656 * void DGMethod<dim>::output_results(const unsigned int cycle) const
1658 * std::string refine_type;
1660 * refine_type = ".aniso";
1662 * refine_type = ".iso";
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);
1671 * grid_out.write_svg(triangulation, svg_output);
1675 * const std::string filename =
1676 * "sol-" + std::to_string(cycle) + refine_type + ".vtu";
1677 * std::cout << " Writing solution to <" << filename << ">..."
1679 * std::ofstream gnuplot_output(filename);
1681 * DataOut<dim> data_out;
1682 * data_out.attach_dof_handler(dof_handler);
1683 * data_out.add_data_vector(solution2, "u");
1685 * data_out.build_patches(degree);
1687 * data_out.write_vtu(gnuplot_output);
1693 * template <int dim>
1694 * void DGMethod<dim>::run()
1696 * for (unsigned int cycle = 0; cycle < 6; ++cycle)
1698 * std::cout << "Cycle " << cycle << ':
' << std::endl;
1704 * Create the rectangular domain.
1707 * Point<dim> p1, p2;
1710 * for (unsigned int i = 0; i < dim; ++i)
1714 * Adjust the number of cells in different directions to obtain
1715 * completely isotropic cells for the original mesh.
1718 * std::vector<unsigned int> repetitions(dim, 1);
1719 * repetitions[0] = 2;
1720 * GridGenerator::subdivided_hyper_rectangle(triangulation,
1725 * triangulation.refine_global(5 - dim);
1731 * std::cout << " Number of active cells: "
1732 * << triangulation.n_active_cells() << std::endl;
1736 * std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
1739 * Timer assemble_timer;
1740 * assemble_system();
1741 * std::cout << " Time of assemble_system: " << assemble_timer.cpu_time()
1745 * output_results(cycle);
1747 * std::cout << std::endl;
1750 * } // namespace Step30
1758 * using namespace Step30;
1762 * If you want to run the program in 3D, simply change the following
1763 * line to <code>const unsigned int dim = 3;</code>.
1766 * const unsigned int dim = 2;
1771 * First, we perform a run with isotropic refinement.
1774 * std::cout << "Performing a " << dim
1775 * << "D run with isotropic refinement..." << std::endl
1776 * << "------------------------------------------------"
1778 * DGMethod<dim> dgmethod_iso(false);
1779 * dgmethod_iso.run();
1785 * Now we do a second run, this time with anisotropic refinement.
1788 * std::cout << std::endl
1789 * << "Performing a " << dim
1790 * << "D run with anisotropic refinement..." << std::endl
1791 * << "--------------------------------------------------"
1793 * DGMethod<dim> dgmethod_aniso(true);
1794 * dgmethod_aniso.run();
1797 * catch (std::exception &exc)
1799 * std::cerr << std::endl
1801 * << "----------------------------------------------------"
1803 * std::cerr << "Exception on processing: " << std::endl
1804 * << exc.what() << std::endl
1805 * << "Aborting!" << std::endl
1806 * << "----------------------------------------------------"
1812 * std::cerr << std::endl
1814 * << "----------------------------------------------------"
1816 * std::cerr << "Unknown exception!" << std::endl
1817 * << "Aborting!" << std::endl
1818 * << "----------------------------------------------------"
1826 <a name="Results"></a><h1>Results</h1>
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.
1833 Performing a 2D run with isotropic refinement...
1834 ------------------------------------------------
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>...
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>...
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>...
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>...
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>...
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>...
1878 Performing a 2D run with anisotropic refinement...
1879 --------------------------------------------------
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>...
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>...
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>...
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>...
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>...
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>...
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.
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):
1933 <table width="80%" align="center">
1936 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-0.iso.9.2.png" alt="">
1939 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-0.aniso.9.2.png" alt="">
1944 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-1.iso.9.2.png" alt="">
1947 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-1.aniso.9.2.png" alt="">
1952 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-2.iso.9.2.png" alt="">
1955 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-2.aniso.9.2.png" alt="">
1960 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-3.iso.9.2.png" alt="">
1963 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-3.aniso.9.2.png" alt="">
1968 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-4.iso.9.2.png" alt="">
1971 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-4.aniso.9.2.png" alt="">
1976 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-5.iso.9.2.png" alt="">
1979 <img src="https://www.dealii.org/images/steps/developer/step-30.grid-5.aniso.9.2.png" alt="">
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
1990 <table width="60%" align="center">
1993 <img src="https://www.dealii.org/images/steps/developer/step-30.sol-1.iso.9.2.png" alt="">
1996 <img src="https://www.dealii.org/images/steps/developer/step-30.sol-1.aniso.9.2.png" alt="">
2001 <img src="https://www.dealii.org/images/steps/developer/step-30.sol-4.iso.9.2.png" alt="">
2004 <img src="https://www.dealii.org/images/steps/developer/step-30.sol-4.aniso.9.2.png" alt="">
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
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.
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.
2038 <a name="PlainProg"></a>
2039 <h1> The plain program</h1>
2040 @include "step-30.cc"