From 6914b1ee04a1db5815cd662ac94cd0ac4e60f2d3 Mon Sep 17 00:00:00 2001 From: Doug Shi-Dong Date: Fri, 11 Oct 2019 15:19:10 -0400 Subject: [PATCH] Fix DoFHandler::n_boundary_dofs() in issue #8897 Only count locally owned dofs in the case of a parallel::distributed::Triangulation. Use suggested syntax on cell iterator. Co-Authored-By: Daniel Arndt Add test case for n_boundary_dofs in parallel Fix n_boundary_dofs() for hp version. Also implement suggested change Start 3212: dofs/n_boundary_dofs_01.mpirun=1.debug 1/8 Test # 3212: dofs/n_boundary_dofs_01.mpirun=1.debug ..... Passed 0.50 sec Start 3213: dofs/n_boundary_dofs_01.mpirun=1.release 2/8 Test # 3213: dofs/n_boundary_dofs_01.mpirun=1.release ... Passed 0.49 sec Start 3214: dofs/n_boundary_dofs_01.mpirun=2.debug 3/8 Test # 3214: dofs/n_boundary_dofs_01.mpirun=2.debug ..... Passed 0.50 sec Start 3215: dofs/n_boundary_dofs_01.mpirun=2.release 4/8 Test # 3215: dofs/n_boundary_dofs_01.mpirun=2.release ... Passed 0.51 sec Start 3216: dofs/n_boundary_dofs_02.debug 5/8 Test # 3216: dofs/n_boundary_dofs_02.debug .............. Passed 0.51 sec Start 3217: dofs/n_boundary_dofs_02.release 6/8 Test # 3217: dofs/n_boundary_dofs_02.release ............ Passed 0.52 sec Start 3218: dofs/n_boundary_dofs_03.debug 7/8 Test # 3218: dofs/n_boundary_dofs_03.debug .............. Passed 0.52 sec Start 3219: dofs/n_boundary_dofs_03.release 8/8 Test # 3219: dofs/n_boundary_dofs_03.release ............ Passed 0.54 sec 100% tests passed, 0 tests failed out of 8 Add hp n_boundary_dofs test. Use unordered_set instead. Fix a || to a && from switching the if-continue statement. Start 3212: dofs/n_boundary_dofs_01.mpirun=1.debug 1/18 Test # 3212: dofs/n_boundary_dofs_01.mpirun=1.debug ..... Passed 51.27 sec Start 3213: dofs/n_boundary_dofs_01.mpirun=1.release 2/18 Test # 3213: dofs/n_boundary_dofs_01.mpirun=1.release ... Passed 17.03 sec Start 3214: dofs/n_boundary_dofs_01.mpirun=2.debug 3/18 Test # 3214: dofs/n_boundary_dofs_01.mpirun=2.debug ..... Passed 1.90 sec Start 3215: dofs/n_boundary_dofs_01.mpirun=2.release 4/18 Test # 3215: dofs/n_boundary_dofs_01.mpirun=2.release ... Passed 1.45 sec Start 3216: dofs/n_boundary_dofs_02.debug 5/18 Test # 3216: dofs/n_boundary_dofs_02.debug .............. Passed 50.17 sec Start 3217: dofs/n_boundary_dofs_02.release 6/18 Test # 3217: dofs/n_boundary_dofs_02.release ............ Passed 16.30 sec Start 3218: dofs/n_boundary_dofs_03.debug 7/18 Test # 3218: dofs/n_boundary_dofs_03.debug .............. Passed 50.67 sec Start 3219: dofs/n_boundary_dofs_03.release 8/18 Test # 3219: dofs/n_boundary_dofs_03.release ............ Passed 16.71 sec Start 5498: hp/n_boundary_dofs.debug 9/18 Test # 5498: hp/n_boundary_dofs.debug ................... Passed 44.68 sec Start 5499: hp/n_boundary_dofs.release 10/18 Test # 5499: hp/n_boundary_dofs.release ................. Passed 11.83 sec Start 5500: hp/n_boundary_dofs_01.mpirun=1.debug 11/18 Test # 5500: hp/n_boundary_dofs_01.mpirun=1.debug ....... Passed 34.48 sec Start 5501: hp/n_boundary_dofs_01.mpirun=1.release 12/18 Test # 5501: hp/n_boundary_dofs_01.mpirun=1.release ..... Passed 12.46 sec Start 5502: hp/n_boundary_dofs_01.mpirun=2.debug 13/18 Test # 5502: hp/n_boundary_dofs_01.mpirun=2.debug ....... Passed 1.75 sec Start 5503: hp/n_boundary_dofs_01.mpirun=2.release 14/18 Test # 5503: hp/n_boundary_dofs_01.mpirun=2.release ..... Passed 1.36 sec Start 5504: hp/n_boundary_dofs_02.debug 15/18 Test # 5504: hp/n_boundary_dofs_02.debug ................ Passed 34.33 sec Start 5505: hp/n_boundary_dofs_02.release 16/18 Test # 5505: hp/n_boundary_dofs_02.release .............. Passed 12.25 sec Start 5506: hp/n_boundary_dofs_03.debug 17/18 Test # 5506: hp/n_boundary_dofs_03.debug ................ Passed 41.88 sec Start 5507: hp/n_boundary_dofs_03.release 18/18 Test # 5507: hp/n_boundary_dofs_03.release .............. Passed 15.45 sec 100% tests passed, 0 tests failed out of 18 --- doc/news/changes/minor/20191011DougShiDong | 7 + include/deal.II/dofs/dof_handler.h | 10 +- include/deal.II/hp/dof_handler.h | 8 +- source/dofs/dof_handler.cc | 138 ++++++++++++------ source/hp/dof_handler.cc | 75 ++++++---- tests/dofs/n_boundary_dofs_01.cc | 65 +++++++++ ...ry_dofs_01.with_p4est=true.mpirun=1.output | 3 + ...ry_dofs_01.with_p4est=true.mpirun=2.output | 3 + tests/hp/n_boundary_dofs_01.cc | 67 +++++++++ ...ry_dofs_01.with_p4est=true.mpirun=1.output | 3 + ...ry_dofs_01.with_p4est=true.mpirun=2.output | 3 + 11 files changed, 302 insertions(+), 80 deletions(-) create mode 100644 doc/news/changes/minor/20191011DougShiDong create mode 100644 tests/dofs/n_boundary_dofs_01.cc create mode 100644 tests/dofs/n_boundary_dofs_01.with_p4est=true.mpirun=1.output create mode 100644 tests/dofs/n_boundary_dofs_01.with_p4est=true.mpirun=2.output create mode 100644 tests/hp/n_boundary_dofs_01.cc create mode 100644 tests/hp/n_boundary_dofs_01.with_p4est=true.mpirun=1.output create mode 100644 tests/hp/n_boundary_dofs_01.with_p4est=true.mpirun=2.output diff --git a/doc/news/changes/minor/20191011DougShiDong b/doc/news/changes/minor/20191011DougShiDong new file mode 100644 index 0000000000..d7c9956ffe --- /dev/null +++ b/doc/news/changes/minor/20191011DougShiDong @@ -0,0 +1,7 @@ +Fixed: Take in account p::d::Triangulation when using +DoFHandler::n_boundary_dofs() and its +hp version. It now returns locally owned dofs instead +of attempting a global count. +
+(Doug Shi-Dong, 2019/10/11) + diff --git a/include/deal.II/dofs/dof_handler.h b/include/deal.II/dofs/dof_handler.h index a29432a77d..fc0a5763d0 100644 --- a/include/deal.II/dofs/dof_handler.h +++ b/include/deal.II/dofs/dof_handler.h @@ -954,15 +954,17 @@ public: n_dofs(const unsigned int level) const; /** - * Return the number of degrees of freedom located on the boundary. + * Return the number of locally owned degrees of freedom located on the + * boundary. */ types::global_dof_index n_boundary_dofs() const; /** - * Return the number of degrees of freedom located on those parts of the - * boundary which have a boundary indicator listed in the given set. The - * reason that a @p map rather than a @p set is used is the same as + * Return the number of locally owned degrees of freedom located on those + * parts of the boundary which have a boundary indicator listed in the given + * set. + * The reason that a @p map rather than a @p set is used is the same as * described in the documentation of that variant of * DoFTools::make_boundary_sparsity_pattern() that takes a map. * diff --git a/include/deal.II/hp/dof_handler.h b/include/deal.II/hp/dof_handler.h index e4850d9fd1..badc4bdf1a 100644 --- a/include/deal.II/hp/dof_handler.h +++ b/include/deal.II/hp/dof_handler.h @@ -778,7 +778,8 @@ namespace hp n_dofs(const unsigned int level) const; /** - * Return the number of degrees of freedom located on the boundary. + * Return the number of locally owned degrees of freedom located on the + * boundary. */ types::global_dof_index n_boundary_dofs() const; @@ -800,8 +801,9 @@ namespace hp &boundary_ids) const; /** - * Return the number of degrees of freedom located on those parts of the - * boundary which have a boundary indicator listed in the given set. The + * Return the number of locally owned degrees of freedom located on those + * parts of the boundary which have a boundary indicator listed in the given + * set. */ types::global_dof_index n_boundary_dofs(const std::set &boundary_ids) const; diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index 21bd935fbf..2145bd5b90 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -35,6 +35,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -1096,36 +1097,51 @@ template types::global_dof_index DoFHandler::n_boundary_dofs() const { - std::set boundary_dofs; - - const unsigned int dofs_per_face = get_fe().dofs_per_face; - std::vector dofs_on_face(dofs_per_face); - - // loop over all faces of all cells - // and see whether they are at a - // boundary. note (i) that we visit - // interior faces twice (which we - // don't care about) but exterior - // faces only once as is - // appropriate, and (ii) that we - // need not take special care of - // single lines (using - // @p{cell->has_boundary_lines}), - // since we do not support - // boundaries of dimension dim-2, - // and so every boundary line is + const FiniteElement &fe = get_fe(); + const unsigned int dofs_per_face = fe.n_dofs_per_face(); + std::vector dofs_on_face(dofs_per_face); + + const IndexSet &owned_dofs = locally_owned_dofs(); + + std::unordered_set boundary_dof_indices; + // Loop over all faces of all cells and see whether they are at a boundary. + // note (i) that we visit interior faces twice (which we don't care about) but + // exterior faces only once as is appropriate, and (ii) that we need not take + // special care of single lines (using @p{cell->has_boundary_lines}), since we + // do not support boundaries of dimension dim-2, and so every boundary line is // also part of a boundary face. - active_cell_iterator cell = begin_active(), endc = end(); - for (; cell != endc; ++cell) - for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) - if (cell->at_boundary(f)) + for (const auto &cell : this->active_cell_iterators()) + { + if (cell->is_locally_owned() && cell->at_boundary()) { - cell->face(f)->get_dof_indices(dofs_on_face); - for (unsigned int i = 0; i < dofs_per_face; ++i) - boundary_dofs.insert(dofs_on_face[i]); + for (unsigned int iface = 0; + iface < dealii::GeometryInfo::faces_per_cell; + ++iface) + { + const auto face = cell->face(iface); + if (face->at_boundary()) + { + cell->face(iface)->get_dof_indices(dofs_on_face); + for (unsigned int idof_face = 0; idof_face < dofs_per_face; + ++idof_face) + { + const unsigned int global_idof_index = + dofs_on_face[idof_face]; + + // If dof is not already in our hashmap, then insert it + // into our set of surface indices. However, do not add if + // it is not locally owned, it will get picked up by + // another processor + if (owned_dofs.is_element(global_idof_index)) + { + boundary_dof_indices.insert(global_idof_index); + } + } + } + } } - - return boundary_dofs.size(); + } + return boundary_dof_indices.size(); } @@ -1139,28 +1155,56 @@ DoFHandler::n_boundary_dofs( boundary_ids.end(), ExcInvalidBoundaryIndicator()); - std::set boundary_dofs; - - const unsigned int dofs_per_face = get_fe().dofs_per_face; - std::vector dofs_on_face(dofs_per_face); - - // same as in the previous - // function, but with a different - // check for the boundary indicator - active_cell_iterator cell = begin_active(), endc = end(); - for (; cell != endc; ++cell) - for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) - if (cell->at_boundary(f) && - (std::find(boundary_ids.begin(), - boundary_ids.end(), - cell->face(f)->boundary_id()) != boundary_ids.end())) + const FiniteElement &fe = get_fe(); + const unsigned int dofs_per_face = fe.n_dofs_per_face(); + std::vector dofs_on_face(dofs_per_face); + + const IndexSet &owned_dofs = locally_owned_dofs(); + + std::unordered_set boundary_dof_indices; + // Loop over all faces of all cells and see whether they are at a boundary. + // note (i) that we visit interior faces twice (which we don't care about) but + // exterior faces only once as is appropriate, and (ii) that we need not take + // special care of single lines (using @p{cell->has_boundary_lines}), since we + // do not support boundaries of dimension dim-2, and so every boundary line is + // also part of a boundary face. + for (const auto &cell : this->active_cell_iterators()) + { + if (cell->is_locally_owned() && cell->at_boundary()) { - cell->face(f)->get_dof_indices(dofs_on_face); - for (unsigned int i = 0; i < dofs_per_face; ++i) - boundary_dofs.insert(dofs_on_face[i]); + for (unsigned int iface = 0; + iface < dealii::GeometryInfo::faces_per_cell; + ++iface) + { + const auto face = cell->face(iface); + const unsigned int boundary_id = face->boundary_id(); + + if (face->at_boundary() && + std::find(boundary_ids.begin(), + boundary_ids.end(), + boundary_id) != boundary_ids.end()) + { + cell->face(iface)->get_dof_indices(dofs_on_face); + for (unsigned int idof_face = 0; idof_face < dofs_per_face; + ++idof_face) + { + const unsigned int global_idof_index = + dofs_on_face[idof_face]; + + // If dof is not already in our hashmap, then insert it + // into our set of surface indices. However, do not add if + // it is not locally owned, it will get picked up by + // another processor + if (owned_dofs.is_element(global_idof_index)) + { + boundary_dof_indices.insert(global_idof_index); + } + } + } + } } - - return boundary_dofs.size(); + } + return boundary_dof_indices.size(); } diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc index 27a0948ab2..c5e1678176 100644 --- a/source/hp/dof_handler.cc +++ b/source/hp/dof_handler.cc @@ -42,6 +42,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -1431,10 +1432,12 @@ namespace hp { Assert(fe_collection.size() > 0, ExcNoFESelected()); - std::set boundary_dofs; - std::vector dofs_on_face; + std::unordered_set boundary_dofs; + std::vector dofs_on_face; dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face()); + const IndexSet &owned_dofs = locally_owned_dofs(); + // loop over all faces to check whether they are at a // boundary. note that we need not take special care of single // lines in 3d (using @p{cell->has_boundary_lines}), since we do @@ -1444,17 +1447,26 @@ namespace hp cell = this->begin_active(), endc = this->end(); for (; cell != endc; ++cell) - for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) - if (cell->at_boundary(f)) - { - const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; - dofs_on_face.resize(dofs_per_face); + if (cell->is_locally_owned() && cell->at_boundary()) + { + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + if (cell->at_boundary(f)) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_face.resize(dofs_per_face); - cell->face(f)->get_dof_indices(dofs_on_face, - cell->active_fe_index()); - for (unsigned int i = 0; i < dofs_per_face; ++i) - boundary_dofs.insert(dofs_on_face[i]); - } + cell->face(f)->get_dof_indices(dofs_on_face, + cell->active_fe_index()); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int global_idof_index = dofs_on_face[i]; + if (owned_dofs.is_element(global_idof_index)) + { + boundary_dofs.insert(global_idof_index); + } + } + } + } return boundary_dofs.size(); } @@ -1472,27 +1484,38 @@ namespace hp // same as above, but with additional checks for set of boundary // indicators - std::set boundary_dofs; - std::vector dofs_on_face; + std::unordered_set boundary_dofs; + std::vector dofs_on_face; dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face()); + const IndexSet &owned_dofs = locally_owned_dofs(); + typename HpDoFHandler::active_cell_iterator cell = this->begin_active(), endc = this->end(); for (; cell != endc; ++cell) - for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) - if (cell->at_boundary(f) && - (boundary_ids.find(cell->face(f)->boundary_id()) != - boundary_ids.end())) - { - const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; - dofs_on_face.resize(dofs_per_face); + if (cell->is_locally_owned() && cell->at_boundary()) + { + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + if (cell->at_boundary(f) && + (boundary_ids.find(cell->face(f)->boundary_id()) != + boundary_ids.end())) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_face.resize(dofs_per_face); - cell->face(f)->get_dof_indices(dofs_on_face, - cell->active_fe_index()); - for (unsigned int i = 0; i < dofs_per_face; ++i) - boundary_dofs.insert(dofs_on_face[i]); - } + cell->face(f)->get_dof_indices(dofs_on_face, + cell->active_fe_index()); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int global_idof_index = dofs_on_face[i]; + if (owned_dofs.is_element(global_idof_index)) + { + boundary_dofs.insert(global_idof_index); + } + } + } + } return boundary_dofs.size(); } diff --git a/tests/dofs/n_boundary_dofs_01.cc b/tests/dofs/n_boundary_dofs_01.cc new file mode 100644 index 0000000000..80b12ed21a --- /dev/null +++ b/tests/dofs/n_boundary_dofs_01.cc @@ -0,0 +1,65 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2008 - 2018 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 that DoFHandler::n_boundary_dofs() yields the correct +// results for parallel::distributed::Triangulations + +#include + +#include + +#include + +#include + +#include "../tests.h" + + +template +void +test() +{ + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + + FE_Q fe(2); + DoFHandler dof_handler(triangulation); + + dof_handler.distribute_dofs(fe); + + const unsigned int local_boundary_dofs = dof_handler.n_boundary_dofs(); + const unsigned int global_boundary_dofs = + Utilities::MPI::sum(local_boundary_dofs, MPI_COMM_WORLD); + deallog << global_boundary_dofs << std::endl; +} + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + deallog.push("2d"); + test<2>(); + deallog.pop(); + + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/dofs/n_boundary_dofs_01.with_p4est=true.mpirun=1.output b/tests/dofs/n_boundary_dofs_01.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..1108787d09 --- /dev/null +++ b/tests/dofs/n_boundary_dofs_01.with_p4est=true.mpirun=1.output @@ -0,0 +1,3 @@ + +DEAL:2d::64 +DEAL:3d::1538 diff --git a/tests/dofs/n_boundary_dofs_01.with_p4est=true.mpirun=2.output b/tests/dofs/n_boundary_dofs_01.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..1108787d09 --- /dev/null +++ b/tests/dofs/n_boundary_dofs_01.with_p4est=true.mpirun=2.output @@ -0,0 +1,3 @@ + +DEAL:2d::64 +DEAL:3d::1538 diff --git a/tests/hp/n_boundary_dofs_01.cc b/tests/hp/n_boundary_dofs_01.cc new file mode 100644 index 0000000000..b7a7b03629 --- /dev/null +++ b/tests/hp/n_boundary_dofs_01.cc @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2008 - 2018 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 that DoFHandler::n_boundary_dofs() yields the correct +// results for parallel::distributed::Triangulations + +#include + +#include + +#include + +#include +#include + +#include "../tests.h" + + +template +void +test() +{ + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + + FE_Q fe(2); + hp::FECollection fe_coll(fe); + hp::DoFHandler dof_handler(triangulation); + + dof_handler.distribute_dofs(fe_coll); + + const unsigned int local_boundary_dofs = dof_handler.n_boundary_dofs(); + const unsigned int global_boundary_dofs = + Utilities::MPI::sum(local_boundary_dofs, MPI_COMM_WORLD); + deallog << global_boundary_dofs << std::endl; +} + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + deallog.push("2d"); + test<2>(); + deallog.pop(); + + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/hp/n_boundary_dofs_01.with_p4est=true.mpirun=1.output b/tests/hp/n_boundary_dofs_01.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..1108787d09 --- /dev/null +++ b/tests/hp/n_boundary_dofs_01.with_p4est=true.mpirun=1.output @@ -0,0 +1,3 @@ + +DEAL:2d::64 +DEAL:3d::1538 diff --git a/tests/hp/n_boundary_dofs_01.with_p4est=true.mpirun=2.output b/tests/hp/n_boundary_dofs_01.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..1108787d09 --- /dev/null +++ b/tests/hp/n_boundary_dofs_01.with_p4est=true.mpirun=2.output @@ -0,0 +1,3 @@ + +DEAL:2d::64 +DEAL:3d::1538 -- 2.39.5