From 9f7645a370a1d8f098105a44e066556dbe85f1c2 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 3 Jan 2017 16:47:53 +0100 Subject: [PATCH] Allow using ConstraintMatrix::shift if initialized with IndexSet --- doc/news/changes/minor/20170103DanielArndt | 4 ++ include/deal.II/lac/constraint_matrix.h | 3 +- source/lac/constraint_matrix.cc | 24 +++++-- tests/lac/constraints_shift_01.cc | 76 ++++++++++++++++++++++ tests/lac/constraints_shift_01.output | 39 +++++++++++ 5 files changed, 140 insertions(+), 6 deletions(-) create mode 100644 doc/news/changes/minor/20170103DanielArndt create mode 100644 tests/lac/constraints_shift_01.cc create mode 100644 tests/lac/constraints_shift_01.output diff --git a/doc/news/changes/minor/20170103DanielArndt b/doc/news/changes/minor/20170103DanielArndt new file mode 100644 index 0000000000..90fb9d735f --- /dev/null +++ b/doc/news/changes/minor/20170103DanielArndt @@ -0,0 +1,4 @@ +Fixed: ConstraintMatrix::shift also works if the object +is initialized with an IndexSet object. +
+(Daniel Arndt, 2017/01/03) diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h index 1c7c99e720..edece7fb80 100644 --- a/include/deal.II/lac/constraint_matrix.h +++ b/include/deal.II/lac/constraint_matrix.h @@ -365,7 +365,8 @@ public: /** * Shift all entries of this matrix down @p offset rows and over @p offset - * columns. + * columns. If this object is initialized with an IndexSet, local_lines are + * shifted as well. * * This function is useful if you are building block matrices, where all * blocks are built by the same DoFHandler object, i.e. the matrix size is diff --git a/source/lac/constraint_matrix.cc b/source/lac/constraint_matrix.cc index 90860a2c69..d41c60868c 100644 --- a/source/lac/constraint_matrix.cc +++ b/source/lac/constraint_matrix.cc @@ -591,11 +591,16 @@ ConstraintMatrix::merge (const ConstraintMatrix &other_constraints, void ConstraintMatrix::shift (const size_type offset) { - //TODO: this doesn't work with IndexSets yet. [TH] - AssertThrow(local_lines.size()==0, ExcNotImplemented()); - - lines_cache.insert (lines_cache.begin(), offset, - numbers::invalid_size_type); + if (local_lines.size() == 0) + lines_cache.insert (lines_cache.begin(), offset, + numbers::invalid_size_type); + else + { + // shift local_lines + IndexSet new_local_lines(local_lines.size()); + new_local_lines.add_indices(local_lines, offset); + std::swap(local_lines, new_local_lines); + } for (std::vector::iterator i = lines.begin(); i != lines.end(); ++i) @@ -606,6 +611,15 @@ void ConstraintMatrix::shift (const size_type offset) j != i->entries.end(); ++j) j->first += offset; } + +#ifdef DEBUG + // make sure that lines, lines_cache and local_lines + // are still linked correctly + for (size_type i=0; i + +#include +#include +#include + +using namespace dealii; + +void test() +{ + const int size = 20; + + dealii::IndexSet index_set(size); + dealii::IndexSet set1(size); + deallog << "Create index set with entries: "; + index_set.add_index(2); + index_set.add_index(5); + index_set.add_index(8); + index_set.print(deallog); + + deallog << "Create ConstraintMatrix with constraints u(2)=.5*u(5), u(5)=.7*u(8)" << std::endl; + dealii::ConstraintMatrix constraints1(index_set); + constraints1.add_line(5); + constraints1.add_entry(5, 8, .7); + constraints1.add_line(2); + constraints1.add_entry(2, 5, .5); + dealii::ConstraintMatrix constraints2(constraints1); + constraints1.print(deallog.get_file_stream()); + + constraints1.shift(size/2); + deallog << "Shifted constraints" << std::endl; + constraints1.print(deallog.get_file_stream()); + + constraints1.merge(constraints2, ConstraintMatrix::no_conflicts_allowed, true); + deallog << "Shifted and merged constraints" << std::endl; + constraints1.print(deallog.get_file_stream()); + + constraints1.close(); + deallog << "Close" << std::endl; + constraints1.print(deallog.get_file_stream()); + + Vector vec(size); + for (unsigned int i=0; i