]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
step-6 uses the new way of dealing with constraints, i.e., applying them on the fly...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 3 May 2014 10:48:10 +0000 (10:48 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 3 May 2014 10:48:10 +0000 (10:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@32879 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-16/step-16.cc
deal.II/examples/step-27/doc/intro.dox
deal.II/examples/step-27/step-27.cc
deal.II/examples/step-6/step-6.cc

index 1eb7c93fac7bfb2b19acb717c7a3bb31fffbf31f..e8a2035f0bd1eb60d83c81fa2ec8571d71e07aa2 100644 (file)
 #include <deal.II/numerics/data_out.h>
 #include <deal.II/numerics/error_estimator.h>
 
-// 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 <deal.II/multigrid/mg_dof_handler.h>
 #include <deal.II/multigrid/mg_constrained_dofs.h>
 #include <deal.II/multigrid/multigrid.h>
 #include <deal.II/multigrid/mg_transfer.h>
@@ -112,7 +112,7 @@ namespace Step16
 
     Triangulation<dim>   triangulation;
     FE_Q<dim>            fe;
-    MGDoFHandler<dim>    mg_dof_handler;
+    DoFHandler<dim>      dof_handler;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
@@ -234,7 +234,7 @@ namespace Step16
     triangulation (Triangulation<dim>::
                    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 <int dim>
   void LaplaceProblem<dim>::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<triangulation.n_levels(); ++l)
       deallog << "   " << 'L' << l << ": "
-              << mg_dof_handler.n_dofs(l);
+              << dof_handler.n_dofs(l);
     deallog  << std::endl;
 
-    sparsity_pattern.reinit (mg_dof_handler.n_dofs(),
-                             mg_dof_handler.n_dofs(),
-                             mg_dof_handler.max_couplings_between_dofs());
-    DoFTools::make_sparsity_pattern (mg_dof_handler, sparsity_pattern);
+    sparsity_pattern.reinit (dof_handler.n_dofs(),
+                             dof_handler.n_dofs(),
+                             dof_handler.max_couplings_between_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
 
-    solution.reinit (mg_dof_handler.n_dofs());
-    system_rhs.reinit (mg_dof_handler.n_dofs());
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
 
     // But it starts to be a wee bit different here, although this still
     // doesn't have anything to do with multigrid methods. step-6 took care of
@@ -280,13 +282,13 @@ namespace Step16
     // for a later clean-up stage:
     constraints.clear ();
     hanging_node_constraints.clear ();
-    DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_constraints);
-    DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
+    DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints);
+    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
 
     typename FunctionMap<dim>::type      dirichlet_boundary_functions;
     ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
     dirichlet_boundary_functions[0] = &homogeneous_dirichlet_bc;
-    VectorTools::interpolate_boundary_values (static_cast<const DoFHandler<dim>&>(mg_dof_handler),
+    VectorTools::interpolate_boundary_values (static_cast<const DoFHandler<dim>&>(dof_handler),
                                               dirichlet_boundary_functions,
                                               constraints);
     constraints.close ();
@@ -299,7 +301,7 @@ namespace Step16
     // about the boundary values as well, so we pass the
     // <code>dirichlet_boundary</code> 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<n_levels; ++level)
       {
         CompressedSparsityPattern csp;
-        csp.reinit(mg_dof_handler.n_dofs(level),
-                   mg_dof_handler.n_dofs(level));
-        MGTools::make_sparsity_pattern(mg_dof_handler, csp, level);
+        csp.reinit(dof_handler.n_dofs(level),
+                   dof_handler.n_dofs(level));
+        MGTools::make_sparsity_pattern(dof_handler, csp, level);
 
         mg_sparsity_patterns[level].copy_from (csp);
 
@@ -352,11 +354,7 @@ namespace Step16
   // @sect4{LaplaceProblem::assemble_system}
 
   // The following function assembles the linear system on the finest level of
-  // the mesh. It is almost exactly the same as in step-6, with the exception
-  // that we don't eliminate hanging nodes and boundary values after
-  // assembling, but while copying local contributions into the global
-  // matrix. This is not only simpler but also more efficient for large
-  // problems.
+  // the mesh. It is almost exactly the same as in step-6.
   //
   // This latter trick is something that only found its way into deal.II over
   // time and wasn't used in the initial version of this tutorial
@@ -382,9 +380,9 @@ namespace Step16
     const Coefficient<dim> coefficient;
     std::vector<double>    coefficient_values (n_q_points);
 
-    typename MGDoFHandler<dim>::active_cell_iterator
-    cell = mg_dof_handler.begin_active(),
-    endc = mg_dof_handler.end();
+    typename DoFHandler<dim>::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 <code>MGConstraints</code> already computed the information
-    // for us when we called initialize in <code>setup_system()</code>.
-    std::vector<std::vector<bool> > interface_dofs
-      = mg_constrained_dofs.get_refinement_edge_indices ();
-    std::vector<std::vector<bool> > 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 <code>MGConstraints</code> knows.
-    // So we simply ask for them by calling <code>get_boundary_indices()</code>.
-    // 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 <code>MGConstraints</code>
+    // already computed the information for us when we called initialize in
+    // <code>setup_system()</code>. So we simply ask for them by calling
+    // <code>get_boundary_indices()</code> and
+    // <code>get_refinement_edge_indices(level)</code>. 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<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
     std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
     for (unsigned int level=0; level<triangulation.n_levels(); ++level)
       {
-        boundary_constraints[level].add_lines (interface_dofs[level]);
         boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices()[level]);
+        boundary_constraints[level].
+        add_lines (mg_constrained_dofs.get_refinement_edge_indices(level));
         boundary_constraints[level].close ();
 
         boundary_interface_constraints[level]
-        .add_lines (boundary_interface_dofs[level]);
+        .add_lines (mg_constrained_dofs.get_refinement_edge_boundary_indices(level));
         boundary_interface_constraints[level].close ();
       }
 
@@ -501,10 +494,10 @@ namespace Step16
     // <code>assemble_system</code>, 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<dim>::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<dim>::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; i<dofs_per_cell; ++i)
           for (unsigned int j=0; j<dofs_per_cell; ++j)
-            if ( !(interface_dofs[cell->level()][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 <int dim>
   void LaplaceProblem<dim>::solve ()
   {
-
-    // Create the object that deals with the transfer between different
-    // refinement levels. We need to pass it the hanging node constraints.
     MGTransferPrebuilt<Vector<double> > 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 <code>MGConstraints</code>
-    // knows about that so pass it as an argument.
-    mg_transfer.build_matrices(mg_dof_handler);
+    mg_transfer.build_matrices(dof_handler);
 
     FullMatrix<double> 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<Vector<double> > mg_matrix(mg_matrices);
+    mg::Matrix<Vector<double> > mg_interface_up(mg_interface_matrices);
+    mg::Matrix<Vector<double> > mg_interface_down(mg_interface_matrices);
 
     // Now, we are ready to set up the V-cycle operator and the multilevel
     // preconditioner.
-    Multigrid<Vector<double> > mg(mg_dof_handler,
+    Multigrid<Vector<double> > 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<dim, Vector<double>, MGTransferPrebuilt<Vector<double> > >
-    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 <int dim>
@@ -719,7 +708,7 @@ namespace Step16
   {
     Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
-    KellyErrorEstimator<dim>::estimate (static_cast<DoFHandler<dim>&>(mg_dof_handler),
+    KellyErrorEstimator<dim>::estimate (dof_handler,
                                         QGauss<dim-1>(3),
                                         typename FunctionMap<dim>::type(),
                                         solution,
@@ -737,7 +726,7 @@ namespace Step16
   {
     DataOut<dim> 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; level<triangulation.n_levels(); ++level)
-          std::cout << mg_dof_handler.n_dofs(level)
+          std::cout << dof_handler.n_dofs(level)
                     << (level == triangulation.n_levels()-1
                         ? ")" : ", ");
         std::cout << std::endl;
index e6a266c7980c9a43d580133dd6d21dda5628ef7d..a1a6e886965a67165ad4afa10d41fe0f924774dc 100644 (file)
@@ -681,101 +681,17 @@ freedom in non-$hp$ mode, but the difference is that each constrained
 hanging node is constrained not only against the two adjacent degrees
 of freedom, but is constrained against many more degrees of freedom.
 
-In previous programs, for example in step-6, we have dealt with eliminating
-these constrained degrees of freedom by first building an object of type
-ConstraintMatrix that stores the constraints, and then "condensing" away these
-degrees of freedom first from the sparsity pattern. The same scheme is used
-for the matrix and right hand side, which are also both first built then
-condensed.
-
-Dealing with sparsity patterns and matrices in this way turns out to be
-inefficient because the effort necessary is 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. The
-solution to this problem is to use the ConstraintMatrix object already at the
-time of creating the sparsity pattern, or while copying local contributions
-into the global matrix object (or global vector).
-
-So, instead of the code snippet (taken from step-6)
-@code
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                          hanging_node_constraints);
-  hanging_node_constraints.close ();
-
-  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
-  hanging_node_constraints.condense (sparsity_pattern);
-@endcode
-we now use the following:
-@code
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                          hanging_node_constraints);
-  hanging_node_constraints.close ();
-
-  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern,
-                                   hanging_node_constraints, false);
-@endcode
-
-Note the last <tt>false</tt> 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 &mdash; 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<dim>::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<dofs_per_cell; ++i)
-       {
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           system_matrix.add (local_dof_indices[i],
-                              local_dof_indices[j],
-                              cell_matrix(i,j));
-         system_rhs(local_dof_indices[i]) += cell_rhs(i);
-       }
-    }
-
-  hanging_node_constraints.condense (system_matrix);
-  hanging_node_constraints.condense (system_rhs);
-@endcode
-We now replace copying and later condensing into one step, as already shown in
-step-17, step-18, and step-22:
-@code
-  typename hp::DoFHandler<dim>::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
-<code>hanging_node_constraints.distribute_local_to_global</code> 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.
 
 
 <h3>The test case</h3>
index 1a4beaca1f84ff99ab0dd070fdbda588e8fc2816..2883bc28419afbd32d92ef17f5bbc6e03d9c4802 100644 (file)
@@ -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 <code>constraints</code> instead of
-  // <code>hanging_node_constraints</code>. 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
-  // <tt>false</tt> 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 <int dim>
   void LaplaceProblem<dim>::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 <code>constraints</code> on the fly.
   }
 
 
index 5dcd8ecee8806a7f3cf891068710260eb5689ee5..07901f9bcbaf077a8621cc07f94c3ac15f8032fe 100644 (file)
@@ -449,11 +449,13 @@ void Step6<dim>::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 <code>solve</code> 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 <code>solve</code>
+  // function.
 }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.