]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Changed DoFTools::extract_constant_modes to work on bool vectors instead of doubles...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 2 Sep 2008 10:33:02 +0000 (10:33 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 2 Sep 2008 10:33:02 +0000 (10:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@16705 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_tools.h
deal.II/deal.II/source/dofs/dof_tools.cc
deal.II/examples/step-31/step-31.cc
deal.II/lac/include/lac/trilinos_precondition_amg.h
deal.II/lac/source/trilinos_precondition_amg.cc
deal.II/lac/source/trilinos_sparse_matrix.cc

index 1a1982eccdde2313e5fa5b6e8c358eb2613ee078..2634881624ed8cd475f95a2c8ae91cca6804a9d5 100644 (file)
@@ -1040,26 +1040,24 @@ class DoFTools
                                      * of as many vectors as there
                                      * are true arguments in 
                                      * <tt>component_select</tt>, 
-                                     * say <tt>d</tt>,
                                      * each of which will be one
                                      * in one component and 
-                                     * zero in all others. We store
-                                     * the null space contiguously
-                                     * in the output vector
-                                     * <tt>constant_modes</tt>
-                                     * to enable the internal 
-                                     * Trilinos structures the 
-                                     * access to it. This means that
-                                     * the length of the vector 
-                                     * will be <tt>d</tt> times
-                                     * the number of degrees of
+                                     * zero in all others. We store 
+                                     * this object in a vector of 
+                                     * vectors, where the outer 
+                                     * vector is of the size of 
+                                     * the number of selected 
+                                     * components, and each inner
+                                     * vector has as many components
+                                     * as there are degrees of 
                                      * freedom in the selected 
                                      * components. Note that any
                                      * matrix associated with this
                                      * null space must have been
-                                     * constructed from the 
-                                     * same argument as in
-                                     * <tt>component_select</tt>.
+                                     * constructed using the 
+                                     * same
+                                     * <tt>component_select</tt>
+                                     * argument.
                                      * 
                                      * The main reason for this
                                      * program is the use of the
@@ -1068,9 +1066,9 @@ class DoFTools
                                      */
     template <class DH>
     static void
-    extract_constant_modes (const DH                &dof_handler,
-                           const std::vector<bool> &component_select,
-                           std::vector<double>     &constant_modes);
+    extract_constant_modes (const DH                        &dof_handler,
+                           const std::vector<bool>         &component_select,
+                           std::vector<std::vector<bool> > &constant_modes);
 
                                     /**
                                      * For each active cell of a DoFHandler
index 19372180ac4ff75be635aa2d7845a686da4d0b52..998f451bedf91f1058ee267a7a4bdcb6ab4536dc 100644 (file)
@@ -3615,9 +3615,9 @@ DoFTools::extract_subdomain_dofs (const DH           &dof_handler,
 
 template <class DH>
 void
-DoFTools::extract_constant_modes (const DH                &dof_handler,
-                                 const std::vector<bool> &component_select,
-                                 std::vector<double>     &constant_modes)
+DoFTools::extract_constant_modes (const DH                        &dof_handler,
+                                 const std::vector<bool>         &component_select,
+                                 std::vector<std::vector<bool> > &constant_modes)
 {
   const unsigned int n_components = dof_handler.get_fe().n_components();
   Assert (n_components == component_select.size(),
@@ -3642,7 +3642,7 @@ DoFTools::extract_constant_modes (const DH                &dof_handler,
   std::vector<bool> temporary_dof_list (dof_handler.n_dofs(), false);
   extract_dofs (dof_handler, component_select, selection_dof_list);
 
-  constant_modes.resize (n_components_selected * n_u, 0.);
+  constant_modes.resize (n_components_selected, std::vector<bool>(n_u, false));
   
   for (unsigned int component=0, component_used=0; 
        component < n_components; ++component, ++component_used)
@@ -3655,12 +3655,15 @@ DoFTools::extract_constant_modes (const DH                &dof_handler,
        unsigned int counter = 0;
        for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
        {
-         if (temporary_dof_list[i])
+         if (selection_dof_list[i])
            {
-             constant_modes [component * n_u + counter] = 1.;
+             if (temporary_dof_list[i])
+               constant_modes [component][counter] = true;
+             else
+               constant_modes [component][counter] = false;
+
+             ++counter;
            }
-         if (selection_dof_list[i])
-           ++counter;
        }
       }
 }
@@ -5608,7 +5611,7 @@ void
 DoFTools::extract_constant_modes<DoFHandler<deal_II_dimension> >
 (const DoFHandler<deal_II_dimension> &dof_handler,
  const std::vector<bool> &selected_components,
- std::vector<double>     &constant_modes);
+ std::vector<std::vector<bool> > &constant_modes);
 
 template
 void
index 170fe0a8d0e523e8d396d5e8ccc04ec6e545891d..cfd21ca5e15299cacd1cee6ee3a14b021618d027 100644 (file)
@@ -1189,13 +1189,13 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
   Amg_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionAMG>
                       (new TrilinosWrappers::PreconditionAMG());
 
-  std::vector<double> null_space;
+  std::vector<std::vector<bool> > null_space;
   std::vector<bool>  velocity_components (dim+1,true);
   velocity_components[dim] = false;
   DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components, 
                                    null_space);
   Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
-                                true, true, null_space, dim, false);
+                                true, true, null_space, false);
 
                                   // TODO: we could throw away the (0,0)
                                   // block here since things have been
index 919fe922f97e1323c3f09b554d51473a5d10f903..93a5ce7be6b0e52fedd8d0c05f3c33607a6d5c1d 100755 (executable)
@@ -97,12 +97,11 @@ namespace TrilinosWrappers
                                        * format specified in
                                        * TrilinosWrappers::SparseMatrix.
                                        */
-      void initialize (const SparseMatrix        &matrix,
-                      const bool                 elliptic = true,
-                      const bool                 higher_order_elements = false,
-                      const std::vector<double> &null_space = std::vector<double>(),
-                      const unsigned int         null_space_dimension = 1,
-                      const bool                 output_details = false);
+      void initialize (const SparseMatrix                    &matrix,
+                      const bool                             elliptic = true,
+                      const bool                             higher_order_elements = false,
+                      const std::vector<std::vector<bool> > &null_space = std::vector<std::vector<bool> > (),
+                      const bool                            output_details = false);
 
                                       /**
                                        * Let Trilinos compute a
@@ -115,12 +114,11 @@ namespace TrilinosWrappers
                                        * can be considered rather 
                                        * inefficient.
                                        */
-      void initialize (const dealii::SparseMatrix<double> &matrix,
-                      const bool                 elliptic = true,
-                      const bool                 higher_order_elements = false,
-                      const std::vector<double> &null_space = std::vector<double>(),
-                      const unsigned int         null_space_dimension = 1,
-                      const bool                 output_details = false);
+      void initialize (const ::dealii::SparseMatrix<double>  &deal_ii_sparse_matrix,
+                      const bool                             elliptic = true,
+                      const bool                             higher_order_elements = false,
+                      const std::vector<std::vector<bool> > &null_space = std::vector<std::vector<bool> > (),
+                      const bool                            output_details = false);
 
                                       /**
                                        * This function can be used 
index 7d287c230d621730529ee01a95276ec11ae7fbb7..2ac0fa1db29e548aca058d910cec313363b5e60c 100755 (executable)
@@ -43,14 +43,14 @@ namespace TrilinosWrappers
 
   void
   PreconditionAMG::
-  initialize (const SparseMatrix        &matrix,
-             const bool                 elliptic,
-             const bool                 higher_order_elements,
-             const std::vector<double> &null_space,
-             const unsigned int         null_space_dimension,
-             const bool                 output_details)
+  initialize (const SparseMatrix                    &matrix,
+             const bool                             elliptic,
+             const bool                             higher_order_elements,
+             const std::vector<std::vector<bool> > &null_space,
+             const bool                             output_details)
   {
     const unsigned int n_rows = matrix.m();
+    const unsigned int null_space_dimension = null_space.size();
 
                                     // Build the AMG preconditioner.
     Teuchos::ParameterList parameter_list;
@@ -68,7 +68,7 @@ namespace TrilinosWrappers
        parameter_list.set("aggregation: block scaling", true);
       }
   
-    parameter_list.set("aggregation: threshold", 1e-12);
+    parameter_list.set("aggregation: threshold", 1e-8);
     
     if (output_details)
       parameter_list.set("ML output", 10);
@@ -77,16 +77,26 @@ namespace TrilinosWrappers
   
     if (higher_order_elements)
       parameter_list.set("aggregation: type", "MIS");
+    
+    std::vector<double> null_space_modes;
   
     if (null_space_dimension > 1)
       {
-       Assert (n_rows * null_space_dimension == null_space.size(),
-               ExcDimensionMismatch(n_rows * null_space_dimension,
-                                   null_space.size()));
+       Assert (n_rows == null_space[0].size(),
+               ExcDimensionMismatch(n_rows,
+                                    null_space[0].size()));
+       
+                                // Reshape null space as a contiguous
+                                // vector of doubles so that Trilinos
+                                // can read from it.
+       null_space_modes.resize (n_rows * null_space_dimension, 0.);
+       for (unsigned int d=0; d<null_space_dimension; ++d)
+         for (unsigned int row=0; row<n_rows; ++row)
+           null_space_modes[d*n_rows + row] = (double)null_space[d][row];
   
        parameter_list.set("null space: type", "pre-computed");
        parameter_list.set("null space: dimension", int(null_space_dimension));
-       parameter_list.set("null space: vectors", (double *)&null_space[0]);
+       parameter_list.set("null space: vectors", &null_space_modes[0]);
       }
 
     multigrid_operator = boost::shared_ptr<ML_Epetra::MultiLevelPreconditioner>
@@ -101,12 +111,11 @@ namespace TrilinosWrappers
 
   void
   PreconditionAMG::
-  initialize (const dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
-             const bool                  elliptic,
-             const bool                  higher_order_elements,
-             const std::vector<double>  &null_space,
-             const unsigned int          null_space_dimension,
-             const bool                  output_details)
+  initialize (const ::dealii::SparseMatrix<double>  &deal_ii_sparse_matrix,
+             const bool                             elliptic,
+             const bool                             higher_order_elements,
+             const std::vector<std::vector<bool> > &null_space,
+             const bool                             output_details)
   {
     const unsigned int n_rows = deal_ii_sparse_matrix.m();
   
@@ -122,7 +131,7 @@ namespace TrilinosWrappers
     Matrix->compress();
 
     initialize (*Matrix, elliptic, higher_order_elements, null_space,
-               null_space_dimension,  output_details);
+               output_details);
   }
   
   
index e04d300dbd35bb6d328818a7f9a0c2cbe00f9390..f83268fe4e9ee83d556ef883720eaffb9441d07c 100755 (executable)
@@ -253,6 +253,18 @@ namespace TrilinosWrappers
 
 
 
+  void
+  SparseMatrix::reinit (const SparseMatrix &sparse_matrix)
+  {
+    matrix.reset();
+    row_map = sparse_matrix.row_map;
+    col_map = sparse_matrix.col_map;
+    matrix = std::auto_ptr<Epetra_FECrsMatrix>(new Epetra_FECrsMatrix(
+                             *sparse_matrix.matrix));
+  }
+
+
+
   void
   SparseMatrix::reinit (const Epetra_Map                     &input_map,
                        const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
@@ -480,10 +492,10 @@ namespace TrilinosWrappers
        int diag_index = (int)(diag_find - col_indices);
 
        for (int j=0; j<num_entries; ++j)
-         if (diag_index != col_indices[j])
+         if (diag_index != j)
            values[j] = 0.;
 
-       if (diag_find && std::fabs(values[diag_index]) > 0.)
+       if (diag_find && std::fabs(values[diag_index]) == 0.)
          values[diag_index] = new_diag_value;
       }
   }
@@ -714,9 +726,14 @@ namespace TrilinosWrappers
 
   void
   SparseMatrix::vmult (Vector       &dst,
-                     const Vector &src) const
+                      const Vector &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
+    
+    Assert (col_map.SameAs(dst.map),
+           ExcMessage ("Column map of matrix does not fit with vector map!"));
+    Assert (row_map.SameAs(src.map),
+           ExcMessage ("Row map of matrix does not fit with vector map!"));
 
     if (!matrix->Filled())
       matrix->FillComplete();
@@ -733,6 +750,11 @@ namespace TrilinosWrappers
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
 
+    Assert (row_map.SameAs(dst.map),
+           ExcMessage ("Row map of matrix does not fit with vector map!"));
+    Assert (col_map.SameAs(src.map),
+           ExcMessage ("Column map of matrix does not fit with vector map!"));
+
     if (!matrix->Filled())
       matrix->FillComplete();
 

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.