]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
When eliminating boundary values from a matrix, we only have to eliminate the respect...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 17 May 2000 12:22:49 +0000 (12:22 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 17 May 2000 12:22:49 +0000 (12:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@2874 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_tools.h
deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/numerics/base.cc
deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc

index b5165e8400bdb89790994741789c0d5bd0cf7af0..af6dc71a7a7487981d5a2d3d8a59b9d4e3009db8 100644 (file)
@@ -83,7 +83,12 @@ class DoFTools
                                      * #make_sparsity_pattern# just
                                      * loops over all cells and
                                      * enters all couplings local to
-                                     * that cell.
+                                     * that cell. As the generation
+                                     * of the sparsity pattern is
+                                     * irrespective of the equation
+                                     * which is solved later on, the
+                                     * resulting sparsity pattern is
+                                     * symmetric.
                                      *
                                      * Since this process is purely
                                      * local, the sparsity pattern
@@ -157,13 +162,18 @@ class DoFTools
                                      * than to the degrees of freedom
                                      * contained in there.
                                      *
-                                     * This function is designed to accept
-                                     * a mask, like the one shown above,
-                                     * through the #mask# parameter, which
-                                     * contains boolean values. It builds
-                                     * the matrix structure just like the
-                                     * previous function, but does not create
-                                     * elements if not specified by the mask.
+                                     * This function is designed to
+                                     * accept a mask, like the one
+                                     * shown above, through the
+                                     * #mask# parameter, which
+                                     * contains boolean values. It
+                                     * builds the matrix structure
+                                     * just like the previous
+                                     * function, but does not create
+                                     * elements if not specified by
+                                     * the mask. If the mask is
+                                     * symmetric, then so will be the
+                                     * resulting sparsity pattern.
                                      *
                                      * The actual type of the
                                      * sparsity pattern may be
index 37c761551aa86562040e72f08e9e3d698151c69e..f56fae4baaab578fea83b3f78027332ebe5673da 100644 (file)
@@ -352,9 +352,9 @@ class MatrixCreator
  *
  * \subsection{Boundary conditions}
  *
- * The #apply_boundar_values# function inserts boundary conditions of
+ * The #apply_boundary_values# function inserts boundary conditions of
  * into a system of equations.  To actually do this you have to specify
- * a list of degree of freedom indices along with the value this degree of
+ * a list of degree of freedom indices along with the values these degrees of
  * freedom shall assume. To see how to get such a list, see the discussion
  * of the #VectorTools::interpolate_boundary_values# function.
  *
@@ -371,6 +371,43 @@ class MatrixCreator
  * eliminating all coupling between this degree of freedom and others. Now
  * also the column consists only of zeroes, apart from the main diagonal entry.
  *
+ * Finding which rows contain an entry in the column for which we are
+ * presently performing a Gauss elimination step is either difficult
+ * or very simple, depending on the circumstances. If the sparsity
+ * pattern is symmetric (whether the matrix is symmetric is irrelevant
+ * here), then we can infer the rows which have a nonzero entry in the
+ * present column by looking at which columns in the present row are
+ * nonempty. In this case, we only need to look into a fixed number of
+ * rows and need not search all rows. On the other hand, if the
+ * sparsity pattern is nonsymmetric, then we need to use an iterative
+ * solver which can handle nonsymmetric matrices in any case, so there
+ * may be no need to do the Gauss elimination anyway. In fact, this is
+ * the way the function works: it takes a parameter
+ * (#elininate_columns#) that specifies whether the sparsity pattern
+ * is symmetric; if so, then the column is eliminated and the right
+ * hand side is also modified accordingly. If not, then only the row
+ * is deleted and the column is not touched at all, and all right hand
+ * side values apart from the one corresponding to the present row
+ * remain unchanged.
+ *
+ * If the sparsity pattern for your matrix is non-symmetric, you must
+ * set the value of this parameter to #false# in any case, since then
+ * we can't eliminate the column without searching all rows, which
+ * would be too expensive (if #N# be the number of rows, and #m# the
+ * number of nonzero elements per row, then eliminating one column is
+ * an #O(N*log(m))# operation, since searching in each row takes
+ * #log(m)# operations). If your sparsity pattern is symmetric, but
+ * your matrix is not, then you might specify #false# as well. If your
+ * sparsity pattern and matrix are both symmetric, you might want to
+ * specify #true# (the complexity of eliminating one row is then
+ * #O(m*log(m))#, since we only have to search #m# rows for the
+ * respective element of the column). Given the fact that #m# is
+ * roughly constant, irrespective of the discretization, and that the
+ * number of boundary nodes is #sqrt(N)# in 2d, the algorithm for
+ * symmetric sparsity patterns is #O(sqrt(N)*m*log(m))#, while it
+ * would be #O(N*sqrt(N)*m*log(m))# for the general case; the latter
+ * is too expensive to be performed.
+ *
  * It seems as if we had to make clear not to overwrite the lines of other
  * boundary nodes when doing the Gauss elimination step. However, since we
  * reset the right hand side when passing such a node, it is not a problem
@@ -381,33 +418,22 @@ class MatrixCreator
  * the right hand side. We need therefore not take special care of other
  * boundary nodes.
  * 
- * To make solving faster, we preset the solution vector with the right boundary
- * values. Since boundary nodes can never be hanging nodes, and since all other
- * entries of the solution vector are zero, we need not condense the solution
- * vector if the condensation process is done in-place. If done by copying
- * matrix and vectors to smaller ones, it would also be necessary to condense
- * the solution vector to preserve the preset boundary values.
- * 
- * It it not clear whether the deletion of coupling between the boundary degree
- * of freedom and other dofs really forces the corresponding entry in the
- * solution vector to have the right value when using iterative solvers,
- * since their search directions may contain components in the direction
- * of the boundary node. For this reason, we perform a very simple line
- * balancing by not setting the main diagonal entry to unity, but rather
- * to the value it had before deleting this line, or to the first nonzero
- * main diagonal entry if it is zero from a previous Gauss elimination
- * step. Of course we have to change
- * the right hand side appropriately. This is not a very good
- * strategy, but it at least should give the main diagonal entry a value
- * in the right order of dimension, which makes the solving process a bit
- * more stable. A refined algorithm would set the entry to the mean of the
- * other diagonal entries, but this seems to be too expensive.
- *
- * Because of the mentioned question, whether or not a preset solution value
- * which does not couple with other degrees of freedom remains its value or
- * not during solving iteratively, it may or may not be necessary to set
- * the correct value after solving again. This question is an open one as of
- * now and may be answered by future experience.
+ * To make solving faster, we preset the solution vector with the
+ * right boundary values. It it not clear whether the deletion of
+ * coupling between the boundary degree of freedom and other dofs
+ * really forces the corresponding entry in the solution vector to
+ * have the right value when using iterative solvers, since their
+ * search directions may contain components in the direction of the
+ * boundary node. For this reason, we perform a very simple line
+ * balancing by not setting the main diagonal entry to unity, but
+ * rather to the value it had before deleting this line, or to the
+ * first nonzero main diagonal entry if it is zero for some reason.
+ * Of course we have to change the right hand side appropriately. This
+ * is not a very good strategy, but it at least should give the main
+ * diagonal entry a value in the right order of dimension, which makes
+ * the solvution process a bit more stable. A refined algorithm would
+ * set the entry to the mean of the other diagonal entries, but this
+ * seems to be too expensive.
  *
  * 
  * @author Wolfgang Bangerth, 1998
@@ -425,7 +451,8 @@ class MatrixTools : public MatrixCreator<dim>
     static void apply_boundary_values (const map<unsigned int,double> &boundary_values,
                                       SparseMatrix<double>  &matrix,
                                       Vector<double>        &solution,
-                                      Vector<double>        &right_hand_side);
+                                      Vector<double>        &right_hand_side,
+                                      const bool             eliminate_columns = true);
 
                                     /**
                                      * Exception
index 7a07c3711cc6e41de81994833f92423ad30e6c53..ddf586f3ae171a41fa3fe1add10224a31eb5f21b 100644 (file)
@@ -140,7 +140,8 @@ void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
                                              boundary_value_list);
   MatrixTools<dim>::apply_boundary_values (boundary_value_list,
                                           system_matrix, solution,
-                                          right_hand_side);  
+                                          right_hand_side,
+                                          true);
 };
 
 
index 8f838033f69b3bf1371194c150561e008d82faf4..4d26689880844606e8d4fba5a9f1b9fe6cc1de41 100644 (file)
@@ -34,7 +34,8 @@ template <int dim>
 void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
                                             const Quadrature<dim>    &q,
                                             SparseMatrix<double>     &matrix,
-                                            const Function<dim> * const a) {
+                                            const Function<dim> * const a)
+{
   Vector<double> dummy;    // no entries, should give an error if accessed
   UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values);
   if (a != 0)
@@ -57,13 +58,15 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 };
 
 
+
 template <int dim>
 void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
                                             const Quadrature<dim>    &q,
                                             SparseMatrix<double>     &matrix,
                                             const Function<dim>      &rhs,
                                             Vector<double>           &rhs_vector,
-                                            const Function<dim> * const a) {
+                                            const Function<dim> * const a)
+{
   UpdateFlags update_flags = UpdateFlags(update_values |
                                         update_q_points |
                                         update_JxW_values);
@@ -85,9 +88,11 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 };
 
 
+
 template <int dim>
 void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
-                                            SparseMatrix<double>     &matrix) {
+                                            SparseMatrix<double>     &matrix)
+{
   const FiniteElement<dim> &fe = dof.get_fe();
 
   const unsigned int dofs_per_cell = fe.dofs_per_cell;
@@ -119,7 +124,8 @@ void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1>    &,
                                                    const FunctionMap      &,
                                                    Vector<double>         &,
                                                    vector<unsigned int>   &,
-                                                   const Function<1>      *) {
+                                                   const Function<1>      *)
+{
   Assert (false, ExcNotImplemented());
 };
 
@@ -133,7 +139,8 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
                                                      const FunctionMap        &rhs,
                                                      Vector<double>           &rhs_vector,
                                                      vector<unsigned int>     &dof_to_boundary_mapping,
-                                                     const Function<dim>      *a) {
+                                                     const Function<dim>      *a)
+{
   const FiniteElement<dim> &fe = dof.get_fe();
   const unsigned int n_components  = fe.n_components();
   const bool         fe_is_system  = (n_components != 1);
@@ -404,11 +411,13 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
 };
 
 
+
 template <int dim>
 void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
                                                const Quadrature<dim>    &q,
                                                SparseMatrix<double>     &matrix,
-                                               const Function<dim> * const a) {
+                                               const Function<dim> * const a)
+{
   const unsigned int n_components  = dof.get_fe().n_components();
   Assert ((n_components==1) || (a==0), ExcNotImplemented());
 
@@ -434,6 +443,8 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
   while ((++assembler).state() == valid);
 };
 
+
+
 /*
 
 template <int dim>
@@ -468,6 +479,8 @@ void MatrixCreator<dim>::create_level_laplace_matrix (unsigned int level,
 
 */
 
+
+
 template <int dim>
 void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
                                                const Quadrature<dim>    &q,
@@ -500,11 +513,15 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 };
 
 
+
 template <int dim>
-void MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &boundary_values,
-                                             SparseMatrix<double>  &matrix,
-                                             Vector<double>   &solution,
-                                             Vector<double>   &right_hand_side) {
+void
+MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &boundary_values,
+                                        SparseMatrix<double>  &matrix,
+                                        Vector<double>   &solution,
+                                        Vector<double>   &right_hand_side,
+                                        const bool        preserve_symmetry)
+{
   Assert (matrix.n() == matrix.m(),
          ExcDimensionsDontMatch(matrix.n(), matrix.m()));
   Assert (matrix.n() == right_hand_side.size(),
@@ -583,47 +600,63 @@ void MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &bo
          new_rhs = right_hand_side(dof_number)
                  = dof->second * first_nonzero_diagonal_entry;
        };
-      
-                                      // store the only nonzero entry
-                                      // of this line for the Gauss
-                                      // elimination step
-      const double diagonal_entry = matrix.diag_element(dof_number);
 
-                                      // do the Gauss step
-      for (unsigned int row=0; row<n_dofs; ++row) 
+
+                                      // if the user wants to have
+                                      // the symmetry of the matrix
+                                      // preserved, and if the
+                                      // sparsity pattern is
+                                      // symmetric, then do a Gauss
+                                      // elimination step with the
+                                      // present row
+      if (preserve_symmetry)
        {
-                                          // we need not handle the
-                                          // row we have already cleared
-         if (row == dof_number)
-           continue;
-
-                                          // check whether the line has
-                                          // an entry in the row corresponding
-                                          // to the dof presently worked with.
-                                          // note again: the first entry is
-                                          // the diagonal entry which we
-                                          // cannot be interested in; following
-                                          // are the other entries in sorted
-                                          // order, so we can use a binary
-                                          // search
-                                          //
-                                          // if this row contains an element
-                                          // for this dof, *p==dof_number
-         const unsigned int * p = lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
-                                              &sparsity_colnums[sparsity_rowstart[row+1]],
-                                              dof_number);
-                                          // check whether this line has
-                                          // an entry in the regarding column
-                                          // (check for ==dof_number and
-                                          // != next_row, since if
-                                          // row==dof_number-1, *p is a
-                                          // past-the-end pointer but points
-                                          // to dof_number anyway...)
-         if ((*p == dof_number) &&
-             (p != &sparsity_colnums[sparsity_rowstart[row+1]]))
-                                            // this line has an entry
-                                            // in the regarding column
+                                          // store the only nonzero entry
+                                          // of this line for the Gauss
+                                          // elimination step
+         const double diagonal_entry = matrix.diag_element(dof_number);
+         
+                                          // we have to loop over all
+                                          // rows of the matrix which
+                                          // have a nonzero entry in
+                                          // the column which we work
+                                          // in presently. if the
+                                          // sparsity pattern is
+                                          // symmetric, then we can
+                                          // get the positions of
+                                          // these rows cheaply by
+                                          // looking at the nonzero
+                                          // column numbers of the
+                                          // present row. we need not
+                                          // look at the first entry,
+                                          // since that is the
+                                          // diagonal element and
+                                          // thus the present row
+         for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
            {
+             const unsigned int row = sparsity_colnums[j];
+
+                                              // find the position of
+                                              // element
+                                              // (row,dof_number)
+             const unsigned int *
+               p = lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
+                               &sparsity_colnums[sparsity_rowstart[row+1]],
+                               dof_number);
+
+                                              // check whether this line has
+                                              // an entry in the regarding column
+                                              // (check for ==dof_number and
+                                              // != next_row, since if
+                                              // row==dof_number-1, *p is a
+                                              // past-the-end pointer but points
+                                              // to dof_number anyway...)
+                                              //
+                                              // there should be such an entry!
+             Assert ((*p == dof_number) &&
+                     (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
+                     ExcInternalError());
+
              const unsigned int global_entry
                = (p - &sparsity_colnums[sparsity_rowstart[0]]);
              
@@ -636,7 +669,6 @@ void MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &bo
            };
        };
 
