From: Peter Munch Date: Sat, 2 Jan 2021 22:25:58 +0000 (+0100) Subject: Add a version of update_neighbors() for simplex meshes X-Git-Tag: v9.3.0-rc1~680^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ac288e3f66633540a69c14aaf8e9888e910d71a0;p=dealii.git Add a version of update_neighbors() for simplex meshes --- diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 1feb29179c..a5b33c7f9c 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -698,187 +698,6 @@ namespace } - - /** - * For a given triangulation: set up the - * neighbor information on all cells. - */ - template - void update_neighbors(Triangulation<1, spacedim> &) - {} - - - template - void - update_neighbors(Triangulation &triangulation) - { - // each face can be neighbored on two sides - // by cells. according to the face's - // intrinsic normal we define the left - // neighbor as the one for which the face - // normal points outward, and store that - // one first; the second one is then - // the right neighbor for which the - // face normal points inward. This - // information depends on the type of cell - // and local number of face for the - // 'standard ordering and orientation' of - // faces and then on the face_orientation - // information for the real mesh. Set up a - // table to have fast access to those - // offsets (0 for left and 1 for - // right). Some of the values are invalid - // as they reference too large face - // numbers, but we just leave them at a - // zero value. - // - // Note, that in 2d for lines as faces the - // normal direction given in the - // GeometryInfo class is not consistent. We - // thus define here that the normal for a - // line points to the right if the line - // points upwards. - // - // There is one more point to - // consider, however: if we have - // dim::cell_iterator dummy; - std::vector::cell_iterator> - adjacent_cells(2 * triangulation.n_raw_faces(), dummy); - - for (const auto &cell : triangulation.cell_iterators()) - for (auto f : cell->face_indices()) - { - const typename Triangulation::face_iterator face = - cell->face(f); - - const unsigned int offset = - (cell->direction_flag() ? - left_right_offset[dim - 2][f][cell->face_orientation(f)] : - 1 - left_right_offset[dim - 2][f][cell->face_orientation(f)]); - - adjacent_cells[2 * face->index() + offset] = cell; - - // if this cell is not refined, but the - // face is, then we'll have to set our - // cell as neighbor for the child faces - // as well. Fortunately the normal - // orientation of children will be just - // the same. - if (dim == 2) - { - if (cell->is_active() && face->has_children()) - { - adjacent_cells[2 * face->child(0)->index() + offset] = cell; - adjacent_cells[2 * face->child(1)->index() + offset] = cell; - } - } - else // -> dim == 3 - { - // We need the same as in 2d - // here. Furthermore, if the face is - // refined with cut_x or cut_y then - // those children again in the other - // direction, and if this cell is - // refined isotropically (along the - // face) then the neighbor will - // (probably) be refined as cut_x or - // cut_y along the face. For those - // neighboring children cells, their - // neighbor will be the current, - // inactive cell, as our children are - // too fine to be neighbors. Catch that - // case by also acting on inactive - // cells with isotropic refinement - // along the face. If the situation - // described is not present, the data - // will be overwritten later on when we - // visit cells on finer levels, so no - // harm will be done. - if (face->has_children() && - (cell->is_active() || - GeometryInfo::face_refinement_case( - cell->refinement_case(), f) == - RefinementCase::isotropic_refinement)) - { - for (unsigned int c = 0; c < face->n_children(); ++c) - adjacent_cells[2 * face->child(c)->index() + offset] = cell; - if (face->child(0)->has_children()) - { - adjacent_cells[2 * face->child(0)->child(0)->index() + - offset] = cell; - adjacent_cells[2 * face->child(0)->child(1)->index() + - offset] = cell; - } - if (face->child(1)->has_children()) - { - adjacent_cells[2 * face->child(1)->child(0)->index() + - offset] = cell; - adjacent_cells[2 * face->child(1)->child(1)->index() + - offset] = cell; - } - } // if cell active and face refined - } // else -> dim==3 - } // for all faces of all cells - - // now loop again over all cells and set the - // corresponding neighbor cell. Note, that we - // have to use the opposite of the - // left_right_offset in this case as we want - // the offset of the neighbor, not our own. - for (const auto &cell : triangulation.cell_iterators()) - for (auto f : cell->face_indices()) - { - const unsigned int offset = - (cell->direction_flag() ? - left_right_offset[dim - 2][f][cell->face_orientation(f)] : - 1 - left_right_offset[dim - 2][f][cell->face_orientation(f)]); - cell->set_neighbor( - f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]); - } - } - - template void update_periodic_face_map_recursively( @@ -1810,6 +1629,12 @@ namespace internal */ virtual ~Policy() = default; + /** + * Update neighbors. + */ + virtual void + update_neighbors(Triangulation &tria) = 0; + /** * Delete children of given cell. */ @@ -1869,6 +1694,12 @@ namespace internal class PolicyWrapper : public Policy { public: + void + update_neighbors(Triangulation &tria) override + { + T::update_neighbors(tria); + } + void delete_children( Triangulation & tria, @@ -2312,6 +2143,188 @@ namespace internal } + + template + static void update_neighbors(Triangulation<1, spacedim> &) + {} + + + template + static void + update_neighbors(Triangulation &triangulation) + { + // each face can be neighbored on two sides + // by cells. according to the face's + // intrinsic normal we define the left + // neighbor as the one for which the face + // normal points outward, and store that + // one first; the second one is then + // the right neighbor for which the + // face normal points inward. This + // information depends on the type of cell + // and local number of face for the + // 'standard ordering and orientation' of + // faces and then on the face_orientation + // information for the real mesh. Set up a + // table to have fast access to those + // offsets (0 for left and 1 for + // right). Some of the values are invalid + // as they reference too large face + // numbers, but we just leave them at a + // zero value. + // + // Note, that in 2d for lines as faces the + // normal direction given in the + // GeometryInfo class is not consistent. We + // thus define here that the normal for a + // line points to the right if the line + // points upwards. + // + // There is one more point to + // consider, however: if we have + // dim::cell_iterator dummy; + std::vector::cell_iterator> + adjacent_cells(2 * triangulation.n_raw_faces(), dummy); + + for (const auto &cell : triangulation.cell_iterators()) + for (auto f : cell->face_indices()) + { + const typename Triangulation::face_iterator face = + cell->face(f); + + const unsigned int offset = + (cell->direction_flag() ? + left_right_offset[dim - 2][f][cell->face_orientation(f)] : + 1 - + left_right_offset[dim - 2][f][cell->face_orientation(f)]); + + adjacent_cells[2 * face->index() + offset] = cell; + + // if this cell is not refined, but the + // face is, then we'll have to set our + // cell as neighbor for the child faces + // as well. Fortunately the normal + // orientation of children will be just + // the same. + if (dim == 2) + { + if (cell->is_active() && face->has_children()) + { + adjacent_cells[2 * face->child(0)->index() + offset] = + cell; + adjacent_cells[2 * face->child(1)->index() + offset] = + cell; + } + } + else // -> dim == 3 + { + // We need the same as in 2d + // here. Furthermore, if the face is + // refined with cut_x or cut_y then + // those children again in the other + // direction, and if this cell is + // refined isotropically (along the + // face) then the neighbor will + // (probably) be refined as cut_x or + // cut_y along the face. For those + // neighboring children cells, their + // neighbor will be the current, + // inactive cell, as our children are + // too fine to be neighbors. Catch that + // case by also acting on inactive + // cells with isotropic refinement + // along the face. If the situation + // described is not present, the data + // will be overwritten later on when we + // visit cells on finer levels, so no + // harm will be done. + if (face->has_children() && + (cell->is_active() || + GeometryInfo::face_refinement_case( + cell->refinement_case(), f) == + RefinementCase::isotropic_refinement)) + { + for (unsigned int c = 0; c < face->n_children(); ++c) + adjacent_cells[2 * face->child(c)->index() + offset] = + cell; + if (face->child(0)->has_children()) + { + adjacent_cells[2 * face->child(0)->child(0)->index() + + offset] = cell; + adjacent_cells[2 * face->child(0)->child(1)->index() + + offset] = cell; + } + if (face->child(1)->has_children()) + { + adjacent_cells[2 * face->child(1)->child(0)->index() + + offset] = cell; + adjacent_cells[2 * face->child(1)->child(1)->index() + + offset] = cell; + } + } // if cell active and face refined + } // else -> dim==3 + } // for all faces of all cells + + // now loop again over all cells and set the + // corresponding neighbor cell. Note, that we + // have to use the opposite of the + // left_right_offset in this case as we want + // the offset of the neighbor, not our own. + for (const auto &cell : triangulation.cell_iterators()) + for (auto f : cell->face_indices()) + { + const unsigned int offset = + (cell->direction_flag() ? + left_right_offset[dim - 2][f][cell->face_orientation(f)] : + 1 - + left_right_offset[dim - 2][f][cell->face_orientation(f)]); + cell->set_neighbor( + f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]); + } + } + + /** * Create a triangulation from given data. */ @@ -9913,6 +9926,75 @@ namespace internal */ struct ImplementationMixedMesh { + template + static void update_neighbors(Triangulation<1, spacedim> &) + {} + + template + void static update_neighbors(Triangulation &triangulation) + { + std::vector> adjacent_cells( + 2 * triangulation.n_raw_faces(), + {numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}); + + const auto set_entry = [&](const auto &face_index, const auto &cell) { + const std::pair cell_pair = { + cell->level(), cell->index()}; + unsigned int index; + + if (adjacent_cells[2 * face_index].first == + numbers::invalid_unsigned_int && + adjacent_cells[2 * face_index].second == + numbers::invalid_unsigned_int) + { + index = 2 * face_index + 0; + } + else + { + Assert(((adjacent_cells[2 * face_index + 1].first == + numbers::invalid_unsigned_int) && + (adjacent_cells[2 * face_index + 1].second == + numbers::invalid_unsigned_int)), + ExcNotImplemented()); + index = 2 * face_index + 1; + } + + adjacent_cells[index] = cell_pair; + }; + + const auto get_entry = + [&](const auto &face_index, + const auto &cell) -> TriaIterator> { + auto test = adjacent_cells[2 * face_index]; + + if (test == std::pair(cell->level(), + cell->index())) + test = adjacent_cells[2 * face_index + 1]; + + if ((test.first != numbers::invalid_unsigned_int) && + (test.second != numbers::invalid_unsigned_int)) + return TriaIterator>(&triangulation, + test.first, + test.second); + else + return typename Triangulation::cell_iterator(); + }; + + for (const auto &cell : triangulation.cell_iterators()) + for (const auto &face : cell->face_iterators()) + { + set_entry(face->index(), cell); + + if (cell->is_active() && face->has_children()) + for (unsigned int c = 0; c < face->n_children(); ++c) + set_entry(face->child(c)->index(), cell); + } + + for (const auto &cell : triangulation.cell_iterators()) + for (auto f : cell->face_indices()) + cell->set_neighbor(f, get_entry(cell->face(f)->index(), cell)); + } + template static void delete_children( @@ -13300,7 +13382,7 @@ Triangulation::execute_coarsening_and_refinement() // finally build up neighbor connectivity information, and set // active cell indices - update_neighbors(*this); + this->policy->update_neighbors(*this); reset_active_cell_indices(); reset_global_cell_indices(); // TODO: better place? diff --git a/tests/simplex/refinement_02.cc b/tests/simplex/refinement_02.cc new file mode 100644 index 0000000000..d74bfe7339 --- /dev/null +++ b/tests/simplex/refinement_02.cc @@ -0,0 +1,77 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// Test the correct setup of neighbors during refinement of simplex mesh. + +#include +#include + +#include + +#include "../tests.h" + +#include "./simplex_grids.h" + +void +test(const unsigned int v) +{ + const unsigned int dim = 2; + + Triangulation tria; + + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + + if (v == 1) + { + tria.begin()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + else if (v == 2) + { + tria.refine_global(1); + } + else if (v == 3) + { + tria.begin()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.begin_active(1)->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + for (const auto &cell : tria.active_cell_iterators()) + { + deallog << cell->level() << ":" << cell->index() << " "; + for (const auto f : cell->face_indices()) + if (cell->at_boundary(f)) + deallog << "-:- "; + else + deallog << cell->neighbor_level(f) << ":" << cell->neighbor_index(f) + << " "; + deallog << std::endl; + } + + deallog << std::endl; +} + +int +main() +{ + initlog(); + test(0); // no refinement + test(1); // refinement of a single cell + test(2); // global refinement +} diff --git a/tests/simplex/refinement_02.with_simplex_support=on.output b/tests/simplex/refinement_02.with_simplex_support=on.output new file mode 100644 index 0000000000..03c1627a4a --- /dev/null +++ b/tests/simplex/refinement_02.with_simplex_support=on.output @@ -0,0 +1,19 @@ + +DEAL::0:0 -:- 0:1 -:- +DEAL::0:1 -:- 0:0 -:- +DEAL:: +DEAL::0:1 -:- 0:0 -:- +DEAL::1:0 -:- 1:3 -:- +DEAL::1:1 -:- 0:1 1:3 +DEAL::1:2 1:3 0:1 -:- +DEAL::1:3 1:1 1:2 1:0 +DEAL:: +DEAL::1:0 -:- 1:3 -:- +DEAL::1:1 -:- 1:6 1:3 +DEAL::1:2 1:3 1:5 -:- +DEAL::1:3 1:1 1:2 1:0 +DEAL::1:4 -:- 1:7 -:- +DEAL::1:5 -:- 1:2 1:7 +DEAL::1:6 1:7 1:1 -:- +DEAL::1:7 1:5 1:6 1:4 +DEAL::