939 *
double product = 1;
940 *
for (
unsigned int d = 0;
d < dim; ++
d)
941 * product *= (p[
d] + 1);
950 * <a name=
"Implementationofthemainclass"></a>
951 * <h3>Implementation of the main
class</h3>
956 * <a name=
"LaplaceProblemLaplaceProblemconstructor"></a>
957 * <h4>LaplaceProblem::LaplaceProblem constructor</h4>
961 * The constructor of
this class is fairly straightforward. It associates
963 * maximal polynomial degree to 7 (in 1
d and 2
d) or 5 (in 3
d and higher). We
964 *
do so because
using higher order polynomial degrees becomes prohibitively
965 * expensive, especially in higher space dimensions.
969 * Following
this, we fill the collections of finite element, and cell and
970 * face quadrature objects. We start with quadratic elements, and each
971 * quadrature formula is chosen so that it is appropriate
for the matching
977 * calculate coefficient in Fourier series as described in the introduction.
983 * In order to resize fourier_coefficients
Table, we use the following
987 *
template <
int dim,
typename T>
991 *
for (
unsigned int d = 0;
d < dim;
d++)
993 * coeff.reinit(size);
997 * LaplaceProblem<dim>::LaplaceProblem()
999 * , max_degree(dim <= 2 ? 7 : 5)
1001 *
for (
unsigned int degree = 2; degree <= max_degree; ++degree)
1003 * fe_collection.push_back(
FE_Q<dim>(degree));
1004 * quadrature_collection.push_back(
QGauss<dim>(degree + 1));
1010 * As described in the introduction, we define the Fourier vectors @f${\bf
1011 * k}@f$
for which we want to compute Fourier coefficients of the solution
1012 * on each cell as follows. In 2
d, we will need coefficients corresponding
1013 * to vectors @f${\bf k}=(2 \pi i, 2\pi j)^
T@f$
for which @f$\
sqrt{i^2+j^2}\le
N@f$,
1014 * with @f$i,j@f$ integers and @f$N@f$ being the maximal polynomial degree we use
1016 * constructor first parameter @f$N@f$ defines the number of coefficients in 1D
1017 * with the total number of coefficients being @f$N^{dim}@f$. Although we will
1018 * not use coefficients corresponding to
1019 * @f$\sqrt{i^2+j^2}> N@f$ and @f$i+j==0@f$, the overhead of their calculation is
1020 * minimal. The transformation matrices for each FiniteElement will be
1021 * calculated only once the first time they are required in the course of
1022 * hp-adaptive refinement. Because we work on the unit cell, we can do all
1023 * this work without a mapping from reference to real cell and consequently
1024 * can precalculate these matrices. The calculation of expansion
1025 * coefficients for a particular set of local degrees of freedom on a given
1026 * cell then follows as a simple matrix-vector product.
1027 * The 3d case is handled analogously.
1030 * const unsigned int N = max_degree;
1034 * We will need to assemble the matrices that do the Fourier transforms
1035 * for each of the finite elements we deal with, i.e. the matrices @f${\cal
1036 * F}_{{\bf k},j}@f$ defined in the introduction. We have to do that for
1037 * each of the finite elements in use. To that end we need a quadrature
1038 * rule. In this example we use the same quadrature formula for each
1039 * finite element, namely that is obtained by iterating a
1040 * 2-point Gauss formula as many times as the maximal exponent we use for
1041 * the term @f$e^{i{\bf k}\cdot{\bf x}}@f$:
1044 * QGauss<1> base_quadrature(2);
1045 * QIterated<dim> quadrature(base_quadrature, N);
1046 * for (unsigned int i = 0; i < fe_collection.size(); i++)
1047 * fourier_q_collection.push_back(quadrature);
1051 * Now we are ready to set-up the FESeries::Fourier object
1054 * const std::vector<unsigned int> n_coefficients_per_direction(
1055 * fe_collection.size(), N);
1056 * fourier = std_cxx14::make_unique<FESeries::Fourier<dim>>(
1057 * n_coefficients_per_direction, fe_collection, fourier_q_collection);
1061 * We need to resize the matrix of fourier coefficients according to the
1062 * number of modes N.
1065 * resize(fourier_coefficients, N);
1072 * <a name="LaplaceProblemLaplaceProblemdestructor"></a>
1073 * <h4>LaplaceProblem::~LaplaceProblem destructor</h4>
1077 * The destructor is unchanged from what we already did in @ref step_6 "step-6":
1080 * template <int dim>
1081 * LaplaceProblem<dim>::~LaplaceProblem()
1083 * dof_handler.clear();
1090 * <a name="LaplaceProblemsetup_system"></a>
1091 * <h4>LaplaceProblem::setup_system</h4>
1095 * This function is again a verbatim copy of what we already did in
1096 * @ref step_6 "step-6". Despite function calls with exactly the same names and arguments,
1097 * the algorithms used internally are different in some aspect since the
1098 * dof_handler variable here is an hp object.
1101 * template <int dim>
1102 * void LaplaceProblem<dim>::setup_system()
1104 * dof_handler.distribute_dofs(fe_collection);
1106 * solution.reinit(dof_handler.n_dofs());
1107 * system_rhs.reinit(dof_handler.n_dofs());
1109 * constraints.clear();
1110 * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1111 * VectorTools::interpolate_boundary_values(dof_handler,
1113 * Functions::ZeroFunction<dim>(),
1115 * constraints.close();
1117 * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
1118 * DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
1119 * sparsity_pattern.copy_from(dsp);
1121 * system_matrix.reinit(sparsity_pattern);
1129 * <a name="LaplaceProblemassemble_system"></a>
1130 * <h4>LaplaceProblem::assemble_system</h4>
1134 * This is the function that assembles the global matrix and right hand side
1135 * vector from the local contributions of each cell. Its main working is as
1136 * has been described in many of the tutorial programs before. The
1137 * significant deviations are the ones necessary for <i>hp</i> finite
1138 * element methods. In particular, that we need to use a collection of
1139 * FEValues object (implemented through the hp::FEValues class), and that we
1140 * have to eliminate constrained degrees of freedom already when copying
1141 * local contributions into global objects. Both of these are explained in
1142 * detail in the introduction of this program.
1146 * One other slight complication is the fact that because we use different
1147 * polynomial degrees on different cells, the matrices and vectors holding
1148 * local contributions do not have the same size on all cells. At the
1149 * beginning of the loop over all cells, we therefore each time have to
1150 * resize them to the correct size (given by
1151 * <code>dofs_per_cell</code>). Because these classes are implement in such
1152 * a way that reducing the size of a matrix or vector does not release the
1153 * currently allocated memory (unless the new size is zero), the process of
1154 * resizing at the beginning of the loop will only require re-allocation of
1155 * memory during the first few iterations. Once we have found in a cell with
1156 * the maximal finite element degree, no more re-allocations will happen
1157 * because all subsequent <code>reinit</code> calls will only set the size
1158 * to something that fits the currently allocated memory. This is important
1159 * since allocating memory is expensive, and doing so every time we visit a
1160 * new cell would take significant compute time.
1163 * template <int dim>
1164 * void LaplaceProblem<dim>::assemble_system()
1166 * hp::FEValues<dim> hp_fe_values(fe_collection,
1167 * quadrature_collection,
1168 * update_values | update_gradients |
1169 * update_quadrature_points |
1170 * update_JxW_values);
1172 * RightHandSide<dim> rhs_function;
1174 * FullMatrix<double> cell_matrix;
1175 * Vector<double> cell_rhs;
1177 * std::vector<types::global_dof_index> local_dof_indices;
1179 * for (const auto &cell : dof_handler.active_cell_iterators())
1181 * const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
1183 * cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1186 * cell_rhs.reinit(dofs_per_cell);
1189 * hp_fe_values.reinit(cell);
1191 * const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
1193 * std::vector<double> rhs_values(fe_values.n_quadrature_points);
1194 * rhs_function.value_list(fe_values.get_quadrature_points(), rhs_values);
1196 * for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
1198 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1200 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1201 * cell_matrix(i, j) +=
1202 * (fe_values.shape_grad(i, q_point) * // grad phi_i(x_q)
1203 * fe_values.shape_grad(j, q_point) * // grad phi_j(x_q)
1204 * fe_values.JxW(q_point)); // dx
1206 * cell_rhs(i) += (fe_values.shape_value(i, q_point) * // phi_i(x_q)
1207 * rhs_values[q_point] * // f(x_q)
1208 * fe_values.JxW(q_point)); // dx
1211 * local_dof_indices.resize(dofs_per_cell);
1212 * cell->get_dof_indices(local_dof_indices);
1214 * constraints.distribute_local_to_global(
1215 * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1224 * <a name="LaplaceProblemsolve"></a>
1225 * <h4>LaplaceProblem::solve</h4>
1229 * The function solving the linear system is entirely unchanged from
1230 * previous examples. We simply try to reduce the initial residual (which
1231 * equals the @f$l_2@f$ norm of the right hand side) by a certain factor:
1234 * template <int dim>
1235 * void LaplaceProblem<dim>::solve()
1237 * SolverControl solver_control(system_rhs.size(),
1238 * 1e-12 * system_rhs.l2_norm());
1239 * SolverCG<Vector<double>> cg(solver_control);
1241 * PreconditionSSOR<SparseMatrix<double>> preconditioner;
1242 * preconditioner.initialize(system_matrix, 1.2);
1244 * cg.solve(system_matrix, solution, system_rhs, preconditioner);
1246 * constraints.distribute(solution);
1254 * <a name="LaplaceProblempostprocess"></a>
1255 * <h4>LaplaceProblem::postprocess</h4>
1259 * After solving the linear system, we will want to postprocess the
1260 * solution. Here, all we do is to estimate the error, estimate the local
1261 * smoothness of the solution as described in the introduction, then write
1262 * graphical output, and finally refine the mesh in both @f$h@f$ and @f$p@f$
1263 * according to the indicators computed before. We do all this in the same
1264 * function because we want the estimated error and smoothness indicators
1265 * not only for refinement, but also include them in the graphical output.
1268 * template <int dim>
1269 * void LaplaceProblem<dim>::postprocess(const unsigned int cycle)
1273 * Let us start with computing estimated error and smoothness indicators,
1274 * which each are one number for each active cell of our
1275 * triangulation. For the error indicator, we use the KellyErrorEstimator
1276 * class as always. Estimating the smoothness is done in the respective
1277 * function of this class; that function is discussed further down below:
1280 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1281 * KellyErrorEstimator<dim>::estimate(
1283 * face_quadrature_collection,
1284 * std::map<types::boundary_id, const Function<dim> *>(),
1286 * estimated_error_per_cell);
1289 * Vector<float> smoothness_indicators(triangulation.n_active_cells());
1290 * estimate_smoothness(smoothness_indicators);
1294 * Next we want to generate graphical output. In addition to the two
1295 * estimated quantities derived above, we would also like to output the
1296 * polynomial degree of the finite elements used on each of the elements
1301 * The way to do that requires that we loop over all cells and poll the
1302 * active finite element index of them using
1303 * <code>cell-@>active_fe_index()</code>. We then use the result of this
1304 * operation and query the finite element collection for the finite
1305 * element with that index, and finally determine the polynomial degree of
1306 * that element. The result we put into a vector with one element per
1307 * cell. The DataOut class requires this to be a vector of
1308 * <code>float</code> or <code>double</code>, even though our values are
1309 * all integers, so that it what we use:
1313 * Vector<float> fe_degrees(triangulation.n_active_cells());
1314 * for (const auto &cell : dof_handler.active_cell_iterators())
1315 * fe_degrees(cell->active_cell_index()) =
1316 * fe_collection[cell->active_fe_index()].degree;
1320 * With now all data vectors available -- solution, estimated errors and
1321 * smoothness indicators, and finite element degrees --, we create a
1322 * DataOut object for graphical output and attach all data. Note that
1323 * the DataOut class has a second template argument (which defaults to
1324 * DoFHandler@<dim@>, which is why we have never seen it in previous
1325 * tutorial programs) that indicates the type of DoF handler to be
1326 * used. Here, we have to use the hp::DoFHandler class:
1329 * DataOut<dim, hp::DoFHandler<dim>> data_out;
1331 * data_out.attach_dof_handler(dof_handler);
1332 * data_out.add_data_vector(solution, "solution");
1333 * data_out.add_data_vector(estimated_error_per_cell, "error");
1334 * data_out.add_data_vector(smoothness_indicators, "smoothness");
1335 * data_out.add_data_vector(fe_degrees, "fe_degree");
1336 * data_out.build_patches();
1340 * The final step in generating output is to determine a file name, open
1341 * the file, and write the data into it (here, we use VTK format):
1344 * const std::string filename =
1345 * "solution-" + Utilities::int_to_string(cycle, 2) + ".vtk";
1346 * std::ofstream output(filename);
1347 * data_out.write_vtk(output);
1352 * After this, we would like to actually refine the mesh, in both @f$h@f$ and
1353 * @f$p@f$. The way we are going to do this is as follows: first, we use the
1354 * estimated error to flag those cells for refinement that have the
1355 * largest error. This is what we have always done:
1359 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
1360 * estimated_error_per_cell,
1366 * Next we would like to figure out which of the cells that have been
1367 * flagged for refinement should actually have @f$p@f$ increased instead of
1368 * @f$h@f$ decreased. The strategy we choose here is that we look at the
1369 * smoothness indicators of those cells that are flagged for refinement,
1370 * and increase @f$p@f$ for those with a smoothness larger than a certain
1371 * threshold. For this, we first have to determine the maximal and
1372 * minimal values of the smoothness indicators of all flagged cells,
1373 * which we do using a loop over all cells and comparing current minimal
1374 * and maximal values. (We start with the minimal and maximal values of
1375 * <i>all</i> cells, a range within which the minimal and maximal values
1376 * on cells flagged for refinement must surely lie.) Absent any better
1377 * strategies, we will then set the threshold above which will increase
1378 * @f$p@f$ instead of reducing @f$h@f$ as the mean value between minimal and
1379 * maximal smoothness indicators on cells flagged for refinement:
1382 * float max_smoothness = *std::min_element(smoothness_indicators.begin(),
1383 * smoothness_indicators.end()),
1384 * min_smoothness = *std::max_element(smoothness_indicators.begin(),
1385 * smoothness_indicators.end());
1386 * for (const auto &cell : dof_handler.active_cell_iterators())
1387 * if (cell->refine_flag_set())
1390 * std::max(max_smoothness,
1391 * smoothness_indicators(cell->active_cell_index()));
1393 * std::min(min_smoothness,
1394 * smoothness_indicators(cell->active_cell_index()));
1396 * const float threshold_smoothness = (max_smoothness + min_smoothness) / 2;
1400 * With this, we can go back, loop over all cells again, and for those
1401 * cells for which (i) the refinement flag is set, (ii) the smoothness
1402 * indicator is larger than the threshold, and (iii) we still have a
1403 * finite element with a polynomial degree higher than the current one
1404 * in the finite element collection, we then increase the polynomial
1405 * degree and in return remove the flag indicating that the cell should
1406 * undergo bisection. For all other cells, the refinement flags remain
1410 * for (const auto &cell : dof_handler.active_cell_iterators())
1411 * if (cell->refine_flag_set() &&
1412 * (smoothness_indicators(cell->active_cell_index()) >
1413 * threshold_smoothness) &&
1414 * (cell->active_fe_index() + 1 < fe_collection.size()))
1416 * cell->clear_refine_flag();
1417 * cell->set_active_fe_index(cell->active_fe_index() + 1);
1422 * At the end of this procedure, we then refine the mesh. During this
1423 * process, children of cells undergoing bisection inherit their mother
1424 * cell's finite element index:
1435 * <a name=
"LaplaceProblemcreate_coarse_grid"></a>
1436 * <h4>LaplaceProblem::create_coarse_grid</h4>
1440 * The following
function is used when creating the
initial grid. It is a
1441 * specialization
for the 2
d case, i.e. a corresponding
function needs to be
1442 * implemented
if the program is
run in anything other then 2
d. The
function
1443 * is actually stolen from @ref step_14
"step-14" and generates the same mesh used already
1444 * there, i.e. the square domain with the square hole in the middle. The
1445 * meaning of the different parts of
this function are explained in the
1446 * documentation of @ref step_14
"step-14":
1450 *
void LaplaceProblem<2>::create_coarse_grid()
1452 *
const unsigned int dim = 2;
1454 *
const std::vector<Point<2>>
vertices = {
1455 * {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0},
1456 * {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5},
1457 * {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0},
1458 * {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5},
1459 * {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}};
1461 *
const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
1462 * cell_vertices = {{{0, 1, 5, 6}},
1468 * {{10, 11, 14, 15}},
1469 * {{12, 13, 17, 18}},
1470 * {{14, 15, 19, 20}},
1471 * {{15, 16, 20, 21}},
1472 * {{16, 17, 21, 22}},
1473 * {{17, 18, 22, 23}}};
1475 *
const unsigned int n_cells = cell_vertices.size();
1478 *
for (
unsigned int i = 0; i <
n_cells; ++i)
1480 *
for (
unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell; ++j)
1481 * cells[i].
vertices[j] = cell_vertices[i][j];
1482 * cells[i].material_id = 0;
1494 * <a name=
"LaplaceProblemrun"></a>
1499 * This
function implements the logic of the program, as did the respective
1500 *
function in most of the previous programs already, see
for example
1501 * @ref step_6
"step-6".
1505 * Basically, it contains the adaptive
loop: in the
first iteration create a
1506 * coarse grid, and then
set up the linear system,
assemble it, solve, and
1507 * postprocess the solution including mesh refinement. Then start over
1508 * again. In the meantime, also output some information
for those staring at
1509 * the screen trying to figure out what the program does:
1512 *
template <
int dim>
1515 *
for (
unsigned int cycle = 0; cycle < 6; ++cycle)
1517 * std::cout <<
"Cycle " << cycle <<
':' << std::endl;
1520 * create_coarse_grid();
1524 * std::cout <<
" Number of active cells : "
1526 * <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1528 * <<
" Number of constraints : "
1529 * << constraints.n_constraints() << std::endl;
1531 * assemble_system();
1533 * postprocess(cycle);
1541 * <a name=
"LaplaceProblemestimate_smoothness"></a>
1542 * <h4>LaplaceProblem::estimate_smoothness</h4>
1546 * As described in the introduction, we will need to take the maximum
1547 * absolute
value of fourier coefficients which correspond to @f$k@f$-vector
1548 * @f$|{\bf k}|=
const@f$. To filter the coefficients
Table we
1550 * to be specified. The predicate should operate on
TableIndices and return
1551 * a pair of <code>
bool</code> and <code>
unsigned int</code>. The latter
1552 * is the
value of the map from TableIndicies to
unsigned int. It is
1553 * used to define subsets of coefficients from which we search for the
one
1554 * with highest absolute
value, i.
e. @f$l^\infty@f$-
norm. The <code>
bool</code>
1555 * parameter defines which indices should be used in processing. In the
1556 * current case we are interested in coefficients which correspond to
1557 * @f$0 < i*i+j*j <
N*
N@f$ and @f$0 < i*i+j*j+k*k <
N*
N@f$ in 2D and 3D, respectively.
1560 * template <
int dim>
1561 * std::pair<
bool,
unsigned int>
1562 * LaplaceProblem<dim>::predicate(const
TableIndices<dim> &ind)
1564 *
unsigned int v = 0;
1565 *
for (
unsigned int i = 0; i < dim; i++)
1566 * v += ind[i] * ind[i];
1567 *
if (v > 0 && v < max_degree * max_degree)
1568 *
return std::make_pair(
true, v);
1570 *
return std::make_pair(
false, v);
1575 * This last
function of significance implements the algorithm to estimate
1576 * the smoothness exponent
using the algorithms explained in detail in the
1577 * introduction. We will therefore only comment on those points that are of
1578 * implementational importance.
1581 *
template <
int dim>
1583 * LaplaceProblem<dim>::estimate_smoothness(
Vector<float> &smoothness_indicators)
1588 * we
set up the
object of
this class in the constructor, what we are left
1589 * to
do here is
apply this class to calculate coefficients and then
1590 * perform linear regression to fit their decay slope.
1597 * First thing to
do is to
loop over all cells and
do our work there, i.e.
1598 * to locally
do the Fourier
transform and estimate the decay coefficient.
1599 * We will use the following array as a scratch array in the
loop to store
1607 * Then here is the
loop:
1610 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1614 * Inside the
loop, we
first need to get the values of the local
1615 * degrees of freedom (which we put into the
1616 * <code>local_dof_values</code> array after setting it to the right
1617 * size) and then need to compute the Fourier
transform by multiplying
1618 *
this vector with the
matrix @f${\cal
F}@f$ corresponding to
this finite
1620 * that has to be provided with the <code>local_dof_values</code>,
1621 * <code>cell->active_fe_index()</code> and a
Table to store
1625 * local_dof_values.
reinit(cell->get_fe().dofs_per_cell);
1626 * cell->get_dof_values(solution, local_dof_values);
1628 * fourier->calculate(local_dof_values,
1629 * cell->active_fe_index(),
1630 * fourier_coefficients);
1634 * The next thing, as explained in the introduction, is that we wanted
1635 * to only fit our exponential decay of Fourier coefficients to the
1636 * largest coefficients
for each possible
value of @f$|{\bf k}|@f$. To
1638 * coefficients into the desired format. We'll only take those Fourier
1639 * coefficients with the largest magnitude for a given
value of @f$|{\bf
1643 * std::pair<std::vector<unsigned int>, std::vector<double>> res =
1644 * FESeries::process_coefficients<dim>(
1645 * fourier_coefficients,
1647 *
return this->predicate(indices);
1655 * The
first vector in the <code>std::pair</code> will store values of
1656 * the predicate, that is @f$i*i+j*j=
const@f$ or @f$i*i+j*j+k*k =
const@f$ in
1657 * 2D or 3D respectively. This vector will be the same
for all the cells
1658 * so we can calculate logarithms of the corresponding Fourier vectors
1659 * @f$|{\bf k}|@f$ only once in the whole
hp-refinement cycle:
1662 *
if (ln_k.size() == 0)
1664 * ln_k.resize(res.first.size(), 0);
1665 *
for (
unsigned int f = 0; f < ln_k.size(); f++)
1667 * std::log(2.0 *
numbers::PI * std::sqrt(1. * res.first[f]));
1672 * We have to calculate the logarithms of absolute values of
1673 * coefficients and use it in a linear regression fit to obtain @f$\mu@f$.
1676 *
for (
double &residual_element : res.second)
1677 * residual_element =
std::log(residual_element);
1679 * std::pair<double, double> fit =
1684 * The
final step is to compute the Sobolev index @f$s=\mu-\frac d2@f$ and
1685 * store it in the vector of estimated values
for each cell:
1688 * smoothness_indicators(cell->active_cell_index()) =
1689 * -fit.first - 1. * dim / 2;
1698 * <a name=
"Themainfunction"></a>
1699 * <h3>The main
function</h3>
1703 * The main
function is again verbatim what we had before: wrap creating and
1704 * running an
object of the main
class into a <code>try</code> block and catch
1705 * whatever exceptions are thrown, thereby producing meaningful output
if
1706 * anything should go wrong:
1713 *
using namespace Step27;
1715 * LaplaceProblem<2> laplace_problem;
1716 * laplace_problem.run();
1718 *
catch (std::exception &exc)
1720 * std::cerr << std::endl
1722 * <<
"----------------------------------------------------"
1724 * std::cerr <<
"Exception on processing: " << std::endl
1725 * << exc.what() << std::endl
1726 * <<
"Aborting!" << std::endl
1727 * <<
"----------------------------------------------------"
1734 * std::cerr << std::endl
1736 * <<
"----------------------------------------------------"
1738 * std::cerr <<
"Unknown exception!" << std::endl
1739 * <<
"Aborting!" << std::endl
1740 * <<
"----------------------------------------------------"
1748 <a name=
"Results"></a><h1>Results</h1>
1751 In
this section, we discuss a few results produced from running the
1752 current tutorial program. More results, in particular the extension to
1753 3
d calculations and determining how much compute time the individual
1754 components of the program take, are given in the @ref hp_paper .
1756 When
run,
this is what the program produces:
1760 [ 66%] Built target @ref step_27
"step-27"
1761 [100%] Run @ref step_27
"step-27" with Release configuration
1763 Number of active cells : 768
1764 Number of degrees of freedom: 3264
1765 Number of constraints : 384
1767 Number of active cells : 966
1768 Number of degrees of freedom: 5245
1769 Number of constraints : 936
1771 Number of active cells : 1143
1772 Number of degrees of freedom: 8441
1773 Number of constraints : 1929
1775 Number of active cells : 1356
1776 Number of degrees of freedom: 12349
1777 Number of constraints : 3046
1779 Number of active cells : 1644
1780 Number of degrees of freedom: 18178
1781 Number of constraints : 4713
1783 Number of active cells : 1728
1784 Number of degrees of freedom: 22591
1785 Number of constraints : 6095
1788 The
first thing we learn from
this is that the number of constrained degrees
1789 of freedom is on the order of 20-25% of the total number of degrees of
1790 freedom, at least on the later grids when we have elements of relatively
1791 high order (in 3
d, the fraction of constrained degrees of freedom can be up
1792 to 30%). This is, in fact, on the same order of magnitude as
for non-@f$hp@f$
1793 discretizations. For example, in the last step of the @ref step_6
"step-6"
1794 program, we have 18353 degrees of freedom, 4432 of which are
1795 constrained. The difference is that in the latter program, each constrained
1796 hanging node is constrained against only the two adjacent degrees of
1797 freedom, whereas in the @f$hp@f$
case, constrained nodes are constrained against
1798 many more degrees of freedom. Note also that the current program also
1799 includes nodes subject to Dirichlet boundary conditions in the list of
1800 constraints. In cycle 0, all the constraints are actually because of
1801 boundary conditions.
1803 Of maybe more interest is to look at the graphical output. First, here is the
1804 solution of the problem:
1806 <img src=
"https://www.dealii.org/images/steps/developer/step-27-solution.png"
1807 alt=
"Elevation plot of the solution, showing the lack of regularity near
1808 the interior (reentrant) corners."
1809 width=
"200" height=
"200">
1811 Secondly, let us look at the sequence of meshes generated:
1813 <div
class=
"threecolumn" style=
"width: 80%">
1815 <img src=
"https://www.dealii.org/images/steps/developer/step-27-mesh-0.svg"
1816 alt=
"Triangulation containing reentrant corners without adaptive refinement."
1817 width=
"200" height=
"200">
1820 <img src=
"https://www.dealii.org/images/steps/developer/step-27-mesh-1.svg"
1821 alt=
"Triangulation containing reentrant corners with one level of
1822 refinement. New cells are placed near the corners."
1823 width=
"200" height=
"200">
1826 <img src=
"https://www.dealii.org/images/steps/developer/step-27-mesh-2.svg"
1827 alt=
"Triangulation containing reentrant corners with two levels of
1828 refinement. New cells are placed near the corners."
1829 width=
"200" height=
"200">
1832 <img src=
"https://www.dealii.org/images/steps/developer/step-27-mesh-3.svg"
1833 alt=
"Triangulation containing reentrant corners with three levels of
1834 refinement. New cells are placed near the corners."
1835 width=
"200" height=
"200">
1838 <img src=
"https://www.dealii.org/images/steps/developer/step-27-mesh-4.svg"
1839 alt=
"Triangulation containing reentrant corners with four levels of
1840 refinement. New cells are placed near the corners."
1841 width=
"200" height=
"200">
1844 <img src=
"https://www.dealii.org/images/steps/developer/step-27-mesh-5.svg"
1845 alt=
"Triangulation containing reentrant corners with five levels of
1846 refinement. New cells are placed near the corners."
1847 width=
"200" height=
"200">
1851 It is clearly visible how the mesh is refined near the corner singularities,
1852 as
one would expect it. More interestingly, we should be curious to see the
1853 distribution of finite element polynomial degrees to these mesh cells, where
1854 grey corresponds to degree two and pink corresponds to degree seven:
1856 <div
class=
"threecolumn" style=
"width: 80%">
1858 <img src=
"https://www.dealii.org/images/steps/developer/step-27-cell-degree-0.png"
1859 alt=
"Initial grid where all cells contain just biquadratic functions."
1860 width=
"200" height=
"200">
1863 <img src=
"https://www.dealii.org/images/steps/developer/step-27-cell-degree-1.png"
1864 alt=
"Depiction of local approximation degrees after one refinement."
1865 width=
"200" height=
"200">
1868 <img src=
"https://www.dealii.org/images/steps/developer/step-27-cell-degree-2.png"
1869 alt=
"Depiction of local approximation degrees after two refinements."
1870 width=
"200" height=
"200">
1873 <img src=
"https://www.dealii.org/images/steps/developer/step-27-cell-degree-3.png"
1874 alt=
"Depiction of local approximation degrees after three refinements."
1875 width=
"200" height=
"200">
1878 <img src=
"https://www.dealii.org/images/steps/developer/step-27-cell-degree-4.png"
1879 alt=
"Depiction of local approximation degrees after four refinements."
1880 width=
"200" height=
"200">
1883 <img src=
"https://www.dealii.org/images/steps/developer/step-27-cell-degree-5.png"
1884 alt=
"Depiction of local approximation degrees after five refinements."
1885 width=
"200" height=
"200">
1889 While
this is certainly not a perfect arrangement, it does make some sense: we
1890 use low order elements close to boundaries and corners where regularity is
1891 low. On the other hand, higher order elements are used where (i) the error was
1892 at
one point fairly large, i.e. mainly in the
general area around the corner
1893 singularities and in the top right corner where the solution is large, and
1894 (ii) where the solution is smooth, i.e. far away from the boundary.
1896 This arrangement of polynomial degrees of course follows from our smoothness
1897 estimator. Here is the estimated smoothness of the solution, with darker colors
1898 indicating least smoothness and lighter indicating the smoothest areas:
1900 <div
class=
"threecolumn" style=
"width: 80%">
1902 <img src=
"https://www.dealii.org/images/steps/developer/step-27-smoothness-0.png"
1903 alt=
"Estimated regularity per cell on the initial grid."
1904 width=
"200" height=
"200">
1907 <img src=
"https://www.dealii.org/images/steps/developer/step-27-smoothness-1.png"
1908 alt=
"Depiction of the estimated regularity per cell after one refinement."
1909 width=
"200" height=
"200">
1912 <img src=
"https://www.dealii.org/images/steps/developer/step-27-smoothness-2.png"
1913 alt=
"Depiction of the estimated regularity per cell after two refinements."
1914 width=
"200" height=
"200">
1917 <img src=
"https://www.dealii.org/images/steps/developer/step-27-smoothness-3.png"
1918 alt=
"Depiction of the estimated regularity per cell after three refinements."
1919 width=
"200" height=
"200">
1922 <img src=
"https://www.dealii.org/images/steps/developer/step-27-smoothness-4.png"
1923 alt=
"Depiction of the estimated regularity per cell after four refinements."
1924 width=
"200" height=
"200">
1927 <img src=
"https://www.dealii.org/images/steps/developer/step-27-smoothness-5.png"
1928 alt=
"Depiction of the estimated regularity per cell after five refinements."
1929 width=
"200" height=
"200">
1933 The primary conclusion
one can draw from
this is that the loss of regularity at
1934 the
internal corners is a highly localized phenomenon; it only seems to impact
1935 the cells adjacent to the corner itself, so when we
refine the mesh the black
1936 coloring is no longer visible. Besides the corners,
this sequence of plots
1937 implies that the smoothness estimates are somewhat independent of the mesh
1938 refinement, particularly when we are far away from boundaries.
1939 It is also obvious that the smoothness estimates are independent of the actual
1940 size of the solution (see the picture of the solution above), as it should be.
1941 A
point of larger concern, however, is that
one realizes on closer inspection
1942 that the estimator we have overestimates the smoothness of the solution on
1943 cells with hanging nodes. This in turn leads to higher polynomial degrees in
1944 these areas, skewing the allocation of finite elements onto cells.
1946 We have no good explanation
for this effect at the moment. One theory is that
1947 the numerical solution on cells with hanging nodes is, of course, constrained
1948 and therefore not entirely
free to explore the
function space to get close to
1949 the exact solution. This lack of degrees of freedom may manifest itself by
1950 yielding numerical solutions on these cells with suppressed oscillation,
1951 meaning a higher degree of smoothness. The estimator picks
this signal up and
1952 the estimated smoothness overestimates the actual
value. However, a definite
1953 answer to what is going on currently eludes the authors of
this program.
1955 The bigger question is, of course, how to avoid
this problem. Possibilities
1956 include estimating the smoothness not on single cells, but cell assemblies or
1957 patches surrounding each cell. It may also be possible to find simple
1958 correction factors
for each cell depending on the number of constrained
1959 degrees of freedom it has. In either
case, there are ample opportunities
for
1960 further research on finding good @f$hp@f$ refinement criteria. On the other hand,
1961 the main
point of the current program was to demonstrate
using the @f$hp@f$
1962 technology in deal.II, which is unaffected by our use of a possible
1963 sub-optimal refinement criterion.
1966 <a name=
"PlainProg"></a>
1967 <h1> The plain program</h1>
1968 @include
"step-27.cc"