From: Martin Kronbichler Date: Tue, 26 Jan 2016 08:02:10 +0000 (+0100) Subject: Fix 2:1 balance on artificial cells over periodic boundaries. X-Git-Tag: v8.4.0-rc2~48^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=64d5ad5efb20bb3f7006d45e0fd7c058b638207f;p=dealii.git Fix 2:1 balance on artificial cells over periodic boundaries. --- diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 8db94c26f9..b0f1149c6a 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -3122,6 +3122,204 @@ namespace parallel + namespace + { + template + void + identify_periodic_vertices_recursively(const GridTools::PeriodicFacePair &periodic, + std::vector &topological_vertex_numbering) + { + // for hanging nodes we will consider all necessary coupling already + // on the parent level, so we just need to consider neighbors of the + // same level + if (periodic.cell[0]->has_children() && + periodic.cell[1]->has_children()) + { + // copy orientations etc. from parent to child + GridTools::PeriodicFacePair periodic_child = periodic; + + // find appropriate pairs of child elements + for (unsigned int cf=0; cfface(periodic.face_idx[0])->n_children(); ++cf) + { + unsigned int c=0; + for (; cn_children(); ++c) + { + if (periodic.cell[0]->child(c)->face(periodic.face_idx[0]) == + periodic.cell[0]->face(periodic.face_idx[0])->child(cf)) + break; + } + Assert(c < periodic.cell[0]->n_children(), + ExcMessage("Face child not found")); + periodic_child.cell[0] = periodic.cell[0]->child(c); + periodic_child.face_idx[0] = periodic.face_idx[0]; + + c=0; + for (; cn_children(); ++c) + { + if (periodic.cell[1]->child(c)->face(periodic.face_idx[1]) == + periodic.cell[1]->face(periodic.face_idx[1])->child(cf)) + break; + } + Assert(c < periodic.cell[1]->n_children(), + ExcMessage("Face child not found")); + periodic_child.cell[1] = periodic.cell[1]->child(c); + periodic_child.face_idx[1] = periodic.face_idx[1]; + + // recursive call into children + identify_periodic_vertices_recursively (periodic_child, + topological_vertex_numbering); + } + } + + // TODO: fix non-standard orientation + Assert(periodic.orientation[0] == true && + periodic.orientation[1] == false && + periodic.orientation[2] == false, + ExcNotImplemented()); + + // TODO: might need to fix the case with corners where several PBC + // meet (in most cases it will get fixed eventually because we visit + // corner dofs more than once, and each time we pass around we connect + // two more vertices; in the worst case of 3 periodicity connects in + // 3D, we are safe after three levels) + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + const unsigned int vi0 = topological_vertex_numbering[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v)]; + const unsigned int vi1 = topological_vertex_numbering[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)]; + const unsigned int min_index = std::min(vi0, vi1); + topological_vertex_numbering[vi0] = topological_vertex_numbering[vi1] = min_index; + } + } + + template + void enforce_mesh_balance_over_periodic_boundaries + (Triangulation &tria, + std::vector::cell_iterator> > periodic_face_pairs_level_0) + { + if (periodic_face_pairs_level_0.empty()) + return; + + std::vector topological_vertex_numbering(tria.n_vertices()); + for (unsigned int i=0; i vertex_level(tria.n_vertices()); + while (continue_iterating) + { + // store highest level one of the cells adjacent to a vertex + // belongs to + std::fill (vertex_level.begin(), vertex_level.end(), 0); + typename Triangulation::active_cell_iterator + cell = tria.begin_active(), endc = tria.end(); + for (; cell!=endc; ++cell) + { + if (cell->refine_flag_set()) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + vertex_level[topological_vertex_numbering[cell->vertex_index(vertex)]] + = std::max (vertex_level[topological_vertex_numbering[cell->vertex_index(vertex)]], + cell->level()+1); + else if (!cell->coarsen_flag_set()) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + vertex_level[topological_vertex_numbering[cell->vertex_index(vertex)]] + = std::max (vertex_level[topological_vertex_numbering[cell->vertex_index(vertex)]], + cell->level()); + else + { + // if coarsen flag is set then tentatively assume + // that the cell will be coarsened. this isn't + // always true (the coarsen flag could be removed + // again) and so we may make an error here. we try + // to correct this by iterating over the entire + // process until we are converged + Assert (cell->coarsen_flag_set(), ExcInternalError()); + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + vertex_level[topological_vertex_numbering[cell->vertex_index(vertex)]] + = std::max (vertex_level[topological_vertex_numbering[cell->vertex_index(vertex)]], + cell->level()-1); + } + } + + continue_iterating = false; + + // loop over all cells in reverse order. do so because we + // can then update the vertex levels on the adjacent + // vertices and maybe already flag additional cells in this + // loop + // + // note that not only may we have to add additional + // refinement flags, but we will also have to remove + // coarsening flags on cells adjacent to vertices that will + // see refinement + for (cell=tria.last_active(); cell != endc; --cell) + if (cell->refine_flag_set() == false) + { + for (unsigned int vertex=0; + vertex::vertices_per_cell; ++vertex) + if (vertex_level[topological_vertex_numbering[cell->vertex_index(vertex)]] >= + cell->level()+1) + { + // remove coarsen flag... + cell->clear_coarsen_flag(); + + // ...and if necessary also refine the current + // cell, at the same time updating the level + // information about vertices + if (vertex_level[topological_vertex_numbering[cell->vertex_index(vertex)]] > + cell->level()+1) + { + cell->set_refine_flag(); + continue_iterating = true; + + for (unsigned int v=0; v::vertices_per_cell; + ++v) + vertex_level[topological_vertex_numbering[cell->vertex_index(v)]] + = std::max (vertex_level[topological_vertex_numbering[cell->vertex_index(v)]], + cell->level()+1); + } + + // continue and see whether we may, for example, + // go into the inner 'if' above based on a + // different vertex + } + } + + // clear coarsen flag if not all children were marked + for (typename Triangulation::cell_iterator cell = tria.begin(); + cell!=tria.end(); ++cell) + { + // nothing to do if we are already on the finest level + if (cell->active()) + continue; + + const unsigned int n_children=cell->n_children(); + unsigned int flagged_children=0; + for (unsigned int child=0; childchild(child)->active() && + cell->child(child)->coarsen_flag_set()) + ++flagged_children; + + // if not all children were flagged for coarsening, remove + // coarsen flags + if (flagged_children < n_children) + for (unsigned int child=0; childchild(child)->active()) + cell->child(child)->clear_coarsen_flag(); + } + } + } + } + + + template void Triangulation::copy_local_forest_to_triangulation () @@ -3272,6 +3470,11 @@ namespace parallel // fix all the flags to make sure we have a consistent mesh this->prepare_coarsening_and_refinement (); + // enforce 2:1 mesh balance over periodic boundaries + if (this->smooth_grid & dealii::Triangulation::limit_level_difference_at_vertices) + enforce_mesh_balance_over_periodic_boundaries(*this, + periodic_face_pairs_level_0); + // see if any flags are still set mesh_changed = false; for (typename Triangulation::active_cell_iterator diff --git a/tests/mpi/distribute_mg_dofs_periodicity_01.cc b/tests/mpi/distribute_mg_dofs_periodicity_01.cc new file mode 100644 index 0000000000..c6aacdc0b1 --- /dev/null +++ b/tests/mpi/distribute_mg_dofs_periodicity_01.cc @@ -0,0 +1,109 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2008 - 2015 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// this tests distribute_mg_dofs on a mesh with periodic boundary +// conditions. we used to forget to enforce the 2:1 level balance over +// periodic boundaries in artificial cells, which lead to the case where the +// coarser MG level did not provide the correct ghost layer (leading to a +// deadlock in the distribute_dofs function for 6 procs in 2D and 28 procs in +// 3D) + +#include "../tests.h" + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include + +// This is C++: +#include +#include + + +template +void +set_owner_recursively(typename Triangulation::cell_iterator cell, + const int owner, + Vector &vec) +{ + if (cell->active()) + { + vec(cell->active_cell_index()) = owner; + } + else + for (unsigned int c=0; cn_children(); ++c) + set_owner_recursively(cell->child(c), owner, vec); +} + +template +void test() +{ + parallel::distributed::Triangulation + tria(MPI_COMM_WORLD, dealii::Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + GridGenerator::subdivided_hyper_cube (tria, 1); + // set periodic boundary conditions in x and z directions + for (typename Triangulation::cell_iterator cell=tria.begin(); + cell != tria.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (f / 2 != 1 && cell->at_boundary(f)) + cell->face(f)->set_all_boundary_ids(f+10); + + std::vector::cell_iterator> > periodic_faces; + GridTools::collect_periodic_faces(tria, 0+10, 1+10, 0, periodic_faces); + if (dim == 3) + GridTools::collect_periodic_faces(tria, 4+10, 5+10, 2, periodic_faces); + tria.add_periodicity(periodic_faces); + tria.refine_global(5); + + FE_Q fe(1); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + deallog << "Number of cells: " << tria.n_global_active_cells() + << std::endl + << "Number of DoFs: " << dof_handler.n_dofs() << std::endl; + + dof_handler.distribute_mg_dofs(fe); + deallog << "Number of DoFs per level: "; + for (unsigned int level=0; level(); + test<3>(); + + return 0; +} diff --git a/tests/mpi/distribute_mg_dofs_periodicity_01.mpirun=28.output b/tests/mpi/distribute_mg_dofs_periodicity_01.mpirun=28.output new file mode 100644 index 0000000000..6180333698 --- /dev/null +++ b/tests/mpi/distribute_mg_dofs_periodicity_01.mpirun=28.output @@ -0,0 +1,7 @@ + +DEAL::Number of cells: 1024 +DEAL::Number of DoFs: 1089 +DEAL::Number of DoFs per level: 4 9 25 81 289 1089 +DEAL::Number of cells: 32768 +DEAL::Number of DoFs: 35937 +DEAL::Number of DoFs per level: 8 27 125 729 4913 35937 diff --git a/tests/mpi/distribute_mg_dofs_periodicity_01.mpirun=6.output b/tests/mpi/distribute_mg_dofs_periodicity_01.mpirun=6.output new file mode 100644 index 0000000000..6180333698 --- /dev/null +++ b/tests/mpi/distribute_mg_dofs_periodicity_01.mpirun=6.output @@ -0,0 +1,7 @@ + +DEAL::Number of cells: 1024 +DEAL::Number of DoFs: 1089 +DEAL::Number of DoFs per level: 4 9 25 81 289 1089 +DEAL::Number of cells: 32768 +DEAL::Number of DoFs: 35937 +DEAL::Number of DoFs per level: 8 27 125 729 4913 35937