From bd99433a0be48d9c2f46e29f73cb4420f12782eb Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 2 Sep 2008 10:33:02 +0000 Subject: [PATCH] Changed DoFTools::extract_constant_modes to work on bool vectors instead of doubles to make it consistent with the other extract_* functions. git-svn-id: https://svn.dealii.org/trunk@16705 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/dofs/dof_tools.h | 32 +++++++------ deal.II/deal.II/source/dofs/dof_tools.cc | 21 +++++---- deal.II/examples/step-31/step-31.cc | 4 +- .../include/lac/trilinos_precondition_amg.h | 22 +++++---- .../lac/source/trilinos_precondition_amg.cc | 45 +++++++++++-------- deal.II/lac/source/trilinos_sparse_matrix.cc | 28 ++++++++++-- 6 files changed, 91 insertions(+), 61 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index 1a1982eccd..2634881624 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -1040,26 +1040,24 @@ class DoFTools * of as many vectors as there * are true arguments in * component_select, - * say d, * each of which will be one * in one component and - * zero in all others. We store - * the null space contiguously - * in the output vector - * constant_modes - * to enable the internal - * Trilinos structures the - * access to it. This means that - * the length of the vector - * will be d 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 - * component_select. + * constructed using the + * same + * component_select + * argument. * * The main reason for this * program is the use of the @@ -1068,9 +1066,9 @@ class DoFTools */ template static void - extract_constant_modes (const DH &dof_handler, - const std::vector &component_select, - std::vector &constant_modes); + extract_constant_modes (const DH &dof_handler, + const std::vector &component_select, + std::vector > &constant_modes); /** * For each active cell of a DoFHandler diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 19372180ac..998f451bed 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -3615,9 +3615,9 @@ DoFTools::extract_subdomain_dofs (const DH &dof_handler, template void -DoFTools::extract_constant_modes (const DH &dof_handler, - const std::vector &component_select, - std::vector &constant_modes) +DoFTools::extract_constant_modes (const DH &dof_handler, + const std::vector &component_select, + std::vector > &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 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(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 > (const DoFHandler &dof_handler, const std::vector &selected_components, - std::vector &constant_modes); + std::vector > &constant_modes); template void diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 170fe0a8d0..cfd21ca5e1 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -1189,13 +1189,13 @@ BoussinesqFlowProblem::build_stokes_preconditioner () Amg_preconditioner = boost::shared_ptr (new TrilinosWrappers::PreconditionAMG()); - std::vector null_space; + std::vector > null_space; std::vector 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 diff --git a/deal.II/lac/include/lac/trilinos_precondition_amg.h b/deal.II/lac/include/lac/trilinos_precondition_amg.h index 919fe922f9..93a5ce7be6 100755 --- a/deal.II/lac/include/lac/trilinos_precondition_amg.h +++ b/deal.II/lac/include/lac/trilinos_precondition_amg.h @@ -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 &null_space = std::vector(), - 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 > &null_space = std::vector > (), + 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 &matrix, - const bool elliptic = true, - const bool higher_order_elements = false, - const std::vector &null_space = std::vector(), - const unsigned int null_space_dimension = 1, - const bool output_details = false); + void initialize (const ::dealii::SparseMatrix &deal_ii_sparse_matrix, + const bool elliptic = true, + const bool higher_order_elements = false, + const std::vector > &null_space = std::vector > (), + const bool output_details = false); /** * This function can be used diff --git a/deal.II/lac/source/trilinos_precondition_amg.cc b/deal.II/lac/source/trilinos_precondition_amg.cc index 7d287c230d..2ac0fa1db2 100755 --- a/deal.II/lac/source/trilinos_precondition_amg.cc +++ b/deal.II/lac/source/trilinos_precondition_amg.cc @@ -43,14 +43,14 @@ namespace TrilinosWrappers void PreconditionAMG:: - initialize (const SparseMatrix &matrix, - const bool elliptic, - const bool higher_order_elements, - const std::vector &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 > &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 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 @@ -101,12 +111,11 @@ namespace TrilinosWrappers void PreconditionAMG:: - initialize (const dealii::SparseMatrix &deal_ii_sparse_matrix, - const bool elliptic, - const bool higher_order_elements, - const std::vector &null_space, - const unsigned int null_space_dimension, - const bool output_details) + initialize (const ::dealii::SparseMatrix &deal_ii_sparse_matrix, + const bool elliptic, + const bool higher_order_elements, + const std::vector > &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); } diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index e04d300dbd..f83268fe4e 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -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(new Epetra_FECrsMatrix( + *sparse_matrix.matrix)); + } + + + void SparseMatrix::reinit (const Epetra_Map &input_map, const ::dealii::SparseMatrix &dealii_sparse_matrix, @@ -480,10 +492,10 @@ namespace TrilinosWrappers int diag_index = (int)(diag_find - col_indices); for (int j=0; j 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(); -- 2.39.5