From 89e3a811793e3d6a9dd037b77d1aec4b085629d4 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 3 May 2014 10:48:10 +0000 Subject: [PATCH] step-6 uses the new way of dealing with constraints, i.e., applying them on the fly, so mention this in step-27 and step-16. In step-16, also avoid using deprecated interfaces to multigrid data structures git-svn-id: https://svn.dealii.org/trunk@32879 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-16/step-16.cc | 175 ++++++++++++------------- deal.II/examples/step-27/doc/intro.dox | 108 ++------------- deal.II/examples/step-27/step-27.cc | 30 +---- deal.II/examples/step-6/step-6.cc | 12 +- 4 files changed, 105 insertions(+), 220 deletions(-) diff --git a/deal.II/examples/step-16/step-16.cc b/deal.II/examples/step-16/step-16.cc index 1eb7c93fac..e8a2035f0b 100644 --- a/deal.II/examples/step-16/step-16.cc +++ b/deal.II/examples/step-16/step-16.cc @@ -62,15 +62,15 @@ #include #include -// These, now, are the include necessary for the multilevel methods. The -// first two declare classes that allow us to enumerate degrees of freedom not -// only on the finest mesh level, but also on intermediate levels (that's what -// the MGDoFHandler class does) as well as allow to access this information -// (iterators and accessors over these cells). +// These, now, are the include necessary for the multilevel methods. The first +// one declares how to handle Dirichlet boundary conditions on each of the +// levels of the multigrid method. For the actual description of the degrees +// of freedom, we do not need any new include file because DoFHandler already +// has all necessary methods implemented. We will only need to distribute the +// DoFs for the levels further down. // // The rest of the include files deals with the mechanics of multigrid as a // linear operator (solver or preconditioner). -#include #include #include #include @@ -112,7 +112,7 @@ namespace Step16 Triangulation triangulation; FE_Q fe; - MGDoFHandler mg_dof_handler; + DoFHandler dof_handler; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; @@ -234,7 +234,7 @@ namespace Step16 triangulation (Triangulation:: limit_level_difference_at_vertices), fe (degree), - mg_dof_handler (triangulation), + dof_handler (triangulation), degree(degree) {} @@ -243,29 +243,31 @@ namespace Step16 // @sect4{LaplaceProblem::setup_system} // The following function extends what the corresponding one in step-6 - // did. The top part, apart from the additional output, does the same: + // did, with the exception of also distributing the degrees of freedom on + // each level of the mesh which is needed for the multigrid hierarchy. template void LaplaceProblem::setup_system () { - mg_dof_handler.distribute_dofs (fe); + dof_handler.distribute_dofs (fe); + dof_handler.distribute_mg_dofs (fe); // Here we output not only the degrees of freedom on the finest level, but // also in the multilevel structure deallog << "Number of degrees of freedom: " - << mg_dof_handler.n_dofs(); + << dof_handler.n_dofs(); for (unsigned int l=0; l::type dirichlet_boundary_functions; ZeroFunction homogeneous_dirichlet_bc (1); dirichlet_boundary_functions[0] = &homogeneous_dirichlet_bc; - VectorTools::interpolate_boundary_values (static_cast&>(mg_dof_handler), + VectorTools::interpolate_boundary_values (static_cast&>(dof_handler), dirichlet_boundary_functions, constraints); constraints.close (); @@ -299,7 +301,7 @@ namespace Step16 // about the boundary values as well, so we pass the // dirichlet_boundary here as well. mg_constrained_dofs.clear(); - mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary_functions); + mg_constrained_dofs.initialize(dof_handler, dirichlet_boundary_functions); // Now for the things that concern the multigrid data structures. First, @@ -337,9 +339,9 @@ namespace Step16 for (unsigned int level=0; level coefficient; std::vector coefficient_values (n_q_points); - typename MGDoFHandler::active_cell_iterator - cell = mg_dof_handler.begin_active(), - endc = mg_dof_handler.end(); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); for (; cell!=endc; ++cell) { cell_matrix = 0; @@ -455,44 +453,39 @@ namespace Step16 // obscure if you're not familiar with the algorithm actually implemented // in deal.II to support multilevel algorithms on adaptive meshes; if some // of the things below seem strange, take a look at the @ref mg_paper. - // - // Our first job is to identify those degrees of freedom on each level - // that are located on interfaces between adaptively refined levels, and - // those that lie on the interface but also on the exterior boundary of - // the domain. As in many other parts of the library, we do this by using - // Boolean masks, i.e. vectors of Booleans each element of which indicates - // whether the corresponding degree of freedom index is an interface DoF - // or not. The MGConstraints already computed the information - // for us when we called initialize in setup_system(). - std::vector > interface_dofs - = mg_constrained_dofs.get_refinement_edge_indices (); - std::vector > boundary_interface_dofs - = mg_constrained_dofs.get_refinement_edge_boundary_indices (); - - // The indices just identified will later be used to decide where the - // assembled value has to be added into on each level. On the other hand, - // we also have to impose zero boundary conditions on the external - // boundary of each level. But this the MGConstraints knows. - // So we simply ask for them by calling get_boundary_indices(). - // The third step is to construct constraints on all those - // degrees of freedom: their value should be zero after each application - // of the level operators. To this end, we construct ConstraintMatrix - // objects for each level, and add to each of these constraints for each - // degree of freedom. Due to the way the ConstraintMatrix stores its data, - // the function to add a constraint on a single degree of freedom and - // force it to be zero is called Constraintmatrix::add_line(); doing so - // for several degrees of freedom at once can be done using + + // Our first job is to identify the boundary conditions for the levels. On + // each level, we impose Dirichlet boundary conditions on the exterior + // boundary of the domain as well as on interfaces between adaptively + // refined levels. As in many other parts of the library, we do this by + // using a mask described by an IndexSet. The MGConstraints + // already computed the information for us when we called initialize in + // setup_system(). So we simply ask for them by calling + // get_boundary_indices() and + // get_refinement_edge_indices(level). Moreover, we have to + // identify the subset of the refinement edge indices which are also + // located on the boundary as they require special treatment in the + // algorithm further down. + + // These three masks are used to fill a ConstraintMatrix objects for each + // level that we use during the assembly of the matrix: the value of the + // associated degrees of freedom should be zero after each application of + // the level operators. Due to the way the ConstraintMatrix stores its + // data, the function to add a constraint on a single degree of freedom + // and force it to be zero is called Constraintmatrix::add_line(); doing + // so for several degrees of freedom at once can be done using // Constraintmatrix::add_lines(): std::vector boundary_constraints (triangulation.n_levels()); std::vector boundary_interface_constraints (triangulation.n_levels()); for (unsigned int level=0; levelassemble_system, with two exceptions: (i) we don't need a // right hand side, and more significantly (ii) we don't just loop over // all active cells, but in fact all cells, active or not. Consequently, - // the correct iterator to use is MGDoFHandler::cell_iterator rather than - // MGDoFHandler::active_cell_iterator. Let's go about it: - typename MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), - endc = mg_dof_handler.end(); + // the correct iterator to use is DoFHandler::cell_iterator rather than + // DoFHandler::active_cell_iterator. Let's go about it: + typename DoFHandler::cell_iterator cell = dof_handler.begin(), + endc = dof_handler.end(); for (; cell!=endc; ++cell) { @@ -526,8 +519,8 @@ namespace Step16 // with a gotcha that is easily forgotten: The indices of global // degrees of freedom we want here are the ones for current level, not // for the global matrix. We therefore need the function - // MGDoFAccessor::get_mg_dof_indices, not - // MGDoFAccessor::get_dof_indices as used in the assembly of the + // DoFAccessor::get_mg_dof_indices, not + // DoFAccessor::get_dof_indices as used in the assembly of the // global system: cell->get_mg_dof_indices (local_dof_indices); @@ -554,7 +547,7 @@ namespace Step16 // we have a symmetric operator, one of these matrices is the // transpose of the other. // - // The way we assemble these matrices is as follows: since the are + // The way we assemble these matrices is as follows: since they are // formed from parts of the local contributions, we first delete all // those parts of the local contributions that we are not interested // in, namely all those elements of the local matrix for which not $i$ @@ -572,8 +565,10 @@ namespace Step16 // necessary. for (unsigned int i=0; ilevel()][local_dof_indices[i]]==true && - interface_dofs[cell->level()][local_dof_indices[j]]==false)) + if ( !(mg_constrained_dofs.get_refinement_edge_indices(cell->level()). + is_element(local_dof_indices[i])==true && + mg_constrained_dofs.get_refinement_edge_indices(cell->level()). + is_element(local_dof_indices[j])==false)) cell_matrix(i,j) = 0; boundary_interface_constraints[cell->level()] @@ -597,8 +592,11 @@ namespace Step16 // the finite element function spaces involved and can often be computed in // a generic way independent of the problem under consideration. In that // case, we can use the MGTransferPrebuilt class that, given the constraints - // on the global level and an MGDoFHandler object computes the matrices - // corresponding to these transfer operators. + // of the final linear system and the MGConstrainedDoFs object that knows + // about the boundary conditions on the each level and the degrees of + // freedom on interfaces between different refinement level can build the + // matrices for those transfer operations from a DoFHandler object with + // level degrees of freedom. // // The second part of the following lines deals with the coarse grid // solver. Since our coarse grid is very coarse indeed, we decide for a @@ -609,15 +607,8 @@ namespace Step16 template void LaplaceProblem::solve () { - - // Create the object that deals with the transfer between different - // refinement levels. We need to pass it the hanging node constraints. MGTransferPrebuilt > mg_transfer(hanging_node_constraints, mg_constrained_dofs); - // Now the prolongation matrix has to be built. This matrix needs to take - // the boundary values on each level into account and needs to know about - // the indices at the refinement edges. The MGConstraints - // knows about that so pass it as an argument. - mg_transfer.build_matrices(mg_dof_handler); + mg_transfer.build_matrices(dof_handler); FullMatrix coarse_matrix; coarse_matrix.copy_from (mg_matrices[0]); @@ -668,13 +659,13 @@ namespace Step16 // the transpose operator for the latter operation, allowing us to // initialize both up and down versions of the operator with the matrices // we already built: - MGMatrix<> mg_matrix(&mg_matrices); - MGMatrix<> mg_interface_up(&mg_interface_matrices); - MGMatrix<> mg_interface_down(&mg_interface_matrices); + mg::Matrix > mg_matrix(mg_matrices); + mg::Matrix > mg_interface_up(mg_interface_matrices); + mg::Matrix > mg_interface_down(mg_interface_matrices); // Now, we are ready to set up the V-cycle operator and the multilevel // preconditioner. - Multigrid > mg(mg_dof_handler, + Multigrid > mg(dof_handler, mg_matrix, coarse_grid_solver, mg_transfer, @@ -683,7 +674,7 @@ namespace Step16 mg.set_edge_matrices(mg_interface_down, mg_interface_up); PreconditionMG, MGTransferPrebuilt > > - preconditioner(mg_dof_handler, mg, mg_transfer); + preconditioner(dof_handler, mg, mg_transfer); // With all this together, we can finally get about solving the linear // system in the usual way: @@ -709,9 +700,7 @@ namespace Step16 // computed. In particular, the first one refines the mesh at the beginning // of each cycle while the second one outputs results at the end of each // such cycle. The functions are almost unchanged from those in step-6, with - // the exception of two minor differences: The KellyErrorEstimator::estimate - // function wants an argument of type DoFHandler, not MGDoFHandler, and so - // we have to cast from derived to base class; and we generate output in VTK + // the exception of one minor difference: we generate output in VTK // format, to use the more modern visualization programs available today // compared to those that were available when step-6 was written. template @@ -719,7 +708,7 @@ namespace Step16 { Vector estimated_error_per_cell (triangulation.n_active_cells()); - KellyErrorEstimator::estimate (static_cast&>(mg_dof_handler), + KellyErrorEstimator::estimate (dof_handler, QGauss(3), typename FunctionMap::type(), solution, @@ -737,7 +726,7 @@ namespace Step16 { DataOut data_out; - data_out.attach_dof_handler (mg_dof_handler); + data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); data_out.build_patches (); @@ -784,10 +773,10 @@ namespace Step16 setup_system (); std::cout << " Number of degrees of freedom: " - << mg_dof_handler.n_dofs() + << dof_handler.n_dofs() << " (by level: "; for (unsigned int level=0; levelfalse argument we pass when creating the sparsity -pattern. It tells the sparsity pattern that constrained entries should not -be inserted into the sparsity pattern. This makes sense for the way we're -going to do assembly, which does the elimination of constraints already -locally. Hence, we will never write the constrained entries. This is the -second key to increase this program's performance, since the system matrix -will have less entries — this saves both memory and computing time -when doing matrix-vector multiplications or applying a preconditioner. - -In a similar vein, we have to slightly modify the way we copy local -contributions into global matrices and vectors. In previous tutorial programs, -we have always used a process like this: -@code - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - ... // assemble local contributions them into global object - - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - ... // assemble local contributions them into global object - - cell->get_dof_indices (local_dof_indices); - - hanging_node_constraints - .distribute_local_to_global (cell_matrix, cell_rhs, - local_dof_indices, - system_matrix, system_rhs); - } -@endcode -Essentially, what the -hanging_node_constraints.distribute_local_to_global call does is -to implement the same loops as before, but whenever we hit a constrained degree -of freedom, the function does the right thing and already condenses it away. - -Using this technique of eliminating constrained nodes already when -transferring local contributions into the global objects, we avoid the problem -of having to go back later and change these objects. Timing these operations -shows that this makes the overall algorithms faster. +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 +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 +by this inefficient algorithm requires at least ${\cal O}(N \log N)$ in the +number of unknowns, whereas an ideal finite element program would of course +only have algorithms that are linear in the number of unknowns. Timing the +sparsity pattern creation as well as the matrix assembly shows that the +algorithm presented in step-6 (and used in the code below) is indeed faster. In our program, we will also treat the boundary conditions as (possibly inhomogeneous) constraints and eliminate the matrix rows and columns to @@ -783,7 +699,7 @@ 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 transfer from local to global data on matrix and vector simultaneously. This -is exactly what we've shown in the step-22 tutorial program. +is exactly what we've shown in step-6.

