From 71c0e410f058221a0a8620cda160be695f1e4b96 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Sun, 1 Jul 2018 23:49:09 +0200 Subject: [PATCH] examples/step-27: Update indenting and modernize --- examples/step-27/doc/intro.dox | 20 ++--- examples/step-27/step-27.cc | 144 +++++++++++++-------------------- 2 files changed, 62 insertions(+), 102 deletions(-) diff --git a/examples/step-27/doc/intro.dox b/examples/step-27/doc/intro.dox index cc1fa99ef8..8e9de94d7c 100644 --- a/examples/step-27/doc/intro.dox +++ b/examples/step-27/doc/intro.dox @@ -149,9 +149,7 @@ inactive on it. The general outline of this reads like this: @code hp::DoFHandler dof_handler (triangulation); - for (typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) + for (auto &cell: dof_handler.active_cell_iterators()) cell->set_active_fe_index (...); dof_handler.distribute_dofs (fe_collection); @endcode @@ -171,9 +169,8 @@ finite element field to ensure that it is continuous. This is conceptually very similar to how we compute hanging node constraints, and in fact the code looks exactly the same: @code - ConstraintMatrix constraints; - DoFTools::make_hanging_node_constraints (dof_handler, - constraints); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints (dof_handler, constraints); @endcode In other words, the DoFTools::make_hanging_node_constraints deals not only with hanging node constraints, but also with $hp$ constraints at @@ -218,10 +215,7 @@ i.e. the interesting part of the loop over all cells would look like this: update_values | update_gradients | update_q_points | update_JxW_values); - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) + for (const auto &cell: dof_handler.active_cell_iterators()) { hp_fe_values.reinit (cell, cell->active_fe_index(), @@ -624,7 +618,7 @@ challenge. Most programs built on deal.II use the DoFTools::make_sparsity_pattern function to allocate the sparsity pattern of a matrix, and later add a few more entries necessary to handle constrained degrees of freedom using -ConstraintMatrix::condense. The sparsity pattern is then compressed using +AffineConstraints::condense. The sparsity pattern is then compressed using SparsityPattern::compress. This method is explained in step-6 and used in most tutorial programs. In order to work, it needs an initial upper estimate for the maximal number of nonzero entries per row, something that can be had @@ -671,7 +665,7 @@ of freedom, but is constrained against many more degrees of freedom. It turns out that the strategy presented first in step-6 to eliminate the constraints while computing the element matrices and vectors with -ConstraintMatrix::distribute_local_to_global is the most efficient approach +AffineConstraints::distribute_local_to_global is the most efficient approach also for this case. The alternative strategy to first build the matrix without constraints and then "condensing" away constrained degrees of freedom is considerably more expensive. It turns out that building the sparsity pattern @@ -685,7 +679,7 @@ In our program, we will also treat the boundary conditions as (possibly inhomogeneous) constraints and eliminate the matrix rows and columns to those as well. All we have to do for this is to call the function that interpolates the Dirichlet boundary conditions already in the setup phase in -order to tell the ConstraintMatrix object about them, and then do the +order to tell the AffineConstraints object about them, and then do the transfer from local to global data on matrix and vector simultaneously. This is exactly what we've shown in step-6. diff --git a/examples/step-27/step-27.cc b/examples/step-27/step-27.cc index c0ea4491b7..81120240d1 100644 --- a/examples/step-27/step-27.cc +++ b/examples/step-27/step-27.cc @@ -115,7 +115,7 @@ namespace Step27 std::vector ln_k; Table> fourier_coefficients; - ConstraintMatrix constraints; + AffineConstraints constraints; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; @@ -328,10 +328,7 @@ namespace Step27 std::vector local_dof_indices; - typename hp::DoFHandler::active_cell_iterator cell = dof_handler - .begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; @@ -354,11 +351,13 @@ namespace Step27 { for (unsigned int j = 0; j < dofs_per_cell; ++j) cell_matrix(i, j) += - (fe_values.shape_grad(i, q_point) * - fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point)); + (fe_values.shape_grad(i, q_point) * // grad phi_i(x_q) + fe_values.shape_grad(j, q_point) * // grad phi_j(x_q) + fe_values.JxW(q_point)); // dx - cell_rhs(i) += (fe_values.shape_value(i, q_point) * - rhs_values[q_point] * fe_values.JxW(q_point)); + cell_rhs(i) += (fe_values.shape_value(i, q_point) * // phi_i(x_q) + rhs_values[q_point] * // f(x_q) + fe_values.JxW(q_point)); // dx } local_dof_indices.resize(dofs_per_cell); @@ -438,14 +437,9 @@ namespace Step27 // all integers, so that it what we use: { Vector fe_degrees(triangulation.n_active_cells()); - { - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) - fe_degrees(cell->active_cell_index()) = - fe_collection[cell->active_fe_index()].degree; - } + for (const auto &cell : dof_handler.active_cell_iterators()) + fe_degrees(cell->active_cell_index()) = + fe_collection[cell->active_fe_index()].degree; // With now all data vectors available -- solution, estimated errors and // smoothness indicators, and finite element degrees --, we create a @@ -499,21 +493,16 @@ namespace Step27 smoothness_indicators.end()), min_smoothness = *std::max_element(smoothness_indicators.begin(), smoothness_indicators.end()); - { - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) - if (cell->refine_flag_set()) - { - max_smoothness = - std::max(max_smoothness, - smoothness_indicators(cell->active_cell_index())); - min_smoothness = - std::min(min_smoothness, - smoothness_indicators(cell->active_cell_index())); - } - } + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->refine_flag_set()) + { + max_smoothness = + std::max(max_smoothness, + smoothness_indicators(cell->active_cell_index())); + min_smoothness = + std::min(min_smoothness, + smoothness_indicators(cell->active_cell_index())); + } const float threshold_smoothness = (max_smoothness + min_smoothness) / 2; // With this, we can go back, loop over all cells again, and for those @@ -524,20 +513,15 @@ namespace Step27 // degree and in return remove the flag indicating that the cell should // undergo bisection. For all other cells, the refinement flags remain // untouched: - { - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) - if (cell->refine_flag_set() && - (smoothness_indicators(cell->active_cell_index()) > - threshold_smoothness) && - (cell->active_fe_index() + 1 < fe_collection.size())) - { - cell->clear_refine_flag(); - cell->set_active_fe_index(cell->active_fe_index() + 1); - } - } + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->refine_flag_set() && + (smoothness_indicators(cell->active_cell_index()) > + threshold_smoothness) && + (cell->active_fe_index() + 1 < fe_collection.size())) + { + cell->clear_refine_flag(); + cell->set_active_fe_index(cell->active_fe_index() + 1); + } // At the end of this procedure, we then refine the mesh. During this // process, children of cells undergoing bisection inherit their mother @@ -561,43 +545,28 @@ namespace Step27 { const unsigned int dim = 2; - static const Point<2> vertices_1[] = { - Point<2>(-1., -1.), Point<2>(-1. / 2, -1.), - Point<2>(0., -1.), Point<2>(+1. / 2, -1.), - Point<2>(+1, -1.), - - Point<2>(-1., -1. / 2.), Point<2>(-1. / 2, -1. / 2.), - Point<2>(0., -1. / 2.), Point<2>(+1. / 2, -1. / 2.), - Point<2>(+1, -1. / 2.), - - Point<2>(-1., 0.), Point<2>(-1. / 2, 0.), - Point<2>(+1. / 2, 0.), Point<2>(+1, 0.), - - Point<2>(-1., 1. / 2.), Point<2>(-1. / 2, 1. / 2.), - Point<2>(0., 1. / 2.), Point<2>(+1. / 2, 1. / 2.), - Point<2>(+1, 1. / 2.), - - Point<2>(-1., 1.), Point<2>(-1. / 2, 1.), - Point<2>(0., 1.), Point<2>(+1. / 2, 1.), - Point<2>(+1, 1.)}; - const unsigned int n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]); - const std::vector> vertices(&vertices_1[0], - &vertices_1[n_vertices]); - static const int cell_vertices[][GeometryInfo::vertices_per_cell] = { - {0, 1, 5, 6}, - {1, 2, 6, 7}, - {2, 3, 7, 8}, - {3, 4, 8, 9}, - {5, 6, 10, 11}, - {8, 9, 12, 13}, - {10, 11, 14, 15}, - {12, 13, 17, 18}, - {14, 15, 19, 20}, - {15, 16, 20, 21}, - {16, 17, 21, 22}, - {17, 18, 22, 23}}; - const unsigned int n_cells = - sizeof(cell_vertices) / sizeof(cell_vertices[0]); + const std::vector> vertices = { + {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0}, // + {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5}, // + {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0}, // + {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5}, // + {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}}; + + const std::vector::vertices_per_cell>> + cell_vertices = {{0, 1, 5, 6}, + {1, 2, 6, 7}, + {2, 3, 7, 8}, + {3, 4, 8, 9}, + {5, 6, 10, 11}, + {8, 9, 12, 13}, + {10, 11, 14, 15}, + {12, 13, 17, 18}, + {14, 15, 19, 20}, + {15, 16, 20, 21}, + {16, 17, 21, 22}, + {17, 18, 22, 23}}; + + const unsigned int n_cells = cell_vertices.size(); std::vector> cells(n_cells, CellData()); for (unsigned int i = 0; i < n_cells; ++i) @@ -698,10 +667,7 @@ namespace Step27 Vector local_dof_values; // Then here is the loop: - typename hp::DoFHandler::active_cell_iterator cell = dof_handler - .begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { // Inside the loop, we first need to get the values of the local // degrees of freedom (which we put into the @@ -752,8 +718,8 @@ namespace Step27 // We have to calculate the logarithms of absolute // values of coefficients and use it in linear regression fit to // obtain $\mu$. - for (unsigned int f = 0; f < res.second.size(); f++) - res.second[f] = std::log(res.second[f]); + for (double &f : res.second) + f = std::log(f); std::pair fit = FESeries::linear_regression(ln_k, res.second); -- 2.39.5