From d3ea260bb98cca4a51a08b90aeaeed58e24e317a Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 9 Mar 2022 15:50:11 +0100 Subject: [PATCH] Introduce internal::MatrixFreeFunctions::ConstraintInfo Add dummy test --- include/deal.II/matrix_free/constraint_info.h | 530 ++++++++++++++++++ tests/matrix_free/constraint_info_01.cc | 223 ++++++++ tests/matrix_free/constraint_info_01.output | 26 + 3 files changed, 779 insertions(+) create mode 100644 include/deal.II/matrix_free/constraint_info.h create mode 100644 tests/matrix_free/constraint_info_01.cc create mode 100644 tests/matrix_free/constraint_info_01.output diff --git a/include/deal.II/matrix_free/constraint_info.h b/include/deal.II/matrix_free/constraint_info.h new file mode 100644 index 0000000000..dbb635a2f6 --- /dev/null +++ b/include/deal.II/matrix_free/constraint_info.h @@ -0,0 +1,530 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_constraint_info_h +#define dealii_matrix_free_constraint_info_h + + +#include + +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + namespace MatrixFreeFunctions + { + /** + * A helper class to apply constraints in matrix-free loops in + * user code. It combines constraint related functionalties from + * MatrixFree and FEEvaluation. + */ + template + class ConstraintInfo + { + public: + void + reinit(const DoFHandler &dof_handler, + const unsigned int n_cells, + const bool use_fast_hanging_node_algorithm = true); + + void + read_dof_indices( + const unsigned int cell_no, + const unsigned int mg_level, + const TriaIterator> &cell, + const dealii::AffineConstraints + & constraints, + const std::shared_ptr &partitioner); + + void + finalize(); + + template + void + read_write_operation(const T & operation, + VectorType & global_vector, + AlignedVector &local_vector, + const unsigned int first_cell, + const unsigned int n_cells, + const unsigned int n_dofs_per_cell, + const bool apply_constraints) const; + + void + apply_hanging_node_constraints( + const unsigned int first_cell, + const unsigned int n_lanes_filled, + const bool transpose, + AlignedVector &evaluation_data_coarse) const; + + private: + // for setup + ConstraintValues constraint_values; + std::vector> dof_indices_per_cell; + std::vector> plain_dof_indices_per_cell; + std::vector>> + constraint_indicator_per_cell; + + std::unique_ptr> hanging_nodes; + std::vector> lexicographic_numbering; + + public: + // for read_write_operation() + std::vector dof_indices; + std::vector> + constraint_indicator; + std::vector> row_starts; + + std::vector plain_dof_indices; + std::vector row_starts_plain_indices; + + // for constraint_pool_begin/end() + std::vector constraint_pool_data; + std::vector constraint_pool_row_index; + + std::vector> shape_infos; + std::vector hanging_node_constraint_masks; + std::vector active_fe_indices; + + private: + inline const typename Number::value_type * + constraint_pool_begin(const unsigned int row) const; + + inline const typename Number::value_type * + constraint_pool_end(const unsigned int row) const; + }; + + + + template + inline void + ConstraintInfo::reinit( + const DoFHandler &dof_handler, + const unsigned int n_cells, + const bool use_fast_hanging_node_algorithm) + { + this->dof_indices_per_cell.resize(n_cells); + this->plain_dof_indices_per_cell.resize(n_cells); + this->constraint_indicator_per_cell.resize(n_cells); + + if (use_fast_hanging_node_algorithm && + dof_handler.get_triangulation().has_hanging_nodes()) + { + hanging_nodes = std::make_unique>( + dof_handler.get_triangulation()); + + hanging_node_constraint_masks.resize(n_cells); + } + + const auto &fes = dof_handler.get_fe_collection(); + lexicographic_numbering.resize(fes.size()); + shape_infos.resize(fes.size()); + + for (unsigned int i = 0; i < fes.size(); ++i) + { + if (fes[i].reference_cell().is_hyper_cube()) + { + const Quadrature<1> dummy_quadrature( + std::vector>(1, Point<1>())); + shape_infos[i].reinit(dummy_quadrature, fes[i], 0); + } + else + { + const auto dummy_quadrature = + fes[i].reference_cell().template get_gauss_type_quadrature( + 1); + shape_infos[i].reinit(dummy_quadrature, fes[i], 0); + } + + lexicographic_numbering[i] = shape_infos[i].lexicographic_numbering; + } + active_fe_indices.resize(n_cells); + } + + + + template + inline void + ConstraintInfo::read_dof_indices( + const unsigned int cell_no, + const unsigned int mg_level, + const TriaIterator> & cell, + const dealii::AffineConstraints &constraints, + const std::shared_ptr & partitioner) + { + std::vector local_dof_indices( + cell->get_fe().n_dofs_per_cell()); + std::vector local_dof_indices_lex( + cell->get_fe().n_dofs_per_cell()); + + if (mg_level == numbers::invalid_unsigned_int) + cell->get_dof_indices(local_dof_indices); + else + cell->get_mg_dof_indices(local_dof_indices); + + { + AssertIndexRange(cell->active_fe_index(), shape_infos.size()); + + const auto &lexicographic_numbering = + shape_infos[cell->active_fe_index()].lexicographic_numbering; + + AssertDimension(lexicographic_numbering.size(), + local_dof_indices.size()); + + for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i) + local_dof_indices_lex[i] = + local_dof_indices[lexicographic_numbering[i]]; + } + + std::pair constraint_iterator(0, 0); + + AssertIndexRange(cell_no, this->constraint_indicator_per_cell.size()); + AssertIndexRange(cell_no, this->dof_indices_per_cell.size()); + AssertIndexRange(cell_no, this->plain_dof_indices_per_cell.size()); + AssertIndexRange(cell_no, this->plain_dof_indices_per_cell.size()); + + auto &constraint_indicator = this->constraint_indicator_per_cell[cell_no]; + auto &dof_indices = this->dof_indices_per_cell[cell_no]; + auto &plain_dof_indices = this->plain_dof_indices_per_cell[cell_no]; + + AssertDimension(constraint_indicator_per_cell[cell_no].size(), 0); + AssertDimension(dof_indices_per_cell[cell_no].size(), 0); + AssertDimension(plain_dof_indices_per_cell[cell_no].size(), 0); + + const auto global_to_local = + [&](const types::global_dof_index global_index) -> unsigned int { + if (partitioner) + return partitioner->global_to_local(global_index); + else + return global_index; + }; + + // plain indices + plain_dof_indices.resize(local_dof_indices_lex.size()); + for (unsigned int i = 0; i < local_dof_indices_lex.size(); ++i) + plain_dof_indices[i] = global_to_local(local_dof_indices_lex[i]); + + if (hanging_nodes) + { + AssertIndexRange(cell_no, this->hanging_node_constraint_masks.size()); + AssertIndexRange(cell_no, this->active_fe_indices.size()); + + std::vector mask(cell->get_fe().n_components()); + hanging_nodes->setup_constraints( + cell, {}, lexicographic_numbering, local_dof_indices_lex, mask); + + hanging_node_constraint_masks[cell_no] = mask[0]; + active_fe_indices[cell_no] = cell->active_fe_index(); + } + + for (auto current_dof : local_dof_indices_lex) + { + const auto *entries_ptr = + constraints.get_constraint_entries(current_dof); + + // dof is constrained + if (entries_ptr != nullptr) + { + const auto & entries = *entries_ptr; + const types::global_dof_index n_entries = entries.size(); + if (n_entries == 1 && + std::abs(entries[0].second - 1.) < + 100 * std::numeric_limits::epsilon()) + { + current_dof = entries[0].first; + goto no_constraint; + } + + constraint_indicator.push_back(constraint_iterator); + constraint_indicator.back().second = + constraint_values.insert_entries(entries); + + // reset constraint iterator for next round + constraint_iterator.first = 0; + + if (n_entries > 0) + { + const std::vector + &constraint_indices = constraint_values.constraint_indices; + for (unsigned int j = 0; j < n_entries; ++j) + { + dof_indices.push_back( + global_to_local(constraint_indices[j])); + } + } + } + else + { + no_constraint: + dof_indices.push_back(global_to_local(current_dof)); + + // make sure constraint_iterator.first is always within the + // bounds of unsigned short + Assert(constraint_iterator.first < + (1 << (8 * sizeof(unsigned short))) - 1, + ExcInternalError()); + constraint_iterator.first++; + } + } + } + + + + template + inline void + ConstraintInfo::finalize() + { + this->dof_indices = {}; + this->plain_dof_indices = {}; + this->constraint_indicator = {}; + + this->row_starts = {}; + this->row_starts.emplace_back(0, 0); + + this->row_starts_plain_indices = {}; + this->row_starts_plain_indices.emplace_back(0); + + for (unsigned int i = 0; i < dof_indices_per_cell.size(); ++i) + { + this->dof_indices.insert(this->dof_indices.end(), + dof_indices_per_cell[i].begin(), + dof_indices_per_cell[i].end()); + this->constraint_indicator.insert( + this->constraint_indicator.end(), + constraint_indicator_per_cell[i].begin(), + constraint_indicator_per_cell[i].end()); + + this->row_starts.emplace_back(this->dof_indices.size(), + this->constraint_indicator.size()); + + this->plain_dof_indices.insert(this->plain_dof_indices.end(), + plain_dof_indices_per_cell[i].begin(), + plain_dof_indices_per_cell[i].end()); + + this->row_starts_plain_indices.emplace_back( + this->plain_dof_indices.size()); + } + + std::vector *> constraints( + constraint_values.constraints.size()); + unsigned int length = 0; + for (const auto &it : constraint_values.constraints) + { + AssertIndexRange(it.second, constraints.size()); + constraints[it.second] = &it.first; + length += it.first.size(); + } + + constraint_pool_data.clear(); + constraint_pool_data.reserve(length); + constraint_pool_row_index.reserve(constraint_values.constraints.size() + + 1); + constraint_pool_row_index.resize(1, 0); + + for (const auto &constraint : constraints) + { + Assert(constraint != nullptr, ExcInternalError()); + constraint_pool_data.insert(constraint_pool_data.end(), + constraint->begin(), + constraint->end()); + constraint_pool_row_index.push_back(constraint_pool_data.size()); + } + + AssertDimension(constraint_pool_data.size(), length); + + dof_indices_per_cell.clear(); + constraint_indicator_per_cell.clear(); + + if (hanging_nodes && std::all_of(hanging_node_constraint_masks.begin(), + hanging_node_constraint_masks.end(), + [](const auto i) { + return i == + ConstraintKinds::unconstrained; + })) + hanging_node_constraint_masks.clear(); + } + + + + template + template + inline void + ConstraintInfo::read_write_operation( + const T & operation, + VectorType & global_vector, + AlignedVector &local_vector, + const unsigned int first_cell, + const unsigned int n_cells, + const unsigned int n_dofs_per_cell, + const bool apply_constraints) const + { + if (apply_constraints == false) + { + for (unsigned int v = 0; v < n_cells; ++v) + { + const unsigned int cell_index = first_cell + v; + const unsigned int *dof_indices = + this->plain_dof_indices.data() + + this->row_starts_plain_indices[cell_index]; + + for (unsigned int i = 0; i < n_dofs_per_cell; ++dof_indices, ++i) + operation.process_dof(*dof_indices, + global_vector, + local_vector[i][v]); + } + + return; + } + + for (unsigned int v = 0; v < n_cells; ++v) + { + const unsigned int cell_index = first_cell + v; + const unsigned int *dof_indices = + this->dof_indices.data() + this->row_starts[cell_index].first; + unsigned int index_indicators = this->row_starts[cell_index].second; + unsigned int next_index_indicators = + this->row_starts[cell_index + 1].second; + + unsigned int ind_local = 0; + for (; index_indicators != next_index_indicators; ++index_indicators) + { + const std::pair indicator = + this->constraint_indicator[index_indicators]; + + // run through values up to next constraint + for (unsigned int j = 0; j < indicator.first; ++j) + operation.process_dof(dof_indices[j], + global_vector, + local_vector[ind_local + j][v]); + + ind_local += indicator.first; + dof_indices += indicator.first; + + // constrained case: build the local value as a linear + // combination of the global value according to constraints + typename Number::value_type value; + operation.pre_constraints(local_vector[ind_local][v], value); + + const typename Number::value_type *data_val = + this->constraint_pool_begin(indicator.second); + const typename Number::value_type *end_pool = + this->constraint_pool_end(indicator.second); + for (; data_val != end_pool; ++data_val, ++dof_indices) + operation.process_constraint(*dof_indices, + *data_val, + global_vector, + value); + + operation.post_constraints(value, local_vector[ind_local][v]); + ind_local++; + } + + AssertIndexRange(ind_local, n_dofs_per_cell + 1); + + for (; ind_local < n_dofs_per_cell; ++dof_indices, ++ind_local) + operation.process_dof(*dof_indices, + global_vector, + local_vector[ind_local][v]); + } + } + + + + template + inline void + ConstraintInfo::apply_hanging_node_constraints( + const unsigned int first_cell, + const unsigned int n_lanes_filled, + const bool transpose, + AlignedVector &evaluation_data_coarse) const + { + if (hanging_node_constraint_masks.size() == 0) + return; + + std::array constraint_mask; + + bool hn_available = false; + + for (unsigned int v = 0; v < n_lanes_filled; ++v) + { + const auto mask = hanging_node_constraint_masks[first_cell + v]; + + constraint_mask[v] = mask; + + hn_available |= (mask != ConstraintKinds::unconstrained); + } + + if (hn_available == true) + { + for (unsigned int v = n_lanes_filled; v < Number::size(); ++v) + constraint_mask[v] = ConstraintKinds::unconstrained; + + for (unsigned int i = 1; i < n_lanes_filled; ++i) + AssertDimension(active_fe_indices[first_cell], + active_fe_indices[first_cell + i]); + + const auto &shape_info = shape_infos[active_fe_indices[first_cell]]; + + dealii::internal::FEEvaluationHangingNodesFactory< + dim, + typename Number::value_type, + Number>::apply(shape_info.n_components, + shape_info.data.front().fe_degree, + shape_info, + transpose, + constraint_mask, + evaluation_data_coarse.begin()); + } + } + + + + template + inline const typename Number::value_type * + ConstraintInfo::constraint_pool_begin( + const unsigned int row) const + { + AssertIndexRange(row, constraint_pool_row_index.size() - 1); + return constraint_pool_data.empty() ? + nullptr : + constraint_pool_data.data() + constraint_pool_row_index[row]; + } + + + + template + inline const typename Number::value_type * + ConstraintInfo::constraint_pool_end( + const unsigned int row) const + { + AssertIndexRange(row, constraint_pool_row_index.size() - 1); + return constraint_pool_data.empty() ? + nullptr : + constraint_pool_data.data() + constraint_pool_row_index[row + 1]; + } + + + } // namespace MatrixFreeFunctions +} // namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/matrix_free/constraint_info_01.cc b/tests/matrix_free/constraint_info_01.cc new file mode 100644 index 0000000000..b283c0db11 --- /dev/null +++ b/tests/matrix_free/constraint_info_01.cc @@ -0,0 +1,223 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 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 internal::MatrixFreeFunctions::ConstraintInfo and compare the +// result with FEEvaluation. + +#include + +#include + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +int +main() +{ + initlog(); + + const int dim = 2; + const int degree = 1; + + using Number = double; + using VectorizedArrayType = VectorizedArray; + using VectorType = LinearAlgebra::distributed::Vector; + + Triangulation tria; + GridGenerator::subdivided_hyper_cube(tria, 2); + tria.begin()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + QGauss quad(degree + 1); + FE_Q fe(degree); + MappingQ1 mapping; + + DoFHandler dof_handler; + dof_handler.reinit(tria); + dof_handler.distribute_dofs(fe); + + AffineConstraints constraints; + + DoFTools::make_zero_boundary_constraints(dof_handler, constraints); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + typename MatrixFree::AdditionalData + additional_data; + additional_data.mapping_update_flags = update_values | update_gradients | + update_JxW_values | + dealii::update_quadrature_points; + + MatrixFree matrix_free; + matrix_free.reinit(mapping, dof_handler, constraints, quad, additional_data); + + internal::MatrixFreeFunctions::ConstraintInfo + constraint_info; + constraint_info.reinit(dof_handler, matrix_free.n_physical_cells(), true); + + for (unsigned int cell = 0, cell_no = 0; cell < matrix_free.n_cell_batches(); + ++cell) + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v, ++cell_no) + constraint_info.read_dof_indices(cell_no, + -1, + matrix_free.get_cell_iterator(cell, v), + constraints, + matrix_free.get_vector_partitioner()); + + constraint_info.finalize(); + + const auto print_stat = [](const auto &dof_info) { + for (const auto i : dof_info.dof_indices) + deallog << i << " "; + deallog << std::endl; + + for (const auto i : dof_info.constraint_indicator) + deallog << "(" << i.first << "," << i.second << ") "; + deallog << std::endl; + + for (const auto i : dof_info.row_starts) + deallog << "(" << i.first << "," << i.second << ") "; + deallog << std::endl; + + deallog << std::endl; + }; + + print_stat(constraint_info); + print_stat(matrix_free.get_dof_info()); + + VectorType src, dst, dst_mf; + + const auto initialize_dof_vectors = [&]() { + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + matrix_free.initialize_dof_vector(dst_mf); + + for (const auto i : src.locally_owned_elements()) + src[i] = i; + }; + + initialize_dof_vectors(); + + AlignedVector local_vector(fe.n_dofs_per_cell()); + + const auto read_dof_values = [&](const unsigned int first_cell, + const unsigned int n_cells, + const unsigned int n_dofs_per_cell) { + internal::VectorReader reader; + constraint_info.read_write_operation( + reader, src, local_vector, first_cell, n_cells, n_dofs_per_cell, true); + constraint_info.apply_hanging_node_constraints(first_cell, + n_cells, + false, + local_vector); + }; + + const auto distribute_local_to_global = + [&](const unsigned int first_cell, + const unsigned int n_cells, + const unsigned int n_dofs_per_cell) { + internal::VectorDistributorLocalToGlobal + writer; + constraint_info.apply_hanging_node_constraints(first_cell, + n_cells, + true, + local_vector); + constraint_info.read_write_operation( + writer, dst, local_vector, first_cell, n_cells, n_dofs_per_cell, true); + }; + + const auto read_dof_values_plain = [&](const unsigned int first_cell, + const unsigned int n_cells, + const unsigned int n_dofs_per_cell) { + internal::VectorReader reader; + constraint_info.read_write_operation( + reader, src, local_vector, first_cell, n_cells, n_dofs_per_cell, false); + }; + + const auto set_dof_values_plain = [&](const unsigned int first_cell, + const unsigned int n_cells, + const unsigned int n_dofs_per_cell) { + internal::VectorSetter writer; + constraint_info.read_write_operation( + writer, dst, local_vector, first_cell, n_cells, n_dofs_per_cell, false); + }; + + unsigned int first_cell = 0; + const unsigned int n_dofs_per_cell = fe.n_dofs_per_cell(); + + FEEvaluation fe_eval(matrix_free); + + for (unsigned int cell = 0; cell < matrix_free.n_cell_batches(); ++cell) + { + { + const unsigned int n_cells = + matrix_free.n_active_entries_per_cell_batch(cell); + + read_dof_values(first_cell, n_cells, n_dofs_per_cell); + distribute_local_to_global(first_cell, n_cells, n_dofs_per_cell); + + first_cell += n_cells; + } + + { + fe_eval.reinit(cell); + + fe_eval.read_dof_values(src); + fe_eval.distribute_local_to_global(dst_mf); + } + } + + dst.print(deallog.get_file_stream()); + dst_mf.print(deallog.get_file_stream()); + deallog << std::endl; + + initialize_dof_vectors(); + first_cell = 0; + + for (unsigned int cell = 0; cell < matrix_free.n_cell_batches(); ++cell) + { + { + const unsigned int n_cells = + matrix_free.n_active_entries_per_cell_batch(cell); + + read_dof_values_plain(first_cell, n_cells, n_dofs_per_cell); + set_dof_values_plain(first_cell, n_cells, n_dofs_per_cell); + + first_cell += n_cells; + } + + { + fe_eval.reinit(cell); + + fe_eval.read_dof_values_plain(src); + fe_eval.set_dof_values_plain(dst_mf); + } + } + + dst.print(deallog.get_file_stream()); + dst_mf.print(deallog.get_file_stream()); +} diff --git a/tests/matrix_free/constraint_info_01.output b/tests/matrix_free/constraint_info_01.output new file mode 100644 index 0000000000..3658842fd7 --- /dev/null +++ b/tests/matrix_free/constraint_info_01.output @@ -0,0 +1,26 @@ + +DEAL::11 11 2 11 2 11 2 2 2 2 +DEAL::(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (1,0) (1,0) (0,0) (0,0) (0,0) (1,0) (0,0) (1,0) (0,0) (1,0) (0,0) (0,0) +DEAL::(0,0) (1,3) (3,5) (5,7) (7,9) (8,12) (9,15) (10,18) +DEAL:: +DEAL::11 11 2 11 2 11 2 2 2 2 +DEAL::(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (1,0) (1,0) (0,0) (0,0) (0,0) (1,0) (0,0) (1,0) (0,0) (1,0) (0,0) (0,0) +DEAL::(0,0) (1,3) (3,5) (5,7) (7,9) (8,12) (9,15) (10,18) (10,18) +DEAL:: +Process #0 +Local range: [0, 14), global size: 14 +Vector data: +0.000e+00 0.000e+00 1.000e+01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 4.400e+01 0.000e+00 0.000e+00 +Process #0 +Local range: [0, 14), global size: 14 +Vector data: +0.000e+00 0.000e+00 1.000e+01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 4.400e+01 0.000e+00 0.000e+00 +DEAL:: +Process #0 +Local range: [0, 14), global size: 14 +Vector data: +0.000e+00 1.000e+00 2.000e+00 3.000e+00 4.000e+00 5.000e+00 6.000e+00 7.000e+00 8.000e+00 9.000e+00 1.000e+01 1.100e+01 1.200e+01 1.300e+01 +Process #0 +Local range: [0, 14), global size: 14 +Vector data: +0.000e+00 1.000e+00 2.000e+00 3.000e+00 4.000e+00 5.000e+00 6.000e+00 7.000e+00 8.000e+00 9.000e+00 1.000e+01 1.100e+01 1.200e+01 1.300e+01 -- 2.39.5