From 4a409fc17164074a04a31fc23f2bf5c90367ca93 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Tue, 19 Nov 2024 09:12:47 -0500 Subject: [PATCH] local sparsity pattern for Q_iso_Q1 --- doc/news/changes/minor/20241120tjhei | 7 + include/deal.II/fe/fe.h | 36 +++ source/dofs/dof_tools_sparsity.cc | 29 ++- source/fe/fe_q_iso_q1.cc | 39 ++++ tests/fe/assemble_q_iso_q1.cc | 255 +++++++++++++++++++++ tests/fe/assemble_q_iso_q1.output | 4 + tests/fe/local_sparsity_q_iso_q1.cc | 132 +++++++++++ tests/fe/local_sparsity_q_iso_q1.output | 293 ++++++++++++++++++++++++ 8 files changed, 791 insertions(+), 4 deletions(-) create mode 100644 doc/news/changes/minor/20241120tjhei create mode 100644 tests/fe/assemble_q_iso_q1.cc create mode 100644 tests/fe/assemble_q_iso_q1.output create mode 100644 tests/fe/local_sparsity_q_iso_q1.cc create mode 100644 tests/fe/local_sparsity_q_iso_q1.output diff --git a/doc/news/changes/minor/20241120tjhei b/doc/news/changes/minor/20241120tjhei new file mode 100644 index 0000000000..457f7f8630 --- /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. +
+(Timo Heister, Luca Heltai, 2024/11/20) diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index 12f972276b..9ba2341046 100644 --- 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::is_primitive(const unsigned int i) const +template +inline const Table<2, bool> & +FiniteElement::get_local_dof_sparsity_pattern() const +{ + return local_dof_sparsity_pattern; +} + + + template inline GeometryPrimitive FiniteElement::get_associated_geometry_primitive( diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc index f7f0b74b65..bfce5678e1 100644 --- 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> 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 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> dof_mask //(fe_collection.size()) = dof_couplings_from_component_couplings(fe_collection, couplings); + std::vector> 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> 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 12ce528b30..3f3ec7e262 100644 --- a/source/fe/fe_q_iso_q1.cc +++ b/source/fe/fe_q_iso_q1.cc @@ -14,10 +14,12 @@ #include +#include #include #include #include +#include #include @@ -50,6 +52,42 @@ FE_Q_iso_Q1::FE_Q_iso_Q1(const unsigned int subdivisions) const QIterated<1> points(trapez, subdivisions); this->initialize(points.get_points()); + + { + const std::vector R = + FETools::lexicographic_to_hierarchic_numbering(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::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 index 0000000000..795c25a970 --- /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 + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include "../tests.h" + + + +template +class Step4 +{ +public: + Step4(); + void + run(); + +private: + void + make_grid(); + void + setup_system(); + void + assemble_preconditioner(); + + Triangulation triangulation; + + FE_Q_iso_Q1 fe_precondition; + DoFHandler dof_handler_precondition; + + AffineConstraints constraints; + + SparsityPattern prec_pattern; + SparseMatrix preconditioner_matrix; + + Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const; +}; + + + +template +double +RightHandSide::value(const Point &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 +double +BoundaryValues::value(const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4() + : fe_precondition(3) + , dof_handler_precondition(triangulation) +{} + + +template +void +Step4::make_grid() +{ + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global(1); +} + + + +template +void +Step4::setup_system() +{ + dof_handler_precondition.distribute_dofs(fe_precondition); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler_precondition, + 0, + BoundaryValues(), + 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 +void +Step4::assemble_preconditioner() +{ + QIterated quadrature_formula(QGauss<1>(2), 3); + + const RightHandSide right_hand_side; + + FEValues 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 cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + typename DoFHandler::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 +void +Step4::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 index 0000000000..895ab16f7b --- /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 index 0000000000..5698e80f5d --- /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 +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +template +void +test(const unsigned int n_subdivisions = 1) +{ + deallog << "dimension: " << dim << ", n_subdivisions = " << n_subdivisions + << std::endl; + Triangulation triangulation; + FE_Q_iso_Q1 fe(n_subdivisions); + DoFHandler dof_handler(triangulation); + SparsityPattern sparsity_pattern; + + GridGenerator::hyper_cube(triangulation); + dof_handler.distribute_dofs(fe); + + Triangulation triangulation_q1; + DoFHandler dof_handler_q1(triangulation_q1); + SparsityPattern sparsity_pattern_q1; + FE_Q fe_q1(1); + GridGenerator::subdivided_hyper_cube(triangulation_q1, n_subdivisions); + dof_handler_q1.distribute_dofs(fe_q1); + + std::vector> support_points(dof_handler.n_dofs()); + DoFTools::map_dofs_to_support_points(StaticMappingQ1::mapping, + dof_handler, + support_points); + + + std::vector> support_points_q1(dof_handler_q1.n_dofs()); + DoFTools::map_dofs_to_support_points(StaticMappingQ1::mapping, + dof_handler_q1, + support_points_q1); + + AssertDimension(support_points.size(), support_points_q1.size()); + + std::vector new_numbers(dof_handler_q1.n_dofs()); + + for (unsigned int i = 0; i < support_points_q1.size(); ++i) + { + const Point &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 index 0000000000..6284eb20e3 --- /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:: -- 2.39.5