+ namespace
+ {
+ template <typename ITERATOR>
+ void
+ identify_periodic_vertices_recursively(const GridTools::PeriodicFacePair<ITERATOR> &periodic,
+ std::vector<unsigned int> &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<ITERATOR> periodic_child = periodic;
+
+ // find appropriate pairs of child elements
+ for (unsigned int cf=0; cf<periodic.cell[0]->face(periodic.face_idx[0])->n_children(); ++cf)
+ {
+ unsigned int c=0;
+ for (; c<periodic.cell[0]->n_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 (; c<periodic.cell[1]->n_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<GeometryInfo<ITERATOR::AccessorType::dimension-1>::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 <int dim, int spacedim>
+ void enforce_mesh_balance_over_periodic_boundaries
+ (Triangulation<dim,spacedim> &tria,
+ std::vector<GridTools::PeriodicFacePair<typename dealii::Triangulation<dim,spacedim>::cell_iterator> > periodic_face_pairs_level_0)
+ {
+ if (periodic_face_pairs_level_0.empty())
+ return;
+
+ std::vector<unsigned int> topological_vertex_numbering(tria.n_vertices());
+ for (unsigned int i=0; i<topological_vertex_numbering.size(); ++i)
+ topological_vertex_numbering[i] = i;
+ for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
+ identify_periodic_vertices_recursively(periodic_face_pairs_level_0[i],
+ topological_vertex_numbering);
+
+ // this code is replicated from grid/tria.cc but using an indirection
+ // for periodic boundary conditions
+ bool continue_iterating = true;
+ std::vector<int> 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<dim,spacedim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<dim,spacedim>::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; child<n_children; ++child)
+ if (cell->child(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; child<n_children; ++child)
+ if (cell->child(child)->active())
+ cell->child(child)->clear_coarsen_flag();
+ }
+ }
+ }
+ }
+
+
+
template <int dim, int spacedim>
void
Triangulation<dim,spacedim>::copy_local_forest_to_triangulation ()
// 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<dim,spacedim>::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<dim,spacedim>::active_cell_iterator
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+
+// This is C++:
+#include <fstream>
+#include <sstream>
+
+
+template <int dim>
+void
+set_owner_recursively(typename Triangulation<dim>::cell_iterator cell,
+ const int owner,
+ Vector<double> &vec)
+{
+ if (cell->active())
+ {
+ vec(cell->active_cell_index()) = owner;
+ }
+ else
+ for (unsigned int c=0; c<cell->n_children(); ++c)
+ set_owner_recursively<dim>(cell->child(c), owner, vec);
+}
+
+template <int dim>
+void test()
+{
+ parallel::distributed::Triangulation<dim>
+ tria(MPI_COMM_WORLD, dealii::Triangulation<dim>::limit_level_difference_at_vertices,
+ parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+ GridGenerator::subdivided_hyper_cube (tria, 1);
+ // set periodic boundary conditions in x and z directions
+ for (typename Triangulation<dim>::cell_iterator cell=tria.begin();
+ cell != tria.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (f / 2 != 1 && cell->at_boundary(f))
+ cell->face(f)->set_all_boundary_ids(f+10);
+
+ std::vector<GridTools::PeriodicFacePair<typename Triangulation<dim>::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<dim> fe(1);
+ DoFHandler<dim> 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<tria.n_global_levels(); ++level)
+ deallog << dof_handler.n_dofs(level) << " ";
+ deallog << std::endl;
+}
+
+
+
+int main (int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+ mpi_initlog();
+
+ test<2>();
+ test<3>();
+
+ return 0;
+}
projects / dealii.git / commitdiff