From: kronbichler Date: Mon, 20 Apr 2009 12:43:51 +0000 (+0000) Subject: Use ConstraintMatrix for boundary conditions here since that is more efficient. Some... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=897b3da9a063674c5d6693cbf20111c29c65a2b3;p=dealii-svn.git Use ConstraintMatrix for boundary conditions here since that is more efficient. Some small updates in comments to be in line with the new step-22. git-svn-id: https://svn.dealii.org/trunk@18656 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-27/Makefile b/deal.II/examples/step-27/Makefile index 4d245e073c..bb6e80f5e9 100644 --- a/deal.II/examples/step-27/Makefile +++ b/deal.II/examples/step-27/Makefile @@ -30,7 +30,7 @@ D = ../../ # shall be deleted when calling `make clean'. Object and backup files, # executables and the like are removed anyway. Here, we give a list of # files in the various output formats that deal.II supports. -clean-up-files = *gmv *gnuplot *gpl *eps *pov +clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk diff --git a/deal.II/examples/step-27/doc/intro.dox b/deal.II/examples/step-27/doc/intro.dox index 54996363bb..90fea7f234 100644 --- a/deal.II/examples/step-27/doc/intro.dox +++ b/deal.II/examples/step-27/doc/intro.dox @@ -719,10 +719,10 @@ Note the last false 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 touch 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. +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, @@ -750,7 +750,7 @@ we have always used a process like this: hanging_node_constraints.condense (system_rhs); @endcode We now replace copying and later condensing into one step, as already shown in -@ref step_17 "step-17" and @ref step_18 "step-18": +@ref step_17 "step-17", @ref step_18 "step-18", and @ref step_22 "step-22": @code typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -762,19 +762,14 @@ We now replace copying and later condensing into one step, as already shown in cell->get_dof_indices (local_dof_indices); hanging_node_constraints - .distribute_local_to_global (cell_matrix, + .distribute_local_to_global (cell_matrix, cell_rhs, local_dof_indices, - system_matrix); - - hanging_node_constraints - .distribute_local_to_global (cell_rhs, - local_dof_indices, - system_rhs); + system_matrix, system_rhs); } @endcode Essentially, what the -hanging_node_constraints.distribute_local_to_global calls do is -to implement the same loop as before, but whenever we hit a constrained degree +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 @@ -782,6 +777,14 @@ 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. +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 +transfer from local to global data on matrix and vector simultaneously. This +is exactly what we've shown in the @ref step_22 "step-22" tutorial program. +

The test case

diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index a5ac374a72..48c044a942 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -118,7 +118,7 @@ class LaplaceProblem hp::QCollection quadrature_collection; hp::QCollection face_quadrature_collection; - ConstraintMatrix hanging_node_constraints; + ConstraintMatrix constraints; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; @@ -212,25 +212,35 @@ LaplaceProblem::~LaplaceProblem () // // This function is again an almost // verbatim copy of what we already did in - // step-6, with the main difference that we - // don't directly build the sparsity - // pattern, but first create an - // intermediate object that we later copy - // into the right data structure. 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 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. Both of these changes are + // 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 second change, maybe hidden in plain + // 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 @@ -244,15 +254,18 @@ void LaplaceProblem::setup_system () solution.reinit (dof_handler.n_dofs()); system_rhs.reinit (dof_handler.n_dofs()); - hanging_node_constraints.clear (); + constraints.clear (); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + constraints); DoFTools::make_hanging_node_constraints (dof_handler, - hanging_node_constraints); - hanging_node_constraints.close (); + constraints); + constraints.close (); CompressedSetSparsityPattern csp (dof_handler.n_dofs(), dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, csp, - hanging_node_constraints, false); + DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false); sparsity_pattern.copy_from (csp); system_matrix.reinit (sparsity_pattern); @@ -319,7 +332,7 @@ void LaplaceProblem::assemble_system () Vector cell_rhs; std::vector local_dof_indices; - + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -359,31 +372,21 @@ void LaplaceProblem::assemble_system () local_dof_indices.resize (dofs_per_cell); cell->get_dof_indices (local_dof_indices); - hanging_node_constraints - .distribute_local_to_global (cell_matrix, - local_dof_indices, - system_matrix); - - hanging_node_constraints - .distribute_local_to_global (cell_rhs, - local_dof_indices, - system_rhs); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); } - // After the steps already discussed above, - // all we still have to do is to treat - // Dirichlet boundary values - // correctly. This works in exactly the - // same way as for non-hp objects: - std::map boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(), - boundary_values); - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - 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. } @@ -409,7 +412,7 @@ void LaplaceProblem::solve () cg.solve (system_matrix, solution, system_rhs, preconditioner); - hanging_node_constraints.distribute (solution); + constraints.distribute (solution); } @@ -515,9 +518,9 @@ void LaplaceProblem::postprocess (const unsigned int cycle) // VTK format): const std::string filename = "solution-" + Utilities::int_to_string (cycle, 2) + - ".gmv"; + ".vtk"; std::ofstream output (filename.c_str()); - data_out.write_gmv (output); + data_out.write_vtk (output); } // After this, we would like to actually @@ -739,7 +742,7 @@ void LaplaceProblem::run () << dof_handler.n_dofs() << std::endl << " Number of constraints : " - << hanging_node_constraints.n_constraints() + << constraints.n_constraints() << std::endl; assemble_system ();