From: Marc Fehling Date: Fri, 3 Jul 2020 14:48:05 +0000 (+0200) Subject: Bugfix: Unify dof indices on ghost interfaces if no finite element dominates. X-Git-Tag: v9.3.0-rc1~1331^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=904369c78ae5e461d9168d8f2f6435c54a27f5f8;p=dealii.git Bugfix: Unify dof indices on ghost interfaces if no finite element dominates. By design, there will be always one dominating element for vertices. --- diff --git a/doc/news/changes/minor/20200703Fehling b/doc/news/changes/minor/20200703Fehling new file mode 100644 index 0000000000..91d334b61a --- /dev/null +++ b/doc/news/changes/minor/20200703Fehling @@ -0,0 +1,6 @@ +Fixed: In parallel hp-adaptive applications, +DoFHandler::distribute_dofs() no longer fails to enumerate degrees of +freedom on ghost interfaces if continuous finite elements do not +dominate each other. +
+(Marc Fehling, 2020/07/03) diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index 00fb849ec2..a1c2a6c279 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -1102,98 +1102,99 @@ namespace internal // find out which is the most dominating finite // element of the ones that are used on this vertex - const unsigned int most_dominating_fe_index = + unsigned int most_dominating_fe_index = dof_handler.get_fe_collection().find_dominating_fe( fe_indices, /*codim=*/dim); - // if we found the most dominating element, then use - // this to eliminate some of the degrees of freedom - // by identification. otherwise, the code that - // computes hanging node constraints will have to - // deal with it by computing appropriate constraints - // along this face/edge - if (most_dominating_fe_index != + // if we haven't found a dominating finite element, + // choose the very first one to be dominant similar + // to compute_vertex_dof_identities() + if (most_dominating_fe_index == numbers::invalid_unsigned_int) - { - // loop over the indices of all the finite - // elements that are not dominating, and - // identify their dofs to the most dominating - // one - for (const auto &other_fe_index : fe_indices) - if (other_fe_index != most_dominating_fe_index) + most_dominating_fe_index = + dealii::internal::DoFAccessorImplementation:: + Implementation::nth_active_fe_index( + dof_handler, + 0, + vertex_index, + 0, + std::integral_constant()); + + // loop over the indices of all the finite + // elements that are not dominating, and + // identify their dofs to the most dominating + // one + for (const auto &other_fe_index : fe_indices) + if (other_fe_index != most_dominating_fe_index) + { + // make sure the entry in the equivalence + // table exists + ensure_existence_of_dof_identities<0>( + dof_handler.get_fe(most_dominating_fe_index), + dof_handler.get_fe(other_fe_index), + vertex_dof_identities[most_dominating_fe_index] + [other_fe_index]); + + // then loop through the identities we + // have. first get the global numbers of the + // dofs we want to identify and make sure they + // are not yet constrained to anything else, + // except for to each other. use the rule that + // we will always constrain the dof with the + // higher fe index to the one with the lower, + // to avoid circular reasoning. + DoFIdentities &identities = + *vertex_dof_identities[most_dominating_fe_index] + [other_fe_index]; + for (const auto &identity : identities) { - // make sure the entry in the equivalence - // table exists - ensure_existence_of_dof_identities<0>( - dof_handler.get_fe(most_dominating_fe_index), - dof_handler.get_fe(other_fe_index), - vertex_dof_identities[most_dominating_fe_index] - [other_fe_index]); + const types::global_dof_index master_dof_index = + dealii::internal::DoFAccessorImplementation:: + Implementation::get_dof_index( + dof_handler, + 0, + vertex_index, + most_dominating_fe_index, + identity.first, + std::integral_constant()); + const types::global_dof_index slave_dof_index = + dealii::internal::DoFAccessorImplementation:: + Implementation::get_dof_index( + dof_handler, + 0, + vertex_index, + other_fe_index, + identity.second, + std::integral_constant()); - // then loop through the identities we - // have. first get the global numbers of the - // dofs we want to identify and make sure they - // are not yet constrained to anything else, - // except for to each other. use the rule that - // we will always constrain the dof with the - // higher fe index to the one with the lower, - // to avoid circular reasoning. - DoFIdentities &identities = - *vertex_dof_identities[most_dominating_fe_index] - [other_fe_index]; - for (const auto &identity : identities) - { - const types::global_dof_index - master_dof_index = dealii::internal:: - DoFAccessorImplementation:: - Implementation::get_dof_index( - dof_handler, - 0, - vertex_index, - most_dominating_fe_index, - identity.first, - std::integral_constant()); - const types::global_dof_index - slave_dof_index = dealii::internal:: - DoFAccessorImplementation:: - Implementation::get_dof_index( - dof_handler, - 0, - vertex_index, - other_fe_index, - identity.second, - std::integral_constant()); - - // check if we are on an interface between - // a locally owned and a ghost cell on which - // we need to work on. - // - // all degrees of freedom belonging to - // dominating fe indices or to a processor - // with a higher rank have been set at this - // point (either in Phase 2, or after the - // first ghost exchange in Phase 5). thus, - // we only have to set the indices of - // degrees of freedom that have been - // previously flagged invalid. - if ((slave_dof_index == - numbers::invalid_dof_index) && - (master_dof_index != - numbers::invalid_dof_index)) - dealii::internal:: - DoFAccessorImplementation:: - Implementation::set_dof_index( - dof_handler, - 0, - vertex_index, - other_fe_index, - identity.second, - std::integral_constant(), - master_dof_index); - } + // check if we are on an interface between + // a locally owned and a ghost cell on which + // we need to work on. + // + // all degrees of freedom belonging to + // dominating fe indices or to a processor + // with a higher rank have been set at this + // point (either in Phase 2, or after the + // first ghost exchange in Phase 5). thus, + // we only have to set the indices of + // degrees of freedom that have been + // previously flagged invalid. + if ((slave_dof_index == + numbers::invalid_dof_index) && + (master_dof_index != + numbers::invalid_dof_index)) + dealii::internal::DoFAccessorImplementation:: + Implementation::set_dof_index( + dof_handler, + 0, + vertex_index, + other_fe_index, + identity.second, + std::integral_constant(), + master_dof_index); } - } + } } } } diff --git a/tests/mpi/hp_unify_dof_indices_11.cc b/tests/mpi/hp_unify_dof_indices_11.cc new file mode 100644 index 0000000000..86405dbe82 --- /dev/null +++ b/tests/mpi/hp_unify_dof_indices_11.cc @@ -0,0 +1,113 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +// --------------------------------------------------------------------- + + + +// have a 2x1 coarse mesh (or 2x1x1) and verify DoF indices in the hp +// case with a FECollection that contains two finite elements that do +// not dominate each other. Here, a (FE_Q(1) x FE_Q(2)) and a +// (FE_Q(2) x FE_Q(1)) element on two separate subdomains face each +// other. the hp code will unify DoF indices on boundaries between all +// subdomains. + + +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +template +void +test() +{ + parallel::distributed::Triangulation triangulation( + MPI_COMM_WORLD, Triangulation::limit_level_difference_at_vertices); + + std::vector reps(dim, 1U); + reps[0] = 2; + Point top_right; + for (unsigned int d = 0; d < dim; ++d) + top_right[d] = (d == 0 ? 2 : 1); + GridGenerator::subdivided_hyper_rectangle(triangulation, + reps, + Point(), + top_right); + Assert(triangulation.n_global_active_cells() == 2, ExcInternalError()); + Assert(triangulation.n_active_cells() == 2, ExcInternalError()); + + hp::FECollection fe; + fe.push_back(FESystem(FE_Q(1), 1, FE_Q(2), 1)); + fe.push_back(FESystem(FE_Q(2), 1, FE_Q(1), 1)); + + hp::DoFHandler dof_handler(triangulation); + for (auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + if (cell->id().to_string() == "0_0:") + cell->set_active_fe_index(0); + if (cell->id().to_string() == "1_0:") + cell->set_active_fe_index(1); + } + } + dof_handler.distribute_dofs(fe); + + deallog << "Processor: " << Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) + << std::endl; + for (auto &cell : dof_handler.active_cell_iterators()) + { + deallog << " Cell: " << cell; + if (cell->is_locally_owned()) + deallog << " is locally owned"; + else if (cell->is_ghost()) + deallog << " is ghost"; + deallog << std::endl; + + std::vector dof_indices( + cell->get_fe().dofs_per_cell); + cell->get_dof_indices(dof_indices); + deallog << " "; + for (auto i : dof_indices) + deallog << i << ' '; + deallog << std::endl; + } + deallog << " n_locally_owned_dofs: " << dof_handler.n_locally_owned_dofs() + << std::endl; + deallog << " n_global_dofs: " << dof_handler.n_dofs() << std::endl; +} + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll log; + + deallog.push("2d"); + test<2>(); + deallog.pop(); + + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/mpi/hp_unify_dof_indices_11.mpirun=1.with_p4est=true.output b/tests/mpi/hp_unify_dof_indices_11.mpirun=1.with_p4est=true.output new file mode 100644 index 0000000000..fbb78d57ac --- /dev/null +++ b/tests/mpi/hp_unify_dof_indices_11.mpirun=1.with_p4est=true.output @@ -0,0 +1,15 @@ + +DEAL:0:2d::Processor: 0 +DEAL:0:2d:: Cell: 0.0 is locally owned +DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 +DEAL:0:2d:: Cell: 0.1 is locally owned +DEAL:0:2d:: 2 3 13 14 6 7 15 16 17 18 19 20 21 +DEAL:0:2d:: n_locally_owned_dofs: 22 +DEAL:0:2d:: n_global_dofs: 22 +DEAL:0:3d::Processor: 0 +DEAL:0:3d:: Cell: 0.0 is locally owned +DEAL:0:3d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 +DEAL:0:3d:: Cell: 0.1 is locally owned +DEAL:0:3d:: 2 3 35 36 6 7 37 38 10 11 39 40 14 15 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 +DEAL:0:3d:: n_locally_owned_dofs: 62 +DEAL:0:3d:: n_global_dofs: 62 diff --git a/tests/mpi/hp_unify_dof_indices_11.mpirun=2.with_p4est=true.output b/tests/mpi/hp_unify_dof_indices_11.mpirun=2.with_p4est=true.output new file mode 100644 index 0000000000..92c03fcd02 --- /dev/null +++ b/tests/mpi/hp_unify_dof_indices_11.mpirun=2.with_p4est=true.output @@ -0,0 +1,31 @@ + +DEAL:0:2d::Processor: 0 +DEAL:0:2d:: Cell: 0.0 is locally owned +DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 +DEAL:0:2d:: Cell: 0.1 is ghost +DEAL:0:2d:: 2 3 13 14 6 7 15 16 17 18 19 20 21 +DEAL:0:2d:: n_locally_owned_dofs: 13 +DEAL:0:2d:: n_global_dofs: 22 +DEAL:0:3d::Processor: 0 +DEAL:0:3d:: Cell: 0.0 is locally owned +DEAL:0:3d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 +DEAL:0:3d:: Cell: 0.1 is ghost +DEAL:0:3d:: 2 3 35 36 6 7 37 38 10 11 39 40 14 15 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 +DEAL:0:3d:: n_locally_owned_dofs: 35 +DEAL:0:3d:: n_global_dofs: 62 + +DEAL:1:2d::Processor: 1 +DEAL:1:2d:: Cell: 0.0 is ghost +DEAL:1:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 +DEAL:1:2d:: Cell: 0.1 is locally owned +DEAL:1:2d:: 2 3 13 14 6 7 15 16 17 18 19 20 21 +DEAL:1:2d:: n_locally_owned_dofs: 9 +DEAL:1:2d:: n_global_dofs: 22 +DEAL:1:3d::Processor: 1 +DEAL:1:3d:: Cell: 0.0 is ghost +DEAL:1:3d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 +DEAL:1:3d:: Cell: 0.1 is locally owned +DEAL:1:3d:: 2 3 35 36 6 7 37 38 10 11 39 40 14 15 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 +DEAL:1:3d:: n_locally_owned_dofs: 27 +DEAL:1:3d:: n_global_dofs: 62 +