-
                                       // preset solution vector
       solution(dof_number) = dof->second;
     };
@@ -648,13 +680,16 @@ MassMatrix<dim>::MassMatrix (const Function<dim> * const rhs,
                             const Function<dim> * const a) :
                Equation<dim> (1),
                right_hand_side (rhs),
-               coefficient (a)   {};
+               coefficient (a)
+{};
+
 
 
 template <int dim>
 void MassMatrix<dim>::assemble (FullMatrix<double>      &cell_matrix,
                                const FEValues<dim>     &fe_values,
-                               const typename DoFHandler<dim>::cell_iterator &) const {
+                               const typename DoFHandler<dim>::cell_iterator &) const
+{
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                     n_q_points    = fe_values.n_quadrature_points;
   const FiniteElement<dim>    &fe  = fe_values.get_fe();
@@ -735,11 +770,13 @@ void MassMatrix<dim>::assemble (FullMatrix<double>      &cell_matrix,
 };
 
 
+
 template <int dim>
 void MassMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
                                Vector<double>      &rhs,
                                const FEValues<dim> &fe_values,
-                               const DoFHandler<dim>::cell_iterator &) const {
+                               const DoFHandler<dim>::cell_iterator &) const
+{
   Assert (right_hand_side != 0, ExcNoRHSSelected());
 
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
@@ -800,10 +837,12 @@ void MassMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
 };
 
 
