]> https://gitweb.dealii.org/ - dealii.git/commitdiff
step-6 now uses ConstraintMatrix::distribute_local_to_global()
authorTimo Heister <timo.heister@gmail.com>
Sat, 3 Nov 2012 02:51:05 +0000 (02:51 +0000)
committerTimo Heister <timo.heister@gmail.com>
Sat, 3 Nov 2012 02:51:05 +0000 (02:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@27325 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/examples/step-6/doc/intro.dox
deal.II/examples/step-6/step-6.cc

index 94f3c942772ec3bd61b266dcaf7d8e49c5265b71..9ff00210fb48471fea2341ce4af69b9fcbe7efbe 100644 (file)
@@ -108,6 +108,12 @@ never working correctly and it is not used.
 
 
 <ol>
+<li> step-6 now uses ConstraintMatrix::distribute_local_to_global()
+instead of condense(), which is the preferred way to use a ConstraintMatrix
+ (and the only sensible way in parallel).
+<br>
+(Timo Heister, 2012/11/02)
+
 <li> Fixed: DoFTools::make_flux_sparsity_pattern wasn't prepared to
 deal with adaptively refined meshes in 1d.
 <br>
index 1d20efef625563449d16ca1dce0386de3c27a7ce..0fc9752982c68c63b0fc54650d57a7fbbd09895a 100644 (file)
@@ -54,3 +54,18 @@ The only other new thing is a method to catch exceptions in the
 <code>main</code> function in order to output some information in case the
 program crashes for some reason.
 
+
+<h3>The ConstraintMatrix</h3>
+
+As explained above, we are going to use an object called ConstraintMatrix
+that will store the constraints at the hanging nodes to insure the solution
+is continuous at these nodes. We could first assemble the system as normal
+and then condense out the degrees of freedom that need to be constrained.
+This is also explained  @ref constraints documentation
+module. Instead we will go the more efficient route and eliminate the
+constrained entries while we are copying from the local to the global system.
+Because boundary conditions can be treated in the same way, we will
+incorporate the them as constraints in the same ConstraintMatrix object.
+This way, we don't need to apply the boundary conditions after assembly 
+(like we did in the earlier steps).
\ No newline at end of file
index caf09197cd07c12e7a57e54d4b2af369b018fc54..90c31101f1eea8b92e110fd30b26d3d9d52d50bf 100644 (file)
@@ -72,7 +72,8 @@
                                  // sure that the degrees of freedom
                                  // on hanging nodes conform to some
                                  // constraints such that the global
-                                 // solution is continuous. The
+                                 // solution is continuous. We are
+// also going to store the boundary conditions in this object. The
                                  // following file contains a class
                                  // which is used to handle these
                                  // constraints:
@@ -114,8 +115,7 @@ using namespace dealii;
                                  // (instead of the global refinement
                                  // in the previous examples), and a
                                  // variable which will hold the
-                                 // constraints associated to the
-                                 // hanging nodes. In addition, we
+                                 // constraints. In addition, we
                                  // have added a destructor to the
                                  // class for reasons that will become
                                  // clear when we discuss its
@@ -144,9 +144,10 @@ class Step6
                                      // This is the new variable in
                                      // the main class. We need an
                                      // object which holds a list of
-                                     // constraints originating from
-                                     // the hanging nodes:
-    ConstraintMatrix     hanging_node_constraints;
+                                     // constraints to hold the
+                                    // hanging nodes and the
+                                    // boundary conditions.    
+    ConstraintMatrix     constraints;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
@@ -421,7 +422,7 @@ void Step6<dim>::setup_system ()
                                    // hanging nodes. In the class
                                    // desclaration, we have already
                                    // allocated space for an object
-                                   // <code>hanging_node_constraints</code>
+                                   // <code>constraints</code>
                                    // that will hold a list of these
                                    // constraints (they form a matrix,
                                    // which is reflected in the name
@@ -435,23 +436,38 @@ void Step6<dim>::setup_system ()
                                    // over from computations on the
                                    // previous mesh before the last
                                    // adaptive refinement):
-  hanging_node_constraints.clear ();
+  constraints.clear ();
   DoFTools::make_hanging_node_constraints (dof_handler,
-                                           hanging_node_constraints);
+                                           constraints);
+
+
+  // Now we are ready to interpolate the ZeroFunction to our
+  // boundary with indicator 0 (the whole boundary) and store
+  // the resulting constraints in our <code>constraints</code>
+  // object. Note that we do not to apply the boundary
+  // conditions after assembly, like we did in earlier steps.
+  // As almost all the stuff,
+  // the interpolation of boundary
+  // values works also for higher
+  // order elements without the need
+  // to change your code for that. We
+  // note that for proper results, it
+  // is important that the
+  // elimination of boundary nodes
+  // from the system of equations
+  // happens *after* the elimination
+  // of hanging nodes. For that reason
+  // we are filling the boundary values into
+  // the ContraintMatrix after the hanging node
+  // constraints.
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            ZeroFunction<dim>(),
+                                            constraints);
+
 
                                    // The next step is <code>closing</code>
