]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Use ConstraintMatrix for boundary conditions here since that is more efficient. Some...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 20 Apr 2009 12:43:51 +0000 (12:43 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 20 Apr 2009 12:43:51 +0000 (12:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@18656 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 4d245e073c89c4d885239ba220e537989462f3ea..bb6e80f5e98f83c795a09b1254da0bf3e7a72437 100644 (file)
@@ -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
 
 
 
index 54996363bb122e7111225cf3c98539f409b90eb7..90fea7f234ee60c5594eca5898d94c3e6b457b13 100644 (file)
@@ -719,10 +719,10 @@ 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 touch 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.
+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,
@@ -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<dim>::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
-<code>hanging_node_constraints.distribute_local_to_global</code> calls do is
-to implement the same loop as before, but whenever we hit a constrained degree
+<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
@@ -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.
+
 
 <h3>The test case</h3>
 
index a5ac374a72d1cc503238ca22b831464affd34d75..48c044a94235615fa602239cd3964514f1ba8489 100644 (file)
@@ -118,7 +118,7 @@ class LaplaceProblem
     hp::QCollection<dim>     quadrature_collection;
     hp::QCollection<dim-1>   face_quadrature_collection;
 
-    ConstraintMatrix     hanging_node_constraints;
+    ConstraintMatrix     constraints;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
@@ -212,25 +212,35 @@ LaplaceProblem<dim>::~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
+                                // <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 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. 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 <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 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<dim>::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<dim>(),
+                                           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<dim>::assemble_system ()
   Vector<double>       cell_rhs;
 
   std::vector<unsigned int> local_dof_indices;
-  
+
   typename hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
@@ -359,31 +372,21 @@ void LaplaceProblem<dim>::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-<i>hp</i> objects:
-  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);
+                                  // 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.
 }
 
 
@@ -409,7 +412,7 @@ void LaplaceProblem<dim>::solve ()
   cg.solve (system_matrix, solution, system_rhs,
            preconditioner);
 
-  hanging_node_constraints.distribute (solution);
+  constraints.distribute (solution);
 }
 
 
@@ -515,9 +518,9 @@ void LaplaceProblem<dim>::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<dim>::run ()
                << dof_handler.n_dofs()
                << std::endl
                << "   Number of constraints       : "
-               << hanging_node_constraints.n_constraints()
+               << constraints.n_constraints()
                << std::endl;
 
       assemble_system ();

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.