]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Apply new function in ConstraintMatrix to step-37.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 30 Sep 2009 15:01:58 +0000 (15:01 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 30 Sep 2009 15:01:58 +0000 (15:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@19619 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-37/doc/results.dox
deal.II/examples/step-37/step-37.cc

index a1ddd0c1c0465e376c74b455a9908d30bb7e4b57..a444ddf67a95b10545687d8a7319ccbdbb965c87 100644 (file)
@@ -9,38 +9,38 @@ output:
 @code
 Cycle 0
 Number of degrees of freedom: 337
-System matrix memory consumption: 0.0224 MBytes.
-Multigrid objects memory consumption: 0.04636 MBytes.
+System matrix memory consumption: 0.0257 MBytes.
+Multigrid objects memory consumption: 0.05071 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 1
 Number of degrees of freedom: 1313
-System matrix memory consumption: 0.08109 MBytes.
-Multigrid objects memory consumption: 0.1635 MBytes.
+System matrix memory consumption: 0.0925 MBytes.
+Multigrid objects memory consumption: 0.1792 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 2
 Number of degrees of freedom: 5185
-System matrix memory consumption: 0.3141 MBytes.
-Multigrid objects memory consumption: 0.6208 MBytes.
+System matrix memory consumption: 0.3551 MBytes.
+Multigrid objects memory consumption: 0.6775 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 3
 Number of degrees of freedom: 20609
-System matrix memory consumption: 1.243 MBytes.
-Multigrid objects memory consumption: 2.434 MBytes.
+System matrix memory consumption: 1.397 MBytes.
+Multigrid objects memory consumption: 2.644 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 4
 Number of degrees of freedom: 82177
-System matrix memory consumption: 4.951 MBytes.
-Multigrid objects memory consumption: 9.658 MBytes.
+System matrix memory consumption: 5.544 MBytes.
+Multigrid objects memory consumption: 10.46 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 5
 Number of degrees of freedom: 328193
-System matrix memory consumption: 19.77 MBytes.
-Multigrid objects memory consumption: 38.5 MBytes.
+System matrix memory consumption: 22.1 MBytes.
+Multigrid objects memory consumption: 41.63 MBytes.
 Convergence in 10 CG iterations.
 @endcode
 
@@ -51,32 +51,32 @@ program in three dimensions:
 @code
 Cycle 0
 Number of degrees of freedom: 517
-System matrix memory consumption: 0.09657 MBytes.
-Multigrid objects memory consumption: 0.1413 MBytes.
+System matrix memory consumption: 0.1 MBytes.
+Multigrid objects memory consumption: 0.1462 MBytes.
 Convergence in 9 CG iterations.
 
 Cycle 1
 Number of degrees of freedom: 3817
-System matrix memory consumption: 0.6334 MBytes.
-Multigrid objects memory consumption: 0.8567 MBytes.
+System matrix memory consumption: 0.6612 MBytes.
+Multigrid objects memory consumption: 0.8894 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 2
 Number of degrees of freedom: 29521
-System matrix memory consumption: 4.882 MBytes.
-Multigrid objects memory consumption: 6.403 MBytes.
+System matrix memory consumption: 5.099 MBytes.
+Multigrid objects memory consumption: 6.653 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 3
 Number of degrees of freedom: 232609
-System matrix memory consumption: 38.68 MBytes.
-Multigrid objects memory consumption: 50.28 MBytes.
+System matrix memory consumption: 40.4 MBytes.
+Multigrid objects memory consumption: 52.24 MBytes.
 Convergence in 11 CG iterations.
 
 Cycle 4
 Number of degrees of freedom: 1847617
-System matrix memory consumption: 308.4 MBytes.
-Multigrid objects memory consumption: 399.6 MBytes.
+System matrix memory consumption: 321.9 MBytes.
+Multigrid objects memory consumption: 415.1 MBytes.
 Convergence in 11 CG iterations.
 @endcode
 
index 038824429967f575d68a2457860711aa9c7f1d03..1d3a3edd8455a9ef8712ed7dbb17e4c5aea89271 100644 (file)
@@ -507,62 +507,69 @@ vmult_on_subrange (const unsigned int    first_cell,
 
                                 // OK, now we are sitting in the loop that
                                 // goes over our chunks of cells. What we
-                                // need to do is five things: First, we
-                                // have to give the full matrices
-                                // containing the solution at cell dofs and
-                                // quadrature points the correct sizes. We
-                                // use the <code>true</code> argument in
-                                // order to specify that this should be
-                                // done fast, i.e., the field will not be
-                                // initialized since we fill them manually
-                                // in a second anyway. Then, we copy the
+                                // need to do is five things: First, we have
+                                // to give the full matrices containing the
+                                // solution at cell dofs and quadrature
+                                // points the correct sizes. We use the
+                                // <code>true</code> argument in order to
+                                // specify that this should be done fast,
+                                // i.e., the field will not be initialized
+                                // since we fill them manually in the very
+                                // next step second anyway. Then, we copy the
                                 // source values from the global vector to
                                 // the local cell range, and we perform a
                                 // matrix-matrix product to tranform the
                                 // values to the quadrature points. It is a
                                 // bit tricky to find out how the matrices
-                                // should be multiplied with each
-                                // other. One way to resolve this is to
+                                // should be multiplied with each other,
+                                // i.e., which matrix needs to be
+                                // transposed. One way to resolve this is to
                                 // look at the matrix dimensions:
                                 // <code>solution_cells</code> has
                                 // <code>current_chunk_size</code> rows and
                                 // <code>matrix_sizes.m</code> columns,
                                 // whereas <code>small_matrix</code> has
                                 // <code>matrix_sizes.m</code> rows and
-                                // <code>matrix_sizes.n</code> columns,
-                                // which is also the size of columns in the
-                                // output matrix
+                                // <code>matrix_sizes.n</code> columns, which
+                                // is also the size of columns in the output
+                                // matrix
                                 // <code>solution_points</code>. Hence, the
-                                // columns of the first matrix are as many
-                                // as there are rows in the second, which
-                                // means that the product is done
-                                // non-transposed for both matrices.
+                                // columns of the first matrix are as many as
+                                // there are rows in the second, which means
+                                // that the product is done non-transposed
+                                // for both matrices.
                                 //
                                 // Once the first product is calculated, we
                                 // apply the derivative information on all
-                                // the cells and all the quadrature points
-                                // by calling the <code>transform</code>
+                                // the cells and all the quadrature points by
+                                // calling the <code>transform</code>
                                 // operation of the
                                 // <code>Transformation</code> class, and
-                                // then use a second matrix-matrix product
-                                // to get back to the solution values at
-                                // the support points. This time, we need
-                                // to transpose the small matrix, indicated
-                                // by a <code>mTmult</code> in the
-                                // operations. The fifth and last step is
-                                // to add the local data into the global
-                                // vector, which is what we did in many
-                                // tutorial programs when assembling right
-                                // hand sides. Just use the
-                                // <code>indices_local_to_global</code>
-                                // field to find out how local dofs and
-                                // global dofs are related to each other.
+                                // then use a second matrix-matrix product to
+                                // get back to the solution values at the
+                                // support points. This time, we need to
+                                // transpose the small matrix, indicated by a
+                                // <code>mTmult</code> in the operations. The
+                                // fifth and last step is to add the local
+                                // data into the global vector, which is what
+                                // we did in many tutorial programs when
+                                // assembling right hand sides. We use the
+                                // <code>indices_local_to_global</code> field
+                                // to find out how local dofs and global dofs
+                                // are related to each other. Since we
+                                // simultaneously apply the constraints, we
+                                // hand this task off to the ConstraintMatrix
+                                // object. Most often, itis used to work on
+                                // one cell at a time, but since we work on a
+                                // whole chunk of dofs, we can do that just
+                                // as easily for all the cells at once.
       solution_cells.reinit (current_chunk_size,matrix_sizes.m, true);
       solution_points.reinit (current_chunk_size,matrix_sizes.n, true);
 
-      for (unsigned int i=0; i<current_chunk_size; ++i)
-       for (unsigned int j=0; j<matrix_sizes.m; ++j)
-         solution_cells(i,j) = (number)src(indices_local_to_global(i+k,j));
+      const unsigned int n_cell_entries = current_chunk_size*matrix_sizes.m;
+      constraints.get_dof_values(src, &indices_local_to_global(k,0),
+                                &solution_cells(0,0),
+                                &solution_cells(0,0)+n_cell_entries);
 
       solution_cells.mmult (solution_points, small_matrix);
 
@@ -574,9 +581,10 @@ vmult_on_subrange (const unsigned int    first_cell,
 
       static Threads::Mutex mutex;
       Threads::Mutex::ScopedLock lock (mutex);
-      for (unsigned int i=0; i<current_chunk_size; ++i)
-       for (unsigned int j=0; j<matrix_sizes.m; ++j)
-         dst(indices_local_to_global(i+k,j)) += (number2)solution_cells(i,j);
+      constraints.distribute_local_to_global (&solution_cells(0,0),
+                                             &solution_cells(0,0)+n_cell_entries,
+                                             &indices_local_to_global(k,0),
+                                             dst);
     }
 }
 
@@ -627,25 +635,13 @@ MatrixFree<number,Transformation>::Tvmult_add (Vector<number2>       &dst,
 
 
 
-                                // The <code>vmult_add</code> function that
-                                // multiplies the matrix with vector
-                                // <code>src</code> and adds the result to
-                                // vector <code>dst</code> first creates a
-                                // copy of the source vector in order to
-                                // apply the constraints. The reason for
-                                // doing this is that constrained dofs are
-                                // zero when used in a solver like CG
-                                // (since they are not real degrees of
-                                // freedom), but the solution at the
-                                // respective nodes might still have
-                                // non-zero values which is necessary to
-                                // represent the field correctly in terms
-                                // of the FE basis functions. Then, we call
+                                // This is the <code>vmult_add</code>
+                                // function that multiplies the matrix with
+                                // vector <code>src</code> and adds the
+                                // result to vector <code>dst</code>. We call
                                 // a %parallel function that applies the
                                 // multiplication on a subrange of cells
-                                // (cf. the @ref threads module), and we
-                                // eventually condense the constraints on
-                                // the resulting vector.
+                                // (cf. the @ref threads module).
                                 //
                                 // TODO: Use WorkStream for parallelization
                                 // instead of apply_to_subranges, once we
@@ -657,33 +653,29 @@ void
 MatrixFree<number,Transformation>::vmult_add (Vector<number2>       &dst,
                                              const Vector<number2> &src) const
 {
-  Vector<number2> src_copy (src);
-  constraints.distribute(src_copy);
-
   parallel::apply_to_subranges (0, matrix_sizes.n_cells,
                                std_cxx1x::bind(&MatrixFree<number,Transformation>::
                                                template vmult_on_subrange<number2>,
                                                this,
                                                _1,_2,
                                                boost::ref(dst),
-                                               boost::cref(src_copy)),
+                                               boost::cref(src)),
                                200);
-  constraints.condense (dst);
 
                                 // One thing to be cautious about: The
-                                // deal.II classes expect that the
-                                // matrix still contains a diagonal
-                                // entry for constrained dofs
-                                // (otherwise, the matrix would be
-                                // singular, which is not what we
+                                // deal.II classes expect that the matrix
+                                // still contains a diagonal entry for
+                                // constrained dofs (otherwise, the matrix
+                                // would be singular, which is not what we
                                 // want). Since the
-                                // <code>condense</code> command of the
-                                // constraint matrix sets those
-                                // constrained elements to zero, we
-                                // have to circumvent that problem by
-                                // setting the diagonal to some
-                                // non-zero value. We simply set it to
-                                // one.
+                                // <code>distribute_local_to_global</code>
+                                // command of the constraint matrix which we
+                                // used for add the local elements into the
+                                // global matrix does not do write anything
+                                // with constrained elements, we have to
+                                // circumvent that problem by setting the
+                                // diagonal to some non-zero value. We simply
+                                // set it to one.
   for (unsigned int i=0; i<matrix_sizes.n_dofs; ++i)
     if (constraints.is_constrained(i) == true)
       dst(i) += 1.0 * src(i);

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.