-                                   // this object. For this note that,
-                                   // in principle, the
-                                   // <code>ConstraintMatrix</code> class can
-                                   // hold other constraints as well,
-                                   // i.e. constraints that do not
-                                   // stem from hanging
-                                   // nodes. Sometimes, it is useful
-                                   // to use such constraints, in
-                                   // which case they may be added to
-                                   // the <code>ConstraintMatrix</code> object
-                                   // after the hanging node
-                                   // constraints were computed. After
+                                   // this object. After
                                    // all constraints have been added,
                                    // they need to be sorted and
                                    // rearranged to perform some
@@ -460,7 +476,7 @@ void Step6<dim>::setup_system ()
                                    // <code>close()</code> function, after which
                                    // no further constraints may be
                                    // added any more:
-  hanging_node_constraints.close ();
+  constraints.close ();
 
                                    // Now we first build our
                                    // compressed sparsity pattern like
@@ -468,26 +484,18 @@ void Step6<dim>::setup_system ()
                                    // examples. Nevertheless, we do
                                    // not copy it to the final
                                    // sparsity pattern immediately.
+  // Note that we call a variant of make_sparsity_pattern that takes
+  // the ConstraintMatrix as the third argument. We are letting
+  // the routine know, the we will never write into the locations
+  // given by <code>constraints</code> by setting the argument
+  // <code>keep_constrained_dofs</code> to false. If we were to
+  // condense the constraints after assembling, we would have
+  // to pass <code>true</code> instead.
   CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
-  DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
-
-                                   // The constrained hanging nodes
-                                   // will later be eliminated from
-                                   // the linear system of
-                                   // equations. When doing so, some
-                                   // additional entries in the global
-                                   // matrix will be set to non-zero
-                                   // values, so we have to reserve
-                                   // some space for them here. Since
-                                   // the process of elimination of
-                                   // these constrained nodes is
-                                   // called <code>condensation</code>, the
-                                   // functions that eliminate them
-                                   // are called <code>condense</code> for both
-                                   // the system matrix and right hand
-                                   // side, as well as for the
-                                   // sparsity pattern.
-  hanging_node_constraints.condense (c_sparsity);
+  DoFTools::make_sparsity_pattern(dof_handler,
+      c_sparsity,
+      constraints,
+      false /*keep_constrained_dofs*/);
 
                                    // Now all non-zero entries of the
                                    // matrix are known (i.e. those
@@ -509,9 +517,10 @@ void Step6<dim>::setup_system ()
                                  // @sect4{Step6::assemble_system}
 
                                  // Next, we have to assemble the
-                                 // matrix again. There are no code
-                                 // changes compared to step-5 except
-                                 // for a single place: We have to use
+                                 // matrix again. There are two code
+                                 // changes compared to step-5:
+//
+// First, we have to use
                                  // a higher-order quadrature formula
                                  // to account for the higher
                                  // polynomial degree in the finite
@@ -523,11 +532,15 @@ void Step6<dim>::setup_system ()
                                  // we had two points for bilinear
                                  // elements. Now we should use three
                                  // points for biquadratic elements.
