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
-
+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,
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> >
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);