From 55fcd6917549c11207b176856b44689512a79468 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 2 Nov 2016 10:02:40 +0100 Subject: [PATCH] instantiate ConstraintMatrix::distribute_local_to_global() for DiagonalMatrix To that end, add a few missing functions to DiagonalMatrix and update documentation. --- include/deal.II/lac/diagonal_matrix.h | 51 +++++ source/lac/constraint_matrix.cc | 1 + source/lac/constraint_matrix.inst.in | 10 + tests/lac/diagonal_matrix_02.cc | 175 ++++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=1.output | 9 + ...h_mpi=true.with_p4est=true.mpirun=3.output | 9 + ...h_mpi=true.with_p4est=true.mpirun=7.output | 9 + 7 files changed, 264 insertions(+) create mode 100644 tests/lac/diagonal_matrix_02.cc create mode 100644 tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output create mode 100644 tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=3.output create mode 100644 tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=7.output diff --git a/include/deal.II/lac/diagonal_matrix.h b/include/deal.II/lac/diagonal_matrix.h index ff9e27b8fd..ad08af7bfc 100644 --- a/include/deal.II/lac/diagonal_matrix.h +++ b/include/deal.II/lac/diagonal_matrix.h @@ -27,6 +27,19 @@ DEAL_II_NAMESPACE_OPEN * VectorType::scale, so the template vector class needs to provide a * @p scale() method. * + * When using this class with ConstraintsMatrix::distribute_local_to_global(), + * the underlying vector needs to provide write access to all entries referenced + * by cells in an assembly process. This means that this class also needs access + * to ghost entries that are owned by other processors than the calling one. + * In practice this requires initialization of the vector as follows + * @code + * DiagonalMatrix > diagonal_matrix; + * LinearAlgebra::distributed::Vector &diagonal_vector = diagonal_matrix.get_vector(); + * diagonal_vector.reinit(locally_owned_dofs, + * locally_relevant_dofs, + * mpi_communicator); + * @code + * * @author Martin Kronbichler, 2016 */ template > @@ -47,6 +60,13 @@ public: */ void reinit (const VectorType &vec); + /** + * Compresses the data structures and allows the resulting matrix to be used + * in all other operations like matrix-vector products. This is a collective + * operation, i.e., it needs to be run on all processors when used in parallel. + */ + void compress (VectorOperation::values operation); + /** * Returns a reference to the underlying vector for manipulation of the * entries on the matrix diagonal. @@ -118,6 +138,16 @@ public: const bool elide_zero_values = true, const bool col_indices_are_sorted = false); + /** + * Add value to the element (i,j). + * + * Due to the storage of this matrix, entries are only added to the diagonal + * of the matrix. All other entries are ignored and no exception is thrown. + */ + void add (const size_type i, + const size_type j, + const value_type value); + /** * Performs a matrix-vector multiplication with the given matrix. */ @@ -199,6 +229,15 @@ DiagonalMatrix::reinit(const VectorType &vec) +template +void +DiagonalMatrix::compress (VectorOperation::values operation) +{ + diagonal.compress(operation); +} + + + template VectorType & DiagonalMatrix::get_vector() @@ -274,6 +313,18 @@ DiagonalMatrix::add (const size_type row, +template +void +DiagonalMatrix::add (const size_type i, + const size_type j, + const value_type value) +{ + if (i == j) + diagonal(i) += value; +} + + + template void DiagonalMatrix::vmult(VectorType &dst, diff --git a/source/lac/constraint_matrix.cc b/source/lac/constraint_matrix.cc index db04b3a4d2..90860a2c69 100644 --- a/source/lac/constraint_matrix.cc +++ b/source/lac/constraint_matrix.cc @@ -39,6 +39,7 @@ #include #include #include +#include #include #include diff --git a/source/lac/constraint_matrix.inst.in b/source/lac/constraint_matrix.inst.in index 4776d349ac..6587bc5606 100644 --- a/source/lac/constraint_matrix.inst.in +++ b/source/lac/constraint_matrix.inst.in @@ -34,6 +34,16 @@ for (S: REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES) template void ConstraintMatrix::distribute_local_to_global > ( const Vector&, const std::vector &, const std::vector &, LinearAlgebra::distributed::T &, const FullMatrix&, bool) const; template void ConstraintMatrix::set_zero >(LinearAlgebra::distributed::T &) const; + template void ConstraintMatrix::distribute_local_to_global > > ( + const FullMatrix &, const std::vector< size_type > &, DiagonalMatrix > &) const; + template void ConstraintMatrix::distribute_local_to_global >, LinearAlgebra::distributed::T > ( + const FullMatrix &, const Vector&, const std::vector< size_type > &, + DiagonalMatrix > &, LinearAlgebra::distributed::T&, + bool, dealii::internal::bool2type) const; + template void ConstraintMatrix::distribute_local_to_global >, T > ( + const FullMatrix &, const Vector&, const std::vector< size_type > &, + DiagonalMatrix > &, T&, + bool, dealii::internal::bool2type) const; } diff --git a/tests/lac/diagonal_matrix_02.cc b/tests/lac/diagonal_matrix_02.cc new file mode 100644 index 0000000000..1e57f6e9d4 --- /dev/null +++ b/tests/lac/diagonal_matrix_02.cc @@ -0,0 +1,175 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// this tests the correctness of distribute_local_to_globa() for DiagonalMatrix + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void test (const bool hanging_nodes = true) +{ + typedef double number; + + parallel::distributed::Triangulation tria (MPI_COMM_WORLD); + GridGenerator::hyper_cube (tria); + tria.refine_global(1); + + if (hanging_nodes) + { + if (tria.begin_active()->is_locally_owned()) + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement (); + } + + FE_Q fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + IndexSet owned_set = dof.locally_owned_dofs(); + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dof, relevant_set); + + ConstraintMatrix constraints (relevant_set); + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction(), + constraints); + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + //std::cout << "Number of cells: " << tria.n_global_active_cells() << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + + // assemble trilinos sparse matrix + TrilinosWrappers::SparseMatrix sparse_matrix; + { + TrilinosWrappers::SparsityPattern csp (owned_set, MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (dof, csp, constraints, true, + Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)); + csp.compress(); + sparse_matrix.reinit (csp); + } + + DiagonalMatrix > diagonal_matrix; + { + LinearAlgebra::distributed::Vector dummy; + dummy.reinit(owned_set, relevant_set, MPI_COMM_WORLD); + diagonal_matrix.reinit(dummy); + } + + { + QGauss quadrature_formula(fe_degree+1); + + FEValues fe_values (dof.get_fe(), quadrature_formula, + update_gradients | update_JxW_values); + + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + cell_matrix = 0; + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + sparse_matrix); + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + diagonal_matrix); + } + } + sparse_matrix.compress(VectorOperation::add); + diagonal_matrix.compress(VectorOperation::add); + + // compare elements: + for (unsigned int i=0; i 1e-10 ) + { + deallog << glob_index << " " << d << " != " << t << std::endl; + } + } + + deallog << "Ok." << std::endl; +} + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads()); + + mpi_initlog(); + + deallog.push("2d"); + test<2,1>(); + test<2,2>(); + deallog.pop(); + + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); +} + diff --git a/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..c6a564409d --- /dev/null +++ b/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,9 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Ok. +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d::Ok. +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Ok. +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Ok. diff --git a/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..c6a564409d --- /dev/null +++ b/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=3.output @@ -0,0 +1,9 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Ok. +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d::Ok. +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Ok. +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Ok. diff --git a/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=7.output b/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=7.output new file mode 100644 index 0000000000..c6a564409d --- /dev/null +++ b/tests/lac/diagonal_matrix_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=7.output @@ -0,0 +1,9 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Ok. +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d::Ok. +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Ok. +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Ok. -- 2.39.5