Allow using ConstraintMatrix::shift if initialized with IndexSet

diff --git a/doc/news/changes/minor/20170103DanielArndt b/doc/news/changes/minor/20170103DanielArndt
new file mode 100644 (file)
index 0000000..90fb9d7
--- /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.
+<br>
+(Daniel Arndt, 2017/01/03)

diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h
index 1c7c99e72032fbdf71688474eebc537f410e73a0..edece7fb80257c51590d0473147dda980f02a4d6 100644 (file)
--- 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 90860a2c69f8e302d868203e4aaaf8ab36be1d8c..d41c60868cae8f77f3c94c2d281c7a3834f686e0 100644 (file)
--- 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<ConstraintLine>::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<lines_cache.size(); ++i)
+    Assert(lines_cache[i] == numbers::invalid_size_type ||
+           calculate_line_index(lines[lines_cache[i]].line) == i,
+           ExcInternalError());
+#endif
 }
 
 

diff --git a/tests/lac/constraints_shift_01.cc b/tests/lac/constraints_shift_01.cc
new file mode 100644 (file)
index 0000000..2031199
--- /dev/null
+++ b/tests/lac/constraints_shift_01.cc
@@ -0,0 +1,76 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test ConstraintMatrix::shift for a ConstraintMatrix object
+// initialized with an IndexSet object
+
+#include "../tests.h"
+#include <fstream>
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/vector.h>
+
+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<double> vec(size);
+  for (unsigned int i=0; i<size; ++i)
+    vec(i) = i;
+  constraints1.distribute(vec);
+  vec.print(deallog.get_file_stream(), 3, true, false);
+}
+
+int main()
+{
+  initlog();
+
+  test();
+
+  return 0;
+}

diff --git a/tests/lac/constraints_shift_01.output b/tests/lac/constraints_shift_01.output
new file mode 100644 (file)
index 0000000..1c15a08
--- /dev/null
+++ b/tests/lac/constraints_shift_01.output
@@ -0,0 +1,39 @@
+
+DEAL::Create index set with entries: {2, 5, 8}
+DEAL::Create ConstraintMatrix with constraints u(2)=.5*u(5), u(5)=.7*u(8)
+    5 8:  0.700000
+    2 5:  0.500000
+DEAL::Shifted constraints
+    15 18:  0.700000
+    12 15:  0.500000
+DEAL::Shifted and merged constraints
+    15 18:  0.700000
+    12 15:  0.500000
+    5 8:  0.700000
+    2 5:  0.500000
+DEAL::Close
+    2 8:  0.350000
+    5 8:  0.700000
+    12 18:  0.350000
+    15 18:  0.700000
+0.000e+00
+1.000e+00
+2.800e+00
+3.000e+00
+4.000e+00
+5.600e+00
+6.000e+00
+7.000e+00
+8.000e+00
+9.000e+00
+1.000e+01
+1.100e+01
+6.300e+00
+1.300e+01
+1.400e+01
+1.260e+01
+1.600e+01
+1.700e+01
+1.800e+01
+1.900e+01
+