From c938272a71a29a7410ac4f730c9b3b6563e37997 Mon Sep 17 00:00:00 2001 From: David Wells Date: Wed, 9 Jun 2021 16:07:26 -0400 Subject: [PATCH] Improve comments. --- source/grid/tria.cc | 140 +++++++++++++++++++++++++++++--------------- 1 file changed, 94 insertions(+), 46 deletions(-) diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 6c2c231e75..27a67d8615 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -5101,22 +5101,20 @@ namespace internal triangulation.vertices_used.size(), ExcInternalError()); - { - typename Triangulation::raw_cell_iterator - cell = triangulation.begin_active(triangulation.levels.size() - 1), - endc = triangulation.end(); - for (; cell != endc; ++cell) - if (cell->used()) - if (cell->refine_flag_set()) - { - triangulation.levels.push_back( - std::make_unique< - internal::TriangulationImplementation::TriaLevel>(dim)); - break; - } - } - + // Check whether a new level is needed. We have to check for + // this on the highest level only + for (const auto &cell : triangulation.active_cell_iterators_on_level( + triangulation.levels.size() - 1)) + if (cell->refine_flag_set()) + { + triangulation.levels.push_back( + std::make_unique< + internal::TriangulationImplementation::TriaLevel>(dim)); + break; + } + // first clear user flags for quads and lines; we're going to + // use them to flag which lines and quads need refinement triangulation.faces->quads.clear_user_data(); for (typename Triangulation::line_iterator line = @@ -5131,6 +5129,18 @@ namespace internal ++quad) quad->clear_user_flag(); + // check how much space is needed on every level we need not + // check the highest level since either + // - on the highest level no cells are flagged for refinement + // - there are, but prepare_refinement added another empty + // level which then is the highest level + + // variables to hold the number of newly to be created + // vertices, lines and quads. as these are stored globally, + // declare them outside the loop over al levels. we need lines + // and quads in pairs for refinement of old ones and lines and + // quads, that can be stored as single ones, as they are newly + // created in the inside of an existing cell unsigned int needed_vertices = 0; unsigned int needed_lines_single = 0; unsigned int needed_quads_single = 0; @@ -5140,22 +5150,22 @@ namespace internal { unsigned int new_cells = 0; - for (const auto &acell : + for (const auto &cell : triangulation.active_cell_iterators_on_level(level)) - if (acell->refine_flag_set()) + if (cell->refine_flag_set()) { - Assert(acell->refine_flag_set() == + Assert(cell->refine_flag_set() == RefinementCase::cut_xyz, ExcInternalError()); - if (acell->reference_cell() == ReferenceCells::Hexahedron) + if (cell->reference_cell() == ReferenceCells::Hexahedron) { ++needed_vertices; needed_lines_single += 6; needed_quads_single += 12; new_cells += 8; } - else if (acell->reference_cell() == + else if (cell->reference_cell() == ReferenceCells::Tetrahedron) { needed_lines_single += 1; @@ -5167,13 +5177,13 @@ namespace internal Assert(false, ExcInternalError()); } - for (const auto face : acell->face_indices()) - if (acell->face(face)->number_of_children() < 4) - acell->face(face)->set_user_flag(); + for (const auto face : cell->face_indices()) + if (cell->face(face)->number_of_children() < 4) + cell->face(face)->set_user_flag(); - for (const auto line : acell->line_indices()) - if (acell->line(line)->has_children() == false) - acell->line(line)->set_user_flag(); + for (const auto line : cell->line_indices()) + if (cell->line(line)->has_children() == false) + cell->line(line)->set_user_flag(); } const unsigned int used_cells = @@ -5189,6 +5199,8 @@ namespace internal reserve_space(triangulation.levels[level + 1]->cells, new_cells); } + // now count the quads and lines which were flagged for + // refinement for (typename Triangulation::quad_iterator quad = triangulation.begin_quad(); quad != triangulation.end_quad(); @@ -5237,6 +5249,7 @@ namespace internal needed_quads_single); + // add to needed vertices how many vertices are already in use needed_vertices += std::count(triangulation.vertices_used.begin(), triangulation.vertices_used.end(), true); @@ -5247,6 +5260,26 @@ namespace internal triangulation.vertices_used.resize(needed_vertices, false); } + /////////////////////////////////////////// + // Before we start with the actual refinement, we do some + // sanity checks if in debug mode. especially, we try to catch + // the notorious problem with lines being twice refined, + // i.e. there are cells adjacent at one line ("around the + // edge", but not at a face), with two cells differing by more + // than one refinement level + // + // this check is very simple to implement here, since we have + // all lines flagged if they shall be refined +#ifdef DEBUG + for (const auto &cell : triangulation.active_cell_iterators()) + if (!cell->refine_flag_set()) + for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n) + if (cell->line(line_n)->has_children()) + for (unsigned int c = 0; c < 2; ++c) + Assert(cell->line(line_n)->child(c)->user_flag_set() == false, + ExcInternalError()); +#endif + unsigned int current_vertex = 0; // helper function - find the next available vertex number and mark it @@ -5287,6 +5320,9 @@ namespace internal Assert(next_unused_line.state() == IteratorState::valid, ExcInternalError()); + // now we found two consecutive unused lines, such + // that the children of a line will be consecutive. + // then set the child pointer of the present line line->set_children(0, next_unused_line->index()); const typename Triangulation::raw_line_iterator @@ -5378,7 +5414,8 @@ namespace internal } } - // 3) create new quads (properties are set below) + // 3) create new quads (properties are set below). Both triangles + // and quads are divided in four. std::array< typename Triangulation::raw_quad_iterator, 4> @@ -5410,11 +5447,10 @@ namespace internal quad->set_refinement_case(RefinementCase<2>::cut_xy); } - std::array vertex_indices; + // Maximum of 9 vertices per refined quad (9 for Quadrilateral, 6 + // for Triangle) + std::array vertex_indices = {}; { - for (auto &i : vertex_indices) - i = 0.0; - unsigned int k = 0; for (const auto i : quad->vertex_indices()) vertex_indices[k++] = quad->vertex_index(i); @@ -5440,10 +5476,10 @@ namespace internal { static constexpr std::array, 2> - index = { - {{{1, 0}}, // child 0, line_orientation=false and true - {{0, - 1}}}}; // child 1, line_orientation=false and true + index = {// child 0, line_orientation=false and true + {{{1, 0}}, + // child 1, line_orientation=false and true + {{0, 1}}}}; lines[k++] = quad->line(l)->child( index[c][quad->line_orientation(l)]); @@ -5670,8 +5706,9 @@ namespace internal hex->center(true, true); } - std::vector< - typename Triangulation::raw_line_iterator> + boost::container::small_vector< + typename Triangulation::raw_line_iterator, + 6> new_lines(n_new_lines); for (unsigned int i = 0; i < n_new_lines; ++i) { @@ -5689,8 +5726,9 @@ namespace internal new_lines[i]->set_manifold_id(hex->manifold_id()); } - std::vector< - typename Triangulation::raw_quad_iterator> + boost::container::small_vector< + typename Triangulation::raw_quad_iterator, + 12> new_quads(n_new_quads); for (unsigned int i = 0; i < n_new_quads; ++i) { @@ -5720,9 +5758,11 @@ namespace internal new_quad->set_line_orientation(j, true); } - std::vector< - typename Triangulation::raw_hex_iterator> - new_hexes(n_new_hexes); + // we always get 8 children per refined cell + std::array< + typename Triangulation::raw_hex_iterator, + 8> + new_hexes; { for (unsigned int i = 0; i < n_new_hexes; ++i) { @@ -5753,6 +5793,17 @@ namespace internal if (i % 2) new_hex->set_parent(hex->index()); + // set the face_orientation flag to true for all + // faces initially, as this is the default value + // which is true for all faces interior to the + // hex. later on go the other way round and + // reset faces that are at the boundary of the + // mother cube + // + // the same is true for the face_flip and + // face_rotation flags. however, the latter two + // are set to false by default as this is the + // standard value for (const auto f : new_hex->face_indices()) { new_hex->set_face_orientation(f, true); @@ -5760,16 +5811,13 @@ namespace internal new_hex->set_face_rotation(f, false); } } - // note these hexes as children to the present cell for (unsigned int i = 0; i < n_new_hexes / 2; ++i) hex->set_children(2 * i, new_hexes[2 * i]->index()); } { // load vertex indices - std::array vertex_indices; - for (auto &i : vertex_indices) - i = 0; + std::array vertex_indices = {}; { unsigned int k = 0; -- 2.39.5