+
 template <int dim>
 void MassMatrix<dim>::assemble (Vector<double>      &rhs,
                                const FEValues<dim> &fe_values,
-                               const DoFHandler<dim>::cell_iterator &) const {
+                               const DoFHandler<dim>::cell_iterator &) const
+{
   Assert (right_hand_side != 0, ExcNoRHSSelected());
 
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
@@ -833,6 +872,7 @@ void MassMatrix<dim>::assemble (Vector<double>      &rhs,
 };
 
 
+
 template <int dim>
 LaplaceMatrix<dim>::LaplaceMatrix (const Function<dim> * const rhs,
                                   const Function<dim> * const a) :
@@ -841,11 +881,13 @@ LaplaceMatrix<dim>::LaplaceMatrix (const Function<dim> * const rhs,
                coefficient (a) {};
 
 
+
 template <int dim>
 void LaplaceMatrix<dim>::assemble (FullMatrix<double>         &cell_matrix,
                                   Vector<double>             &rhs,
                                   const FEValues<dim>        &fe_values,
-                                  const DoFHandler<dim>::cell_iterator &) const {
+                                  const DoFHandler<dim>::cell_iterator &) const
+{
   Assert (right_hand_side != 0, ExcNoRHSSelected());
   
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
@@ -909,10 +951,12 @@ void LaplaceMatrix<dim>::assemble (FullMatrix<double>         &cell_matrix,
 };
 
 
+
 template <int dim>
 void LaplaceMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
                                   const FEValues<dim> &fe_values,
-                                  const DoFHandler<dim>::cell_iterator &) const {
+                                  const DoFHandler<dim>::cell_iterator &) const
+{
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                     n_q_points    = fe_values.n_quadrature_points;
 
@@ -963,10 +1007,12 @@ void LaplaceMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
 };
 
 
+
 template <int dim>
 void LaplaceMatrix<dim>::assemble (Vector<double>      &rhs,
                                   const FEValues<dim> &fe_values,
-                                  const DoFHandler<dim>::cell_iterator &) const {
+                                  const DoFHandler<dim>::cell_iterator &) const
+{
   Assert (right_hand_side != 0, ExcNoRHSSelected());
 
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
@@ -997,6 +1043,7 @@ void LaplaceMatrix<dim>::assemble (Vector<double>      &rhs,
 };
 
 
+
 template<int dim>
 void
 MatrixCreator<dim>::create_interpolation_matrix(const FiniteElement<dim> &high,
@@ -1032,6 +1079,9 @@ MatrixCreator<dim>::create_interpolation_matrix(const FiniteElement<dim> &high,
 }
 
 
+
+// explicit instantiations
+
 template class MatrixCreator<deal_II_dimension>;
 template class MatrixTools<deal_II_dimension>;
 template class MassMatrix<deal_II_dimension>;
index c2587cb9791593e2ade2bd17156460dc156e73d2..8c3fb2a2043033cc342c255004739674fec95086 100644 (file)
@@ -437,7 +437,8 @@ void VectorTools::project (const DoFHandler<dim>    &dof,
   constraints.condense (tmp);
   if (boundary_values.size() != 0)
     MatrixTools<dim>::apply_boundary_values (boundary_values,
-                                            mass_matrix, vec, tmp);
+                                            mass_matrix, vec, tmp,
+                                            true);
 
   SolverControl           control(1000,1e-16);
   PrimitiveVectorMemory<> memory;

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.