From: Luca Heltai 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()). +
+(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 + +#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]}