local sparsity pattern for Q_iso_Q1

diff --git a/doc/news/changes/minor/20241120tjhei b/doc/news/changes/minor/20241120tjhei
new file mode 100644 (file)
index 0000000..457f7f8
--- /dev/null
+++ b/doc/news/changes/minor/20241120tjhei
@@ -0,0 +1,7 @@
+New: FiniteElement::get_local_dof_sparsity_pattern() can be used to
+provide the coupling between DoFs within a cell and is used in
+make_sparsity_pattern() to generate matrices with fewer nonzero
+entries depending on the element. FE_Q_iso_Q1 now provides this
+coupling information.
+<br>
+(Timo Heister, Luca Heltai, 2024/11/20)

diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h
index 12f972276bd48d7b0cccea398d4c8a2b23ab8b7d..9ba23410461c77605f5970b9d0e082e216bb0af1 100644 (file)
--- a/include/deal.II/fe/fe.h
+++ b/include/deal.II/fe/fe.h
@@ -1610,6 +1610,25 @@ public:
   bool
   is_primitive(const unsigned int i) const;
 
+  /**
+   * Return a table (or sparsity pattern) with information which of the
+   * local DoFs per cell couple with each other.
+   *
+   * If the FiniteElement returns an empty (0 by 0) table, which is the
+   * default implementation for any class that does not override this
+   * function, all DoFs are assumed to couple.
+   *
+   * Otherwise, this function returns a square table with n_dofs_per_cell
+   * rows and columns. The entry (i,j) then denotes if the DoF i and j
+   * are coupled. This should be true if the support of the shape functions
+   * have a non-empty intersection.
+   *
+   * An example for an element that does make use of this feature is
+   * FE_Q_iso_Q1.
+   */
+  virtual const Table<2, bool> &
+  get_local_dof_sparsity_pattern() const;
+
   /**
    * Number of base elements in a mixed discretization.
    *
@@ -2621,6 +2640,14 @@ protected:
    */
   const bool cached_primitivity;
 