The test case

diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index 1a4beaca1f..2883bc2841 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -192,26 +192,10 @@ namespace Step27 // @sect4{LaplaceProblem::setup_system} // - // This function is again an almost verbatim copy of what we already did in - // step-6. The first change is that we append the Dirichlet boundary - // conditions to the ConstraintMatrix object, which we consequently call - // just constraints instead of - // hanging_node_constraints. The second difference is that we - // don't directly build the sparsity pattern, but first create an - // intermediate object that we later copy into the usual SparsityPattern - // data structure, since this is more efficient for the problem with many - // entries per row (and different number of entries in different rows). In - // another slight deviation, we do not first build the sparsity pattern and - // then condense away constrained degrees of freedom, but pass the - // constraint matrix object directly to the function that builds the - // sparsity pattern. We disable the insertion of constrained entries with - // false as fourth argument in the DoFTools::make_sparsity_pattern - // function. All of these changes are explained in the introduction of this - // program. - // - // The last change, maybe hidden in plain sight, is that the dof_handler - // variable here is an hp object -- nevertheless all the function calls we - // had before still work in exactly the same way as they always did. + // This function is again a verbatim copy of what we already did in + // step-6. Despite function calls with exactly the same names and arguments, + // the algorithms used internally are different in some aspect since the + // dof_handler variable here is an hp object. template void LaplaceProblem::setup_system () { @@ -324,12 +308,6 @@ namespace Step27 local_dof_indices, system_matrix, system_rhs); } - - // Now with the loop over all cells finished, we are done for this - // function. The steps we still had to do at this point in earlier - // tutorial programs, namely condensing hanging node constraints and - // applying Dirichlet boundary conditions, have been taken care of by the - // ConstraintMatrix object constraints on the fly. } diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc index 5dcd8ecee8..07901f9bcb 100644 --- a/deal.II/examples/step-6/step-6.cc +++ b/deal.II/examples/step-6/step-6.cc @@ -449,11 +449,13 @@ void Step6::assemble_system () system_matrix, system_rhs); } - // Now we are done assembling the linear system. The constrained nodes are - // still in the linear system (there is a one on the diagonal of the matrix - // and all other entries for this line are set to zero) but the computed - // values are invalid. We compute the correct values for these nodes at the - // end of the solve function. + // Now we are done assembling the linear system. The constraint matrix took + // care of applying the boundary conditions and also eliminated hanging node + // constraints. The constrained nodes are still in the linear system (there + // is a one on the diagonal of the matrix and all other entries for this + // line are set to zero) but the computed values are invalid. We compute the + // correct values for these nodes at the end of the solve + // function. } -- 2.39.5