From: Wolfgang Bangerth Date: Fri, 17 May 2019 14:43:33 +0000 (-0600) Subject: Updates to the step-63 program. X-Git-Tag: v9.2.0-rc1~1447^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d485e60608978e554e7d8f58505f669f03c8ad58;p=dealii.git Updates to the step-63 program. --- diff --git a/examples/step-63/step-63.cc b/examples/step-63/step-63.cc index d963fa1aba..a97292ca7e 100644 --- a/examples/step-63/step-63.cc +++ b/examples/step-63/step-63.cc @@ -74,19 +74,23 @@ #include +// @sect3{MeshWorker data} + +// As always, we will be putting everything related to this program +// into a namespace of its own. +// +// Since we will be using the MeshWorker framework, the first step is +// to define the following structures needed by the assemble_cell() +// function used by MeshWorker::mesh_loop(): `ScratchData` +// contains an FEValues object which is needed for assembling +// a cell's local contribution, while `CopyData` contains the +// output from a cell's local contribution and necessary information +// to copy that to the global system. (Their purpose is also explained +// in the documentation of the WorkStream class.) namespace Step63 { using namespace dealii; - // @sect3{MeshWorker Data} - - // The following are structures needed by the assemble_cell() - // function used by Meshworker::mesh_loop(). ScratchData - // contains an FeValues object which is needed for assembling - // a cell's local contribution, while CopyData contains the - // output from a cell's local contribution and necessary information - // to copy that to the global system. - template struct ScratchData { @@ -126,10 +130,12 @@ namespace Step63 // @sect3{Problem parameters} + // The second step is to define the classes that deal with run-time + // parameters to be read from an input file. + // // We will use ParameterHandler to pass in parameters at runtime. The - // structure Settings parses and stores these parameters to be queried + // structure `Settings` parses and stores the parameters to be queried // throughout the program. - struct Settings { enum DoFRenumberingStrategy @@ -151,8 +157,11 @@ namespace Step63 bool output; }; + + void Settings::get_parameters(const std::string &prm_filename) { + /* First declare the parameters... */ ParameterHandler prm; prm.declare_entry("Epsilon", @@ -186,11 +195,12 @@ namespace Step63 Patterns::Bool(), "Generate graphical output: true|false"); + /* ...and then try to read their values from the input file: */ if (prm_filename.empty()) { prm.print_parameters(std::cout, ParameterHandler::Text); AssertThrow( - false, ExcMessage("please pass a .prm file as the first argument!")); + false, ExcMessage("Please pass a .prm file as the first argument!")); } prm.parse_input(prm_filename); @@ -209,6 +219,10 @@ namespace Step63 dof_renumbering = DoFRenumberingStrategy::upstream; else if (renumbering == "random") dof_renumbering = DoFRenumberingStrategy::random; + else + AssertThrow(false, + ExcMessage("The parameter has " + "an invalid value.")); with_streamline_diffusion = prm.get_bool("With streamline diffusion"); output = prm.get_bool("Output"); @@ -221,12 +235,19 @@ namespace Step63 // will play a role in the speed of convergence for multiplicative // methods. Here we define functions which return a specific ordering // of cells to be used by the block smoothers. - - // For each type of cell ordering, we define a function for the active - // mesh and one for a level mesh. While the only reordering necessary - // for solving the system will be on the level meshes, we include the - // active reordering for visualization purposes in output_results(). - + // + // For each type of cell ordering, we define a function for the + // active mesh and one for a level mesh (i.e., for the cells at one + // level of a multigrid hierarchy). While the only reordering + // necessary for solving the system will be on the level meshes, we + // include the active reordering for visualization purposes in + // output_results(). + // + // For the two downstream ordering functions, we first create an + // array with all of the relevant cells that we then sort in + // downstream direction using a "comparator" object. The output of + // the functions is then simply an array of the indices of the cells + // in the just computed order. template std::vector create_downstream_cell_ordering(const DoFHandler &dof_handler, @@ -235,13 +256,12 @@ namespace Step63 { std::vector::level_cell_iterator> ordered_cells; ordered_cells.reserve(dof_handler.get_triangulation().n_cells(level)); - const DoFRenumbering:: - CompareDownstream::level_cell_iterator, dim> - comparator(direction); - for (const auto &cell : dof_handler.cell_iterators_on_level(level)) ordered_cells.push_back(cell); + const DoFRenumbering:: + CompareDownstream::level_cell_iterator, dim> + comparator(direction); std::sort(ordered_cells.begin(), ordered_cells.end(), comparator); std::vector ordered_indices; @@ -253,6 +273,8 @@ namespace Step63 return ordered_indices; } + + template std::vector create_downstream_cell_ordering(const DoFHandler &dof_handler, @@ -260,13 +282,12 @@ namespace Step63 { std::vector::active_cell_iterator> ordered_cells; ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells()); - const DoFRenumbering:: - CompareDownstream::active_cell_iterator, dim> - comparator(direction); - for (const auto &cell : dof_handler.active_cell_iterators()) ordered_cells.push_back(cell); + const DoFRenumbering:: + CompareDownstream::active_cell_iterator, dim> + comparator(direction); std::sort(ordered_cells.begin(), ordered_cells.end(), comparator); std::vector ordered_indices; @@ -278,28 +299,30 @@ namespace Step63 return ordered_indices; } + + // The functions that produce a random ordering are similar in + // spirit in that they first put information about all cells into an + // array. But then, instead of sorting them, they shuffle the + // elements randomly using the facilities C++ offers to generate + // random numbers. The way this is done is by iterating over all + // elements of the array, drawing a random number for another + // element before that, and then exchanging these elements. The + // result is a random shuffle of the elements of the array. template std::vector create_random_cell_ordering(const DoFHandler &dof_handler, const unsigned int level) { - const unsigned int n_cells = dof_handler.get_triangulation().n_cells(level); - std::vector ordered_cells; - ordered_cells.reserve(n_cells); - + ordered_cells.reserve(dof_handler.get_triangulation().n_cells(level)); for (const auto &cell : dof_handler.cell_iterators_on_level(level)) ordered_cells.push_back(cell->index()); - // Shuffle the elements: std::mt19937 random_number_generator; - for (unsigned int i = 1; i < n_cells; ++i) + for (unsigned int i = 1; i < ordered_cells.size(); ++i) { - // Get a random number between 0 and i (inclusive): const unsigned int j = std::uniform_int_distribution<>(0, i)(random_number_generator); - - // If possible, swap the elements: if (i != j) std::swap(ordered_cells[i], ordered_cells[j]); } @@ -307,28 +330,22 @@ namespace Step63 return ordered_cells; } + + template std::vector create_random_cell_ordering(const DoFHandler &dof_handler) { - const unsigned int n_cells = - dof_handler.get_triangulation().n_active_cells(); - std::vector ordered_cells; - ordered_cells.reserve(n_cells); - + ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells()); for (const auto &cell : dof_handler.active_cell_iterators()) ordered_cells.push_back(cell->index()); - // Shuffle the elements: std::mt19937 random_number_generator; - for (unsigned int i = 1; i < n_cells; ++i) + for (unsigned int i = 1; i < ordered_cells.size(); ++i) { - // Get a random number between 0 and i (inclusive): const unsigned int j = std::uniform_int_distribution<>(0, i)(random_number_generator); - - // If possible, swap the elements: if (i != j) std::swap(ordered_cells[i], ordered_cells[j]); } @@ -337,7 +354,7 @@ namespace Step63 } - // @sect3{Right-hand Side and Boundary Values} + // @sect3{Right-hand side and boundary values} // The problem solved in this tutorial is an adaptation of Ex. 3.1.3 found // on pg. 118 of . The // main difference being that we add a hole in the center of our domain with - // zero Dirichlet boundary. - - // We have a zero right-hand side: + // zero Dirichlet boundary conditions. + // + // For a complete description, we need classes that implement the + // zero right-hand side first (we could of course have just used + // Functions::ZeroFunction): template class RightHandSide : public Function { @@ -364,6 +383,8 @@ namespace Step63 const unsigned int component = 0) const override; }; + + template double RightHandSide::value(const Point &, const unsigned int component) const @@ -389,9 +410,9 @@ namespace Step63 } - // We have Dirichlet boundary conditions. On a connected portion of the + // We also have Dirichlet boundary conditions. On a connected portion of the // outer, square boundary we set the value to 1, and we set the value to 0 - // everywhere else (including the inner, circular boundary). + // everywhere else (including the inner, circular boundary): template class BoundaryValues : public Function { @@ -409,6 +430,7 @@ namespace Step63 }; + template double BoundaryValues::value(const Point & p, const unsigned int component) const @@ -444,13 +466,14 @@ namespace Step63 - // @sect3{Streamline Diffusion} + // @sect3{Streamline diffusion} - // Streamline diffusion stabilization constant. Parameter design is taken - // from On - // Discontinuity-Capturing Methods for Convection-Diffusion Equations by - // Volker John and Petr Knobloch. + // The streamline diffusion method has a stabilization constant that + // we need to be able to compute. The choice of how this parameter + // is computed is taken from On + // Discontinuity-Capturing Methods for Convection-Diffusion + // Equations by Volker John and Petr Knobloch. template double compute_stabilization_delta(const double hk, const double eps, @@ -469,10 +492,10 @@ namespace Step63 // This is the main class of the program, and should look very similar to // step-16. The major difference is that, since we are defining our multigrid - // smoother at runtime, we choose to define a function create_smoother() and a - // class object mg_smoother which is a std::unique_ptr to a smoother that is - // derived from MGSmoother. Note that for smoothers derived from - // RelaxationBlock, we must include a smoother_data object for each level. + // smoother at runtime, we choose to define a function `create_smoother()` and + // a class object `mg_smoother` which is a `std::unique_ptr` to a smoother + // that is derived from MGSmoother. Note that for smoothers derived from + // RelaxationBlock, we must include a `smoother_data` object for each level. // This will contain information about the cell ordering and the method of // inverting cell matrices. @@ -548,18 +571,22 @@ namespace Step63 , settings(settings) { advection_direction[0] = -std::sin(numbers::PI / 6.0); - if (dim > 1) + if (dim >= 2) advection_direction[1] = std::cos(numbers::PI / 6.0); - if (dim > 2) + if (dim >= 3) AssertThrow(false, ExcNotImplemented()); } - // @sect4{AdvectionProblem::setup_system} - - // Here we set up the DoFHandler, ConstraintMatrix, and sparsity patterns for - // both active and multigrid level meshes. + // @sect4{AdvectionProblem::setup_system()} + // Here we first set up the DoFHandler, AffineConstraints, and + // SparsityPattern objects for both active and multigrid level meshes. + // + // We could renumber the active DoFs with the DoFRenumbering class, + // but the smoothers only act on multigrid levels and as such, this + // would not matter for the computations. Instead, we will renumber the + // DoFs on each multigrid level below. template void AdvectionProblem::setup_system() { @@ -567,11 +594,6 @@ namespace Step63 dof_handler.distribute_dofs(fe); - // We could renumber the active DoFs with the DoFRenumbering class - // here, but the smoothers only act on multigrid levels and as such, this - // would not matter for the computations. Instead, we will renumber the - // DoFs on each multigrid level below. - solution.reinit(dof_handler.n_dofs()); system_rhs.reinit(dof_handler.n_dofs()); @@ -595,9 +617,15 @@ namespace Step63 dof_handler.distribute_mg_dofs(); - // Renumber DoFs on each level in downstream or upstream direction if - // needed. This is only necessary for point smoothers (SOR and Jacobi) as - // the block smoothers operate on cells (see create_smoother()): + // Having enumerated the global degrees of freedom as well as (in + // the last line above) the level degrees of freedom, let us + // renumber the level degrees of freedom to get a better smoother + // as explained in the introduction. The first block below + // renumbers DoFs on each level in downstream or upstream + // direction if needed. This is only necessary for point smoothers + // (SOR and Jacobi) as the block smoothers operate on cells (see + // `create_smoother()`). The blocks below then also implement + // random numbering. if (settings.smoother_type == "SOR" || settings.smoother_type == "Jacobi") { if (settings.dof_renumbering == @@ -628,6 +656,10 @@ namespace Step63 Assert(false, ExcNotImplemented()); } + // The rest of the function just sets up data structures. The last + // lines of the code below is unlike the other GMG tutorials, as + // it sets up both the interface in and out matrices. We need this + // since our problem is non-symmetric. mg_constrained_dofs.clear(); mg_constrained_dofs.initialize(dof_handler); @@ -660,8 +692,6 @@ namespace Step63 level); mg_interface_sparsity_patterns[level].copy_from(dsp); - // Unlike the other GMG tutorials, we need both interface in and out - // matrices since our problem is non-symmetric. mg_interface_in[level].reinit(mg_interface_sparsity_patterns[level]); mg_interface_out[level].reinit(mg_interface_sparsity_patterns[level]); } @@ -669,12 +699,13 @@ namespace Step63 } - // @sect4{AdvectionProblem::assemble_cell} + // @sect4{AdvectionProblem::assemble_cell()} - // Here we define the assembly of the linear system on each cell to be used by - // the mesh_loop() function below. This one function assembles the cell matrix - // for both and active and a level cell, and only assembles a right-hand side - // if called for an active cell. + // Here we define the assembly of the linear system on each cell to + // be used by the mesh_loop() function below. This one function + // assembles the cell matrix for either an active or a level cell + // (whatever it is passed as its first argument), and only assembles + // a right-hand side if called with an active cell. template template @@ -711,20 +742,21 @@ namespace Step63 // using streamline diffusion, setting $\delta=0$ negates this contribution // below and we are left with the standard, Galerkin finite element // assembly. - const double delta = settings.with_streamline_diffusion ? - compute_stabilization_delta(cell->diameter(), - settings.epsilon, - advection_direction, - settings.fe_degree) : - 0.0; + const double delta = (settings.with_streamline_diffusion ? + compute_stabilization_delta(cell->diameter(), + settings.epsilon, + advection_direction, + settings.fe_degree) : + 0.0); for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { + // The assembly of the local matrix has two parts. First + // the Galerkin contribution: copy_data.cell_matrix(i, j) += - // Galerkin contribution: (settings.epsilon * scratch_data.fe_values.shape_grad(i, q_point) * scratch_data.fe_values.shape_grad(j, q_point) * @@ -732,14 +764,14 @@ namespace Step63 (scratch_data.fe_values.shape_value(i, q_point) * (advection_direction * scratch_data.fe_values.shape_grad(j, q_point)) * - scratch_data.fe_values.JxW(q_point)) + - // Streamline diffusion contribution: - delta * - (advection_direction * - scratch_data.fe_values.shape_grad(j, q_point)) * - (advection_direction * - scratch_data.fe_values.shape_grad(i, q_point)) * - scratch_data.fe_values.JxW(q_point) - + scratch_data.fe_values.JxW(q_point)) + // and then the streamline diffusion contribution: + + delta * + (advection_direction * + scratch_data.fe_values.shape_grad(j, q_point)) * + (advection_direction * + scratch_data.fe_values.shape_grad(i, q_point)) * + scratch_data.fe_values.JxW(q_point) - delta * settings.epsilon * trace(scratch_data.fe_values.shape_hessian(j, q_point)) * (advection_direction * @@ -748,20 +780,21 @@ namespace Step63 } if (cell->is_level_cell() == false) { + // The same applies to the right hand side. First the + // Galerkin contribution: copy_data.cell_rhs(i) += - // Galerkin contribution: scratch_data.fe_values.shape_value(i, q_point) * - rhs_values[q_point] * scratch_data.fe_values.JxW(q_point) + - // Streamline diffusion contribution: - delta * rhs_values[q_point] * advection_direction * - scratch_data.fe_values.shape_grad(i, q_point) * - scratch_data.fe_values.JxW(q_point); + rhs_values[q_point] * scratch_data.fe_values.JxW(q_point) + // and then the streamline diffusion contribution: + + delta * rhs_values[q_point] * advection_direction * + scratch_data.fe_values.shape_grad(i, q_point) * + scratch_data.fe_values.JxW(q_point); } } } - // @sect4{AdvectionProblem::assemble_system_and_multigrid} + // @sect4{AdvectionProblem::assemble_system_and_multigrid()} // Here we employ MeshWorker::mesh_loop() to go over cells and assemble the // system_matrix, system_rhs, and all mg_matrices for us. @@ -769,16 +802,14 @@ namespace Step63 template void AdvectionProblem::assemble_system_and_multigrid() { - auto cell_worker_active = + const auto cell_worker_active = [&](const decltype(dof_handler.begin_active()) &cell, ScratchData & scratch_data, CopyData & copy_data) { this->assemble_cell(cell, scratch_data, copy_data); }; - - - auto copier_active = [&](const CopyData ©_data) { + const auto copier_active = [&](const CopyData ©_data) { constraints.distribute_local_to_global(copy_data.cell_matrix, copy_data.cell_rhs, copy_data.local_dof_indices, @@ -814,25 +845,28 @@ namespace Step63 boundary_constraints[level].close(); } - auto cell_worker_mg = [&](const decltype(dof_handler.begin_mg()) &cell, - ScratchData &scratch_data, - CopyData & copy_data) { - this->assemble_cell(cell, scratch_data, copy_data); - }; + const auto cell_worker_mg = + [&](const decltype(dof_handler.begin_mg()) &cell, + ScratchData & scratch_data, + CopyData & copy_data) { + this->assemble_cell(cell, scratch_data, copy_data); + }; - auto copier_mg = [&](const CopyData ©_data) { + const auto copier_mg = [&](const CopyData ©_data) { boundary_constraints[copy_data.level].distribute_local_to_global( copy_data.cell_matrix, copy_data.local_dof_indices, mg_matrices[copy_data.level]); - // If (i,j) is an interface_out dof pair, then (j,i) is an interface_in - // dof pair. Note: for interface_in, we load the transpose of the - // interface entries, i.e., the entry for dof pair (j,i) is stored in - // interface_in(i,j). This is an optimization for the symmetric case - // which allows only one matrix to be used when setting the edge_matrices - // in solve(). Here, however, since our problem is non-symmetric, we must - // store both interface_in and interface_out matrices. + // If $(i,j)$ is an `interface_out` dof pair, then $(j,i)$ is an + // `interface_in` dof pair. Note: For `interface_in`, we load + // the transpose of the interface entries, i.e., the entry for + // dof pair $(j,i)$ is stored in `interface_in(i,j)`. This is an + // optimization for the symmetric case which allows only one + // matrix to be used when setting the edge_matrices in + // solve(). Here, however, since our problem is non-symmetric, + // we must store both `interface_in` and `interface_out` + // matrices. for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i) for (unsigned int j = 0; j < copy_data.dofs_per_cell; ++j) if (mg_constrained_dofs.is_interface_matrix_entry( @@ -861,11 +895,12 @@ namespace Step63 } - // @sect4{AdvectionProblem::setup_smoother} + // @sect4{AdvectionProblem::setup_smoother()} - // Here we set up the smoother based on the settings in the .prm. The two - // options that are of significance is the number of pre- and post-smoothing - // steps on each level of the multigrid v-cycle and the relaxation parameter. + // Next, we set up the smoother based on the settings in the `.prm` file. The + // two options that are of significance is the number of pre- and + // post-smoothing steps on each level of the multigrid v-cycle and the + // relaxation parameter. // Since multiplicative methods tend to be more powerful than additive method, // fewer smoothing steps are required to see convergence indepedent of mesh @@ -874,21 +909,21 @@ namespace Step63 // smoother below. // The relaxation parameter for point smoothers is chosen based on trial and - // error, and they reflect values necessary to keep the iteration counts in + // error, and reflects values necessary to keep the iteration counts in // the GMRES solve constant (or as close as possible) as we refine the mesh. - // The two values given for both "Jacobi" and "SOR" in the .prm files are for - // degree 1 and degree 3 finite elements. If the user wants to change to + // The two values given for both "Jacobi" and "SOR" in the `.prm` files are + // for degree 1 and degree 3 finite elements. If the user wants to change to // another degree, they may need to adjust these numbers. For block smoothers, // this parameter has a more straightforward interpretation, namely that for // additive methods in 2D, a DoF can have a repeated contribution from up to 4 // cells, therefore we must relax these methods by 0.25 to compensate. This is - // not an issue for multiplicative methods as each cell inverse application + // not an issue for multiplicative methods as each cell's inverse application // carries new information to all its DoFs. // Finally, as mentioned above, the point smoothers only operate on DoFs, and // the block smoothers on cells, so only the block smoothers need to be given // information regarding cell orderings. DoF ordering for point smoothers has - // already been taken care of in setup_system(). + // already been taken care of in `setup_system()`. template void AdvectionProblem::setup_smoother() @@ -938,8 +973,6 @@ namespace Step63 std::vector ordered_indices; switch (settings.dof_renumbering) { - // Order the cells downstream with respect - // to the advection direction. case Settings::DoFRenumberingStrategy::downstream: ordered_indices = create_downstream_cell_ordering(dof_handler, @@ -947,9 +980,6 @@ namespace Step63 level); break; - // Order the cells upstream with respect to the advection - // direction, i.e., downstream with respect to the negative - // of the advection direction. case Settings::DoFRenumberingStrategy::upstream: ordered_indices = create_downstream_cell_ordering(dof_handler, @@ -957,13 +987,11 @@ namespace Step63 level); break; - // Order the cells randomly. case Settings::DoFRenumberingStrategy::random: ordered_indices = create_random_cell_ordering(dof_handler, level); break; - // Keep the default cell ordering (z-order, see Glossary). case Settings::DoFRenumberingStrategy::none: break; @@ -1004,34 +1032,38 @@ namespace Step63 } - // @sect4{AdvectionProblem::solve} + // @sect4{AdvectionProblem::solve()} // Before we can solve the system, we must first set up the multigrid // preconditioner. This requires the setup of the transfer between levels, // the coarse matrix solver, and the smoother. This setup follows almost // identically to Step-16, the main difference being the various smoothers // defined above and the fact that we need different interface edge matrices - // for in and out since our problem is non-symetric. (In reality, for this + // for in and out since our problem is non-symmetric. (In reality, for this // tutorial these interface matrices are empty since we are only using global // refinement, and thus have no refinement edges. However, we have still // included both here since if one made the simple switch to an adaptively // refined method, the program would still run correctly.) - // The last thing to note is that since our problem is non-symetric, we must + // The last thing to note is that since our problem is non-symmetric, we must // use an appropriate Krylov subspace method. We choose here to // use GMRES since it offers the guarantee of residual reduction in each - // iteration. The major disatvantage to GMRES is that, for each iteration, we - // must store an additional temporary vector as well as compute an additional - // scalar product. This requirement is relaxed by using the restarted GMRES + // iteration. The major disavantage of GMRES is that, for each iteration, we + // the number of stored temporary vectors increases by one, and one also needs + // to compute a scalar product with all previously stored vectors. This is + // rather expensive. This requirement is relaxed by using the restarted GMRES // method which puts a cap on the number of vectors we are required to store - // at any one time (here we resart after 50 temporary vectors, or 48 - // iterations). This then has the disatvantage that we lose information we + // at any one time (here we restart after 50 temporary vectors, or 48 + // iterations). This then has the disadvantage that we lose information we // have gathered throughout the iteration and therefore we could see slower - // convergence. However, the goal of this tutorial is to have very low + // convergence. As a consequence, where to restart is a question of balancing + // memory consumption, CPU effort, and convergence speed. + // However, the goal of this tutorial is to have very low // iteration counts by using a powerful GMG preconditioner, so we have picked // the restart length such that all of the results shown below converge prior - // and thus we have a standard GMRES method. If the user is interested, - // another sutaible method offered in deal.II would be BiCGStab. + // to restart happening, and thus we have a standard GMRES method. If the user + // is interested, another sutaible method offered in deal.II would be + // BiCGStab. template void AdvectionProblem::solve() @@ -1084,24 +1116,27 @@ namespace Step63 } - // @sect4{AdvectionProblem::output_results} + // @sect4{AdvectionProblem::output_results()} + // The final function of interest generates graphical output. // Here we output the solution and cell ordering in a .vtu format. + // At the top of the function, we generate an index for each cell to + // visualize the ordering used by the smoothers. Note that we do + // this only for the active cells instead of the levels, where the + // smoothers are actually used. For the point smoothers we renumber + // DoFs instead of cells, so this is only an approximation of what + // happens in reality. Finally, the random ordering is not the + // random ordering we actually use (see `create_smoother()` for that). + // + // The (integer) ordering of cells is then copied into a (floating + // point) vector for graphical output. template void AdvectionProblem::output_results(const unsigned int cycle) const { - // We generate an index for each cell to visualize the ordering used - // by the smoothers. Note that we do this only for the active cells - // instead of the levels, where the smoothers are actually used. For the - // point smoothers we renumber DoFs instead of cells, so this is only an - // approximation of what happens in reality. Finally, the random ordering - // is not the random ordering we actually use (see create_smoother() for - // that). const unsigned int n_active_cells = triangulation.n_active_cells(); Vector cell_indices(n_active_cells); { - // First generate a permutation vector for the cell indices: std::vector ordered_indices; switch (settings.dof_renumbering) { @@ -1131,11 +1166,12 @@ namespace Step63 break; } - // Then copy the permutation in ordered_indices into an output vector: for (unsigned int i = 0; i < n_active_cells; ++i) cell_indices(ordered_indices[i]) = static_cast(i); } + // The remainder of the function is then straightforward, given + // previous tutorial programs: DataOut data_out; data_out.attach_dof_handler(dof_handler); data_out.add_data_vector(solution, "solution"); @@ -1149,12 +1185,18 @@ namespace Step63 } - // @sect4{AdvectionProblem::run} + // @sect4{AdvectionProblem::run()} // As in most tutorials, this function creates/refines the mesh and calls - // the various functions defined above to setup, assemble, solve, and output + // the various functions defined above to set up, assemble, solve, and output // the results. + // In cycle zero, we generate the mesh for the on the square + // [-1,1]^dim with a hole of radius 3/10 units centered + // at the origin. For objects with `manifold_id` equal to one + // (namely, the faces adjacent to the hole), we assign a spherical + // manifold. + template void AdvectionProblem::run() { @@ -1165,15 +1207,11 @@ namespace Step63 if (cycle == 0) { - // We are solving on the square [-1,1]^dim with a hole - // of radius 3/10 units centered at the origin. GridGenerator::hyper_cube_with_cylindrical_hole(triangulation, 0.3, 1.0); - // Set manifold for the inner (curved) boundary. - static const SphericalManifold manifold_description( - Point(0, 0)); + const SphericalManifold manifold_description(Point(0, 0)); triangulation.set_manifold(1, manifold_description); } @@ -1202,12 +1240,12 @@ namespace Step63 // @sect4{The main function} -// Here the main function is like most tutorials. The only interesting bit -// is that we require the user to pass a .prm file as a sole command line -// argument (see Step-19 for a complete discussion of parameter files). If no -// parameter file is given, the program will output the contents of a sample -// parameter file with all default values to the screen that the user can then -// copy and paste into their own .prm file. +// Finally, the main function is like most tutorials. The only +// interesting bit is that we require the user to pass a `.prm` file +// as a sole command line argument. If no parameter file is given, the +// program will output the contents of a sample parameter file with +// all default values to the screen that the user can then copy and +// paste into their own `.prm` file. int main(int argc, char *argv[]) {