From cbd85a1089c25636acd2cfe9a5396ae17063fa72 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 29 Jun 2022 23:08:53 +0200 Subject: [PATCH] Restructure execture_refinement(Triangulation<3, spacedim>) --- source/grid/tria.cc | 762 +++++++++++++++++++++----------------------- 1 file changed, 358 insertions(+), 404 deletions(-) diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 5fc10e7bf5..91f08c776e 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -5048,6 +5048,10 @@ namespace internal { static const int dim = 3; static const unsigned int X = numbers::invalid_unsigned_int; + using raw_line_iterator = + typename Triangulation::raw_line_iterator; + using raw_quad_iterator = + typename Triangulation::raw_quad_iterator; Assert(spacedim == 3, ExcNotImplemented()); @@ -5271,19 +5275,13 @@ namespace internal typename Triangulation::active_line_iterator line = triangulation.begin_active_line(), endl = triangulation.end_line(); - typename Triangulation::raw_line_iterator - next_unused_line = triangulation.begin_raw_line(); + raw_line_iterator next_unused_line = triangulation.begin_raw_line(); for (; line != endl; ++line) { if (line->user_flag_set() == false) continue; - current_vertex = - get_next_unused_vertex(current_vertex, - triangulation.vertices_used); - triangulation.vertices[current_vertex] = line->center(true); - next_unused_line = triangulation.faces->lines.template next_free_pair_object<1>( triangulation); @@ -5295,31 +5293,33 @@ namespace internal // then set the child pointer of the present line line->set_children(0, next_unused_line->index()); - const typename Triangulation::raw_line_iterator - children[2] = {next_unused_line, ++next_unused_line}; + const std::array children{ + {next_unused_line, ++next_unused_line}}; AssertIsNotUsed(children[0]); AssertIsNotUsed(children[1]); + current_vertex = + get_next_unused_vertex(current_vertex, + triangulation.vertices_used); + triangulation.vertices[current_vertex] = line->center(true); + children[0]->set_bounding_object_indices( {line->vertex_index(0), current_vertex}); children[1]->set_bounding_object_indices( {current_vertex, line->vertex_index(1)}); - children[0]->set_used_flag(); - children[1]->set_used_flag(); - children[0]->clear_children(); - children[1]->clear_children(); - children[0]->clear_user_data(); - children[1]->clear_user_data(); - children[0]->clear_user_flag(); - children[1]->clear_user_flag(); - - children[0]->set_boundary_id_internal(line->boundary_id()); - children[1]->set_boundary_id_internal(line->boundary_id()); - - children[0]->set_manifold_id(line->manifold_id()); - children[1]->set_manifold_id(line->manifold_id()); + const auto manifold_id = line->manifold_id(); + const auto boundary_id = line->boundary_id(); + for (const auto &child : children) + { + child->set_used_flag(); + child->clear_children(); + child->clear_user_data(); + child->clear_user_flag(); + child->set_boundary_id_internal(boundary_id); + child->set_manifold_id(manifold_id); + } line->clear_user_flag(); } @@ -5330,10 +5330,6 @@ namespace internal typename Triangulation::quad_iterator quad = triangulation.begin_quad(), endq = triangulation.end_quad(); - typename Triangulation::raw_line_iterator - next_unused_line = triangulation.begin_raw_line(); - typename Triangulation::raw_quad_iterator - next_unused_quad = triangulation.begin_raw_quad(); for (; quad != endq; ++quad) { @@ -5342,129 +5338,107 @@ namespace internal const auto reference_face_type = quad->reference_cell(); - // 1) create new vertex (at the center of the face) + // 1) create new lines (property is set later) + // maximum of 4 new lines (4 quadrilateral, 3 triangle) + std::array new_lines; if (reference_face_type == ReferenceCells::Quadrilateral) { - current_vertex = - get_next_unused_vertex(current_vertex, - triangulation.vertices_used); - triangulation.vertices[current_vertex] = - quad->center(true, true); + for (unsigned int l = 0; l < 2; ++l) + { + auto next_unused_line = + triangulation.faces->lines + .template next_free_pair_object<1>(triangulation); + new_lines[2 * l] = next_unused_line; + new_lines[2 * l + 1] = ++next_unused_line; + } + } + else if (reference_face_type == ReferenceCells::Triangle) + { + for (unsigned int l = 0; l < 3; ++l) + new_lines[l] = + triangulation.faces->lines + .template next_free_single_object<1>(triangulation); + } + else + { + Assert(false, ExcNotImplemented()); } - // 2) create new lines (property is set later) - boost::container::small_vector< - typename Triangulation::raw_line_iterator, - GeometryInfo::lines_per_cell> - new_lines(quad->n_lines()); - { - for (unsigned int i = 0; i < new_lines.size(); ++i) - { - if (reference_face_type == ReferenceCells::Quadrilateral) - { - if (i % 2 == 0) - next_unused_line = - triangulation.faces->lines - .template next_free_pair_object<1>(triangulation); - } - else if (reference_face_type == ReferenceCells::Triangle) - { - next_unused_line = - triangulation.faces->lines - .template next_free_single_object<1>(triangulation); - } - else - { - Assert(false, ExcNotImplemented()); - } - - new_lines[i] = next_unused_line; - ++next_unused_line; - AssertIsNotUsed(new_lines[i]); - } - } + for (const unsigned int line : quad->line_indices()) + { + AssertIsNotUsed(new_lines[line]); + (void)line; + } - // 3) create new quads (properties are set below). Both triangles + // 2) create new quads (properties are set below). Both triangles // and quads are divided in four. - std::array< - typename Triangulation::raw_quad_iterator, - 4> - new_quads; - { - next_unused_quad = - triangulation.faces->quads.template next_free_pair_object<2>( - triangulation); - - new_quads[0] = next_unused_quad; - AssertIsNotUsed(new_quads[0]); + std::array new_quads; + for (unsigned int q = 0; q < 2; ++q) + { + auto next_unused_quad = + triangulation.faces->quads + .template next_free_pair_object<2>(triangulation); - ++next_unused_quad; - new_quads[1] = next_unused_quad; - AssertIsNotUsed(new_quads[1]); + new_quads[2 * q] = next_unused_quad; + new_quads[2 * q + 1] = ++next_unused_quad; - next_unused_quad = - triangulation.faces->quads.template next_free_pair_object<2>( - triangulation); - new_quads[2] = next_unused_quad; - AssertIsNotUsed(new_quads[2]); + quad->set_children(2 * q, new_quads[2 * q]->index()); + } + quad->set_refinement_case(RefinementCase<2>::cut_xy); - ++next_unused_quad; - new_quads[3] = next_unused_quad; - AssertIsNotUsed(new_quads[3]); + for (const auto &quad : new_quads) + { + AssertIsNotUsed(quad); + (void)quad; + } - quad->set_children(0, new_quads[0]->index()); - quad->set_children(2, new_quads[2]->index()); - quad->set_refinement_case(RefinementCase<2>::cut_xy); - } + // 3) set vertex indices and set new vertex // Maximum of 9 vertices per refined quad (9 for Quadrilateral, 6 // for Triangle) std::array vertex_indices = {}; - { - unsigned int k = 0; - for (const auto i : quad->vertex_indices()) - vertex_indices[k++] = quad->vertex_index(i); + unsigned int k = 0; + for (const auto i : quad->vertex_indices()) + vertex_indices[k++] = quad->vertex_index(i); - for (const auto i : quad->line_indices()) - vertex_indices[k++] = - quad->line(i)->child(0)->vertex_index(1); + for (const auto i : quad->line_indices()) + vertex_indices[k++] = quad->line(i)->child(0)->vertex_index(1); - vertex_indices[k++] = current_vertex; - } + if (reference_face_type == ReferenceCells::Quadrilateral) + { + current_vertex = + get_next_unused_vertex(current_vertex, + triangulation.vertices_used); + vertex_indices[k++] = current_vertex; - boost::container::small_vector< - typename Triangulation::raw_line_iterator, - 12> - lines(reference_face_type == ReferenceCells::Quadrilateral ? - 12 : - 9); - { - unsigned int k = 0; + triangulation.vertices[current_vertex] = + quad->center(true, true); + } - for (unsigned int l = 0; l < quad->n_lines(); ++l) - for (unsigned int c = 0; c < 2; ++c) - { - static constexpr std::array, - 2> - 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)]); - } + // 4) set new lines on quads and their properties + std::array lines; + unsigned int n_lines = 0; + for (unsigned int l = 0; l < quad->n_lines(); ++l) + for (unsigned int c = 0; c < 2; ++c) + { + static constexpr dealii::ndarray index = + {{// child 0, line_orientation=false and true + {{1, 0}}, + // child 1, line_orientation=false and true + {{0, 1}}}}; + + lines[n_lines++] = + quad->line(l)->child(index[c][quad->line_orientation(l)]); + } - for (unsigned int l = 0; l < new_lines.size(); ++l) - lines[k++] = new_lines[l]; - } + for (unsigned int l = 0; l < quad->n_lines(); ++l) + lines[n_lines++] = new_lines[l]; - boost::container::small_vector line_indices( - lines.size()); - for (unsigned int i = 0; i < line_indices.size(); ++i) + std::array line_indices; + for (unsigned int i = 0; i < n_lines; ++i) line_indices[i] = lines[i]->index(); - static constexpr std::array, 12> + static constexpr dealii::ndarray line_vertices_quad{{{{0, 4}}, {{4, 2}}, {{1, 5}}, @@ -5478,21 +5452,13 @@ namespace internal {{4, 8}}, {{8, 5}}}}; - static constexpr std::array, 4> + static constexpr dealii::ndarray quad_lines_quad{{{{0, 8, 4, 10}}, {{8, 2, 5, 11}}, {{1, 9, 10, 6}}, {{9, 3, 11, 7}}}}; - static constexpr std:: - array, 4>, 4> - quad_line_vertices_quad{ - {{{{{0, 4}}, {{6, 8}}, {{0, 6}}, {{4, 8}}}}, - {{{{6, 8}}, {{1, 5}}, {{6, 1}}, {{8, 5}}}}, - {{{{4, 2}}, {{8, 7}}, {{4, 8}}, {{2, 7}}}}, - {{{{8, 7}}, {{5, 3}}, {{8, 5}}, {{7, 3}}}}}}; - - static constexpr std::array, 12> + static constexpr dealii::ndarray line_vertices_tri{{{{0, 3}}, {{3, 1}}, {{1, 4}}, @@ -5506,19 +5472,18 @@ namespace internal {{X, X}}, {{X, X}}}}; - static constexpr std::array, 4> + static constexpr dealii::ndarray quad_lines_tri{{{{0, 8, 5, X}}, {{1, 2, 6, X}}, {{7, 3, 4, X}}, {{6, 7, 8, X}}}}; - static constexpr std:: - array, 4>, 4> - quad_line_vertices_tri{ - {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}}, - {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}}, - {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}}, - {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}}; + static constexpr dealii::ndarray + quad_line_vertices_tri{ + {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}}, + {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}}, + {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}}, + {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}}; const auto &line_vertices = (reference_face_type == ReferenceCells::Quadrilateral) ? @@ -5528,14 +5493,9 @@ namespace internal (reference_face_type == ReferenceCells::Quadrilateral) ? quad_lines_quad : quad_lines_tri; - const auto &quad_line_vertices = - (reference_face_type == ReferenceCells::Quadrilateral) ? - quad_line_vertices_quad : - quad_line_vertices_tri; - // 4) set properties of lines - for (unsigned int i = 0, j = lines.size() - new_lines.size(); - i < new_lines.size(); + for (unsigned int i = 0, j = 2 * quad->n_lines(); + i < quad->n_lines(); ++i, ++j) { auto &new_line = new_lines[i]; @@ -5560,12 +5520,12 @@ namespace internal triangulation.faces->quad_reference_cell[new_quad->index()] = reference_face_type; - if (new_quad->n_lines() == 3) + if (reference_face_type == ReferenceCells::Triangle) new_quad->set_bounding_object_indices( {line_indices[quad_lines[i][0]], line_indices[quad_lines[i][1]], line_indices[quad_lines[i][2]]}); - else if (new_quad->n_lines() == 4) + else if (reference_face_type == ReferenceCells::Quadrilateral) new_quad->set_bounding_object_indices( {line_indices[quad_lines[i][0]], line_indices[quad_lines[i][1]], @@ -5585,38 +5545,54 @@ namespace internal std::set s; #endif - // ... and fix orientation of faces (lines) of quad - for (const auto f : new_quad->line_indices()) + // ... and fix orientation of lines of face for triangles, + // using an expensive algorithm, quadrilaterals are treated + // a few lines below by a cheaper algorithm + if (reference_face_type == ReferenceCells::Triangle) { - std::array vertices_0, vertices_1; + for (const auto f : new_quad->line_indices()) + { + std::array vertices_0, vertices_1; - for (unsigned int v = 0; v < 2; ++v) - vertices_0[v] = - lines[quad_lines[i][f]]->vertex_index(v); + for (unsigned int v = 0; v < 2; ++v) + vertices_0[v] = + lines[quad_lines[i][f]]->vertex_index(v); - for (unsigned int v = 0; v < 2; ++v) - vertices_1[v] = - vertex_indices[quad_line_vertices[i][f][v]]; + for (unsigned int v = 0; v < 2; ++v) + vertices_1[v] = + vertex_indices[quad_line_vertices_tri[i][f][v]]; - const auto orientation = - ReferenceCells::Line.compute_orientation(vertices_0, - vertices_1); + const auto orientation = + ReferenceCells::Line.compute_orientation( + vertices_0, vertices_1); #ifdef DEBUG - for (const auto i : vertices_0) - s.insert(i); - for (const auto i : vertices_1) - s.insert(i); + for (const auto i : vertices_0) + s.insert(i); + for (const auto i : vertices_1) + s.insert(i); #endif - new_quad->set_line_orientation(f, orientation); - } + new_quad->set_line_orientation(f, orientation); + } #ifdef DEBUG - AssertDimension( - s.size(), - (reference_face_type == ReferenceCells::Quadrilateral ? 4 : - 3)); + AssertDimension(s.size(), 3); #endif + } + } + + // fix orientation of lines of faces for quadrilaterals with + // cheap algorithm + if (reference_face_type == ReferenceCells::Quadrilateral) + { + static constexpr dealii::ndarray + quad_child_boundary_lines{ + {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}}; + + for (unsigned int i = 0; i < 4; ++i) + for (unsigned int j = 0; j < 2; ++j) + new_quads[quad_child_boundary_lines[i][j]] + ->set_line_orientation(i, quad->line_orientation(i)); } quad->clear_user_flag(); @@ -5626,16 +5602,20 @@ namespace internal typename Triangulation<3, spacedim>::DistortedCellList cells_with_distorted_children; + typename Triangulation::active_hex_iterator hex = + triangulation.begin_active_hex(0); for (unsigned int level = 0; level != triangulation.levels.size() - 1; ++level) { - typename Triangulation::active_hex_iterator - hex = triangulation.begin_active_hex(level), - endh = triangulation.begin_active_hex(level + 1); typename Triangulation::raw_hex_iterator next_unused_hex = triangulation.begin_raw_hex(level + 1); + Assert(hex == triangulation.end() || + hex->level() >= static_cast(level), + ExcInternalError()); - for (; hex != endh; ++hex) + for (; hex != triangulation.end() && + hex->level() == static_cast(level); + ++hex) { if (hex->refine_flag_set() == RefinementCase::no_refinement) @@ -5666,20 +5646,7 @@ namespace internal else Assert(false, ExcNotImplemented()); - // Hexes add a single new internal vertex - if (reference_cell_type == ReferenceCells::Hexahedron) - { - current_vertex = - get_next_unused_vertex(current_vertex, - triangulation.vertices_used); - triangulation.vertices[current_vertex] = - hex->center(true, true); - } - - boost::container::small_vector< - typename Triangulation::raw_line_iterator, - 6> - new_lines(n_new_lines); + std::array new_lines; for (unsigned int i = 0; i < n_new_lines; ++i) { new_lines[i] = @@ -5696,10 +5663,7 @@ namespace internal new_lines[i]->set_manifold_id(hex->manifold_id()); } - boost::container::small_vector< - typename Triangulation::raw_quad_iterator, - 12> - new_quads(n_new_quads); + std::array new_quads; for (unsigned int i = 0; i < n_new_quads; ++i) { new_quads[i] = @@ -5747,7 +5711,7 @@ namespace internal auto &new_hex = new_hexes[i]; - // TODO: children have the same type as the parent + // children have the same type as the parent triangulation.levels[new_hex->level()] ->reference_cell[new_hex->index()] = reference_cell_type; @@ -5763,6 +5727,7 @@ 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 @@ -5803,15 +5768,22 @@ namespace internal { for (const unsigned int i : hex->face_indices()) vertex_indices[k++] = - middle_vertex_index(hex->face(i)); + hex->face(i)->child(0)->vertex_index(3); + // Set single new vertex in the center + current_vertex = + get_next_unused_vertex(current_vertex, + triangulation.vertices_used); vertex_indices[k++] = current_vertex; + + triangulation.vertices[current_vertex] = + hex->center(true, true); } } // set up new lines { - static constexpr std::array, 6> + static constexpr dealii::ndarray new_line_vertices_hex = {{{{22, 26}}, {{26, 23}}, {{20, 26}}, @@ -5819,7 +5791,7 @@ namespace internal {{24, 26}}, {{26, 25}}}}; - static constexpr std::array, 6> + static constexpr dealii::ndarray new_line_vertices_tet = {{{{6, 8}}, {{X, X}}, {{X, X}}, @@ -5832,7 +5804,7 @@ namespace internal new_line_vertices_hex : new_line_vertices_tet; - for (unsigned int i = 0; i < new_lines.size(); ++i) + for (unsigned int i = 0; i < n_new_lines; ++i) new_lines[i]->set_bounding_object_indices( {vertex_indices[new_line_vertices[i][0]], vertex_indices[new_line_vertices[i][1]]}); @@ -5840,10 +5812,8 @@ namespace internal // set up new quads { - boost::container::small_vector< - typename Triangulation::raw_line_iterator, - 30> - relevant_lines(0); + boost::container::small_vector + relevant_lines; if (reference_cell_type == ReferenceCells::Hexahedron) { @@ -5851,8 +5821,8 @@ namespace internal for (unsigned int f = 0, k = 0; f < 6; ++f) for (unsigned int c = 0; c < 4; ++c, ++k) { - static constexpr std:: - array, 4> + static constexpr dealii:: + ndarray temp = { {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}}; @@ -5919,7 +5889,7 @@ namespace internal ++i) relevant_line_indices[i] = relevant_lines[i]->index(); - static constexpr std::array, 12> + static constexpr dealii::ndarray new_quad_lines_hex = {{{{10, 28, 16, 24}}, {{28, 14, 17, 25}}, {{11, 29, 24, 20}}, @@ -5933,7 +5903,7 @@ namespace internal {{3, 25, 26, 12}}, {{25, 7, 27, 13}}}}; - static constexpr std::array, 12> + static constexpr dealii::ndarray new_quad_lines_tet = {{{{2, 3, 8, X}}, {{0, 9, 5, X}}, {{1, 6, 11, X}}, @@ -5947,37 +5917,35 @@ namespace internal {{X, X, X, X}}, {{X, X, X, X}}}}; - static constexpr std:: - array, 4>, 12> - table_hex = { - {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}}, - {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}}, - {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}}, - {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}}, - {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}}, - {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}}, - {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}}, - {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}}, - {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}}, - {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}}, - {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}}, - {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}}; - - static constexpr std:: - array, 4>, 12> - table_tet = { - {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}}, - {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}}, - {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}}, - {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}}, - {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}}, - {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}}, - {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}}, - {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}}, - {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, - {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, - {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, - {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}}; + static constexpr dealii::ndarray + table_hex = { + {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}}, + {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}}, + {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}}, + {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}}, + {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}}, + {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}}, + {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}}, + {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}}, + {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}}, + {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}}, + {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}}, + {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}}; + + static constexpr dealii::ndarray + table_tet = { + {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}}, + {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}}, + {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}}, + {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}}, + {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}}, + {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}}, + {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}}, + {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}}, + {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, + {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, + {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, + {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}}; const auto &new_quad_lines = (reference_cell_type == ReferenceCells::Hexahedron) ? @@ -5989,23 +5957,11 @@ namespace internal table_hex : table_tet; - for (unsigned int q = 0; q < new_quads.size(); ++q) - { - for (unsigned int l = 0; l < 3; ++l) - { - std::array vertices_0, vertices_1; + static constexpr dealii::ndarray + representative_lines{ + {{{0, 2}}, {{2, 0}}, {{3, 3}}, {{1, 1}}}}; - for (unsigned int v = 0; v < 2; ++v) - vertices_0[v] = - relevant_lines[new_quad_lines[q][l]] - ->vertex_index(v); - - for (unsigned int v = 0; v < 2; ++v) - vertices_1[v] = vertex_indices[table[q][l][v]]; - } - } - - for (unsigned int q = 0; q < new_quads.size(); ++q) + for (unsigned int q = 0; q < n_new_quads; ++q) { auto &new_quad = new_quads[q]; @@ -6023,8 +5979,22 @@ namespace internal else Assert(false, ExcNotImplemented()); - for (const auto l : new_quad->line_indices()) + // On hexes, we must only determine a single line + // according to the representative_lines array above + // (this saves expensive operations), for tets we do + // all lines manually + const unsigned int n_compute_lines = + reference_cell_type == ReferenceCells::Hexahedron ? + 1 : + new_quad->n_lines(); + for (unsigned int line = 0; line < n_compute_lines; + ++line) { + const unsigned int l = + (reference_cell_type == + ReferenceCells::Hexahedron) ? + representative_lines[q % 4][0] : + line; std::array vertices_0, vertices_1; for (unsigned int v = 0; v < 2; ++v) @@ -6040,6 +6010,15 @@ namespace internal vertices_0, vertices_1); new_quad->set_line_orientation(l, orientation); + + // on a hex, inject the status of the current line + // also to the line on the other quad along the + // same direction + if (reference_cell_type == + ReferenceCells::Hexahedron) + new_quads[representative_lines[q % 4][1] + q - + (q % 4)] + ->set_line_orientation(l, orientation); } } } @@ -6050,11 +6029,10 @@ namespace internal if (reference_cell_type == ReferenceCells::Hexahedron) { - for (unsigned int i = 0; i < new_quads.size(); ++i) + for (unsigned int i = 0; i < n_new_quads; ++i) quad_indices[i] = new_quads[i]->index(); - for (unsigned int f = 0, k = new_quads.size(); f < 6; - ++f) + for (unsigned int f = 0, k = n_new_quads; f < 6; ++f) for (unsigned int c = 0; c < 4; ++c, ++k) quad_indices[k] = hex->face(f)->isotropic_child_index( @@ -6066,11 +6044,10 @@ namespace internal } else if (reference_cell_type == ReferenceCells::Tetrahedron) { - for (unsigned int i = 0; i < new_quads.size(); ++i) + for (unsigned int i = 0; i < n_new_quads; ++i) quad_indices[i] = new_quads[i]->index(); - for (unsigned int f = 0, k = new_quads.size(); f < 4; - ++f) + for (unsigned int f = 0, k = n_new_quads; f < 4; ++f) for (unsigned int c = 0; c < 4; ++c, ++k) { quad_indices[k] = hex->face(f)->child_index( @@ -6092,7 +6069,7 @@ namespace internal Assert(false, ExcNotImplemented()); } - static constexpr std::array, 8> + static constexpr dealii::ndarray cell_quads_hex = {{ {{12, 0, 20, 4, 28, 8}}, // bottom children {{0, 16, 22, 6, 29, 9}}, // @@ -6104,7 +6081,7 @@ namespace internal {{3, 19, 7, 27, 11, 35}} // }}; - static constexpr std::array, 8> + static constexpr dealii::ndarray cell_quads_tet{{{{8, 13, 16, 0, X, X}}, {{9, 12, 1, 21, X, X}}, {{10, 2, 17, 20, X, X}}, @@ -6114,118 +6091,61 @@ namespace internal {{19, 7, 6, 3, X, X}}, {{23, 5, 2, 7, X, X}}}}; - static constexpr std:: - array, 6>, 8> - cell_face_vertices_hex{{{{{{0, 8, 16, 20}}, - {{10, 24, 22, 26}}, - {{0, 16, 10, 22}}, - {{8, 20, 24, 26}}, - {{0, 10, 8, 24}}, - {{16, 22, 20, 26}}}}, - {{{{10, 24, 22, 26}}, - {{1, 9, 17, 21}}, - {{10, 22, 1, 17}}, - {{24, 26, 9, 21}}, - {{10, 1, 24, 9}}, - {{22, 17, 26, 21}}}}, - {{{{8, 2, 20, 18}}, - {{24, 11, 26, 23}}, - {{8, 20, 24, 26}}, - {{2, 18, 11, 23}}, - {{8, 24, 2, 11}}, - {{20, 26, 18, 23}}}}, - {{{{24, 11, 26, 23}}, - {{9, 3, 21, 19}}, - {{24, 26, 9, 21}}, - {{11, 23, 3, 19}}, - {{24, 9, 11, 3}}, - {{26, 21, 23, 19}}}}, - {{{{16, 20, 4, 12}}, - {{22, 26, 14, 25}}, - {{16, 4, 22, 14}}, - {{20, 12, 26, 25}}, - {{16, 22, 20, 26}}, - {{4, 14, 12, 25}}}}, - {{{{22, 26, 14, 25}}, - {{17, 21, 5, 13}}, - {{22, 14, 17, 5}}, - {{26, 25, 21, 13}}, - {{22, 17, 26, 21}}, - {{14, 5, 25, 13}}}}, - {{{{20, 18, 12, 6}}, - {{26, 23, 25, 15}}, - {{20, 12, 26, 25}}, - {{18, 6, 23, 15}}, - {{20, 26, 18, 23}}, - {{12, 25, 6, 15}}}}, - {{{{26, 23, 25, 15}}, - {{21, 19, 13, 7}}, - {{26, 25, 21, 13}}, - {{23, 15, 19, 7}}, - {{26, 21, 23, 19}}, - {{25, 13, 15, 7}}}}}}; - - static constexpr std:: - array, 6>, 8> - cell_face_vertices_tet{{{{{{0, 4, 6, X}}, - {{4, 0, 7, X}}, - {{0, 6, 7, X}}, - {{6, 4, 7, X}}, - {{X, X, X, X}}, - {{X, X, X, X}}}}, - {{{{4, 1, 5, X}}, - {{1, 4, 8, X}}, - {{4, 5, 8, X}}, - {{5, 1, 8, X}}, - {{X, X, X, X}}, - {{X, X, X, X}}}}, - {{{{6, 5, 2, X}}, - {{5, 6, 9, X}}, - {{6, 2, 9, X}}, - {{2, 5, 9, X}}, - {{X, X, X, X}}, - {{X, X, X, X}}}}, - {{{{7, 8, 9, X}}, - {{8, 7, 3, X}}, - {{7, 9, 3, X}}, - {{9, 8, 3, X}}, - {{X, X, X, X}}, - {{X, X, X, X}}}}, - {{{{4, 5, 6, X}}, - {{5, 4, 8, X}}, - {{4, 6, 8, X}}, - {{6, 5, 8, X}}, - {{X, X, X, X}}, - {{X, X, X, X}}}}, - {{{{4, 7, 8, X}}, - {{7, 4, 6, X}}, - {{4, 8, 6, X}}, - {{8, 7, 6, X}}, - {{X, X, X, X}}, - {{X, X, X, X}}}}, - {{{{6, 9, 7, X}}, - {{9, 6, 8, X}}, - {{6, 7, 8, X}}, - {{7, 9, 8, X}}, - {{X, X, X, X}}, - {{X, X, X, X}}}}, - {{{{5, 8, 9, X}}, - {{8, 5, 6, X}}, - {{5, 9, 6, X}}, - {{9, 8, 6, X}}, - {{X, X, X, X}}, - {{X, X, X, X}}}}}}; + static constexpr dealii::ndarray + cell_face_vertices_tet{{{{{{0, 4, 6, X}}, + {{4, 0, 7, X}}, + {{0, 6, 7, X}}, + {{6, 4, 7, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{4, 1, 5, X}}, + {{1, 4, 8, X}}, + {{4, 5, 8, X}}, + {{5, 1, 8, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{6, 5, 2, X}}, + {{5, 6, 9, X}}, + {{6, 2, 9, X}}, + {{2, 5, 9, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{7, 8, 9, X}}, + {{8, 7, 3, X}}, + {{7, 9, 3, X}}, + {{9, 8, 3, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{4, 5, 6, X}}, + {{5, 4, 8, X}}, + {{4, 6, 8, X}}, + {{6, 5, 8, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{4, 7, 8, X}}, + {{7, 4, 6, X}}, + {{4, 8, 6, X}}, + {{8, 7, 6, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{6, 9, 7, X}}, + {{9, 6, 8, X}}, + {{6, 7, 8, X}}, + {{7, 9, 8, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{5, 8, 9, X}}, + {{8, 5, 6, X}}, + {{5, 9, 6, X}}, + {{9, 8, 6, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}}}; const auto &cell_quads = (reference_cell_type == ReferenceCells::Hexahedron) ? cell_quads_hex : cell_quads_tet; - const auto &cell_face_vertices = - (reference_cell_type == ReferenceCells::Hexahedron) ? - cell_face_vertices_hex : - cell_face_vertices_tet; - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; ++c) @@ -6233,11 +6153,42 @@ namespace internal auto &new_hex = new_hexes[c]; if (new_hex->n_faces() == 4) - new_hex->set_bounding_object_indices( - {quad_indices[cell_quads[c][0]], - quad_indices[cell_quads[c][1]], - quad_indices[cell_quads[c][2]], - quad_indices[cell_quads[c][3]]}); + { + new_hex->set_bounding_object_indices( + {quad_indices[cell_quads[c][0]], + quad_indices[cell_quads[c][1]], + quad_indices[cell_quads[c][2]], + quad_indices[cell_quads[c][3]]}); + + // for tets, we need to go through the faces and + // figure the orientation out the hard way + for (const auto f : new_hex->face_indices()) + { + std::array vertices_0, + vertices_1; + + const auto &face = new_hex->face(f); + + for (const auto i : face->vertex_indices()) + vertices_0[i] = face->vertex_index(i); + + for (const auto i : face->vertex_indices()) + vertices_1[i] = + vertex_indices[cell_face_vertices_tet[c][f] + [i]]; + + const auto orientation = + face->reference_cell().compute_orientation( + vertices_1, vertices_0); + + new_hex->set_face_orientation( + f, Utilities::get_bit(orientation, 0)); + new_hex->set_face_flip( + f, Utilities::get_bit(orientation, 2)); + new_hex->set_face_rotation( + f, Utilities::get_bit(orientation, 1)); + } + } else if (new_hex->n_faces() == 6) new_hex->set_bounding_object_indices( {quad_indices[cell_quads[c][0]], @@ -6248,32 +6199,35 @@ namespace internal quad_indices[cell_quads[c][5]]}); else Assert(false, ExcNotImplemented()); - - for (const auto f : new_hex->face_indices()) - { - std::array vertices_0, vertices_1; - - const auto &face = new_hex->face(f); - - for (const auto i : face->vertex_indices()) - vertices_0[i] = face->vertex_index(i); - - for (const auto i : face->vertex_indices()) - vertices_1[i] = - vertex_indices[cell_face_vertices[c][f][i]]; - - const auto orientation = - face->reference_cell().compute_orientation( - vertices_1, vertices_0); - - new_hex->set_face_orientation( - f, Utilities::get_bit(orientation, 0)); - new_hex->set_face_flip( - f, Utilities::get_bit(orientation, 2)); - new_hex->set_face_rotation( - f, Utilities::get_bit(orientation, 1)); - } } + + // for hexes, we can simply inherit the orientation values + // from the parent on the outer faces; the inner faces can + // be skipped as their orientation is always the default + // one set above + static constexpr dealii::ndarray + face_to_child_indices_hex{{{{0, 2, 4, 6}}, + {{1, 3, 5, 7}}, + {{0, 1, 4, 5}}, + {{2, 3, 6, 7}}, + {{0, 1, 2, 3}}, + {{4, 5, 6, 7}}}}; + if (hex->n_faces() == 6) + for (const auto f : hex->face_indices()) + { + const unsigned int orientation = + hex->face_orientation(f); + const unsigned int flip = hex->face_flip(f); + const unsigned int rotation = hex->face_rotation(f); + for (unsigned int c = 0; c < 4; ++c) + { + auto &new_hex = + new_hexes[face_to_child_indices_hex[f][c]]; + new_hex->set_face_orientation(f, orientation); + new_hex->set_face_flip(f, flip); + new_hex->set_face_rotation(f, rotation); + } + } } } -- 2.39.5