From: Luca Heltai <luca.heltai@sissa.it>
Date: Fri, 13 Dec 2019 12:24:19 +0000 (+0100)
Subject: Contiguous tensor product IndexSet.
X-Git-Tag: v9.2.0-rc1~784^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=38d080f7020456db2f59570ff47be4bb9e61b26c;p=dealii.git

Contiguous tensor product IndexSet.
---

diff --git a/doc/news/changes/minor/20191212LucaHeltai b/doc/news/changes/minor/20191212LucaHeltai
new file mode 100644
index 0000000000..08443f5814
--- /dev/null
+++ b/doc/news/changes/minor/20191212LucaHeltai
@@ -0,0 +1,5 @@
+New: IndexSet::tensor_product() creates a new IndexSet with global size equal to
+this->size()*other.size(), containing for every element n of this IndexSet, the indices
+of the other IndexSet, contained in the interval [n*other.size(), (n+1)*other.size()).
+<br>
+(Luca Heltai, 2019/12/12)
diff --git a/include/deal.II/base/index_set.h b/include/deal.II/base/index_set.h
index 82f6f60d57..5c8224e5b1 100644
--- a/include/deal.II/base/index_set.h
+++ b/include/deal.II/base/index_set.h
@@ -352,6 +352,25 @@ public:
   void
   subtract_set(const IndexSet &other);
 
+  /**
+   * Return a new IndexSet, with global size equal to
+   * `this->size()*other.size()`, containing for every element `n` of this
+   * IndexSet, the entries in the half open range `[n*other.size(),
+   * (n+1)*other.size())` of the @p other IndexSet.
+   *
+   * The name results from the perspective that one starts with an IndexSet and
+   * takes the tensor product with another IndexSet with `other.size()`
+   * elements; this results in a matrix of size `this->size()` times
+   * `other.size()` that has ones in exactly the rows for which this IndexSet
+   * contained an index and in the columns for which the @p other IndexSet
+   * contained an index. This matrix is then "unrolled" again by going through
+   * each row one by one and reindexing the entries of the matrix in consecutive
+   * order. A one in the matrix then corresponds to an entry in the reindexed
+   * IndexSet that is returned by this function.
+   */
+  IndexSet
+  tensor_product(const IndexSet &other) const;
+
   /**
    * Remove and return the last element of the last range.
    * This function throws an exception if the IndexSet is empty.
diff --git a/source/base/index_set.cc b/source/base/index_set.cc
index 069d231f5b..0ab8372083 100644
--- a/source/base/index_set.cc
+++ b/source/base/index_set.cc
@@ -306,6 +306,18 @@ IndexSet::subtract_set(const IndexSet &other)
 
 
 
+IndexSet
+IndexSet::tensor_product(const IndexSet &other) const
+{
+  IndexSet set(this->size() * other.size());
+  for (const auto el : *this)
+    set.add_indices(other, el * other.size());
+  set.compress();
+  return set;
+}
+
+
+
 IndexSet::size_type
 IndexSet::pop_back()
 {
diff --git a/tests/base/index_set_contiguous_tensor_product_01.cc b/tests/base/index_set_contiguous_tensor_product_01.cc
new file mode 100644
index 0000000000..6568a29765
--- /dev/null
+++ b/tests/base/index_set_contiguous_tensor_product_01.cc
@@ -0,0 +1,48 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// test IndexSet::contiguous_tensor_product()
+
+#include <deal.II/base/index_set.h>
+
+#include "../tests.h"
+
+
+int
+main()
+{
+  initlog();
+  IndexSet set(10);
+  set.add_index(2);
+  set.add_range(4, 7);
+
+  set.print(deallog);
+
+  const auto tensor2 = set.tensor_product(complete_index_set(2));
+  deallog << "Size: " << tensor2.size() << ": ";
+  tensor2.print(deallog);
+
+
+  const auto tensor3 = set.tensor_product(complete_index_set(3));
+  deallog << "Size: " << tensor3.size() << ": ";
+  tensor3.print(deallog);
+
+  IndexSet other(5);
+  other.add_range(0, 3);
+  const auto tensor4 = set.tensor_product(other);
+  deallog << "Size: " << tensor4.size() << ": ";
+  tensor4.print(deallog);
+}
diff --git a/tests/base/index_set_contiguous_tensor_product_01.output b/tests/base/index_set_contiguous_tensor_product_01.output
new file mode 100644
index 0000000000..9a3edb9170
--- /dev/null
+++ b/tests/base/index_set_contiguous_tensor_product_01.output
@@ -0,0 +1,5 @@
+
+DEAL::{2, [4,6]}
+DEAL::Size: 20: {[4,5], [8,13]}
+DEAL::Size: 30: {[6,8], [12,20]}
+DEAL::Size: 50: {[10,12], [20,22], [25,27], [30,32]}