From 8f02dbf178bf791fc00387b33ec6a6800048aad8 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Thu, 29 Nov 2012 16:05:54 +0000 Subject: [PATCH] Allow ConstraintMatrix::distribute_local_to_global with ChunkSparseMatrix. git-svn-id: https://svn.dealii.org/trunk@27703 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/deal.II/lac/chunk_sparse_matrix.h | 55 +++++++- deal.II/source/lac/constraint_matrix.cc | 4 + .../constraints_local_to_global_chunk.cc | 125 ++++++++++++++++++ .../cmp/generic | 4 + 4 files changed, 182 insertions(+), 6 deletions(-) create mode 100644 tests/deal.II/constraints_local_to_global_chunk.cc create mode 100644 tests/deal.II/constraints_local_to_global_chunk/cmp/generic diff --git a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h index 4784340b14..7f47346c8e 100644 --- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h @@ -378,6 +378,31 @@ public: const unsigned int j, const number value); + /** + * Add an array of values given by + * values in the given + * global matrix row at columns + * specified by col_indices in the + * sparse matrix. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + template + void add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number2 *values, + const bool elide_zero_values = true, + const bool col_indices_are_sorted = false); + /** * Multiply the entire matrix by a * fixed factor. @@ -1277,13 +1302,31 @@ void ChunkSparseMatrix::add (const unsigned int i, Assert (cols != 0, ExcNotInitialized()); - const unsigned int index = compute_location(i,j); - Assert ((index != ChunkSparsityPattern::invalid_entry) || - (value == 0.), - ExcInvalidIndex(i,j)); - if (value != 0.) - val[index] += value; + { + const unsigned int index = compute_location(i,j); + Assert ((index != ChunkSparsityPattern::invalid_entry), + ExcInvalidIndex(i,j)); + + val[index] += value; + } +} + + + +template +template +inline +void ChunkSparseMatrix::add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number2 *values, + const bool elide_zero_values, + const bool col_indices_are_sorted) +{ + // TODO: could be done more efficiently... + for (unsigned int col=0; col(values[col])); } diff --git a/deal.II/source/lac/constraint_matrix.cc b/deal.II/source/lac/constraint_matrix.cc index b652bbe942..aa0e4f72fa 100644 --- a/deal.II/source/lac/constraint_matrix.cc +++ b/deal.II/source/lac/constraint_matrix.cc @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -2373,7 +2374,10 @@ BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix, BlockVector); MATRIX_FUNCTIONS(SparseMatrixEZ); MATRIX_FUNCTIONS(SparseMatrixEZ); +MATRIX_FUNCTIONS(ChunkSparseMatrix); +MATRIX_FUNCTIONS(ChunkSparseMatrix); MATRIX_VECTOR_FUNCTIONS(SparseMatrixEZ, Vector); +MATRIX_VECTOR_FUNCTIONS(ChunkSparseMatrix, Vector); // BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrixEZ); // BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrixEZ, Vector); diff --git a/tests/deal.II/constraints_local_to_global_chunk.cc b/tests/deal.II/constraints_local_to_global_chunk.cc new file mode 100644 index 0000000000..5638055142 --- /dev/null +++ b/tests/deal.II/constraints_local_to_global_chunk.cc @@ -0,0 +1,125 @@ +//------------------ constraints_local_to_global_chunk.cc ------------------ +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2012 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//------------------ constraints_local_to_global_chunk.cc ------------------ + + +// this function tests the correctness of the implementation of +// ConstraintMatrix::distribute_local_to_global for ChunkSparseMatrix by +// comparing the results with a sparse matrix. As a test case, we use a square +// mesh that is refined once globally and then the first cell is refined +// adaptively. + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +std::ofstream logfile("constraints_local_to_global_chunk/output"); + +template +void test (unsigned int chunk_size) +{ + Triangulation tria; + GridGenerator::hyper_cube (tria); + tria.begin()->face(0)->set_boundary_indicator(1); + tria.refine_global(1); + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + FE_Q fe (1); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values (dof, 1, ZeroFunction(), + constraints); + constraints.close(); + + SparsityPattern sparsity; + ChunkSparsityPattern chunk_sparsity; + { + CompressedSimpleSparsityPattern csp (dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, csp, constraints, false); + sparsity.copy_from (csp); + chunk_sparsity.copy_from (csp, chunk_size); + } + SparseMatrix sparse (sparsity); + ChunkSparseMatrix chunk_sparse (chunk_sparsity); + + FullMatrix local_mat (fe.dofs_per_cell, fe.dofs_per_cell); + std::vector local_dof_indices (fe.dofs_per_cell); + + // loop over cells, fill local matrix with + // random values, insert both into sparse and + // full matrix. Make some random entries equal + // to zero + typename DoFHandler::active_cell_iterator + cell = dof.begin_active(), endc = dof.end(); + unsigned int counter = 0; + for ( ; cell != endc; ++cell) + { + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_mat, local_dof_indices, + sparse); + constraints.distribute_local_to_global (local_mat, local_dof_indices, + chunk_sparse); + } + + // now check that the entries are indeed the + // same + double frobenius = 0.; + for (unsigned int i=0; i::abs_square(sparse.el(i,j) - + chunk_sparse.el(i,j)); + deallog << "Difference between chunk and sparse matrix: " + << std::sqrt(frobenius) << std::endl; +} + + +int main () +{ + deallog << std::setprecision (2); + logfile << std::setprecision (2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-14); + + test<2>(1); + test<2>(2); + test<2>(5); +} + diff --git a/tests/deal.II/constraints_local_to_global_chunk/cmp/generic b/tests/deal.II/constraints_local_to_global_chunk/cmp/generic new file mode 100644 index 0000000000..52403f3e9d --- /dev/null +++ b/tests/deal.II/constraints_local_to_global_chunk/cmp/generic @@ -0,0 +1,4 @@ + +DEAL::Difference between chunk and sparse matrix: 0 +DEAL::Difference between chunk and sparse matrix: 0 +DEAL::Difference between chunk and sparse matrix: 0 -- 2.39.5