--- /dev/null
+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.
+<br>
+(Denis Davydov, 2019/01/30)
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<unsigned int> &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
+ void
+ DoFInfo::get_dof_indices_on_cell_batch(std::vector<unsigned int> &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 <typename number>
void
DoFInfo ::read_dof_indices(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/distributed/tria.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/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, int fe_degree, typename number = double>
+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<dim> 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<dim> fe(fe_degree);
+ DoFHandler<dim> 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<double> 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<dim>(),
+ constraints);
+ constraints.close();
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+ std::shared_ptr<MatrixFree<dim, number>> mf_data(
+ new MatrixFree<dim, number>());
+ {
+ const QGauss<1> quad(fe_degree + 1);
+ typename MatrixFree<dim, number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, number>::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<number>::n_array_elements;
+
+ std::vector<unsigned int> 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<types::global_dof_index, Point<dim>> support_points;
+ MappingQ1<dim> mapping;
+ DoFTools::map_dofs_to_support_points(mapping, dof, support_points);
+
+ const std::string prefix =
+ std::is_same<float, number>::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();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/distributed/tria.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_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, int fe_degree, typename number = double>
+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<dim> 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<dim> feq(fe_degree);
+ FESystem<dim> fe(feq, dim);
+ DoFHandler<dim> 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<double> 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<dim>(2),
+ constraints);
+ constraints.close();
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+ std::shared_ptr<MatrixFree<dim, number>> mf_data(
+ new MatrixFree<dim, number>());
+ {
+ const QGauss<1> quad(fe_degree + 1);
+ typename MatrixFree<dim, number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, number>::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<number>::n_array_elements;
+
+ std::vector<unsigned int> 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<types::global_dof_index, Point<dim>> support_points;
+ MappingQ1<dim> mapping;
+ DoFTools::map_dofs_to_support_points(mapping, dof, support_points);
+
+ const std::string prefix =
+ std::is_same<float, number>::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();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff
