]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added the function DoFTools::extract_constant_modes(). This simplifies the preconditi...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 1 Sep 2008 12:46:00 +0000 (12:46 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 1 Sep 2008 12:46:00 +0000 (12:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@16701 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

index cd67ee3180a6ee90e1797afb4708ccf77d7939db..1a1982eccdde2313e5fa5b6e8c358eb2613ee078 100644 (file)
@@ -1008,6 +1008,70 @@ class DoFTools
                            const unsigned int  subdomain_id,
                            std::vector<bool>  &selected_dofs);
 
+                                    /**
+                                     * Extract a vector that 
+                                     * represents the constant
+                                     * modes of the DoFHandler
+                                     * for the components 
+                                     * chosen by <tt>component_select</tt>.
+                                     * The constant modes on a
+                                     * discretization are the 
+                                     * null space of a Laplace
+                                     * operator on the selected
+                                     * components with Neumann
+                                     * boundary conditions 
+                                     * applied. The null space is
+                                     * a necessary ingredient for
+                                     * obtaining a good AMG
+                                     * preconditioner when using
+                                     * the class 
+                                     * TrilinosWrappers::PreconditionAMG.
+                                     * Since the ML AMG package
+                                     * only works on algebraic
+                                     * properties of the respective
+                                     * matrix, it has no chance 
+                                     * to detect whether the matrix
+                                     * comes from a scalar or a
+                                     * vector valued problem. However,
+                                     * a near null space supplies
+                                     * exactly the needed information
+                                     * about these components.
+                                     * The null space will consist
+                                     * 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
+                                     * 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>.
+                                     * 
+                                     * The main reason for this
+                                     * program is the use of the
+                                     * null space with the
+                                     * AMG preconditioner.
+                                     */
+    template <class DH>
+    static void
+    extract_constant_modes (const DH                &dof_handler,
+                           const std::vector<bool> &component_select,
+                           std::vector<double>     &constant_modes);
+
                                     /**
                                      * For each active cell of a DoFHandler
                                      * or hp::DoFHandler, extract the active
index ffb981acd2d692036962468a9875bbdc30446ffe..19372180ac4ff75be635aa2d7845a686da4d0b52 100644 (file)
@@ -3613,7 +3613,60 @@ 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)
+{
+  const unsigned int n_components = dof_handler.get_fe().n_components();
+  Assert (n_components == component_select.size(),
+         ExcDimensionMismatch(n_components,
+                              component_select.size()));
+
+                                // First count the number of dofs
+                                // in the current component.
+  unsigned int n_components_selected = 0;
+  std::vector<unsigned int> component_list (n_components, 0);
+  for (unsigned int d=0; d<n_components; ++d)
+    {
+      component_list[d] = component_select[d];
+      n_components_selected += component_select[d];
+    }
+
+  std::vector<unsigned int> dofs_per_block(2);
+  count_dofs_per_block(dof_handler, dofs_per_block, component_list);
+
+  const unsigned int n_u = dofs_per_block[1];
+  std::vector<bool> selection_dof_list (dof_handler.n_dofs(), false);
+  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.);
+  
+  for (unsigned int component=0, component_used=0; 
+       component < n_components; ++component, ++component_used)
+    if (component_select[component])
+      {
+       std::vector<bool> selection_mask (n_components, false);
+       selection_mask[component] = true;
+       extract_dofs (dof_handler, selection_mask, temporary_dof_list);
+  
+       unsigned int counter = 0;
+       for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+       {
+         if (temporary_dof_list[i])
+           {
+             constant_modes [component * n_u + counter] = 1.;
+           }
+         if (selection_dof_list[i])
+           ++counter;
+       }
+      }
+}
+
+
+
 template <class DH>
 void
 DoFTools::get_active_fe_indices (const DH                  &dof_handler,
@@ -5550,6 +5603,13 @@ DoFTools::extract_subdomain_dofs<hp::DoFHandler<deal_II_dimension> >
  const unsigned int     subdomain_id,
  std::vector<bool>     &selected_dofs);
 
+template
+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);
+
 template
 void
 DoFTools::get_active_fe_indices<DoFHandler<deal_II_dimension> >
index 49c44bc3258878ffe8e75cf662e70000cc4a3467..c2ac126e028c9b3ed57a3a99e031081c84c2b9ba 100644 (file)
@@ -1188,42 +1188,12 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
       
   Amg_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionAMG>
                       (new TrilinosWrappers::PreconditionAMG());
-      
-  const unsigned int n_u = stokes_preconditioner_matrix.block(0,0).m();
-  std::vector<double> null_space (dim * n_u, 0.);
-      
-  std::vector<bool> precondition_dof_list (stokes_dof_handler.n_dofs(), false);
-      
-  for (unsigned int component=0; component < dim; ++component)
-    {
-      std::vector<bool> precondition_mask (dim + 1, false);
-      precondition_mask[component] = true;
-      DoFTools::extract_dofs (stokes_dof_handler, precondition_mask, 
-                             precondition_dof_list);
-
-                                      // TODO: The current implementation 
-                                      // assumes that we are working on 
-                                      // the first components of a system when
-                                      // writing into the null vector.
-                                      // Change this to the general case,
-                                      // probably use something similar as
-                                      // for block vectors.
-      unsigned int counter = 0;
-      for (unsigned int i=0; i<stokes_dof_handler.n_dofs(); ++i)
-       {
-         if (precondition_dof_list[i])
-           {
-             Assert(i < n_u,
-                    ExcMessage("Could not correctly locate "
-                               "preconditioner dofs in system!"));
-             null_space [component * n_u + i] = 1.;
-             ++counter;
-           }
-       }
-      Assert (counter == n_u / dim,
-             ExcDimensionMismatch(counter, n_u / dim));
-    }
-       
+
+  std::vector<double> null_space;
+  std::vector<bool>  velocity_components (dim+1,true);
+  velocity[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);
 

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.