* 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<LinearAlgebra::distributed::Vector<double> > diagonal_matrix;
+ * LinearAlgebra::distributed::Vector<double> &diagonal_vector = diagonal_matrix.get_vector();
+ * diagonal_vector.reinit(locally_owned_dofs,
+ * locally_relevant_dofs,
+ * mpi_communicator);
+ * @code
+ *
* @author Martin Kronbichler, 2016
*/
template <typename VectorType=Vector<double> >
*/
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.
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.
*/
+template <typename VectorType>
+void
+DiagonalMatrix<VectorType>::compress (VectorOperation::values operation)
+{
+ diagonal.compress(operation);
+}
+
+
+
template <typename VectorType>
VectorType &
DiagonalMatrix<VectorType>::get_vector()
+template <typename VectorType>
+void
+DiagonalMatrix<VectorType>::add (const size_type i,
+ const size_type j,
+ const value_type value)
+{
+ if (i == j)
+ diagonal(i) += value;
+}
+
+
+
template <typename VectorType>
void
DiagonalMatrix<VectorType>::vmult(VectorType &dst,
template void ConstraintMatrix::distribute_local_to_global<LinearAlgebra::distributed::T<S> > (
const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, LinearAlgebra::distributed::T<S> &, const FullMatrix<S>&, bool) const;
template void ConstraintMatrix::set_zero<LinearAlgebra::distributed::T<S> >(LinearAlgebra::distributed::T<S> &) const;
+ template void ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> > > (
+ const FullMatrix<S> &, const std::vector< size_type > &, DiagonalMatrix<LinearAlgebra::distributed::T<S> > &) const;
+ template void ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> >, LinearAlgebra::distributed::T<S> > (
+ const FullMatrix<S> &, const Vector<S>&, const std::vector< size_type > &,
+ DiagonalMatrix<LinearAlgebra::distributed::T<S> > &, LinearAlgebra::distributed::T<S>&,
+ bool, dealii::internal::bool2type<false>) const;
+ template void ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> >, T<S> > (
+ const FullMatrix<S> &, const Vector<S>&, const std::vector< size_type > &,
+ DiagonalMatrix<LinearAlgebra::distributed::T<S> > &, T<S>&,
+ bool, dealii::internal::bool2type<false>) const;
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/diagonal_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+
+
+
+template <int dim, int fe_degree>
+void test (const bool hanging_nodes = true)
+{
+ typedef double number;
+
+ parallel::distributed::Triangulation<dim> 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<dim> fe (fe_degree);
+ DoFHandler<dim> 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<dim>(),
+ 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<LinearAlgebra::distributed::Vector<double> > diagonal_matrix;
+ {
+ LinearAlgebra::distributed::Vector<double> dummy;
+ dummy.reinit(owned_set, relevant_set, MPI_COMM_WORLD);
+ diagonal_matrix.reinit(dummy);
+ }
+
+ {
+ QGauss<dim> quadrature_formula(fe_degree+1);
+
+ FEValues<dim> 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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::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_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point)) *
+ fe_values.JxW(q_point);
+ }
+
+ cell->get_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<owned_set.n_elements(); ++i)
+ {
+ const unsigned int glob_index =
+ owned_set.nth_index_in_set (i);
+ if (constraints.is_constrained(glob_index))
+ continue;
+
+ const double d = diagonal_matrix(glob_index, glob_index);
+ const double t = sparse_matrix.diag_element(glob_index);
+ if ( std::abs(d - t)/t > 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();
+}
+
projects / dealii.git / commitdiff