From: Denis Davydov Date: Fri, 25 Jan 2019 21:18:17 +0000 (+0100) Subject: Add internal::MatrixFreeFunctions::DoFInfo::get_dof_indices_on_cell_batch() X-Git-Tag: v9.1.0-rc1~377^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fc4d48b4ca8a7bd66098c0ef8c589dbb4c4eb5b8;p=dealii.git Add internal::MatrixFreeFunctions::DoFInfo::get_dof_indices_on_cell_batch() --- diff --git a/doc/doxygen/images/dofinfo_get_dof_indices.png b/doc/doxygen/images/dofinfo_get_dof_indices.png new file mode 100644 index 0000000000..8366cacc63 Binary files /dev/null and b/doc/doxygen/images/dofinfo_get_dof_indices.png differ diff --git a/doc/news/changes/minor/20190130DenisDavydov b/doc/news/changes/minor/20190130DenisDavydov new file mode 100644 index 0000000000..4cf4ec962c --- /dev/null +++ b/doc/news/changes/minor/20190130DenisDavydov @@ -0,0 +1,4 @@ +New: Add internal::MatrixFreeFunctions::DoFInfo::get_dof_indices_on_cell_batch() to +return locally owned DoFs used by matrix-free framework on a given cell. +
+(Denis Davydov, 2019/01/30) diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index 6839cd1d35..8f80cf2334 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -102,6 +102,30 @@ namespace internal fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const; + /** + * Populate the vector @p locall_indices with locally owned degrees of freedom + * stored on the cell block @p cell. + * If @p with_constraints is `true`, then the returned vector will contain indices + * required to resolve constraints. + * + * The image below illustrates the output of this function for cell blocks + * zero and one with zero Dirichlet boundary conditions at the bottom of + * the domain. Note that due to the presence of constraints, the DoFs + * returned by this function for the case `with_constraints = true` are + * not a simple union + * of per cell DoFs on the cell block @p cell. + * + * @image html dofinfo_get_dof_indices.png + * + * @note The returned indices may contain duplicates. The unique set can be + * obtain using `std::sort()` followed by `std::unique()` and + * `std::vector::erase()`. + */ + void + get_dof_indices_on_cell_batch(std::vector &locall_indices, + const unsigned int cell, + const bool with_constraints = true) const; + /** * This internal method takes the local indices on a cell and fills them * into this class. It resolves the constraints and distributes the diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 153e6a8236..a2a7c6ae68 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -183,6 +183,73 @@ namespace internal + void + DoFInfo::get_dof_indices_on_cell_batch(std::vector &my_rows, + const unsigned int cell, + const bool apply_constraints) const + { + const unsigned int n_fe_components = start_components.back(); + const unsigned int fe_index = + dofs_per_cell.size() == 1 ? 0 : cell_active_fe_index[cell]; + const unsigned int dofs_this_cell = dofs_per_cell[fe_index]; + + const unsigned int n_vectorization = vectorization_length; + constexpr auto dof_access_index = dof_access_cell; + AssertIndexRange(cell, + n_vectorization_lanes_filled[dof_access_index].size()); + const unsigned int n_vectorization_actual = + n_vectorization_lanes_filled[dof_access_index][cell]; + + // we might have constraints, so the final number + // of indices is not known a priori. + // conservatively reserve the maximum without constraints + my_rows.reserve(n_vectorization * dofs_this_cell); + my_rows.resize(0); + unsigned int total_size = 0; + for (unsigned int v = 0; v < n_vectorization_actual; ++v) + { + const unsigned int ib = + (cell * n_vectorization + v) * n_fe_components; + const unsigned int ie = + (cell * n_vectorization + v + 1) * n_fe_components; + + // figure out constraints by comparing constraint_indicator row + // shift for this cell within the block as compared to the next + // one + const bool has_constraints = + row_starts[ib].second != row_starts[ib + n_fe_components].second; + + auto do_copy = [&](const unsigned int *begin, + const unsigned int *end) { + const unsigned int shift = total_size; + total_size += (end - begin); + my_rows.resize(total_size); + std::copy(begin, end, my_rows.begin() + shift); + }; + + if (!has_constraints || apply_constraints) + { + const unsigned int *begin = + dof_indices.data() + row_starts[ib].first; + const unsigned int *end = + dof_indices.data() + row_starts[ie].first; + do_copy(begin, end); + } + else + { + Assert(row_starts_plain_indices[cell * n_vectorization + v] != + numbers::invalid_unsigned_int, + ExcNotInitialized()); + const unsigned int *begin = + plain_dof_indices.data() + + row_starts_plain_indices[cell * n_vectorization + v]; + const unsigned int *end = begin + dofs_this_cell; + do_copy(begin, end); + } + } + } + + template void DoFInfo ::read_dof_indices( diff --git a/tests/matrix_free/dof_info_01.cc b/tests/matrix_free/dof_info_01.cc new file mode 100644 index 0000000000..462abcac78 --- /dev/null +++ b/tests/matrix_free/dof_info_01.cc @@ -0,0 +1,203 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// Tests DoFInfo::get_dof_indices() + +#include + +#include +#include + +#include + +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + + + +template +void +test(const bool adaptive_ref = true) +{ + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + const unsigned int this_mpi_core = + dealii::Utilities::MPI::this_mpi_process(mpi_communicator); + + parallel::distributed::Triangulation tria(mpi_communicator); + GridGenerator::hyper_cube(tria, 0, 1, true); + tria.refine_global(1); + if (adaptive_ref) + { + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + if (cell->center().norm() < 0.5) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + if (cell->center()[0] < 0.2) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + else + { + tria.refine_global(1); + } + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + IndexSet owned_set = dof.locally_owned_dofs(); + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs(dof, relevant_set); + + AffineConstraints constraints(relevant_set); + DoFTools::make_hanging_node_constraints(dof, constraints); + // constrain bottom part of the boundary (lower in y direction) + VectorTools::interpolate_boundary_values(dof, + 2, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + std::shared_ptr> mf_data( + new MatrixFree()); + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.tasks_block_size = 7; + mf_data->reinit(dof, constraints, quad, data); + } + + const unsigned int n_cells = mf_data->n_macro_cells(); + const auto & dof_info = mf_data->get_dof_info(); + constexpr unsigned int n_vectorization = + VectorizedArray::n_array_elements; + + std::vector my_rows; + my_rows.reserve(fe.dofs_per_cell * n_vectorization); + + constexpr auto dof_access_index = + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; + + bool checked_interleaved = false; + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + deallog << "Cell: " << cell << std::endl; + checked_interleaved |= + (dof_info.index_storage_variants[dof_access_index][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved); + auto get_and_log = [&](const bool apply_constraints) { + dof_info.get_dof_indices_on_cell_batch(my_rows, + cell, + apply_constraints); + // sort and make unique to make visual inspection easier: + std::sort(my_rows.begin(), my_rows.end()); + my_rows.erase(std::unique(my_rows.begin(), my_rows.end()), + my_rows.end()); + for (auto el : my_rows) + deallog << " " << el; + deallog << std::endl; + }; + get_and_log(true); + get_and_log(false); + } + + // if we don't have adaptive refinement, we should have + // some cells interleaved (at least that is the intention of this + // part of the test) + Assert(checked_interleaved || adaptive_ref, ExcInternalError()); + + // output in Gnuplot + if (dim == 2) + { + std::map> support_points; + MappingQ1 mapping; + DoFTools::map_dofs_to_support_points(mapping, dof, support_points); + + const std::string prefix = + std::is_same::value ? "float_" : "double_"; + const std::string href = (adaptive_ref ? "" : "global_"); + const std::string base_filename = + prefix + href + "grid" + dealii::Utilities::int_to_string(dim) + "_" + + dealii::Utilities::int_to_string(fe_degree) + "_p" + + dealii::Utilities::int_to_string(this_mpi_core); + + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" + << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "unset grid" << std::endl + << "unset border" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels tc rgb 'red' nopoint notitle, '-' with labels point pt 4 offset 1,1 notitle" + << std::endl; + GridOut().write_gnuplot(tria, f); + f << "e" << std::endl; + + // output cell blocks: + for (unsigned int cell = 0; cell < n_cells; ++cell) + for (unsigned int c = 0; c < mf_data->n_components_filled(cell); ++c) + { + const auto dof_cell = mf_data->get_cell_iterator(cell, c); + f << dof_cell->center() << " \"" << cell << "\"\n"; + } + + f << std::flush; + f << "e" << std::endl << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, support_points); + + f << "e" << std::endl; + } +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + MPILogInitAll mpi_init_log; + + deallog.push("2d"); + test<2, 1>(); + test<2, 1, float>(); + test<2, 1>(false); + deallog.pop(); +} diff --git a/tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..226a54f496 --- /dev/null +++ b/tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,34 @@ + +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 7 18 19 21 22 +DEAL:0:2d:: 6 7 16 17 18 19 20 21 22 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 2 7 9 10 21 22 23 24 +DEAL:0:2d:: 0 6 7 8 10 21 22 23 24 25 26 +DEAL:0:2d::Cell: 2 +DEAL:0:2d:: 2 3 7 9 10 11 12 24 +DEAL:0:2d:: 0 1 2 3 7 8 9 10 11 12 24 25 26 +DEAL:0:2d::Cell: 3 +DEAL:0:2d:: 2 3 4 5 9 11 12 14 15 +DEAL:0:2d:: 2 3 4 5 9 11 12 13 14 15 +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 2 7 9 10 18 19 21 22 23 24 +DEAL:0:2d:: 0 6 7 8 10 16 17 18 19 20 21 22 23 24 25 26 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 2 3 4 5 7 9 10 11 12 14 15 24 +DEAL:0:2d:: 0 1 2 3 4 5 7 8 9 10 11 12 13 14 15 24 25 26 +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 2 3 5 6 7 8 +DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 5 8 10 12 13 14 +DEAL:0:2d:: 4 5 8 9 10 11 12 13 14 +DEAL:0:2d::Cell: 2 +DEAL:0:2d:: 6 7 8 15 16 17 18 19 20 +DEAL:0:2d:: 6 7 8 15 16 17 18 19 20 +DEAL:0:2d::Cell: 3 +DEAL:0:2d:: 8 13 14 17 20 21 22 23 24 +DEAL:0:2d:: 8 13 14 17 20 21 22 23 24 diff --git a/tests/matrix_free/dof_info_02.cc b/tests/matrix_free/dof_info_02.cc new file mode 100644 index 0000000000..a980641a68 --- /dev/null +++ b/tests/matrix_free/dof_info_02.cc @@ -0,0 +1,205 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// Tests DoFInfo::get_dof_indices() for multi-component system + +#include + +#include +#include + +#include +#include + +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + + + +template +void +test(const bool adaptive_ref = true) +{ + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + const unsigned int this_mpi_core = + dealii::Utilities::MPI::this_mpi_process(mpi_communicator); + + parallel::distributed::Triangulation tria(mpi_communicator); + GridGenerator::hyper_cube(tria, 0, 1, true); + tria.refine_global(1); + if (adaptive_ref) + { + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + if (cell->center().norm() < 0.5) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + if (cell->center()[0] < 0.2) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + else + { + tria.refine_global(1); + } + + FE_Q feq(fe_degree); + FESystem fe(feq, dim); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + IndexSet owned_set = dof.locally_owned_dofs(); + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs(dof, relevant_set); + + AffineConstraints constraints(relevant_set); + DoFTools::make_hanging_node_constraints(dof, constraints); + // constrain bottom part of the boundary (lower in y direction) + VectorTools::interpolate_boundary_values(dof, + 2, + Functions::ZeroFunction(2), + constraints); + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + std::shared_ptr> mf_data( + new MatrixFree()); + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.tasks_block_size = 7; + mf_data->reinit(dof, constraints, quad, data); + } + + const unsigned int n_cells = mf_data->n_macro_cells(); + const auto & dof_info = mf_data->get_dof_info(); + constexpr unsigned int n_vectorization = + VectorizedArray::n_array_elements; + + std::vector my_rows; + my_rows.reserve(fe.dofs_per_cell * n_vectorization); + + constexpr auto dof_access_index = + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; + + bool checked_interleaved = false; + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + deallog << "Cell: " << cell << std::endl; + checked_interleaved |= + (dof_info.index_storage_variants[dof_access_index][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved); + auto get_and_log = [&](const bool apply_constraints) { + dof_info.get_dof_indices_on_cell_batch(my_rows, + cell, + apply_constraints); + // sort and make unique to make visual inspection easier: + std::sort(my_rows.begin(), my_rows.end()); + my_rows.erase(std::unique(my_rows.begin(), my_rows.end()), + my_rows.end()); + for (auto el : my_rows) + deallog << " " << el; + deallog << std::endl; + }; + get_and_log(true); + get_and_log(false); + } + + // if we don't have adaptive refinement, we should have + // some cells interleaved (at least that is the intention of this + // part of the test) + Assert(checked_interleaved || adaptive_ref, ExcInternalError()); + + // output in Gnuplot + if (dim == 2) + { + std::map> support_points; + MappingQ1 mapping; + DoFTools::map_dofs_to_support_points(mapping, dof, support_points); + + const std::string prefix = + std::is_same::value ? "float_" : "double_"; + const std::string href = (adaptive_ref ? "" : "global_"); + const std::string base_filename = + prefix + href + "grid" + dealii::Utilities::int_to_string(dim) + "_" + + dealii::Utilities::int_to_string(fe_degree) + "_p" + + dealii::Utilities::int_to_string(this_mpi_core); + + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" + << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "unset grid" << std::endl + << "unset border" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels tc rgb 'red' nopoint notitle, '-' with labels point pt 4 offset 1,1 notitle" + << std::endl; + GridOut().write_gnuplot(tria, f); + f << "e" << std::endl; + + // output cell blocks: + for (unsigned int cell = 0; cell < n_cells; ++cell) + for (unsigned int c = 0; c < mf_data->n_components_filled(cell); ++c) + { + const auto dof_cell = mf_data->get_cell_iterator(cell, c); + f << dof_cell->center() << " \"" << cell << "\"\n"; + } + + f << std::flush; + f << "e" << std::endl << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, support_points); + + f << "e" << std::endl; + } +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + MPILogInitAll mpi_init_log; + + deallog.push("2d"); + test<2, 1>(); + test<2, 1, float>(); + test<2, 1>(false); + deallog.pop(); +} diff --git a/tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..1da9dc63ea --- /dev/null +++ b/tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,34 @@ + +DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2] +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 14 15 36 37 38 39 42 43 44 45 +DEAL:0:2d:: 12 13 14 15 32 33 34 35 36 37 38 39 40 41 42 43 44 45 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 4 5 14 15 18 19 20 21 42 43 44 45 46 47 48 49 +DEAL:0:2d:: 0 1 12 13 14 15 16 17 20 21 42 43 44 45 46 47 48 49 50 51 52 53 +DEAL:0:2d::Cell: 2 +DEAL:0:2d:: 4 5 6 7 14 15 18 19 20 21 22 23 24 25 48 49 +DEAL:0:2d:: 0 1 2 3 4 5 6 7 14 15 16 17 18 19 20 21 22 23 24 25 48 49 50 51 52 53 +DEAL:0:2d::Cell: 3 +DEAL:0:2d:: 4 5 6 7 8 9 10 11 18 19 22 23 24 25 28 29 30 31 +DEAL:0:2d:: 4 5 6 7 8 9 10 11 18 19 22 23 24 25 26 27 28 29 30 31 +DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2] +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 4 5 14 15 18 19 20 21 36 37 38 39 42 43 44 45 46 47 48 49 +DEAL:0:2d:: 0 1 12 13 14 15 16 17 20 21 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 4 5 6 7 8 9 10 11 14 15 18 19 20 21 22 23 24 25 28 29 30 31 48 49 +DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 48 49 50 51 52 53 +DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2] +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 4 5 6 7 10 11 12 13 14 15 16 17 +DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 10 11 16 17 20 21 24 25 26 27 28 29 +DEAL:0:2d:: 8 9 10 11 16 17 18 19 20 21 22 23 24 25 26 27 28 29 +DEAL:0:2d::Cell: 2 +DEAL:0:2d:: 12 13 14 15 16 17 30 31 32 33 34 35 36 37 38 39 40 41 +DEAL:0:2d:: 12 13 14 15 16 17 30 31 32 33 34 35 36 37 38 39 40 41 +DEAL:0:2d::Cell: 3 +DEAL:0:2d:: 16 17 26 27 28 29 34 35 40 41 42 43 44 45 46 47 48 49 +DEAL:0:2d:: 16 17 26 27 28 29 34 35 40 41 42 43 44 45 46 47 48 49