* 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
*/
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
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(),
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)
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;
}
}
}
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
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
* 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
* 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
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;
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);
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>
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();
Matrix->compress();
initialize (*Matrix, elliptic, higher_order_elements, null_space,
- null_space_dimension, output_details);
+ output_details);
}
+ 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,
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;
}
}
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();
{
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();