From 35f131b4720779a0979a69ce3070e2d08ad74154 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 3 Sep 2021 10:10:47 +0200 Subject: [PATCH] Fix MGTransferMatrixFree::interpolate_to_mg() for PBC --- .../deal.II/multigrid/mg_transfer_internal.h | 13 ++ .../multigrid/mg_transfer_matrix_free.h | 5 + source/multigrid/mg_transfer_internal.cc | 66 ++++----- tests/multigrid/interpolate_to_mg_01.cc | 134 ++++++++++++++++++ ...e_to_mg_01.with_p4est=true.mpirun=2.output | 11 ++ 5 files changed, 194 insertions(+), 35 deletions(-) create mode 100644 tests/multigrid/interpolate_to_mg_01.cc create mode 100644 tests/multigrid/interpolate_to_mg_01.with_p4est=true.mpirun=2.output diff --git a/include/deal.II/multigrid/mg_transfer_internal.h b/include/deal.II/multigrid/mg_transfer_internal.h index 7367914349..246abdfd16 100644 --- a/include/deal.II/multigrid/mg_transfer_internal.h +++ b/include/deal.II/multigrid/mg_transfer_internal.h @@ -137,6 +137,19 @@ namespace internal MGLevelObject> &vector_partitioners); + + + /** + * Helper function for setup_transfer. Checks for identity constrained + * dofs and replace with the indices of the dofs to which they are + * constrained + */ + void + resolve_identity_constraints( + const MGConstrainedDoFs * mg_constrained_dofs, + const unsigned int level, + std::vector &dof_indices); + } // namespace MGTransfer } // namespace internal diff --git a/include/deal.II/multigrid/mg_transfer_matrix_free.h b/include/deal.II/multigrid/mg_transfer_matrix_free.h index 5d6324be13..f87fff7d48 100644 --- a/include/deal.II/multigrid/mg_transfer_matrix_free.h +++ b/include/deal.II/multigrid/mg_transfer_matrix_free.h @@ -536,6 +536,7 @@ MGTransferMatrixFree::interpolate_to_mg( // do the transfer from level to level-1: dst[max_level].update_ghost_values(); + for (unsigned int level = max_level; level > min_level; --level) { // auxiliary vector which always has ghost elements @@ -569,6 +570,10 @@ MGTransferMatrixFree::interpolate_to_mg( for (unsigned int child = 0; child < cell->n_children(); ++child) { cell->child(child)->get_mg_dof_indices(dof_indices); + + internal::MGTransfer::resolve_identity_constraints( + this->mg_constrained_dofs, level, dof_indices); + for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) dof_values_fine(i) = (*input)(dof_indices[i]); fe.get_restriction_matrix(child, cell->refinement_case()) diff --git a/source/multigrid/mg_transfer_internal.cc b/source/multigrid/mg_transfer_internal.cc index 3bed763b93..d3636be4a1 100644 --- a/source/multigrid/mg_transfer_internal.cc +++ b/source/multigrid/mg_transfer_internal.cc @@ -641,39 +641,6 @@ namespace internal - namespace - { - /** - * Helper function for setup_transfer. Checks for identity constrained - * dofs and replace with the indices of the dofs to which they are - * constrained - */ - void - replace(const MGConstrainedDoFs * mg_constrained_dofs, - const unsigned int level, - std::vector &dof_indices) - { - if (mg_constrained_dofs != nullptr && - mg_constrained_dofs->get_level_constraints(level).n_constraints() > - 0) - for (auto &ind : dof_indices) - if (mg_constrained_dofs->get_level_constraints(level) - .is_identity_constrained(ind)) - { - Assert(mg_constrained_dofs->get_level_constraints(level) - .get_constraint_entries(ind) - ->size() == 1, - ExcInternalError()); - ind = mg_constrained_dofs->get_level_constraints(level) - .get_constraint_entries(ind) - ->front() - .first; - } - } - } // namespace - - - // Sets up most of the internal data structures of the MGTransferMatrixFree // class template @@ -809,7 +776,9 @@ namespace internal continue; cell->child(c)->get_mg_dof_indices(local_dof_indices); - replace(mg_constrained_dofs, level, local_dof_indices); + resolve_identity_constraints(mg_constrained_dofs, + level, + local_dof_indices); const IndexSet &owned_level_dofs = dof_handler.locally_owned_mg_dofs(level); @@ -880,7 +849,9 @@ namespace internal { cell->get_mg_dof_indices(local_dof_indices); - replace(mg_constrained_dofs, level - 1, local_dof_indices); + resolve_identity_constraints(mg_constrained_dofs, + level - 1, + local_dof_indices); const IndexSet &owned_level_dofs_l0 = dof_handler.locally_owned_mg_dofs(0); @@ -1030,6 +1001,31 @@ namespace internal } } + + + void + resolve_identity_constraints( + const MGConstrainedDoFs * mg_constrained_dofs, + const unsigned int level, + std::vector &dof_indices) + { + if (mg_constrained_dofs != nullptr && + mg_constrained_dofs->get_level_constraints(level).n_constraints() > 0) + for (auto &ind : dof_indices) + if (mg_constrained_dofs->get_level_constraints(level) + .is_identity_constrained(ind)) + { + Assert(mg_constrained_dofs->get_level_constraints(level) + .get_constraint_entries(ind) + ->size() == 1, + ExcInternalError()); + ind = mg_constrained_dofs->get_level_constraints(level) + .get_constraint_entries(ind) + ->front() + .first; + } + } + } // namespace MGTransfer } // namespace internal diff --git a/tests/multigrid/interpolate_to_mg_01.cc b/tests/multigrid/interpolate_to_mg_01.cc new file mode 100644 index 0000000000..9d13b93004 --- /dev/null +++ b/tests/multigrid/interpolate_to_mg_01.cc @@ -0,0 +1,134 @@ +// --------------------------------------------------------------------- +// +// 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 function MGTransferMatrixFree::interpolate_to_mg() for periodic +// boundaries + +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +template +class RightHandSideFunction : public Function +{ +public: + RightHandSideFunction(const unsigned int n_components) + : Function(n_components) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const + { + return p[component % dim] * p[component % dim]; + } +}; + +template +void +test() +{ + const unsigned int fe_degree = 1; + + parallel::distributed::Triangulation tria( + MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + + GridGenerator::hyper_cube(tria, -1., 1., true); + + std::vector>>> + tria_matched_pairs; + GridTools::collect_periodic_faces(tria, 0, 1, 0, tria_matched_pairs); + tria.add_periodicity(tria_matched_pairs); + tria.refine_global(2); + + const FE_Q fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + AffineConstraints constraints; + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + std::vector< + GridTools::PeriodicFacePair::cell_iterator>> + dof_matched_pairs; + GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, dof_matched_pairs); + DoFTools::make_periodicity_constraints(dof_matched_pairs, + constraints); + constraints.close(); + + LinearAlgebra::distributed::Vector qsol( + dof_handler.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD); + + VectorTools::interpolate(MappingQ1(), + dof_handler, + RightHandSideFunction(1), + qsol); + + MGLevelObject> mg_qsol; + MGConstrainedDoFs mg_constrained_dofs; + MGTransferMatrixFree mg_transfer; + + unsigned int n_tria_levels = tria.n_global_levels(); + mg_qsol.resize(0, n_tria_levels - 1); + + mg_constrained_dofs.initialize(dof_handler); + mg_transfer.initialize_constraints(mg_constrained_dofs); + mg_transfer.build(dof_handler); + + mg_transfer.interpolate_to_mg(dof_handler, mg_qsol, qsol); + + for (unsigned int i = mg_qsol.min_level(); i <= mg_qsol.max_level(); ++i) + deallog << mg_qsol[i].l2_norm() << std::endl; + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + test<2>(); +} diff --git a/tests/multigrid/interpolate_to_mg_01.with_p4est=true.mpirun=2.output b/tests/multigrid/interpolate_to_mg_01.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..5bc629fc1d --- /dev/null +++ b/tests/multigrid/interpolate_to_mg_01.with_p4est=true.mpirun=2.output @@ -0,0 +1,11 @@ + +DEAL:0::2.00000 +DEAL:0::2.44949 +DEAL:0::3.25960 +DEAL:0::OK + +DEAL:1::2.00000 +DEAL:1::2.44949 +DEAL:1::3.25960 +DEAL:1::OK + -- 2.39.5