-                                 //
+
+// Second, to copy the local matrix and vector on each cell
+// into the global system, we are no longer using a hand-written
+// loop. Instead, we use <code>ConstraintMatrix::distribute_local_to_global</code>
+// that internally executes this loop and eliminates all the constraints
+// at the same time.
+//
                                  // The rest of the code that forms
-                                 // the local contributions and
-                                 // transfers them into the global
-                                 // objects remains unchanged. It is
+                                 // the local contributions remains unchanged. It is
                                  // worth noting, however, that under
                                  // the hood several things are
                                  // different than before. First, the
@@ -596,67 +609,19 @@ void Step6<dim>::assemble_system ()
           }
 
       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);
-        }
+      // transfer the contributions from @p cell_matrix and @cell_rhs into the global objects.
+      constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
     }
-
-                                   // After the system of equations
-                                   // has been assembled just as for
-                                   // the previous examples, we still
-                                   // have to eliminate the
-                                   // constraints due to hanging
-                                   // nodes. This is done using the
-                                   // following two function calls:
-  hanging_node_constraints.condense (system_matrix);
-  hanging_node_constraints.condense (system_rhs);
-                                   // Using them, degrees of freedom
-                                   // associated to hanging nodes have
-                                   // been removed from the linear
-                                   // system and the independent
-                                   // variables are only the regular
-                                   // nodes. The constrained nodes are
+// 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 (the <code>condense</code>
-                                   // function modifies the system so
-                                   // that the values in the solution
-                                   // corresponding to constrained
-                                   // nodes are invalid, but that the
-                                   // system still has a well-defined
-                                   // solution; we compute the correct
+                                   // are invalid. We compute the correct
                                    // values for these nodes at the
-                                   // end of the <code>solve</code> function).
-
-                                   // As almost all the stuff before,
-                                   // the interpolation of boundary
-                                   // values works also for higher
-                                   // order elements without the need
-                                   // to change your code for that. We
-                                   // note that for proper results, it
-                                   // is important that the
-                                   // elimination of boundary nodes
-                                   // from the system of equations
-                                   // happens *after* the elimination
-                                   // of hanging nodes.
-  std::map<unsigned int,double> boundary_values;
-  VectorTools::interpolate_boundary_values (dof_handler,
-                                            0,
-                                            ZeroFunction<dim>(),
-                                            boundary_values);
-  MatrixTools::apply_boundary_values (boundary_values,
-                                      system_matrix,
-                                      solution,
-                                      system_rhs);
+                                   // end of the <code>solve</code> function.
 }
 
 
@@ -671,8 +636,8 @@ void Step6<dim>::assemble_system ()
                                  // have to incorporate hanging node
                                  // constraints. As mentioned above,
                                  // the degrees of freedom
-                                 // corresponding to hanging node
-                                 // constraints have been removed from
+                                 // from the ConstraintMatrix corresponding to hanging node
+                                 // constraints and boundary values have been removed from
                                  // the linear system by giving the
                                  // rows and columns of the matrix a
                                  // special treatment. This way, the
@@ -684,7 +649,7 @@ void Step6<dim>::assemble_system ()
                                  // constraints to assign to them the
                                  // values that they should have. This
                                  // process, called <code>distributing</code>
-                                 // hanging nodes, computes the values
+                                 // constraints, computes the values
                                  // of constrained nodes from the
                                  // values of the unconstrained ones,
                                  // and requires only a single
@@ -703,7 +668,7 @@ void Step6<dim>::solve ()
   solver.solve (system_matrix, solution, system_rhs,
                 preconditioner);
 
-  hanging_node_constraints.distribute (solution);
+  constraints.distribute (solution);
 }
 
 
@@ -1068,7 +1033,7 @@ void Step6<dim>::run ()
     }
 
                                    // After we have finished computing
-                                   // the solution on the finesh mesh,
+                                   // the solution on the finest mesh,
                                    // and writing all the grids to
                                    // disk, we want to also write the
                                    // actual solution on this final

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.