From: Peter Munch Date: Tue, 12 Jul 2022 21:02:38 +0000 (+0200) Subject: Introduce inverse of MF::get_cell_iterator() X-Git-Tag: v9.5.0-rc1~1023^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b619fd3910910253e31766027b2563c16e99ecce;p=dealii.git Introduce inverse of MF::get_cell_iterator() --- diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index d8e519ab6d..60334d1be2 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -1656,6 +1656,16 @@ public: get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const; + /** + * Get MatrixFree index associated to a deal.II @p cell. To get + * the actual cell batch index and lane, do the postprocessing + * `index / VectorizedArrayType::size()` and `index % + * VectorizedArrayType::size()`. + */ + unsigned int + get_matrix_free_cell_index( + const typename Triangulation::cell_iterator &cell) const; + /** * Return the cell iterator in deal.II speak to an interior/exterior cell of * a face in a pair of a face batch and lane index. The second element of @@ -2043,6 +2053,11 @@ private: */ std::vector> cell_level_index; + /** + * Conversion from deal.II index (active or level index) to MatrixFree index + * (inverse of cell_level_index). + */ + std::vector mf_cell_indices; /** * For discontinuous Galerkin, the cell_level_index includes cells that are diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 90417b5f06..2e42dcbc39 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -270,6 +270,19 @@ MatrixFree::get_cell_level_and_index( +template +unsigned int +MatrixFree::get_matrix_free_cell_index( + const typename Triangulation::cell_iterator &cell) const +{ + return mf_cell_indices[(this->get_mg_level() == + numbers::invalid_unsigned_int) ? + cell->active_cell_index() : + cell->index()]; +} + + + template std::pair::cell_iterator, unsigned int> MatrixFree::get_face_iterator( @@ -802,6 +815,29 @@ MatrixFree::internal_reinit( mapping_is_initialized = true; } + + // set up map: deal.II index -> MatrixFree index + { + const auto &tria = dof_handler[0]->get_triangulation(); + const auto mg_level = this->get_mg_level(); + + mf_cell_indices.resize((mg_level == numbers::invalid_unsigned_int) ? + tria.n_active_cells() : + tria.n_cells(mg_level), + numbers::invalid_unsigned_int); + + for (unsigned int cell = 0; + cell < n_cell_batches() + n_ghost_cell_batches(); + ++cell) + for (unsigned int v = 0; v < n_active_entries_per_cell_batch(cell); ++v) + { + const auto tria_cell = get_cell_iterator(cell, v); + mf_cell_indices[(mg_level == numbers::invalid_unsigned_int) ? + tria_cell->active_cell_index() : + tria_cell->index()] = + cell * VectorizedArrayType::size() + v; + } + } } diff --git a/tests/matrix_free/iterators.cc b/tests/matrix_free/iterators_01.cc similarity index 100% rename from tests/matrix_free/iterators.cc rename to tests/matrix_free/iterators_01.cc diff --git a/tests/matrix_free/iterators.output b/tests/matrix_free/iterators_01.output similarity index 100% rename from tests/matrix_free/iterators.output rename to tests/matrix_free/iterators_01.output diff --git a/tests/matrix_free/iterators_02.cc b/tests/matrix_free/iterators_02.cc new file mode 100644 index 0000000000..ead60661fc --- /dev/null +++ b/tests/matrix_free/iterators_02.cc @@ -0,0 +1,96 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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 MatrixFree::get_matrix_free_cell_index() and +// MatrixFree::get_cell_iterator(), whether they indeed return +// inverse information. + + +#include +#include + +#include + +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + const unsigned int fe_degree = 1; + + using Number = double; + using VectorizedArrayType = VectorizedArray; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(3); + + FE_Q fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingQ mapping(1); + + QGauss<1> quad(fe_degree + 1); + + AffineConstraints constraints; + + MatrixFree matrix_free; + matrix_free.reinit(mapping, dof_handler, constraints, quad); + + for (const auto &cell : tria.active_cell_iterators()) + { + const auto mf_index = matrix_free.get_matrix_free_cell_index(cell); + + AssertThrow(cell == matrix_free.get_cell_iterator( + mf_index / VectorizedArrayType::size(), + mf_index % VectorizedArrayType::size()), + ExcInternalError()); + } + + deallog << "OK!" << std::endl; + + for (unsigned int i = 0; i < matrix_free.n_cell_batches(); ++i) + for (unsigned int v = 0; v < matrix_free.n_active_entries_per_cell_batch(i); + ++v) + AssertDimension(matrix_free.get_matrix_free_cell_index( + matrix_free.get_cell_iterator(i, v)), + i * VectorizedArrayType::size() + v); + + deallog << "OK!" << std::endl; +} + + +int +main() +{ + initlog(); + + test<2>(); +} diff --git a/tests/matrix_free/iterators_02.output b/tests/matrix_free/iterators_02.output new file mode 100644 index 0000000000..043986ad9a --- /dev/null +++ b/tests/matrix_free/iterators_02.output @@ -0,0 +1,3 @@ + +DEAL::OK! +DEAL::OK!