+  /**
+   * This table is returned by get_local_dof_sparsity_pattern(). Its meaning
+   * is described there in detail. This variable has to be
+   * filled by any class derived from FiniteElement that does not want to keep
+   * the default empty table (meaning all DoFs couple within the cell).
+   */
+  Table<2, bool> local_dof_sparsity_pattern;
+
   /**
    * Return the size of interface constraint matrices. Since this is
    * needed in every derived finite element class when initializing
@@ -3345,6 +3372,15 @@ FiniteElement<dim, spacedim>::is_primitive(const unsigned int i) const
 
 
 
+template <int dim, int spacedim>
+inline const Table<2, bool> &
+FiniteElement<dim, spacedim>::get_local_dof_sparsity_pattern() const
+{
+  return local_dof_sparsity_pattern;
+}
+
+
+
 template <int dim, int spacedim>
 inline GeometryPrimitive
 FiniteElement<dim, spacedim>::get_associated_geometry_primitive(

diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc
index f7f0b74b652bd075811a14fdaa1f9de942eb948c..bfce5678e1216763a5cba70df5933167c7d82c2c 100644 (file)
--- a/source/dofs/dof_tools_sparsity.cc
+++ b/source/dofs/dof_tools_sparsity.cc
@@ -82,6 +82,13 @@ namespace DoFTools
                  "locally owned one does not make sense."));
       }
 
+    const auto                 &fe_collection = dof.get_fe_collection();
+    std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
+    for (unsigned int f = 0; f < fe_collection.size(); ++f)
+      {
+        fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
+      }
+
     std::vector<types::global_dof_index> dofs_on_this_cell;
     dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
 
@@ -100,9 +107,16 @@ namespace DoFTools
           // make sparsity pattern for this cell. if no constraints pattern
           // was given, then the following call acts as if simply no
           // constraints existed
-          constraints.add_entries_local_to_global(dofs_on_this_cell,
-                                                  sparsity,
-                                                  keep_constrained_dofs);
+          const types::fe_index fe_index = cell->active_fe_index();
+          if (fe_dof_mask[fe_index].empty())
+            constraints.add_entries_local_to_global(dofs_on_this_cell,
+                                                    sparsity,
+                                                    keep_constrained_dofs);
+          else
+            constraints.add_entries_local_to_global(dofs_on_this_cell,
+                                                    sparsity,
+                                                    keep_constrained_dofs,
+                                                    fe_dof_mask[fe_index]);
         }
   }
 
@@ -152,6 +166,12 @@ namespace DoFTools
     const std::vector<Table<2, Coupling>> dof_mask //(fe_collection.size())
       = dof_couplings_from_component_couplings(fe_collection, couplings);
 
+    std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
+    for (unsigned int f = 0; f < fe_collection.size(); ++f)
+      {
+        fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
+      }
+
     // Convert the dof_mask to bool_dof_mask so we can pass it
     // to constraints.add_entries_local_to_global()
     std::vector<Table<2, bool>> bool_dof_mask(fe_collection.size());
@@ -163,7 +183,8 @@ namespace DoFTools
         bool_dof_mask[f].fill(false);
         for (unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
           for (unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
-            if (dof_mask[f](i, j) != none)
+            if (dof_mask[f](i, j) != none &&
+                (fe_dof_mask[f].empty() || fe_dof_mask[f](i, j)))
               bool_dof_mask[f](i, j) = true;
       }
 

diff --git a/source/fe/fe_q_iso_q1.cc b/source/fe/fe_q_iso_q1.cc
index 12ce528b305737830196609f9cb5a8aa526cf759..3f3ec7e2628831cbfde0e98be0499f92b1fd766b 100644 (file)
--- a/source/fe/fe_q_iso_q1.cc
+++ b/source/fe/fe_q_iso_q1.cc
@@ -14,10 +14,12 @@
 
 
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
 
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/fe_tools.h>
 
 #include <deal.II/lac/vector.h>
 
@@ -50,6 +52,42 @@ FE_Q_iso_Q1<dim, spacedim>::FE_Q_iso_Q1(const unsigned int subdivisions)
   const QIterated<1>  points(trapez, subdivisions);
 
   this->initialize(points.get_points());
+
+  {
+    const std::vector<unsigned int> R =
+      FETools::lexicographic_to_hierarchic_numbering<dim>(subdivisions);
+
+    this->local_dof_sparsity_pattern.reinit(this->n_dofs_per_cell(),
+                                            this->n_dofs_per_cell());
+
+    {
+      this->local_dof_sparsity_pattern.fill(false);
+
+      const unsigned int N   = Utilities::pow(points.size(), dim);
+      const int          N1d = points.size();
+      for (unsigned int i = 0; i < N; ++i)
+        for (unsigned int j = 0; j < N; ++j)
+          {
+            // compute l1 distance:
+            int distance = 0;
+
+            int xi = i;
+            int xj = j;
+            for (unsigned int d = 0; d < dim; ++d)
+              {
+                int current_distance = std::abs((xi % N1d) - (xj % N1d));
+                xi /= N1d;
+                xj /= N1d;
+                distance = std::max(distance, current_distance);
+                if (distance > 1)
+                  break;
+              }
+
+            if (distance <= 1)
+              this->local_dof_sparsity_pattern(R[i], R[j]) = true;
+          }
+    }
+  }
 }
 
 
@@ -178,6 +216,7 @@ FE_Q_iso_Q1<dim, spacedim>::compare_for_domination(
 }
 
 
+
 // explicit instantiations
 #include "fe_q_iso_q1.inst"
 

diff --git a/tests/fe/assemble_q_iso_q1.cc b/tests/fe/assemble_q_iso_q1.cc
new file mode 100644 (file)
index 0000000..795c25a
--- /dev/null
+++ b/tests/fe/assemble_q_iso_q1.cc
@@ -0,0 +1,255 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2013 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+// Assemble a 1d,2d,3d Poisson problem with FE_Q_iso_Q1 elements to
+// check that the sparsity pattern is computed correctly.  The problem
+// is taken from step-4
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+class Step4
+{
+public:
+  Step4();
+  void
+  run();
+
+private:
+  void
+  make_grid();
+  void
+  setup_system();
+  void
+  assemble_preconditioner();
+
+  Triangulation<dim> triangulation;
+
+  FE_Q_iso_Q1<dim> fe_precondition;
+  DoFHandler<dim>  dof_handler_precondition;
+
+  AffineConstraints<double> constraints;
+
+  SparsityPattern      prec_pattern;
+  SparseMatrix<double> preconditioner_matrix;
+
+  Vector<double> system_rhs;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+double
+RightHandSide<dim>::value(const Point<dim> &p,
+                          const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i = 0; i < dim; ++i)
+    return_value += 4 * std::pow(p[i], 4);
+
+  return return_value;
+}
+
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value(const Point<dim> &p,
+                           const unsigned int /*component*/) const
+{
+  return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4()
+  : fe_precondition(3)
+  , dof_handler_precondition(triangulation)
+{}
+
+
+template <int dim>
+void
+Step4<dim>::make_grid()
+{
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+  triangulation.refine_global(1);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::setup_system()
+{
+  dof_handler_precondition.distribute_dofs(fe_precondition);
+
+  constraints.clear();
+  std::map<unsigned int, double> boundary_values;
+  VectorTools::interpolate_boundary_values(dof_handler_precondition,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           constraints);
+  constraints.close();
+
+  {
+    DynamicSparsityPattern c_sparsity(dof_handler_precondition.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler_precondition,
+                                    c_sparsity,
+                                    constraints,
+                                    false);
+    prec_pattern.copy_from(c_sparsity);
+    preconditioner_matrix.reinit(prec_pattern);
+  }
+
+  system_rhs.reinit(dof_handler_precondition.n_dofs());
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::assemble_preconditioner()
+{
+  QIterated<dim> quadrature_formula(QGauss<1>(2), 3);
+
+  const RightHandSide<dim> right_hand_side;
+
+  FEValues<dim> fe_values(fe_precondition,
+                          quadrature_formula,
+                          update_values | update_gradients |
+                            update_quadrature_points | update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe_precondition.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs(dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler_precondition.begin_active(),
+    endc = dof_handler_precondition.end();
+
+  for (; cell != endc; ++cell)
+    {
+      fe_values.reinit(cell);
+      cell_matrix = 0;
+      cell_rhs    = 0;
+
+      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              cell_matrix(i, j) +=
+                (fe_values.shape_grad(i, q_point) *
+                 fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
+          }
+
+      cell->get_dof_indices(local_dof_indices);
+      constraints.distribute_local_to_global(cell_matrix,
+                                             local_dof_indices,
+                                             preconditioner_matrix);
+    }
+  preconditioner_matrix.compress(VectorOperation::add);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::run()
+{
+  deallog.push("dim " + std::to_string(dim));
+
+  make_grid();
+
+  setup_system();
+  assemble_preconditioner();
+  deallog << "nnz: " << preconditioner_matrix.n_nonzero_elements() << std::endl;
+
+  deallog.pop();
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog(true);
+
+  {
+    Step4<1> test;
+    test.run();
+  }
+  {
+    Step4<2> test;
+    test.run();
+  }
+  {
+    Step4<3> test;
+    test.run();
+  }
+}

diff --git a/tests/fe/assemble_q_iso_q1.output b/tests/fe/assemble_q_iso_q1.output
new file mode 100644 (file)
index 0000000..895ab16
--- /dev/null
+++ b/tests/fe/assemble_q_iso_q1.output
@@ -0,0 +1,4 @@
+
+DEAL:dim 1::nnz: 17
+DEAL:dim 2::nnz: 193
+DEAL:dim 3::nnz: 2415

diff --git a/tests/fe/local_sparsity_q_iso_q1.cc b/tests/fe/local_sparsity_q_iso_q1.cc
new file mode 100644 (file)
index 0000000..5698e80
--- /dev/null
+++ b/tests/fe/local_sparsity_q_iso_q1.cc
@@ -0,0 +1,132 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2015 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+// Assemble the sparsity pattern for a Q1 finite element space and a an iso Q
+// element space with the same number of subdivisions, and check we get the same
+// sparsity patterns.
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(const unsigned int n_subdivisions = 1)
+{
+  deallog << "dimension: " << dim << ", n_subdivisions = " << n_subdivisions
+          << std::endl;
+  Triangulation<dim> triangulation;
+  FE_Q_iso_Q1<dim>   fe(n_subdivisions);
+  DoFHandler<dim>    dof_handler(triangulation);
+  SparsityPattern    sparsity_pattern;
+
+  GridGenerator::hyper_cube(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  Triangulation<dim> triangulation_q1;
+  DoFHandler<dim>    dof_handler_q1(triangulation_q1);
+  SparsityPattern    sparsity_pattern_q1;
+  FE_Q<dim>          fe_q1(1);
+  GridGenerator::subdivided_hyper_cube(triangulation_q1, n_subdivisions);
+  dof_handler_q1.distribute_dofs(fe_q1);
+
+  std::vector<Point<dim>> support_points(dof_handler.n_dofs());
+  DoFTools::map_dofs_to_support_points(StaticMappingQ1<dim>::mapping,
+                                       dof_handler,
+                                       support_points);
+
+
+  std::vector<Point<dim>> support_points_q1(dof_handler_q1.n_dofs());
+  DoFTools::map_dofs_to_support_points(StaticMappingQ1<dim>::mapping,
+                                       dof_handler_q1,
+                                       support_points_q1);
+
+  AssertDimension(support_points.size(), support_points_q1.size());
+
+  std::vector<types::global_dof_index> new_numbers(dof_handler_q1.n_dofs());
+
+  for (unsigned int i = 0; i < support_points_q1.size(); ++i)
+    {
+      const Point<dim> &p = support_points_q1[i];
+      for (unsigned int j = 0; j < support_points.size(); ++j)
+        {
+          if (p.distance(support_points[j]) < 1e-10)
+            {
+              new_numbers[i] = j;
+              break;
+            }
+        }
+    }
+  dof_handler_q1.renumber_dofs(new_numbers);
+
+  {
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler, dsp);
+    dsp.print(deallog.get_file_stream());
+    sparsity_pattern.copy_from(dsp);
+    sparsity_pattern.print(deallog.get_file_stream());
+    deallog << std::endl;
+  }
+
+  {
+    DynamicSparsityPattern dsp(dof_handler_q1.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler_q1, dsp);
+    sparsity_pattern_q1.copy_from(dsp);
+  }
+
+  // Check that the two sparsity patterns are the same
+  auto el = sparsity_pattern.begin();
+  for (const auto &el_q1 : sparsity_pattern_q1)
+    {
+      if ((el->row() != el_q1.row()) || (el->column() != el_q1.column()))
+        {
+          deallog << "el_q1    : " << el_q1.row() << ", " << el_q1.column()
+                  << std::endl;
+          deallog << "el_q_iso : " << el->row() << ", " << el->column()
+                  << std::endl;
+        }
+      ++el;
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  test<1>(1);
+  test<1>(2);
+  test<1>(3);
+
+  test<2>(1);
+  test<2>(2);
+  test<2>(3);
+
+  test<3>(1);
+  test<3>(2);
+  test<3>(3);
+}

diff --git a/tests/fe/local_sparsity_q_iso_q1.output b/tests/fe/local_sparsity_q_iso_q1.output
new file mode 100644 (file)
index 0000000..6284eb2
--- /dev/null
+++ b/tests/fe/local_sparsity_q_iso_q1.output
@@ -0,0 +1,293 @@
+
+DEAL::dimension: 1, n_subdivisions = 1
+[0,0,1]
+[1,0,1]
+[0,0,1]
+[1,1,0]
+DEAL::
+DEAL::dimension: 1, n_subdivisions = 2
+[0,0,2]
+[1,1,2]
+[2,0,1,2]
+[0,0,2]
+[1,1,2]
+[2,2,0,1]
+DEAL::
+DEAL::dimension: 1, n_subdivisions = 3
+[0,0,2]
+[1,1,3]
+[2,0,2,3]
+[3,1,2,3]
+[0,0,2]
+[1,1,3]
+[2,2,0,3]
+[3,3,1,2]
+DEAL::
+DEAL::dimension: 2, n_subdivisions = 1
+[0,0,1,2,3]
+[1,0,1,2,3]
+[2,0,1,2,3]
+[3,0,1,2,3]
+[0,0,1,2,3]
+[1,1,0,2,3]
+[2,2,0,1,3]
+[3,3,0,1,2]
+DEAL::
+DEAL::dimension: 2, n_subdivisions = 2
+[0,0,4,6,8]
+[1,1,5,6,8]
+[2,2,4,7,8]
+[3,3,5,7,8]
+[4,0,2,4,6,7,8]
+[5,1,3,5,6,7,8]
+[6,0,1,4,5,6,8]
+[7,2,3,4,5,7,8]
+[8,0,1,2,3,4,5,6,7,8]
+[0,0,4,6,8]
+[1,1,5,6,8]
+[2,2,4,7,8]
+[3,3,5,7,8]
+[4,4,0,2,6,7,8]
+[5,5,1,3,6,7,8]
+[6,6,0,1,4,5,8]
+[7,7,2,3,4,5,8]
+[8,8,0,1,2,3,4,5,6,7]
+DEAL::
+DEAL::dimension: 2, n_subdivisions = 3
+[0,0,4,8,12]
+[1,1,6,9,13]
+[2,2,5,10,14]
+[3,3,7,11,15]
+[4,0,4,5,8,12,14]
+[5,2,4,5,10,12,14]
+[6,1,6,7,9,13,15]
+[7,3,6,7,11,13,15]
+[8,0,4,8,9,12,13]
+[9,1,6,8,9,12,13]
+[10,2,5,10,11,14,15]
+[11,3,7,10,11,14,15]
+[12,0,4,5,8,9,12,13,14,15]
+[13,1,6,7,8,9,12,13,14,15]
+[14,2,4,5,10,11,12,13,14,15]
+[15,3,6,7,10,11,12,13,14,15]
+[0,0,4,8,12]
+[1,1,6,9,13]
+[2,2,5,10,14]
+[3,3,7,11,15]
+[4,4,0,5,8,12,14]
+[5,5,2,4,10,12,14]
+[6,6,1,7,9,13,15]
+[7,7,3,6,11,13,15]
+[8,8,0,4,9,12,13]
+[9,9,1,6,8,12,13]
+[10,10,2,5,11,14,15]
+[11,11,3,7,10,14,15]
+[12,12,0,4,5,8,9,13,14,15]
+[13,13,1,6,7,8,9,12,14,15]
+[14,14,2,4,5,10,11,12,13,15]
+[15,15,3,6,7,10,11,12,13,14]
+DEAL::
+DEAL::dimension: 3, n_subdivisions = 1
+[0,0,1,2,3,4,5,6,7]
+[1,0,1,2,3,4,5,6,7]
+[2,0,1,2,3,4,5,6,7]
+[3,0,1,2,3,4,5,6,7]
+[4,0,1,2,3,4,5,6,7]
+[5,0,1,2,3,4,5,6,7]
+[6,0,1,2,3,4,5,6,7]
+[7,0,1,2,3,4,5,6,7]
+[0,0,1,2,3,4,5,6,7]
+[1,1,0,2,3,4,5,6,7]
+[2,2,0,1,3,4,5,6,7]
+[3,3,0,1,2,4,5,6,7]
+[4,4,0,1,2,3,5,6,7]
+[5,5,0,1,2,3,4,6,7]
+[6,6,0,1,2,3,4,5,7]
+[7,7,0,1,2,3,4,5,6]
+DEAL::
+DEAL::dimension: 3, n_subdivisions = 2
+[0,0,8,10,16,20,22,24,26]
+[1,1,9,10,17,21,22,24,26]
+[2,2,8,11,18,20,23,24,26]
+[3,3,9,11,19,21,23,24,26]
+[4,4,12,14,16,20,22,25,26]
+[5,5,13,14,17,21,22,25,26]
+[6,6,12,15,18,20,23,25,26]
+[7,7,13,15,19,21,23,25,26]
+[8,0,2,8,10,11,16,18,20,22,23,24,26]
+[9,1,3,9,10,11,17,19,21,22,23,24,26]
+[10,0,1,8,9,10,16,17,20,21,22,24,26]
+[11,2,3,8,9,11,18,19,20,21,23,24,26]
+[12,4,6,12,14,15,16,18,20,22,23,25,26]
+[13,5,7,13,14,15,17,19,21,22,23,25,26]
+[14,4,5,12,13,14,16,17,20,21,22,25,26]
+[15,6,7,12,13,15,18,19,20,21,23,25,26]
+[16,0,4,8,10,12,14,16,20,22,24,25,26]
+[17,1,5,9,10,13,14,17,21,22,24,25,26]
+[18,2,6,8,11,12,15,18,20,23,24,25,26]
+[19,3,7,9,11,13,15,19,21,23,24,25,26]
+[20,0,2,4,6,8,10,11,12,14,15,16,18,20,22,23,24,25,26]
+[21,1,3,5,7,9,10,11,13,14,15,17,19,21,22,23,24,25,26]
+[22,0,1,4,5,8,9,10,12,13,14,16,17,20,21,22,24,25,26]
+[23,2,3,6,7,8,9,11,12,13,15,18,19,20,21,23,24,25,26]
+[24,0,1,2,3,8,9,10,11,16,17,18,19,20,21,22,23,24,26]
+[25,4,5,6,7,12,13,14,15,16,17,18,19,20,21,22,23,25,26]
+[26,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]
+[0,0,8,10,16,20,22,24,26]
+[1,1,9,10,17,21,22,24,26]
+[2,2,8,11,18,20,23,24,26]
+[3,3,9,11,19,21,23,24,26]
+[4,4,12,14,16,20,22,25,26]
+[5,5,13,14,17,21,22,25,26]
+[6,6,12,15,18,20,23,25,26]
+[7,7,13,15,19,21,23,25,26]
+[8,8,0,2,10,11,16,18,20,22,23,24,26]
+[9,9,1,3,10,11,17,19,21,22,23,24,26]
+[10,10,0,1,8,9,16,17,20,21,22,24,26]
+[11,11,2,3,8,9,18,19,20,21,23,24,26]
+[12,12,4,6,14,15,16,18,20,22,23,25,26]
+[13,13,5,7,14,15,17,19,21,22,23,25,26]
+[14,14,4,5,12,13,16,17,20,21,22,25,26]
+[15,15,6,7,12,13,18,19,20,21,23,25,26]
+[16,16,0,4,8,10,12,14,20,22,24,25,26]
+[17,17,1,5,9,10,13,14,21,22,24,25,26]
+[18,18,2,6,8,11,12,15,20,23,24,25,26]
+[19,19,3,7,9,11,13,15,21,23,24,25,26]
+[20,20,0,2,4,6,8,10,11,12,14,15,16,18,22,23,24,25,26]
+[21,21,1,3,5,7,9,10,11,13,14,15,17,19,22,23,24,25,26]
+[22,22,0,1,4,5,8,9,10,12,13,14,16,17,20,21,24,25,26]
+[23,23,2,3,6,7,8,9,11,12,13,15,18,19,20,21,24,25,26]
+[24,24,0,1,2,3,8,9,10,11,16,17,18,19,20,21,22,23,26]
+[25,25,4,5,6,7,12,13,14,15,16,17,18,19,20,21,22,23,26]
+[26,26,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]
+DEAL::
+DEAL::dimension: 3, n_subdivisions = 3
+[0,0,8,12,24,32,40,48,56]
+[1,1,10,13,26,36,42,49,57]
+[2,2,9,14,28,33,44,50,58]
+[3,3,11,15,30,37,46,51,59]
+[4,4,16,20,25,34,41,52,60]
+[5,5,18,21,27,38,43,53,61]
+[6,6,17,22,29,35,45,54,62]
+[7,7,19,23,31,39,47,55,63]
+[8,0,8,9,12,24,32,33,40,48,50,56,58]
+[9,2,8,9,14,28,32,33,44,48,50,56,58]
+[10,1,10,11,13,26,36,37,42,49,51,57,59]
+[11,3,10,11,15,30,36,37,46,49,51,57,59]
+[12,0,8,12,13,24,32,40,42,48,49,56,57]
+[13,1,10,12,13,26,36,40,42,48,49,56,57]
+[14,2,9,14,15,28,33,44,46,50,51,58,59]
+[15,3,11,14,15,30,37,44,46,50,51,58,59]
+[16,4,16,17,20,25,34,35,41,52,54,60,62]
+[17,6,16,17,22,29,34,35,45,52,54,60,62]
+[18,5,18,19,21,27,38,39,43,53,55,61,63]
+[19,7,18,19,23,31,38,39,47,53,55,61,63]
+[20,4,16,20,21,25,34,41,43,52,53,60,61]
+[21,5,18,20,21,27,38,41,43,52,53,60,61]
+[22,6,17,22,23,29,35,45,47,54,55,62,63]
+[23,7,19,22,23,31,39,45,47,54,55,62,63]
+[24,0,8,12,24,25,32,34,40,41,48,56,60]
+[25,4,16,20,24,25,32,34,40,41,52,56,60]
+[26,1,10,13,26,27,36,38,42,43,49,57,61]
+[27,5,18,21,26,27,36,38,42,43,53,57,61]
+[28,2,9,14,28,29,33,35,44,45,50,58,62]
+[29,6,17,22,28,29,33,35,44,45,54,58,62]
+[30,3,11,15,30,31,37,39,46,47,51,59,63]
+[31,7,19,23,30,31,37,39,46,47,55,59,63]
+[32,0,8,9,12,24,25,32,33,34,35,40,41,48,50,56,58,60,62]
+[33,2,8,9,14,28,29,32,33,34,35,44,45,48,50,56,58,60,62]
+[34,4,16,17,20,24,25,32,33,34,35,40,41,52,54,56,58,60,62]
+[35,6,16,17,22,28,29,32,33,34,35,44,45,52,54,56,58,60,62]
+[36,1,10,11,13,26,27,36,37,38,39,42,43,49,51,57,59,61,63]
+[37,3,10,11,15,30,31,36,37,38,39,46,47,49,51,57,59,61,63]
+[38,5,18,19,21,26,27,36,37,38,39,42,43,53,55,57,59,61,63]
+[39,7,18,19,23,30,31,36,37,38,39,46,47,53,55,57,59,61,63]
+[40,0,8,12,13,24,25,32,34,40,41,42,43,48,49,56,57,60,61]
+[41,4,16,20,21,24,25,32,34,40,41,42,43,52,53,56,57,60,61]
+[42,1,10,12,13,26,27,36,38,40,41,42,43,48,49,56,57,60,61]
+[43,5,18,20,21,26,27,36,38,40,41,42,43,52,53,56,57,60,61]
+[44,2,9,14,15,28,29,33,35,44,45,46,47,50,51,58,59,62,63]
+[45,6,17,22,23,28,29,33,35,44,45,46,47,54,55,58,59,62,63]
+[46,3,11,14,15,30,31,37,39,44,45,46,47,50,51,58,59,62,63]
+[47,7,19,22,23,30,31,37,39,44,45,46,47,54,55,58,59,62,63]
+[48,0,8,9,12,13,24,32,33,40,42,48,49,50,51,56,57,58,59]
+[49,1,10,11,12,13,26,36,37,40,42,48,49,50,51,56,57,58,59]
+[50,2,8,9,14,15,28,32,33,44,46,48,49,50,51,56,57,58,59]
+[51,3,10,11,14,15,30,36,37,44,46,48,49,50,51,56,57,58,59]
+[52,4,16,17,20,21,25,34,35,41,43,52,53,54,55,60,61,62,63]
+[53,5,18,19,20,21,27,38,39,41,43,52,53,54,55,60,61,62,63]
+[54,6,16,17,22,23,29,34,35,45,47,52,53,54,55,60,61,62,63]
+[55,7,18,19,22,23,31,38,39,45,47,52,53,54,55,60,61,62,63]
+[56,0,8,9,12,13,24,25,32,33,34,35,40,41,42,43,48,49,50,51,56,57,58,59,60,61,62,63]
+[57,1,10,11,12,13,26,27,36,37,38,39,40,41,42,43,48,49,50,51,56,57,58,59,60,61,62,63]
+[58,2,8,9,14,15,28,29,32,33,34,35,44,45,46,47,48,49,50,51,56,57,58,59,60,61,62,63]
+[59,3,10,11,14,15,30,31,36,37,38,39,44,45,46,47,48,49,50,51,56,57,58,59,60,61,62,63]
+[60,4,16,17,20,21,24,25,32,33,34,35,40,41,42,43,52,53,54,55,56,57,58,59,60,61,62,63]
+[61,5,18,19,20,21,26,27,36,37,38,39,40,41,42,43,52,53,54,55,56,57,58,59,60,61,62,63]
+[62,6,16,17,22,23,28,29,32,33,34,35,44,45,46,47,52,53,54,55,56,57,58,59,60,61,62,63]
+[63,7,18,19,22,23,30,31,36,37,38,39,44,45,46,47,52,53,54,55,56,57,58,59,60,61,62,63]
+[0,0,8,12,24,32,40,48,56]
+[1,1,10,13,26,36,42,49,57]
+[2,2,9,14,28,33,44,50,58]
+[3,3,11,15,30,37,46,51,59]
+[4,4,16,20,25,34,41,52,60]
+[5,5,18,21,27,38,43,53,61]
+[6,6,17,22,29,35,45,54,62]
+[7,7,19,23,31,39,47,55,63]
+[8,8,0,9,12,24,32,33,40,48,50,56,58]
+[9,9,2,8,14,28,32,33,44,48,50,56,58]
+[10,10,1,11,13,26,36,37,42,49,51,57,59]
+[11,11,3,10,15,30,36,37,46,49,51,57,59]
+[12,12,0,8,13,24,32,40,42,48,49,56,57]
+[13,13,1,10,12,26,36,40,42,48,49,56,57]
+[14,14,2,9,15,28,33,44,46,50,51,58,59]
+[15,15,3,11,14,30,37,44,46,50,51,58,59]
+[16,16,4,17,20,25,34,35,41,52,54,60,62]
+[17,17,6,16,22,29,34,35,45,52,54,60,62]
+[18,18,5,19,21,27,38,39,43,53,55,61,63]
+[19,19,7,18,23,31,38,39,47,53,55,61,63]
+[20,20,4,16,21,25,34,41,43,52,53,60,61]
+[21,21,5,18,20,27,38,41,43,52,53,60,61]
+[22,22,6,17,23,29,35,45,47,54,55,62,63]
+[23,23,7,19,22,31,39,45,47,54,55,62,63]
+[24,24,0,8,12,25,32,34,40,41,48,56,60]
+[25,25,4,16,20,24,32,34,40,41,52,56,60]
+[26,26,1,10,13,27,36,38,42,43,49,57,61]
+[27,27,5,18,21,26,36,38,42,43,53,57,61]
+[28,28,2,9,14,29,33,35,44,45,50,58,62]
+[29,29,6,17,22,28,33,35,44,45,54,58,62]
+[30,30,3,11,15,31,37,39,46,47,51,59,63]
+[31,31,7,19,23,30,37,39,46,47,55,59,63]
+[32,32,0,8,9,12,24,25,33,34,35,40,41,48,50,56,58,60,62]
+[33,33,2,8,9,14,28,29,32,34,35,44,45,48,50,56,58,60,62]
+[34,34,4,16,17,20,24,25,32,33,35,40,41,52,54,56,58,60,62]
+[35,35,6,16,17,22,28,29,32,33,34,44,45,52,54,56,58,60,62]
+[36,36,1,10,11,13,26,27,37,38,39,42,43,49,51,57,59,61,63]
+[37,37,3,10,11,15,30,31,36,38,39,46,47,49,51,57,59,61,63]
+[38,38,5,18,19,21,26,27,36,37,39,42,43,53,55,57,59,61,63]
+[39,39,7,18,19,23,30,31,36,37,38,46,47,53,55,57,59,61,63]
+[40,40,0,8,12,13,24,25,32,34,41,42,43,48,49,56,57,60,61]
+[41,41,4,16,20,21,24,25,32,34,40,42,43,52,53,56,57,60,61]
+[42,42,1,10,12,13,26,27,36,38,40,41,43,48,49,56,57,60,61]
+[43,43,5,18,20,21,26,27,36,38,40,41,42,52,53,56,57,60,61]
+[44,44,2,9,14,15,28,29,33,35,45,46,47,50,51,58,59,62,63]
+[45,45,6,17,22,23,28,29,33,35,44,46,47,54,55,58,59,62,63]
+[46,46,3,11,14,15,30,31,37,39,44,45,47,50,51,58,59,62,63]
+[47,47,7,19,22,23,30,31,37,39,44,45,46,54,55,58,59,62,63]
+[48,48,0,8,9,12,13,24,32,33,40,42,49,50,51,56,57,58,59]
+[49,49,1,10,11,12,13,26,36,37,40,42,48,50,51,56,57,58,59]
+[50,50,2,8,9,14,15,28,32,33,44,46,48,49,51,56,57,58,59]
+[51,51,3,10,11,14,15,30,36,37,44,46,48,49,50,56,57,58,59]
+[52,52,4,16,17,20,21,25,34,35,41,43,53,54,55,60,61,62,63]
+[53,53,5,18,19,20,21,27,38,39,41,43,52,54,55,60,61,62,63]
+[54,54,6,16,17,22,23,29,34,35,45,47,52,53,55,60,61,62,63]
+[55,55,7,18,19,22,23,31,38,39,45,47,52,53,54,60,61,62,63]
+[56,56,0,8,9,12,13,24,25,32,33,34,35,40,41,42,43,48,49,50,51,57,58,59,60,61,62,63]
+[57,57,1,10,11,12,13,26,27,36,37,38,39,40,41,42,43,48,49,50,51,56,58,59,60,61,62,63]
+[58,58,2,8,9,14,15,28,29,32,33,34,35,44,45,46,47,48,49,50,51,56,57,59,60,61,62,63]
+[59,59,3,10,11,14,15,30,31,36,37,38,39,44,45,46,47,48,49,50,51,56,57,58,60,61,62,63]
+[60,60,4,16,17,20,21,24,25,32,33,34,35,40,41,42,43,52,53,54,55,56,57,58,59,61,62,63]
+[61,61,5,18,19,20,21,26,27,36,37,38,39,40,41,42,43,52,53,54,55,56,57,58,59,60,62,63]
+[62,62,6,16,17,22,23,28,29,32,33,34,35,44,45,46,47,52,53,54,55,56,57,58,59,60,61,63]
+[63,63,7,18,19,22,23,30,31,36,37,38,39,44,45,46,47,52,53,54,55,56,57,58,59,60,61,62]
+DEAL::