5908 const bool check_for_distorted_cells)
5910 static const int dim = 3;
5912 using raw_line_iterator =
5914 using raw_quad_iterator =
5925 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
5927 if (cell->refine_flag_set())
5953 unsigned int needed_vertices = 0;
5954 unsigned int needed_lines_single = 0;
5955 unsigned int needed_quads_single = 0;
5956 unsigned int needed_lines_pair = 0;
5957 unsigned int needed_quads_pair = 0;
5960 unsigned int new_cells = 0;
5962 for (
const auto &cell :
5964 if (cell->refine_flag_set())
5967 Assert(cell->refine_flag_set() ==
5973 new_cells += cell->reference_cell().n_isotropic_children();
5978 needed_lines_single += 6;
5979 needed_quads_single += 12;
5981 else if (cell->reference_cell() ==
5984 needed_lines_single += 1;
5985 needed_quads_single += 8;
5996 for (
const auto face : cell->face_indices())
5997 if (cell->face(face)->n_children() == 0)
5998 cell->face(face)->set_user_flag();
6000 Assert(cell->face(face)->n_children() ==
6001 cell->reference_cell()
6002 .face_reference_cell(face)
6003 .n_isotropic_children(),
6006 for (
const auto line : cell->line_indices())
6007 if (cell->line(line)->has_children() ==
false)
6008 cell->line(line)->set_user_flag();
6010 Assert(cell->line(line)->n_children() == 2,
6014 const unsigned int used_cells =
6020 used_cells + new_cells,
6034 if (quad->user_flag_set() ==
false)
6039 needed_quads_pair += 4;
6040 needed_lines_pair += 4;
6041 needed_vertices += 1;
6045 needed_quads_pair += 4;
6046 needed_lines_single += 3;
6059 if (line->user_flag_set() ==
false)
6062 needed_lines_pair += 2;
6063 needed_vertices += 1;
6068 needed_lines_single);
6071 needed_quads_single);
6074 needed_quads_single);
6078 needed_vertices += std::count(
triangulation.vertices_used.begin(),
6099 for (
const auto &cell :
triangulation.active_cell_iterators())
6100 if (!cell->refine_flag_set())
6101 for (
unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
6102 if (cell->line(line_n)->has_children())
6103 for (
unsigned int c = 0; c < 2; ++c)
6104 Assert(cell->line(line_n)->child(c)->user_flag_set() ==
false,
6108 unsigned int current_vertex = 0;
6112 auto get_next_unused_vertex = [](
const unsigned int current_vertex,
6113 std::vector<bool> &vertices_used) {
6114 unsigned int next_vertex = current_vertex;
6115 while (next_vertex < vertices_used.size() &&
6116 vertices_used[next_vertex] ==
true)
6119 vertices_used[next_vertex] =
true;
6129 raw_line_iterator next_unused_line =
triangulation.begin_raw_line();
6131 for (; line != endl; ++line)
6133 if (line->user_flag_set() ==
false)
6137 triangulation.faces->lines.template next_free_pair_object<1>(
6145 line->set_children(0, next_unused_line->index());
6147 const std::array<raw_line_iterator, 2> children{
6148 {next_unused_line, ++next_unused_line}};
6154 get_next_unused_vertex(current_vertex,
6156 triangulation.vertices[current_vertex] = line->center(
true);
6158 children[0]->set_bounding_object_indices(
6159 {line->vertex_index(0), current_vertex});
6160 children[1]->set_bounding_object_indices(
6161 {current_vertex, line->vertex_index(1)});
6163 const auto manifold_id = line->manifold_id();
6164 const auto boundary_id = line->boundary_id();
6165 for (
const auto &child : children)
6167 child->set_used_flag();
6168 child->clear_children();
6170 child->clear_user_flag();
6171 child->set_boundary_id_internal(boundary_id);
6172 child->set_manifold_id(manifold_id);
6175 line->clear_user_flag();
6185 for (; quad != endq; ++quad)
6187 if (quad->user_flag_set() ==
false)
6190 const auto reference_face_type = quad->reference_cell();
6194 std::array<raw_line_iterator, 4> new_lines;
6197 for (
unsigned int l = 0; l < 2; ++l)
6199 auto next_unused_line =
6202 new_lines[2 * l] = next_unused_line;
6203 new_lines[2 * l + 1] = ++next_unused_line;
6208 for (
unsigned int l = 0; l < 3; ++l)
6218 for (
const unsigned int line : quad->line_indices())
6226 std::array<raw_quad_iterator, 4> new_quads;
6227 for (
unsigned int q = 0; q < 2; ++q)
6229 auto next_unused_quad =
6233 new_quads[2 * q] = next_unused_quad;
6234 new_quads[2 * q + 1] = ++next_unused_quad;
6236 quad->set_children(2 * q, new_quads[2 * q]->
index());
6240 for (
const auto &quad : new_quads)
6252 for (
const auto i : quad->vertex_indices())
6255 for (
const auto i : quad->line_indices())
6261 get_next_unused_vertex(current_vertex,
6266 quad->center(
true,
true);
6270 std::array<raw_line_iterator, 12> lines;
6271 unsigned int n_lines = 0;
6272 for (
unsigned int l = 0; l < quad->
n_lines(); ++l)
6273 for (
unsigned int c = 0; c < 2; ++c)
6275 static constexpr ::ndarray<unsigned int, 2, 2>
index =
6282 quad->line(l)->child(
index[c][quad->line_orientation(l)]);
6285 for (
unsigned int l = 0; l < quad->
n_lines(); ++l)
6286 lines[n_lines++] = new_lines[l];
6288 std::array<int, 12> line_indices;
6289 for (
unsigned int i = 0; i < n_lines; ++i)
6290 line_indices[i] = lines[i]->
index();
6292 static constexpr ::ndarray<unsigned int, 12, 2>
6293 line_vertices_quad{{{{0, 4}},
6306 static constexpr ::ndarray<unsigned int, 4, 4>
6307 quad_lines_quad{{{{0, 8, 4, 10}},
6312 static constexpr ::ndarray<unsigned int, 12, 2>
6313 line_vertices_tri{{{{0, 3}},
6326 static constexpr ::ndarray<unsigned int, 4, 4>
6327 quad_lines_tri{{{{0, 8, 5, X}},
6332 static constexpr ::ndarray<unsigned int, 4, 4, 2>
6333 quad_line_vertices_tri{
6334 {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
6335 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
6336 {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
6337 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
6339 const auto &line_vertices =
6341 line_vertices_quad :
6343 const auto &quad_lines =
6348 for (
unsigned int i = 0, j = 2 * quad->
n_lines();
6352 auto &new_line = new_lines[i];
6353 new_line->set_bounding_object_indices(
6356 new_line->set_used_flag();
6357 new_line->clear_user_flag();
6358 new_line->clear_user_data();
6359 new_line->clear_children();
6360 new_line->set_boundary_id_internal(quad->boundary_id());
6361 new_line->set_manifold_id(quad->manifold_id());
6365 for (
unsigned int i = 0; i < new_quads.size(); ++i)
6367 auto &new_quad = new_quads[i];
6372 reference_face_type);
6375 new_quad->set_bounding_object_indices(
6376 {line_indices[quad_lines[i][0]],
6377 line_indices[quad_lines[i][1]],
6378 line_indices[quad_lines[i][2]]});
6380 new_quad->set_bounding_object_indices(
6381 {line_indices[quad_lines[i][0]],
6382 line_indices[quad_lines[i][1]],
6383 line_indices[quad_lines[i][2]],
6384 line_indices[quad_lines[i][3]]});
6388 new_quad->set_used_flag();
6389 new_quad->clear_user_flag();
6390 new_quad->clear_user_data();
6391 new_quad->clear_children();
6392 new_quad->set_boundary_id_internal(quad->boundary_id());
6393 new_quad->set_manifold_id(quad->manifold_id());
6396 std::set<unsigned int> s;
6404 for (
const auto f : new_quad->line_indices())
6406 const std::array<unsigned int, 2> vertices_0 = {
6407 {lines[quad_lines[i][f]]->vertex_index(0),
6408 lines[quad_lines[i][f]]->vertex_index(1)}};
6410 const std::array<unsigned int, 2> vertices_1 = {
6414 const auto orientation =
6420 for (
const auto i : vertices_0)
6422 for (
const auto i : vertices_1)
6426 new_quad->set_line_orientation(f, orientation);
6438 static constexpr ::ndarray<unsigned int, 4, 2>
6439 quad_child_boundary_lines{
6440 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
6442 for (
unsigned int i = 0; i < 4; ++i)
6443 for (
unsigned int j = 0; j < 2; ++j)
6444 new_quads[quad_child_boundary_lines[i][j]]
6445 ->set_line_orientation(i, quad->line_orientation(i));
6448 quad->clear_user_flag();
6453 cells_with_distorted_children;
6463 hex->level() >=
static_cast<int>(
level),
6467 hex->level() ==
static_cast<int>(
level);
6470 if (hex->refine_flag_set() ==
6474 const auto &reference_cell_type = hex->reference_cell();
6477 hex->clear_refine_flag();
6478 hex->set_refinement_case(ref_case);
6480 unsigned int n_new_lines = 0;
6481 unsigned int n_new_quads = 0;
6482 unsigned int n_new_hexes = 0;
6499 std::array<raw_line_iterator, 6> new_lines;
6500 for (
unsigned int i = 0; i < n_new_lines; ++i)
6507 new_lines[i]->set_used_flag();
6508 new_lines[i]->clear_user_flag();
6509 new_lines[i]->clear_user_data();
6510 new_lines[i]->clear_children();
6511 new_lines[i]->set_boundary_id_internal(
6513 new_lines[i]->set_manifold_id(hex->manifold_id());
6516 std::array<raw_quad_iterator, 12> new_quads;
6517 for (
unsigned int i = 0; i < n_new_quads; ++i)
6523 auto &new_quad = new_quads[i];
6529 reference_cell_type.face_reference_cell(0));
6532 new_quad->set_used_flag();
6533 new_quad->clear_user_flag();
6534 new_quad->clear_user_data();
6535 new_quad->clear_children();
6536 new_quad->set_boundary_id_internal(
6538 new_quad->set_manifold_id(hex->manifold_id());
6539 for (
const auto j : new_quads[i]->line_indices())
6540 new_quad->set_line_orientation(j,
true);
6549 for (
unsigned int i = 0; i < n_new_hexes; ++i)
6558 new_hexes[i] = next_unused_hex;
6560 auto &new_hex = new_hexes[i];
6564 ->reference_cell[new_hex->index()] =
6565 reference_cell_type;
6568 new_hex->set_used_flag();
6569 new_hex->clear_user_flag();
6570 new_hex->clear_user_data();
6571 new_hex->clear_children();
6572 new_hex->set_material_id(hex->material_id());
6573 new_hex->set_manifold_id(hex->manifold_id());
6574 new_hex->set_subdomain_id(hex->subdomain_id());
6577 new_hex->set_parent(hex->index());
6582 for (
const auto f : new_hex->face_indices())
6583 new_hex->set_combined_face_orientation(
6587 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
6588 hex->set_children(2 * i, new_hexes[2 * i]->index());
6600 const unsigned int n_vertices =
6602 for (
unsigned int i = 0; i < n_vertices; ++i)
6605 const std::array<unsigned int, 12> line_indices =
6606 TriaAccessorImplementation::Implementation::
6607 get_line_indices_of_cell(*hex);
6619 for (
unsigned int l = 0; l < n_lines; ++l)
6629 for (
const unsigned int i : hex->face_indices())
6631 hex->face(i)->child(0)->vertex_index(3);
6635 get_next_unused_vertex(current_vertex,
6640 hex->center(
true,
true);
6644 unsigned int chosen_line_tetrahedron = 0;
6648 static constexpr ::ndarray<unsigned int, 6, 2>
6649 new_line_vertices = {{{{22, 26}},
6655 for (
unsigned int i = 0; i < n_new_lines; ++i)
6656 new_lines[i]->set_bounding_object_indices(
6666 static constexpr ::ndarray<unsigned int, 3, 2>
6667 new_line_vertices = {{{{6, 8}}, {{5, 7}}, {{4, 9}}}};
6670 double min_distance =
6671 std::numeric_limits<double>::infinity();
6672 for (
unsigned int i = 0; i < new_line_vertices.size();
6675 const double current_distance =
6680 if (current_distance < min_distance)
6682 chosen_line_tetrahedron = i;
6683 min_distance = current_distance;
6687 new_lines[0]->set_bounding_object_indices(
6689 [new_line_vertices[chosen_line_tetrahedron][0]],
6691 [new_line_vertices[chosen_line_tetrahedron][1]]});
6696 boost::container::small_vector<raw_line_iterator, 30>
6701 relevant_lines.resize(30);
6702 for (
unsigned int f = 0, k = 0; f < 6; ++f)
6703 for (
unsigned int c = 0; c < 4; ++c, ++k)
6706 ndarray<unsigned int, 4, 2>
6708 {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}};
6714 standard_to_real_face_vertex(
6716 hex->face_orientation(f),
6718 hex->face_rotation(f)))
6720 standard_to_real_face_line(
6722 hex->face_orientation(f),
6724 hex->face_rotation(f)));
6727 for (
unsigned int i = 0, k = 24; i < 6; ++i, ++k)
6728 relevant_lines[k] = new_lines[i];
6743 relevant_lines.resize(13);
6746 for (
unsigned int f = 0; f < 4; ++f)
6747 for (
unsigned int l = 0; l < 3; ++l, ++k)
6751 array<std::array<unsigned int, 3>, 6>
6752 table = {{{{1, 0, 2}},
6765 .get_combined_orientation(
6771 relevant_lines[k++] = new_lines[0];
6777 boost::container::small_vector<unsigned int, 30>
6778 relevant_line_indices(relevant_lines.size());
6779 for (
unsigned int i = 0; i < relevant_line_indices.size();
6781 relevant_line_indices[i] = relevant_lines[i]->
index();
6783 static constexpr ::ndarray<unsigned int, 12, 4>
6784 new_quad_lines_hex = {{{{10, 28, 16, 24}},
6795 {{25, 7, 27, 13}}}};
6812 static constexpr ::ndarray<unsigned int, 3, 12, 4>
6813 new_quad_lines_tet = {{
6853 static constexpr ::ndarray<unsigned int, 12, 4, 2>
6855 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
6856 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
6857 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
6858 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
6859 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
6860 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
6861 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
6862 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
6863 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
6864 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
6865 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
6866 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
6871 static constexpr ::ndarray<unsigned int, 3, 12, 4, 2>
6874 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
6875 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
6876 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
6877 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
6878 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
6879 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
6880 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
6881 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
6882 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6883 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6884 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6885 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}},
6887 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
6888 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
6889 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
6890 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
6891 {{{{5, 4}}, {{4, 7}}, {{7, 5}}, {{X, X}}}},
6892 {{{{6, 5}}, {{5, 7}}, {{7, 6}}, {{X, X}}}},
6893 {{{{8, 7}}, {{7, 5}}, {{5, 8}}, {{X, X}}}},
6894 {{{{7, 9}}, {{9, 5}}, {{5, 7}}, {{X, X}}}},
6895 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6896 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6897 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6898 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}},
6900 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
6901 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
6902 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
6903 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
6904 {{{{5, 4}}, {{4, 9}}, {{9, 5}}, {{X, X}}}},
6905 {{{{4, 6}}, {{6, 9}}, {{9, 4}}, {{X, X}}}},
6906 {{{{7, 4}}, {{4, 9}}, {{9, 7}}, {{X, X}}}},
6907 {{{{4, 8}}, {{8, 9}}, {{9, 4}}, {{X, X}}}},
6908 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6909 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6910 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
6911 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}}}};
6913 const auto &new_quad_lines =
6915 new_quad_lines_hex :
6916 new_quad_lines_tet[chosen_line_tetrahedron];
6921 table_tet[chosen_line_tetrahedron];
6923 static constexpr ::ndarray<unsigned int, 4, 2>
6924 representative_lines{
6925 {{{0, 2}}, {{2, 0}}, {{3, 3}}, {{1, 1}}}};
6927 for (
unsigned int q = 0; q < n_new_quads; ++q)
6929 auto &new_quad = new_quads[q];
6931 if (new_quad->n_lines() == 3)
6932 new_quad->set_bounding_object_indices(
6933 {relevant_line_indices[new_quad_lines[q][0]],
6934 relevant_line_indices[new_quad_lines[q][1]],
6935 relevant_line_indices[new_quad_lines[q][2]]});
6936 else if (new_quad->n_lines() == 4)
6937 new_quad->set_bounding_object_indices(
6938 {relevant_line_indices[new_quad_lines[q][0]],
6939 relevant_line_indices[new_quad_lines[q][1]],
6940 relevant_line_indices[new_quad_lines[q][2]],
6941 relevant_line_indices[new_quad_lines[q][3]]});
6949 const unsigned int n_compute_lines =
6953 for (
unsigned int line = 0; line < n_compute_lines;
6956 const unsigned int l =
6957 (reference_cell_type ==
6959 representative_lines[q % 4][0] :
6962 const std::array<unsigned int, 2> vertices_0 = {
6963 {relevant_lines[new_quad_lines[q][l]]
6965 relevant_lines[new_quad_lines[q][l]]
6966 ->vertex_index(1)}};
6968 const std::array<unsigned int, 2> vertices_1 = {
6972 const auto orientation =
6977 new_quad->set_line_orientation(l, orientation);
6982 if (reference_cell_type ==
6984 new_quads[representative_lines[q % 4][1] + q -
6986 ->set_line_orientation(l, orientation);
6993 std::array<int, 36> quad_indices;
6997 for (
unsigned int i = 0; i < n_new_quads; ++i)
6998 quad_indices[i] = new_quads[i]->
index();
7000 for (
unsigned int f = 0, k = n_new_quads; f < 6; ++f)
7001 for (
unsigned int c = 0; c < 4; ++c, ++k)
7003 hex->face(f)->isotropic_child_index(
7006 hex->face_orientation(f),
7008 hex->face_rotation(f)));
7018 for (
unsigned int i = 0; i < n_new_quads; ++i)
7019 quad_indices[i] = new_quads[i]->
index();
7021 for (
unsigned int f = 0, k = n_new_quads; f < 4; ++f)
7022 for (
unsigned int c = 0; c < 4; ++c, ++k)
7024 quad_indices[k] = hex->face(f)->child_index(
7028 .standard_to_real_face_vertex(
7033 .get_combined_orientation(
7044 static constexpr ::ndarray<unsigned int, 8, 6>
7046 {{12, 0, 20, 4, 28, 8}},
7047 {{0, 16, 22, 6, 29, 9}},
7048 {{13, 1, 4, 24, 30, 10}},
7049 {{1, 17, 6, 26, 31, 11}},
7050 {{14, 2, 21, 5, 8, 32}},
7051 {{2, 18, 23, 7, 9, 33}},
7052 {{15, 3, 5, 25, 10, 34}},
7053 {{3, 19, 7, 27, 11, 35}}
7064 static constexpr ::ndarray<unsigned int, 3, 8, 6>
7066 {{{{8, 13, 16, 0, X, X}},
7067 {{9, 12, 1, 21, X, X}},
7068 {{10, 2, 17, 20, X, X}},
7069 {{3, 14, 18, 22, X, X}},
7070 {{11, 1, 4, 5, X, X}},
7071 {{15, 0, 4, 6, X, X}},
7072 {{19, 7, 6, 3, X, X}},
7073 {{23, 5, 2, 7, X, X}}}},
7076 {{8, 13, 16, 0, X, X}},
7077 {{9, 12, 1, 21, X, X}},
7078 {{10, 2, 17, 20, X, X}},
7079 {{3, 14, 18, 22, X, X}},
7080 {{11, 4, 0, 5, X, X}},
7081 {{15, 4, 1, 6, X, X}},
7082 {{19, 2, 5, 7, X, X}},
7083 {{23, 6, 7, 3, X, X}},
7087 {{8, 13, 16, 0, X, X}},
7088 {{9, 12, 1, 21, X, X}},
7089 {{10, 2, 17, 20, X, X}},
7090 {{3, 14, 18, 22, X, X}},
7091 {{11, 4, 5, 2, X, X}},
7092 {{15, 6, 7, 3, X, X}},
7093 {{19, 5, 0, 6, X, X}},
7094 {{23, 1, 4, 7, X, X}},
7099 static constexpr ::ndarray<unsigned int, 4, 3>
7100 reference_cell_face_vertices = {
7101 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
7109 static constexpr ::ndarray<unsigned int, 3, 8, 4>
7110 cell_vertices_tet = {{
7139 const auto &cell_quads =
7142 cell_quads_tet[chosen_line_tetrahedron];
7144 for (
unsigned int c = 0;
7145 c < GeometryInfo<dim>::max_children_per_cell;
7148 auto &new_hex = new_hexes[c];
7150 if (new_hex->n_faces() == 4)
7152 new_hex->set_bounding_object_indices(
7153 {quad_indices[cell_quads[c][0]],
7154 quad_indices[cell_quads[c][1]],
7155 quad_indices[cell_quads[c][2]],
7156 quad_indices[cell_quads[c][3]]});
7161 for (
const auto f : new_hex->face_indices())
7163 const auto &face = new_hex->face(f);
7165 Assert(face->n_vertices() == 3,
7168 const std::array<unsigned int, 3> vertices_0 = {
7169 {face->vertex_index(0),
7170 face->vertex_index(1),
7171 face->vertex_index(2)}};
7173 const auto &new_hex_vertices =
7174 cell_vertices_tet[chosen_line_tetrahedron][c];
7175 const auto face_f_vertices =
7176 reference_cell_face_vertices[f];
7178 const std::array<unsigned int, 3> vertices_1 = {
7181 [new_hex_vertices[face_f_vertices[0]]],
7183 [new_hex_vertices[face_f_vertices[1]]],
7185 [new_hex_vertices[face_f_vertices[2]]],
7188 new_hex->set_combined_face_orientation(
7190 face->reference_cell()
7191 .get_combined_orientation(
7196 else if (new_hex->n_faces() == 6)
7197 new_hex->set_bounding_object_indices(
7198 {quad_indices[cell_quads[c][0]],
7199 quad_indices[cell_quads[c][1]],
7200 quad_indices[cell_quads[c][2]],
7201 quad_indices[cell_quads[c][3]],
7202 quad_indices[cell_quads[c][4]],
7203 quad_indices[cell_quads[c][5]]});
7212 static constexpr ::ndarray<unsigned int, 6, 4>
7213 face_to_child_indices_hex{{{{0, 2, 4, 6}},
7220 for (
const auto f : hex->face_indices())
7222 const unsigned char combined_orientation =
7223 hex->combined_face_orientation(f);
7224 for (
unsigned int c = 0; c < 4; ++c)
7225 new_hexes[face_to_child_indices_hex[f][c]]
7226 ->set_combined_face_orientation(
7227 f, combined_orientation);
7232 if (check_for_distorted_cells &&
7233 has_distorted_children<dim, spacedim>(hex))
7242 return cells_with_distorted_children;
7252 const bool check_for_distorted_cells)
7254 const unsigned int dim = 3;
7257 bool flag_isotropic_mesh =
true;
7261 for (; cell != endc; ++cell)
7271 flag_isotropic_mesh =
false;
7275 if (flag_isotropic_mesh)
7277 check_for_distorted_cells);
7288 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
7290 if (cell->refine_flag_set())
7307 line->clear_user_flag();
7312 quad->clear_user_flag();
7335 unsigned int needed_vertices = 0;
7336 unsigned int needed_lines_single = 0;
7337 unsigned int needed_quads_single = 0;
7338 unsigned int needed_lines_pair = 0;
7339 unsigned int needed_quads_pair = 0;
7344 unsigned int new_cells = 0;
7346 for (
const auto &acell :
7348 if (acell->refine_flag_set())
7358 ++needed_quads_single;
7366 ++needed_lines_single;
7367 needed_quads_single += 4;
7374 needed_lines_single += 6;
7375 needed_quads_single += 12;
7389 for (
const unsigned int face :
7393 aface = acell->face(face);
7400 acell->face_orientation(face),
7401 acell->face_flip(face),
7402 acell->face_rotation(face));
7407 if (face_ref_case ==
7410 if (aface->n_active_descendants() < 4)
7413 aface->set_user_flag();
7415 else if (aface->refinement_case() != face_ref_case)
7426 Assert(aface->refinement_case() ==
7428 dim - 1>::isotropic_refinement ||
7429 aface->refinement_case() ==
7432 aface->set_user_index(face_ref_case);
7438 for (
unsigned int line = 0;
7439 line < GeometryInfo<dim>::lines_per_cell;
7443 !acell->line(line)->has_children())
7444 acell->line(line)->set_user_flag();
7450 const unsigned int used_cells =
7460 used_cells + new_cells,
7474 if (quad->user_flag_set())
7480 needed_quads_pair += 4;
7481 needed_lines_pair += 4;
7482 needed_vertices += 1;
7484 if (quad->user_index())
7488 needed_quads_pair += 2;
7489 needed_lines_single += 1;
7502 if (quad->has_children())
7504 Assert(quad->refinement_case() ==
7507 if ((face_refinement_cases[quad->user_index()] ==
7509 (quad->child(0)->line_index(1) + 1 !=
7510 quad->child(2)->line_index(1))) ||
7511 (face_refinement_cases[quad->user_index()] ==
7513 (quad->child(0)->line_index(3) + 1 !=
7514 quad->child(1)->line_index(3))))
7515 needed_lines_pair += 2;
7524 if (line->user_flag_set())
7526 needed_lines_pair += 2;
7527 needed_vertices += 1;
7533 needed_lines_single);
7537 needed_quads_single);
7540 needed_quads_single);
7544 needed_vertices += std::count(
triangulation.vertices_used.begin(),
7568 for (
const auto &cell :
triangulation.active_cell_iterators())
7569 if (!cell->refine_flag_set())
7570 for (
unsigned int line = 0;
7571 line < GeometryInfo<dim>::lines_per_cell;
7573 if (cell->line(line)->has_children())
7574 for (
unsigned int c = 0; c < 2; ++c)
7575 Assert(cell->line(line)->child(c)->user_flag_set() ==
false,
7587 unsigned int next_unused_vertex = 0;
7598 for (; line != endl; ++line)
7599 if (line->user_flag_set())
7605 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
7606 ++next_unused_vertex;
7610 "Internal error: During refinement, the triangulation wants to access an element of the 'vertices' array but it turns out that the array is not large enough."));
7613 triangulation.vertices[next_unused_vertex] = line->center(
true);
7627 line->set_children(0, next_unused_line->index());
7631 children[2] = {next_unused_line, ++next_unused_line};
7638 children[0]->set_bounding_object_indices(
7639 {line->vertex_index(0), next_unused_vertex});
7640 children[1]->set_bounding_object_indices(
7641 {next_unused_vertex, line->vertex_index(1)});
7643 children[0]->set_used_flag();
7644 children[1]->set_used_flag();
7645 children[0]->clear_children();
7646 children[1]->clear_children();
7649 children[0]->clear_user_flag();
7650 children[1]->clear_user_flag();
7652 children[0]->set_boundary_id_internal(line->boundary_id());
7653 children[1]->set_boundary_id_internal(line->boundary_id());
7655 children[0]->set_manifold_id(line->manifold_id());
7656 children[1]->set_manifold_id(line->manifold_id());
7661 line->clear_user_flag();
7691 for (
unsigned int loop = 0; loop < 2; ++loop)
7704 for (; quad != endq; ++quad)
7706 if (quad->user_index())
7709 face_refinement_cases[quad->user_index()];
7718 if (aniso_quad_ref_case == quad->refinement_case())
7721 Assert(quad->refinement_case() ==
7723 quad->refinement_case() ==
7728 Assert(quad->user_index() ==
7730 quad->user_index() ==
7760 quad->line(2)->child(0)->vertex_index(1);
7762 quad->line(3)->child(0)->vertex_index(1);
7767 quad->line(0)->child(0)->vertex_index(1);
7769 quad->line(1)->child(0)->vertex_index(1);
7772 new_line->set_bounding_object_indices(
7774 new_line->set_used_flag();
7775 new_line->clear_user_flag();
7777 new_line->clear_children();
7778 new_line->set_boundary_id_internal(quad->boundary_id());
7779 new_line->set_manifold_id(quad->manifold_id());
7787 const unsigned int index[2][2] = {
7799 new_quads[0] = next_unused_quad;
7803 new_quads[1] = next_unused_quad;
7808 new_quads[0]->set_bounding_object_indices(
7809 {
static_cast<int>(quad->line_index(0)),
7812 ->child(
index[0][quad->line_orientation(2)])
7815 ->child(
index[0][quad->line_orientation(3)])
7817 new_quads[1]->set_bounding_object_indices(
7819 static_cast<int>(quad->line_index(1)),
7821 ->child(
index[1][quad->line_orientation(2)])
7824 ->child(
index[1][quad->line_orientation(3)])
7829 new_quads[0]->set_bounding_object_indices(
7831 ->child(
index[0][quad->line_orientation(0)])
7834 ->child(
index[0][quad->line_orientation(1)])
7836 static_cast<int>(quad->line_index(2)),
7837 new_line->index()});
7838 new_quads[1]->set_bounding_object_indices(
7840 ->child(
index[1][quad->line_orientation(0)])
7843 ->child(
index[1][quad->line_orientation(1)])
7846 static_cast<int>(quad->line_index(3))});
7849 for (
const auto &new_quad : new_quads)
7851 new_quad->set_used_flag();
7852 new_quad->clear_user_flag();
7854 new_quad->clear_children();
7855 new_quad->set_boundary_id_internal(quad->boundary_id());
7856 new_quad->set_manifold_id(quad->manifold_id());
7860 for (
unsigned int j = 0;
7861 j < GeometryInfo<dim>::lines_per_face;
7863 new_quad->set_line_orientation(j,
true);
7869 new_quads[0]->set_line_orientation(
7870 0, quad->line_orientation(0));
7871 new_quads[0]->set_line_orientation(
7872 2, quad->line_orientation(2));
7873 new_quads[1]->set_line_orientation(
7874 1, quad->line_orientation(1));
7875 new_quads[1]->set_line_orientation(
7876 3, quad->line_orientation(3));
7879 new_quads[0]->set_line_orientation(
7880 3, quad->line_orientation(3));
7881 new_quads[1]->set_line_orientation(
7882 2, quad->line_orientation(2));
7886 new_quads[0]->set_line_orientation(
7887 1, quad->line_orientation(1));
7888 new_quads[1]->set_line_orientation(
7889 0, quad->line_orientation(0));
7895 if (quad->refinement_case() ==
7918 if (aniso_quad_ref_case ==
7921 old_child[0] = quad->child(0)->line(1);
7922 old_child[1] = quad->child(2)->line(1);
7926 Assert(aniso_quad_ref_case ==
7930 old_child[0] = quad->child(0)->line(3);
7931 old_child[1] = quad->child(1)->line(3);
7934 if (old_child[0]->
index() + 1 != old_child[1]->
index())
7940 spacedim>::raw_line_iterator
7943 new_child[0] = new_child[1] =
7945 .template next_free_pair_object<1>(
7949 new_child[0]->set_used_flag();
7950 new_child[1]->set_used_flag();
7952 const int old_index_0 = old_child[0]->index(),
7953 old_index_1 = old_child[1]->index(),
7954 new_index_0 = new_child[0]->index(),
7955 new_index_1 = new_child[1]->index();
7959 for (
unsigned int q = 0;
7962 for (
unsigned int l = 0;
7963 l < GeometryInfo<dim>::lines_per_face;
7966 const int this_index =
7968 .get_bounding_object_indices(q)[l];
7969 if (this_index == old_index_0)
7971 .get_bounding_object_indices(q)[l] =
7973 else if (this_index == old_index_1)
7975 .get_bounding_object_indices(q)[l] =
7980 for (
unsigned int i = 0; i < 2; ++i)
7982 Assert(!old_child[i]->has_children(),
7985 new_child[i]->set_bounding_object_indices(
7986 {old_child[i]->vertex_index(0),
7987 old_child[i]->vertex_index(1)});
7988 new_child[i]->set_boundary_id_internal(
7989 old_child[i]->boundary_id());
7990 new_child[i]->set_manifold_id(
7991 old_child[i]->manifold_id());
7992 new_child[i]->set_user_index(
7993 old_child[i]->user_index());
7994 if (old_child[i]->user_flag_set())
7995 new_child[i]->set_user_flag();
7997 new_child[i]->clear_user_flag();
7999 new_child[i]->clear_children();
8001 old_child[i]->clear_user_flag();
8002 old_child[i]->clear_user_index();
8003 old_child[i]->clear_used_flag();
8009 if (aniso_quad_ref_case ==
8012 new_line->set_children(
8013 0, quad->child(0)->line_index(1));
8014 Assert(new_line->child(1) ==
8015 quad->child(2)->line(1),
8052 quad_iterator switch_1 = quad->child(1),
8053 switch_2 = quad->child(2);
8054 const int switch_1_index = switch_1->index();
8055 const int switch_2_index = switch_2->index();
8056 for (
unsigned int l = 0;
8059 for (
unsigned int h = 0;
8063 for (
const unsigned int q :
8066 const int face_index =
8068 ->cells.get_bounding_object_indices(
8070 if (face_index == switch_1_index)
8072 ->cells.get_bounding_object_indices(
8073 h)[q] = switch_2_index;
8074 else if (face_index == switch_2_index)
8076 ->cells.get_bounding_object_indices(
8077 h)[q] = switch_1_index;
8081 const unsigned int switch_1_lines[4] = {
8082 switch_1->line_index(0),
8083 switch_1->line_index(1),
8084 switch_1->line_index(2),
8085 switch_1->line_index(3)};
8086 const bool switch_1_line_orientations[4] = {
8087 switch_1->line_orientation(0),
8088 switch_1->line_orientation(1),
8089 switch_1->line_orientation(2),
8090 switch_1->line_orientation(3)};
8092 switch_1->boundary_id();
8093 const unsigned int switch_1_user_index =
8094 switch_1->user_index();
8095 const bool switch_1_user_flag =
8096 switch_1->user_flag_set();
8098 switch_1_refinement_case =
8099 switch_1->refinement_case();
8100 const int switch_1_first_child_pair =
8101 (switch_1_refinement_case ?
8102 switch_1->child_index(0) :
8104 const int switch_1_second_child_pair =
8105 (switch_1_refinement_case ==
8107 switch_1->child_index(2) :
8110 switch_1->set_bounding_object_indices(
8111 {switch_2->line_index(0),
8112 switch_2->line_index(1),
8113 switch_2->line_index(2),
8114 switch_2->line_index(3)});
8115 switch_1->set_line_orientation(
8116 0, switch_2->line_orientation(0));
8117 switch_1->set_line_orientation(
8118 1, switch_2->line_orientation(1));
8119 switch_1->set_line_orientation(
8120 2, switch_2->line_orientation(2));
8121 switch_1->set_line_orientation(
8122 3, switch_2->line_orientation(3));
8123 switch_1->set_boundary_id_internal(
8124 switch_2->boundary_id());
8125 switch_1->set_manifold_id(switch_2->manifold_id());
8126 switch_1->set_user_index(switch_2->user_index());
8127 if (switch_2->user_flag_set())
8128 switch_1->set_user_flag();
8130 switch_1->clear_user_flag();
8131 switch_1->clear_refinement_case();
8132 switch_1->set_refinement_case(
8133 switch_2->refinement_case());
8134 switch_1->clear_children();
8135 if (switch_2->refinement_case())
8136 switch_1->set_children(0,
8137 switch_2->child_index(0));
8138 if (switch_2->refinement_case() ==
8140 switch_1->set_children(2,
8141 switch_2->child_index(2));
8143 switch_2->set_bounding_object_indices(
8147 switch_1_lines[3]});
8148 switch_2->set_line_orientation(
8149 0, switch_1_line_orientations[0]);
8150 switch_2->set_line_orientation(
8151 1, switch_1_line_orientations[1]);
8152 switch_2->set_line_orientation(
8153 2, switch_1_line_orientations[2]);
8154 switch_2->set_line_orientation(
8155 3, switch_1_line_orientations[3]);
8156 switch_2->set_boundary_id_internal(
8157 switch_1_boundary_id);
8158 switch_2->set_manifold_id(switch_1->manifold_id());
8159 switch_2->set_user_index(switch_1_user_index);
8160 if (switch_1_user_flag)
8161 switch_2->set_user_flag();
8163 switch_2->clear_user_flag();
8164 switch_2->clear_refinement_case();
8165 switch_2->set_refinement_case(
8166 switch_1_refinement_case);
8167 switch_2->clear_children();
8168 switch_2->set_children(0,
8169 switch_1_first_child_pair);
8170 switch_2->set_children(2,
8171 switch_1_second_child_pair);
8173 new_quads[0]->set_refinement_case(
8175 new_quads[0]->set_children(0, quad->child_index(0));
8176 new_quads[1]->set_refinement_case(
8178 new_quads[1]->set_children(0, quad->child_index(2));
8182 new_quads[0]->set_refinement_case(
8184 new_quads[0]->set_children(0, quad->child_index(0));
8185 new_quads[1]->set_refinement_case(
8187 new_quads[1]->set_children(0, quad->child_index(2));
8188 new_line->set_children(
8189 0, quad->child(0)->line_index(3));
8190 Assert(new_line->child(1) ==
8191 quad->child(1)->line(3),
8194 quad->clear_children();
8198 quad->set_children(0, new_quads[0]->
index());
8200 quad->set_refinement_case(aniso_quad_ref_case);
8207 if (quad->user_flag_set())
8219 ++next_unused_vertex;
8223 "Internal error: During refinement, the triangulation wants to access an element of the 'vertices' array but it turns out that the array is not large enough."));
8231 quad->refinement_case();
8237 quad->child(0)->set_user_index(
8239 quad->child(1)->set_user_index(
8245 middle_line = quad->child(0)->line(1);
8247 middle_line = quad->child(0)->line(3);
8254 if (!middle_line->has_children())
8262 middle_line->center(
true);
8270 .template next_free_pair_object<1>(
8275 middle_line->set_children(
8276 0, next_unused_line->index());
8280 raw_line_iterator children[2] = {
8281 next_unused_line, ++next_unused_line};
8289 children[0]->set_bounding_object_indices(
8290 {middle_line->vertex_index(0),
8291 next_unused_vertex});
8292 children[1]->set_bounding_object_indices(
8293 {next_unused_vertex,
8294 middle_line->vertex_index(1)});
8296 children[0]->set_used_flag();
8297 children[1]->set_used_flag();
8298 children[0]->clear_children();
8299 children[1]->clear_children();
8302 children[0]->clear_user_flag();
8303 children[1]->clear_user_flag();
8305 children[0]->set_boundary_id_internal(
8306 middle_line->boundary_id());
8307 children[1]->set_boundary_id_internal(
8308 middle_line->boundary_id());
8310 children[0]->set_manifold_id(
8311 middle_line->manifold_id());
8312 children[1]->set_manifold_id(
8313 middle_line->manifold_id());
8319 quad->clear_user_flag();
8349 quad->center(
true,
true);
8358 for (
unsigned int i = 0; i < 4; ++i)
8371 new_lines[i] = next_unused_line;
8395 quad->line(0)->child(0)->vertex_index(1),
8396 quad->line(1)->child(0)->vertex_index(1),
8397 quad->line(2)->child(0)->vertex_index(1),
8398 quad->line(3)->child(0)->vertex_index(1),
8399 next_unused_vertex};
8401 new_lines[0]->set_bounding_object_indices(
8403 new_lines[1]->set_bounding_object_indices(
8405 new_lines[2]->set_bounding_object_indices(
8407 new_lines[3]->set_bounding_object_indices(
8410 for (
const auto &new_line : new_lines)
8412 new_line->set_used_flag();
8413 new_line->clear_user_flag();
8415 new_line->clear_children();
8416 new_line->set_boundary_id_internal(quad->boundary_id());
8417 new_line->set_manifold_id(quad->manifold_id());
8436 const unsigned int index[2][2] = {
8440 const int line_indices[12] = {
8442 ->child(
index[0][quad->line_orientation(0)])
8445 ->child(
index[1][quad->line_orientation(0)])
8448 ->child(
index[0][quad->line_orientation(1)])
8451 ->child(
index[1][quad->line_orientation(1)])
8454 ->child(
index[0][quad->line_orientation(2)])
8457 ->child(
index[1][quad->line_orientation(2)])
8460 ->child(
index[0][quad->line_orientation(3)])
8463 ->child(
index[1][quad->line_orientation(3)])
8465 new_lines[0]->index(),
8466 new_lines[1]->index(),
8467 new_lines[2]->index(),
8468 new_lines[3]->index()};
8480 new_quads[0] = next_unused_quad;
8484 new_quads[1] = next_unused_quad;
8490 new_quads[2] = next_unused_quad;
8494 new_quads[3] = next_unused_quad;
8498 quad->set_children(0, new_quads[0]->
index());
8499 quad->set_children(2, new_quads[2]->
index());
8502 new_quads[0]->set_bounding_object_indices(
8507 new_quads[1]->set_bounding_object_indices(
8512 new_quads[2]->set_bounding_object_indices(
8517 new_quads[3]->set_bounding_object_indices(
8522 for (
const auto &new_quad : new_quads)
8524 new_quad->set_used_flag();
8525 new_quad->clear_user_flag();
8527 new_quad->clear_children();
8528 new_quad->set_boundary_id_internal(quad->boundary_id());
8529 new_quad->set_manifold_id(quad->manifold_id());
8533 for (
unsigned int j = 0;
8534 j < GeometryInfo<dim>::lines_per_face;
8536 new_quad->set_line_orientation(j,
true);
8542 new_quads[0]->set_line_orientation(
8543 0, quad->line_orientation(0));
8544 new_quads[0]->set_line_orientation(
8545 2, quad->line_orientation(2));
8546 new_quads[1]->set_line_orientation(
8547 1, quad->line_orientation(1));
8548 new_quads[1]->set_line_orientation(
8549 2, quad->line_orientation(2));
8550 new_quads[2]->set_line_orientation(
8551 0, quad->line_orientation(0));
8552 new_quads[2]->set_line_orientation(
8553 3, quad->line_orientation(3));
8554 new_quads[3]->set_line_orientation(
8555 1, quad->line_orientation(1));
8556 new_quads[3]->set_line_orientation(
8557 3, quad->line_orientation(3));
8561 quad->clear_user_flag();
8572 cells_with_distorted_children;
8586 for (; hex != endh; ++hex)
8587 if (hex->refine_flag_set())
8595 hex->clear_refine_flag();
8596 hex->set_refinement_case(ref_case);
8607 unsigned int n_new_lines = 0;
8608 unsigned int n_new_quads = 0;
8609 unsigned int n_new_hexes = 0;
8640 new_lines(n_new_lines);
8641 for (
unsigned int i = 0; i < n_new_lines; ++i)
8648 new_lines[i]->set_used_flag();
8649 new_lines[i]->clear_user_flag();
8650 new_lines[i]->clear_user_data();
8651 new_lines[i]->clear_children();
8653 new_lines[i]->set_boundary_id_internal(
8657 new_lines[i]->set_manifold_id(hex->manifold_id());
8664 new_quads(n_new_quads);
8665 for (
unsigned int i = 0; i < n_new_quads; ++i)
8672 new_quads[i]->set_used_flag();
8673 new_quads[i]->clear_user_flag();
8674 new_quads[i]->clear_user_data();
8675 new_quads[i]->clear_children();
8677 new_quads[i]->set_boundary_id_internal(
8681 new_quads[i]->set_manifold_id(hex->manifold_id());
8684 for (
unsigned int j = 0;
8685 j < GeometryInfo<dim>::lines_per_face;
8687 new_quads[i]->set_line_orientation(j,
true);
8696 new_hexes(n_new_hexes);
8697 for (
unsigned int i = 0; i < n_new_hexes; ++i)
8706 new_hexes[i] = next_unused_hex;
8709 new_hexes[i]->set_used_flag();
8710 new_hexes[i]->clear_user_flag();
8711 new_hexes[i]->clear_user_data();
8712 new_hexes[i]->clear_children();
8715 new_hexes[i]->set_material_id(hex->material_id());
8716 new_hexes[i]->set_manifold_id(hex->manifold_id());
8717 new_hexes[i]->set_subdomain_id(subdomainid);
8720 new_hexes[i]->set_parent(hex->index());
8732 for (
const unsigned int f :
8734 new_hexes[i]->set_combined_face_orientation(
8739 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
8740 hex->set_children(2 * i, new_hexes[2 * i]->index());
8747 const bool f_or[6] = {hex->face_orientation(0),
8748 hex->face_orientation(1),
8749 hex->face_orientation(2),
8750 hex->face_orientation(3),
8751 hex->face_orientation(4),
8752 hex->face_orientation(5)};
8755 const bool f_fl[6] = {hex->face_flip(0),
8763 const bool f_ro[6] = {hex->face_rotation(0),
8764 hex->face_rotation(1),
8765 hex->face_rotation(2),
8766 hex->face_rotation(3),
8767 hex->face_rotation(4),
8768 hex->face_rotation(5)};
8771 const unsigned char f_co[6] = {
8772 hex->combined_face_orientation(0),
8773 hex->combined_face_orientation(1),
8774 hex->combined_face_orientation(2),
8775 hex->combined_face_orientation(3),
8776 hex->combined_face_orientation(4),
8777 hex->combined_face_orientation(5)};
8787 const unsigned int child_at_origin[2][2][2] = {
8874 raw_line_iterator lines[4] = {
8875 hex->face(2)->child(0)->line(
8876 (hex->face(2)->refinement_case() ==
8880 hex->face(3)->child(0)->line(
8881 (hex->face(3)->refinement_case() ==
8885 hex->face(4)->child(0)->line(
8886 (hex->face(4)->refinement_case() ==
8890 hex->face(5)->child(0)->line(
8891 (hex->face(5)->refinement_case() ==
8897 unsigned int line_indices[4];
8898 for (
unsigned int i = 0; i < 4; ++i)
8899 line_indices[i] = lines[i]->
index();
8906 bool line_orientation[4];
8912 const unsigned int middle_vertices[2] = {
8913 hex->line(2)->child(0)->vertex_index(1),
8914 hex->line(7)->child(0)->vertex_index(1)};
8916 for (
unsigned int i = 0; i < 4; ++i)
8917 if (lines[i]->vertex_index(i % 2) ==
8918 middle_vertices[i % 2])
8919 line_orientation[i] =
true;
8924 Assert(lines[i]->vertex_index((i + 1) % 2) ==
8925 middle_vertices[i % 2],
8927 line_orientation[i] =
false;
8932 new_quads[0]->set_bounding_object_indices(
8938 new_quads[0]->set_line_orientation(
8939 0, line_orientation[0]);
8940 new_quads[0]->set_line_orientation(
8941 1, line_orientation[1]);
8942 new_quads[0]->set_line_orientation(
8943 2, line_orientation[2]);
8944 new_quads[0]->set_line_orientation(
8945 3, line_orientation[3]);
8978 const int quad_indices[11] = {
8979 new_quads[0]->index(),
8981 hex->face(0)->index(),
8983 hex->face(1)->index(),
8985 hex->face(2)->child_index(
8986 child_at_origin[hex->face(2)->refinement_case() -
8987 1][f_fl[2]][f_ro[2]]),
8988 hex->face(2)->child_index(
8990 child_at_origin[hex->face(2)->refinement_case() -
8991 1][f_fl[2]][f_ro[2]]),
8993 hex->face(3)->child_index(
8994 child_at_origin[hex->face(3)->refinement_case() -
8995 1][f_fl[3]][f_ro[3]]),
8996 hex->face(3)->child_index(
8998 child_at_origin[hex->face(3)->refinement_case() -
8999 1][f_fl[3]][f_ro[3]]),
9001 hex->face(4)->child_index(
9002 child_at_origin[hex->face(4)->refinement_case() -
9003 1][f_fl[4]][f_ro[4]]),
9004 hex->face(4)->child_index(
9006 child_at_origin[hex->face(4)->refinement_case() -
9007 1][f_fl[4]][f_ro[4]]),
9009 hex->face(5)->child_index(
9010 child_at_origin[hex->face(5)->refinement_case() -
9011 1][f_fl[5]][f_ro[5]]),
9012 hex->face(5)->child_index(
9014 child_at_origin[hex->face(5)->refinement_case() -
9015 1][f_fl[5]][f_ro[5]])
9019 new_hexes[0]->set_bounding_object_indices(
9026 new_hexes[1]->set_bounding_object_indices(
9101 raw_line_iterator lines[4] = {
9102 hex->face(0)->child(0)->line(
9103 (hex->face(0)->refinement_case() ==
9107 hex->face(1)->child(0)->line(
9108 (hex->face(1)->refinement_case() ==
9112 hex->face(4)->child(0)->line(
9113 (hex->face(4)->refinement_case() ==
9117 hex->face(5)->child(0)->line(
9118 (hex->face(5)->refinement_case() ==
9124 unsigned int line_indices[4];
9125 for (
unsigned int i = 0; i < 4; ++i)
9126 line_indices[i] = lines[i]->
index();
9133 bool line_orientation[4];
9139 const unsigned int middle_vertices[2] = {
9140 hex->line(0)->child(0)->vertex_index(1),
9141 hex->line(5)->child(0)->vertex_index(1)};
9143 for (
unsigned int i = 0; i < 4; ++i)
9144 if (lines[i]->vertex_index(i % 2) ==
9145 middle_vertices[i % 2])
9146 line_orientation[i] =
true;
9150 Assert(lines[i]->vertex_index((i + 1) % 2) ==
9151 middle_vertices[i % 2],
9153 line_orientation[i] =
false;
9158 new_quads[0]->set_bounding_object_indices(
9164 new_quads[0]->set_line_orientation(
9165 0, line_orientation[2]);
9166 new_quads[0]->set_line_orientation(
9167 1, line_orientation[3]);
9168 new_quads[0]->set_line_orientation(
9169 2, line_orientation[0]);
9170 new_quads[0]->set_line_orientation(
9171 3, line_orientation[1]);
9204 const int quad_indices[11] = {
9205 new_quads[0]->index(),
9207 hex->face(0)->child_index(
9208 child_at_origin[hex->face(0)->refinement_case() -
9209 1][f_fl[0]][f_ro[0]]),
9210 hex->face(0)->child_index(
9212 child_at_origin[hex->face(0)->refinement_case() -
9213 1][f_fl[0]][f_ro[0]]),
9215 hex->face(1)->child_index(
9216 child_at_origin[hex->face(1)->refinement_case() -
9217 1][f_fl[1]][f_ro[1]]),
9218 hex->face(1)->child_index(
9220 child_at_origin[hex->face(1)->refinement_case() -
9221 1][f_fl[1]][f_ro[1]]),
9223 hex->face(2)->index(),
9225 hex->face(3)->index(),
9227 hex->face(4)->child_index(
9228 child_at_origin[hex->face(4)->refinement_case() -
9229 1][f_fl[4]][f_ro[4]]),
9230 hex->face(4)->child_index(
9232 child_at_origin[hex->face(4)->refinement_case() -
9233 1][f_fl[4]][f_ro[4]]),
9235 hex->face(5)->child_index(
9236 child_at_origin[hex->face(5)->refinement_case() -
9237 1][f_fl[5]][f_ro[5]]),
9238 hex->face(5)->child_index(
9240 child_at_origin[hex->face(5)->refinement_case() -
9241 1][f_fl[5]][f_ro[5]])
9245 new_hexes[0]->set_bounding_object_indices(
9252 new_hexes[1]->set_bounding_object_indices(
9329 raw_line_iterator lines[4] = {
9330 hex->face(0)->child(0)->line(
9331 (hex->face(0)->refinement_case() ==
9335 hex->face(1)->child(0)->line(
9336 (hex->face(1)->refinement_case() ==
9340 hex->face(2)->child(0)->line(
9341 (hex->face(2)->refinement_case() ==
9345 hex->face(3)->child(0)->line(
9346 (hex->face(3)->refinement_case() ==
9352 unsigned int line_indices[4];
9353 for (
unsigned int i = 0; i < 4; ++i)
9354 line_indices[i] = lines[i]->
index();
9361 bool line_orientation[4];
9367 const unsigned int middle_vertices[2] = {
9368 middle_vertex_index<dim, spacedim>(hex->line(8)),
9369 middle_vertex_index<dim, spacedim>(hex->line(11))};
9371 for (
unsigned int i = 0; i < 4; ++i)
9372 if (lines[i]->vertex_index(i % 2) ==
9373 middle_vertices[i % 2])
9374 line_orientation[i] =
true;
9378 Assert(lines[i]->vertex_index((i + 1) % 2) ==
9379 middle_vertices[i % 2],
9381 line_orientation[i] =
false;
9386 new_quads[0]->set_bounding_object_indices(
9392 new_quads[0]->set_line_orientation(
9393 0, line_orientation[0]);
9394 new_quads[0]->set_line_orientation(
9395 1, line_orientation[1]);
9396 new_quads[0]->set_line_orientation(
9397 2, line_orientation[2]);
9398 new_quads[0]->set_line_orientation(
9399 3, line_orientation[3]);
9433 const int quad_indices[11] = {
9434 new_quads[0]->index(),
9436 hex->face(0)->child_index(
9437 child_at_origin[hex->face(0)->refinement_case() -
9438 1][f_fl[0]][f_ro[0]]),
9439 hex->face(0)->child_index(
9441 child_at_origin[hex->face(0)->refinement_case() -
9442 1][f_fl[0]][f_ro[0]]),
9444 hex->face(1)->child_index(
9445 child_at_origin[hex->face(1)->refinement_case() -
9446 1][f_fl[1]][f_ro[1]]),
9447 hex->face(1)->child_index(
9449 child_at_origin[hex->face(1)->refinement_case() -
9450 1][f_fl[1]][f_ro[1]]),
9452 hex->face(2)->child_index(
9453 child_at_origin[hex->face(2)->refinement_case() -
9454 1][f_fl[2]][f_ro[2]]),
9455 hex->face(2)->child_index(
9457 child_at_origin[hex->face(2)->refinement_case() -
9458 1][f_fl[2]][f_ro[2]]),
9460 hex->face(3)->child_index(
9461 child_at_origin[hex->face(3)->refinement_case() -
9462 1][f_fl[3]][f_ro[3]]),
9463 hex->face(3)->child_index(
9465 child_at_origin[hex->face(3)->refinement_case() -
9466 1][f_fl[3]][f_ro[3]]),
9468 hex->face(4)->index(),
9470 hex->face(5)->index()
9473 new_hexes[0]->set_bounding_object_indices(
9480 new_hexes[1]->set_bounding_object_indices(
9512 new_lines[0]->set_bounding_object_indices(
9513 {middle_vertex_index<dim, spacedim>(hex->face(4)),
9514 middle_vertex_index<dim, spacedim>(hex->face(5))});
9582 spacedim>::raw_line_iterator lines[13] = {
9583 hex->face(0)->child(0)->line(
9584 (hex->face(0)->refinement_case() ==
9588 hex->face(1)->child(0)->line(
9589 (hex->face(1)->refinement_case() ==
9593 hex->face(2)->child(0)->line(
9594 (hex->face(2)->refinement_case() ==
9598 hex->face(3)->child(0)->line(
9599 (hex->face(3)->refinement_case() ==
9607 0, f_or[4], f_fl[4], f_ro[4]))
9610 1, f_or[4], f_fl[4], f_ro[4])),
9614 3, f_or[4], f_fl[4], f_ro[4]))
9617 0, f_or[4], f_fl[4], f_ro[4])),
9621 0, f_or[4], f_fl[4], f_ro[4]))
9624 3, f_or[4], f_fl[4], f_ro[4])),
9628 3, f_or[4], f_fl[4], f_ro[4]))
9631 2, f_or[4], f_fl[4], f_ro[4])),
9636 0, f_or[5], f_fl[5], f_ro[5]))
9639 1, f_or[5], f_fl[5], f_ro[5])),
9643 3, f_or[5], f_fl[5], f_ro[5]))
9646 0, f_or[5], f_fl[5], f_ro[5])),
9650 0, f_or[5], f_fl[5], f_ro[5]))
9653 3, f_or[5], f_fl[5], f_ro[5])),
9657 3, f_or[5], f_fl[5], f_ro[5]))
9660 2, f_or[5], f_fl[5], f_ro[5])),
9665 unsigned int line_indices[13];
9666 for (
unsigned int i = 0; i < 13; ++i)
9667 line_indices[i] = lines[i]->
index();
9674 bool line_orientation[13];
9678 const unsigned int middle_vertices[4] = {
9679 hex->line(0)->child(0)->vertex_index(1),
9680 hex->line(1)->child(0)->vertex_index(1),
9681 hex->line(2)->child(0)->vertex_index(1),
9682 hex->line(3)->child(0)->vertex_index(1),
9688 for (
unsigned int i = 0; i < 4; ++i)
9689 if (lines[i]->vertex_index(0) == middle_vertices[i])
9690 line_orientation[i] =
true;
9694 Assert(lines[i]->vertex_index(1) ==
9697 line_orientation[i] =
false;
9706 for (
unsigned int i = 4; i < 12; ++i)
9707 if (lines[i]->vertex_index((i + 1) % 2) ==
9708 middle_vertex_index<dim, spacedim>(
9709 hex->face(3 + i / 4)))
9710 line_orientation[i] =
true;
9715 Assert(lines[i]->vertex_index(i % 2) ==
9716 (middle_vertex_index<dim, spacedim>(
9717 hex->face(3 + i / 4))),
9719 line_orientation[i] =
false;
9724 line_orientation[12] =
true;
9749 new_quads[0]->set_bounding_object_indices(
9754 new_quads[1]->set_bounding_object_indices(
9759 new_quads[2]->set_bounding_object_indices(
9764 new_quads[3]->set_bounding_object_indices(
9770 new_quads[0]->set_line_orientation(
9771 0, line_orientation[2]);
9772 new_quads[0]->set_line_orientation(
9773 2, line_orientation[4]);
9774 new_quads[0]->set_line_orientation(
9775 3, line_orientation[8]);
9777 new_quads[1]->set_line_orientation(
9778 1, line_orientation[3]);
9779 new_quads[1]->set_line_orientation(
9780 2, line_orientation[5]);
9781 new_quads[1]->set_line_orientation(
9782 3, line_orientation[9]);
9784 new_quads[2]->set_line_orientation(
9785 0, line_orientation[6]);
9786 new_quads[2]->set_line_orientation(
9787 1, line_orientation[10]);
9788 new_quads[2]->set_line_orientation(
9789 2, line_orientation[0]);
9791 new_quads[3]->set_line_orientation(
9792 0, line_orientation[7]);
9793 new_quads[3]->set_line_orientation(
9794 1, line_orientation[11]);
9795 new_quads[3]->set_line_orientation(
9796 3, line_orientation[1]);
9829 const int quad_indices[20] = {
9830 new_quads[0]->index(),
9831 new_quads[1]->index(),
9832 new_quads[2]->index(),
9833 new_quads[3]->index(),
9835 hex->face(0)->child_index(
9836 child_at_origin[hex->face(0)->refinement_case() -
9837 1][f_fl[0]][f_ro[0]]),
9838 hex->face(0)->child_index(
9840 child_at_origin[hex->face(0)->refinement_case() -
9841 1][f_fl[0]][f_ro[0]]),
9843 hex->face(1)->child_index(
9844 child_at_origin[hex->face(1)->refinement_case() -
9845 1][f_fl[1]][f_ro[1]]),
9846 hex->face(1)->child_index(
9848 child_at_origin[hex->face(1)->refinement_case() -
9849 1][f_fl[1]][f_ro[1]]),
9851 hex->face(2)->child_index(
9852 child_at_origin[hex->face(2)->refinement_case() -
9853 1][f_fl[2]][f_ro[2]]),
9854 hex->face(2)->child_index(
9856 child_at_origin[hex->face(2)->refinement_case() -
9857 1][f_fl[2]][f_ro[2]]),
9859 hex->face(3)->child_index(
9860 child_at_origin[hex->face(3)->refinement_case() -
9861 1][f_fl[3]][f_ro[3]]),
9862 hex->face(3)->child_index(
9864 child_at_origin[hex->face(3)->refinement_case() -
9865 1][f_fl[3]][f_ro[3]]),
9867 hex->face(4)->isotropic_child_index(
9869 0, f_or[4], f_fl[4], f_ro[4])),
9870 hex->face(4)->isotropic_child_index(
9872 1, f_or[4], f_fl[4], f_ro[4])),
9873 hex->face(4)->isotropic_child_index(
9875 2, f_or[4], f_fl[4], f_ro[4])),
9876 hex->face(4)->isotropic_child_index(
9878 3, f_or[4], f_fl[4], f_ro[4])),
9880 hex->face(5)->isotropic_child_index(
9882 0, f_or[5], f_fl[5], f_ro[5])),
9883 hex->face(5)->isotropic_child_index(
9885 1, f_or[5], f_fl[5], f_ro[5])),
9886 hex->face(5)->isotropic_child_index(
9888 2, f_or[5], f_fl[5], f_ro[5])),
9889 hex->face(5)->isotropic_child_index(
9891 3, f_or[5], f_fl[5], f_ro[5]))};
9893 new_hexes[0]->set_bounding_object_indices(
9900 new_hexes[1]->set_bounding_object_indices(
9907 new_hexes[2]->set_bounding_object_indices(
9914 new_hexes[3]->set_bounding_object_indices(
9946 new_lines[0]->set_bounding_object_indices(
9947 {middle_vertex_index<dim, spacedim>(hex->face(2)),
9948 middle_vertex_index<dim, spacedim>(hex->face(3))});
10016 spacedim>::raw_line_iterator lines[13] = {
10017 hex->face(0)->child(0)->line(
10018 (hex->face(0)->refinement_case() ==
10022 hex->face(1)->child(0)->line(
10023 (hex->face(1)->refinement_case() ==
10027 hex->face(4)->child(0)->line(
10028 (hex->face(4)->refinement_case() ==
10032 hex->face(5)->child(0)->line(
10033 (hex->face(5)->refinement_case() ==
10041 0, f_or[2], f_fl[2], f_ro[2]))
10044 3, f_or[2], f_fl[2], f_ro[2])),
10048 3, f_or[2], f_fl[2], f_ro[2]))
10051 2, f_or[2], f_fl[2], f_ro[2])),
10055 0, f_or[2], f_fl[2], f_ro[2]))
10058 1, f_or[2], f_fl[2], f_ro[2])),
10062 3, f_or[2], f_fl[2], f_ro[2]))
10065 0, f_or[2], f_fl[2], f_ro[2])),
10070 0, f_or[3], f_fl[3], f_ro[3]))
10073 3, f_or[3], f_fl[3], f_ro[3])),
10077 3, f_or[3], f_fl[3], f_ro[3]))
10080 2, f_or[3], f_fl[3], f_ro[3])),
10084 0, f_or[3], f_fl[3], f_ro[3]))
10087 1, f_or[3], f_fl[3], f_ro[3])),
10091 3, f_or[3], f_fl[3], f_ro[3]))
10094 0, f_or[3], f_fl[3], f_ro[3])),
10099 unsigned int line_indices[13];
10100 for (
unsigned int i = 0; i < 13; ++i)
10101 line_indices[i] = lines[i]->
index();
10108 bool line_orientation[13];
10112 const unsigned int middle_vertices[4] = {
10113 hex->line(8)->child(0)->vertex_index(1),
10114 hex->line(9)->child(0)->vertex_index(1),
10115 hex->line(2)->child(0)->vertex_index(1),
10116 hex->line(6)->child(0)->vertex_index(1),
10121 for (
unsigned int i = 0; i < 4; ++i)
10122 if (lines[i]->vertex_index(0) == middle_vertices[i])
10123 line_orientation[i] =
true;
10127 Assert(lines[i]->vertex_index(1) ==
10128 middle_vertices[i],
10130 line_orientation[i] =
false;
10139 for (
unsigned int i = 4; i < 12; ++i)
10140 if (lines[i]->vertex_index((i + 1) % 2) ==
10141 middle_vertex_index<dim, spacedim>(
10142 hex->face(1 + i / 4)))
10143 line_orientation[i] =
true;
10148 Assert(lines[i]->vertex_index(i % 2) ==
10149 (middle_vertex_index<dim, spacedim>(
10150 hex->face(1 + i / 4))),
10152 line_orientation[i] =
false;
10157 line_orientation[12] =
true;
10183 new_quads[0]->set_bounding_object_indices(
10187 line_indices[10]});
10188 new_quads[1]->set_bounding_object_indices(
10192 line_indices[11]});
10193 new_quads[2]->set_bounding_object_indices(
10197 line_indices[12]});
10198 new_quads[3]->set_bounding_object_indices(
10204 new_quads[0]->set_line_orientation(
10205 0, line_orientation[0]);
10206 new_quads[0]->set_line_orientation(
10207 2, line_orientation[6]);
10208 new_quads[0]->set_line_orientation(
10209 3, line_orientation[10]);
10211 new_quads[1]->set_line_orientation(
10212 1, line_orientation[1]);
10213 new_quads[1]->set_line_orientation(
10214 2, line_orientation[7]);
10215 new_quads[1]->set_line_orientation(
10216 3, line_orientation[11]);
10218 new_quads[2]->set_line_orientation(
10219 0, line_orientation[4]);
10220 new_quads[2]->set_line_orientation(
10221 1, line_orientation[8]);
10222 new_quads[2]->set_line_orientation(
10223 2, line_orientation[2]);
10225 new_quads[3]->set_line_orientation(
10226 0, line_orientation[5]);
10227 new_quads[3]->set_line_orientation(
10228 1, line_orientation[9]);
10229 new_quads[3]->set_line_orientation(
10230 3, line_orientation[3]);
10264 const int quad_indices[20] = {
10265 new_quads[0]->index(),
10266 new_quads[1]->index(),
10267 new_quads[2]->index(),
10268 new_quads[3]->index(),
10270 hex->face(0)->child_index(
10271 child_at_origin[hex->face(0)->refinement_case() -
10272 1][f_fl[0]][f_ro[0]]),
10273 hex->face(0)->child_index(
10275 child_at_origin[hex->face(0)->refinement_case() -
10276 1][f_fl[0]][f_ro[0]]),
10278 hex->face(1)->child_index(
10279 child_at_origin[hex->face(1)->refinement_case() -
10280 1][f_fl[1]][f_ro[1]]),
10281 hex->face(1)->child_index(
10283 child_at_origin[hex->face(1)->refinement_case() -
10284 1][f_fl[1]][f_ro[1]]),
10286 hex->face(2)->isotropic_child_index(
10288 0, f_or[2], f_fl[2], f_ro[2])),
10289 hex->face(2)->isotropic_child_index(
10291 1, f_or[2], f_fl[2], f_ro[2])),
10292 hex->face(2)->isotropic_child_index(
10294 2, f_or[2], f_fl[2], f_ro[2])),
10295 hex->face(2)->isotropic_child_index(
10297 3, f_or[2], f_fl[2], f_ro[2])),
10299 hex->face(3)->isotropic_child_index(
10301 0, f_or[3], f_fl[3], f_ro[3])),
10302 hex->face(3)->isotropic_child_index(
10304 1, f_or[3], f_fl[3], f_ro[3])),
10305 hex->face(3)->isotropic_child_index(
10307 2, f_or[3], f_fl[3], f_ro[3])),
10308 hex->face(3)->isotropic_child_index(
10310 3, f_or[3], f_fl[3], f_ro[3])),
10312 hex->face(4)->child_index(
10313 child_at_origin[hex->face(4)->refinement_case() -
10314 1][f_fl[4]][f_ro[4]]),
10315 hex->face(4)->child_index(
10317 child_at_origin[hex->face(4)->refinement_case() -
10318 1][f_fl[4]][f_ro[4]]),
10320 hex->face(5)->child_index(
10321 child_at_origin[hex->face(5)->refinement_case() -
10322 1][f_fl[5]][f_ro[5]]),
10323 hex->face(5)->child_index(
10325 child_at_origin[hex->face(5)->refinement_case() -
10326 1][f_fl[5]][f_ro[5]])};
10337 new_hexes[0]->set_bounding_object_indices(
10344 new_hexes[1]->set_bounding_object_indices(
10350 quad_indices[18]});
10351 new_hexes[2]->set_bounding_object_indices(
10358 new_hexes[3]->set_bounding_object_indices(
10364 quad_indices[19]});
10391 new_lines[0]->set_bounding_object_indices(
10393 {middle_vertex_index<dim, spacedim>(hex->face(0)),
10394 middle_vertex_index<dim, spacedim>(hex->face(1))});
10463 spacedim>::raw_line_iterator lines[13] = {
10464 hex->face(2)->child(0)->line(
10465 (hex->face(2)->refinement_case() ==
10469 hex->face(3)->child(0)->line(
10470 (hex->face(3)->refinement_case() ==
10474 hex->face(4)->child(0)->line(
10475 (hex->face(4)->refinement_case() ==
10479 hex->face(5)->child(0)->line(
10480 (hex->face(5)->refinement_case() ==
10488 0, f_or[0], f_fl[0], f_ro[0]))
10491 1, f_or[0], f_fl[0], f_ro[0])),
10495 3, f_or[0], f_fl[0], f_ro[0]))
10498 0, f_or[0], f_fl[0], f_ro[0])),
10502 0, f_or[0], f_fl[0], f_ro[0]))
10505 3, f_or[0], f_fl[0], f_ro[0])),
10509 3, f_or[0], f_fl[0], f_ro[0]))
10512 2, f_or[0], f_fl[0], f_ro[0])),
10517 0, f_or[1], f_fl[1], f_ro[1]))
10520 1, f_or[1], f_fl[1], f_ro[1])),
10524 3, f_or[1], f_fl[1], f_ro[1]))
10527 0, f_or[1], f_fl[1], f_ro[1])),
10531 0, f_or[1], f_fl[1], f_ro[1]))
10534 3, f_or[1], f_fl[1], f_ro[1])),
10538 3, f_or[1], f_fl[1], f_ro[1]))
10541 2, f_or[1], f_fl[1], f_ro[1])),
10546 unsigned int line_indices[13];
10548 for (
unsigned int i = 0; i < 13; ++i)
10549 line_indices[i] = lines[i]->
index();
10556 bool line_orientation[13];
10560 const unsigned int middle_vertices[4] = {
10561 hex->line(8)->child(0)->vertex_index(1),
10562 hex->line(10)->child(0)->vertex_index(1),
10563 hex->line(0)->child(0)->vertex_index(1),
10564 hex->line(4)->child(0)->vertex_index(1),
10569 for (
unsigned int i = 0; i < 4; ++i)
10570 if (lines[i]->vertex_index(0) == middle_vertices[i])
10571 line_orientation[i] =
true;
10575 Assert(lines[i]->vertex_index(1) ==
10576 middle_vertices[i],
10578 line_orientation[i] =
false;
10587 for (
unsigned int i = 4; i < 12; ++i)
10588 if (lines[i]->vertex_index((i + 1) % 2) ==
10589 middle_vertex_index<dim, spacedim>(
10590 hex->face(i / 4 - 1)))
10591 line_orientation[i] =
true;
10596 Assert(lines[i]->vertex_index(i % 2) ==
10597 (middle_vertex_index<dim, spacedim>(
10598 hex->face(i / 4 - 1))),
10600 line_orientation[i] =
false;
10605 line_orientation[12] =
true;
10625 new_quads[0]->set_bounding_object_indices(
10629 line_indices[12]});
10630 new_quads[1]->set_bounding_object_indices(
10635 new_quads[2]->set_bounding_object_indices(
10640 new_quads[3]->set_bounding_object_indices(
10646 new_quads[0]->set_line_orientation(
10647 0, line_orientation[6]);
10648 new_quads[0]->set_line_orientation(
10649 1, line_orientation[10]);
10650 new_quads[0]->set_line_orientation(
10651 2, line_orientation[0]);
10653 new_quads[1]->set_line_orientation(
10654 0, line_orientation[7]);
10655 new_quads[1]->set_line_orientation(
10656 1, line_orientation[11]);
10657 new_quads[1]->set_line_orientation(
10658 3, line_orientation[1]);
10660 new_quads[2]->set_line_orientation(
10661 0, line_orientation[2]);
10662 new_quads[2]->set_line_orientation(
10663 2, line_orientation[4]);
10664 new_quads[2]->set_line_orientation(
10665 3, line_orientation[8]);
10667 new_quads[3]->set_line_orientation(
10668 1, line_orientation[3]);
10669 new_quads[3]->set_line_orientation(
10670 2, line_orientation[5]);
10671 new_quads[3]->set_line_orientation(
10672 3, line_orientation[9]);
10706 const int quad_indices[20] = {
10707 new_quads[0]->index(),
10708 new_quads[1]->index(),
10709 new_quads[2]->index(),
10710 new_quads[3]->index(),
10712 hex->face(0)->isotropic_child_index(
10714 0, f_or[0], f_fl[0], f_ro[0])),
10715 hex->face(0)->isotropic_child_index(
10717 1, f_or[0], f_fl[0], f_ro[0])),
10718 hex->face(0)->isotropic_child_index(
10720 2, f_or[0], f_fl[0], f_ro[0])),
10721 hex->face(0)->isotropic_child_index(
10723 3, f_or[0], f_fl[0], f_ro[0])),
10725 hex->face(1)->isotropic_child_index(
10727 0, f_or[1], f_fl[1], f_ro[1])),
10728 hex->face(1)->isotropic_child_index(
10730 1, f_or[1], f_fl[1], f_ro[1])),
10731 hex->face(1)->isotropic_child_index(
10733 2, f_or[1], f_fl[1], f_ro[1])),
10734 hex->face(1)->isotropic_child_index(
10736 3, f_or[1], f_fl[1], f_ro[1])),
10738 hex->face(2)->child_index(
10739 child_at_origin[hex->face(2)->refinement_case() -
10740 1][f_fl[2]][f_ro[2]]),
10741 hex->face(2)->child_index(
10743 child_at_origin[hex->face(2)->refinement_case() -
10744 1][f_fl[2]][f_ro[2]]),
10746 hex->face(3)->child_index(
10747 child_at_origin[hex->face(3)->refinement_case() -
10748 1][f_fl[3]][f_ro[3]]),
10749 hex->face(3)->child_index(
10751 child_at_origin[hex->face(3)->refinement_case() -
10752 1][f_fl[3]][f_ro[3]]),
10754 hex->face(4)->child_index(
10755 child_at_origin[hex->face(4)->refinement_case() -
10756 1][f_fl[4]][f_ro[4]]),
10757 hex->face(4)->child_index(
10759 child_at_origin[hex->face(4)->refinement_case() -
10760 1][f_fl[4]][f_ro[4]]),
10762 hex->face(5)->child_index(
10763 child_at_origin[hex->face(5)->refinement_case() -
10764 1][f_fl[5]][f_ro[5]]),
10765 hex->face(5)->child_index(
10767 child_at_origin[hex->face(5)->refinement_case() -
10768 1][f_fl[5]][f_ro[5]])};
10770 new_hexes[0]->set_bounding_object_indices(
10777 new_hexes[1]->set_bounding_object_indices(
10784 new_hexes[2]->set_bounding_object_indices(
10790 quad_indices[18]});
10791 new_hexes[3]->set_bounding_object_indices(
10797 quad_indices[19]});
10829 ++next_unused_vertex;
10833 "Internal error: During refinement, the triangulation wants to access an element of the 'vertices' array but it turns out that the array is not large enough."));
10844 hex->center(
true,
true);
10866 middle_vertex_index<dim, spacedim>(hex->face(0)),
10867 middle_vertex_index<dim, spacedim>(hex->face(1)),
10868 middle_vertex_index<dim, spacedim>(hex->face(2)),
10869 middle_vertex_index<dim, spacedim>(hex->face(3)),
10870 middle_vertex_index<dim, spacedim>(hex->face(4)),
10871 middle_vertex_index<dim, spacedim>(hex->face(5)),
10872 next_unused_vertex};
10874 new_lines[0]->set_bounding_object_indices(
10876 new_lines[1]->set_bounding_object_indices(
10878 new_lines[2]->set_bounding_object_indices(
10880 new_lines[3]->set_bounding_object_indices(
10882 new_lines[4]->set_bounding_object_indices(
10884 new_lines[5]->set_bounding_object_indices(
10954 spacedim>::raw_line_iterator lines[30] = {
10958 0, f_or[0], f_fl[0], f_ro[0]))
10961 1, f_or[0], f_fl[0], f_ro[0])),
10965 3, f_or[0], f_fl[0], f_ro[0]))
10968 0, f_or[0], f_fl[0], f_ro[0])),
10972 0, f_or[0], f_fl[0], f_ro[0]))
10975 3, f_or[0], f_fl[0], f_ro[0])),
10979 3, f_or[0], f_fl[0], f_ro[0]))
10982 2, f_or[0], f_fl[0], f_ro[0])),
10987 0, f_or[1], f_fl[1], f_ro[1]))
10990 1, f_or[1], f_fl[1], f_ro[1])),
10994 3, f_or[1], f_fl[1], f_ro[1]))
10997 0, f_or[1], f_fl[1], f_ro[1])),
11001 0, f_or[1], f_fl[1], f_ro[1]))
11004 3, f_or[1], f_fl[1], f_ro[1])),
11008 3, f_or[1], f_fl[1], f_ro[1]))
11011 2, f_or[1], f_fl[1], f_ro[1])),
11016 0, f_or[2], f_fl[2], f_ro[2]))
11019 1, f_or[2], f_fl[2], f_ro[2])),
11023 3, f_or[2], f_fl[2], f_ro[2]))
11026 0, f_or[2], f_fl[2], f_ro[2])),
11030 0, f_or[2], f_fl[2], f_ro[2]))
11033 3, f_or[2], f_fl[2], f_ro[2])),
11037 3, f_or[2], f_fl[2], f_ro[2]))
11040 2, f_or[2], f_fl[2], f_ro[2])),
11045 0, f_or[3], f_fl[3], f_ro[3]))
11048 1, f_or[3], f_fl[3], f_ro[3])),
11052 3, f_or[3], f_fl[3], f_ro[3]))
11055 0, f_or[3], f_fl[3], f_ro[3])),
11059 0, f_or[3], f_fl[3], f_ro[3]))
11062 3, f_or[3], f_fl[3], f_ro[3])),
11066 3, f_or[3], f_fl[3], f_ro[3]))
11069 2, f_or[3], f_fl[3], f_ro[3])),
11074 0, f_or[4], f_fl[4], f_ro[4]))
11077 1, f_or[4], f_fl[4], f_ro[4])),
11081 3, f_or[4], f_fl[4], f_ro[4]))
11084 0, f_or[4], f_fl[4], f_ro[4])),
11088 0, f_or[4], f_fl[4], f_ro[4]))
11091 3, f_or[4], f_fl[4], f_ro[4])),
11095 3, f_or[4], f_fl[4], f_ro[4]))
11098 2, f_or[4], f_fl[4], f_ro[4])),
11103 0, f_or[5], f_fl[5], f_ro[5]))
11106 1, f_or[5], f_fl[5], f_ro[5])),
11110 3, f_or[5], f_fl[5], f_ro[5]))
11113 0, f_or[5], f_fl[5], f_ro[5])),
11117 0, f_or[5], f_fl[5], f_ro[5]))
11120 3, f_or[5], f_fl[5], f_ro[5])),
11124 3, f_or[5], f_fl[5], f_ro[5]))
11127 2, f_or[5], f_fl[5], f_ro[5])),
11137 unsigned int line_indices[30];
11138 for (
unsigned int i = 0; i < 30; ++i)
11139 line_indices[i] = lines[i]->
index();
11146 bool line_orientation[30];
11154 for (
unsigned int i = 0; i < 24; ++i)
11155 if (lines[i]->vertex_index((i + 1) % 2) ==
11157 line_orientation[i] =
true;
11162 Assert(lines[i]->vertex_index(i % 2) ==
11165 line_orientation[i] =
false;
11170 for (
unsigned int i = 24; i < 30; ++i)
11171 line_orientation[i] =
true;
11203 new_quads[0]->set_bounding_object_indices(
11207 line_indices[24]});
11208 new_quads[1]->set_bounding_object_indices(
11212 line_indices[25]});
11213 new_quads[2]->set_bounding_object_indices(
11217 line_indices[20]});
11218 new_quads[3]->set_bounding_object_indices(
11222 line_indices[21]});
11223 new_quads[4]->set_bounding_object_indices(
11227 line_indices[28]});
11228 new_quads[5]->set_bounding_object_indices(
11232 line_indices[29]});
11233 new_quads[6]->set_bounding_object_indices(
11238 new_quads[7]->set_bounding_object_indices(
11243 new_quads[8]->set_bounding_object_indices(
11247 line_indices[26]});
11248 new_quads[9]->set_bounding_object_indices(
11252 line_indices[27]});
11253 new_quads[10]->set_bounding_object_indices(
11257 line_indices[12]});
11258 new_quads[11]->set_bounding_object_indices(
11262 line_indices[13]});
11267 new_quads[0]->set_line_orientation(
11268 0, line_orientation[10]);
11269 new_quads[0]->set_line_orientation(
11270 2, line_orientation[16]);
11272 new_quads[1]->set_line_orientation(
11273 1, line_orientation[14]);
11274 new_quads[1]->set_line_orientation(
11275 2, line_orientation[17]);
11277 new_quads[2]->set_line_orientation(
11278 0, line_orientation[11]);
11279 new_quads[2]->set_line_orientation(
11280 3, line_orientation[20]);
11282 new_quads[3]->set_line_orientation(
11283 1, line_orientation[15]);
11284 new_quads[3]->set_line_orientation(
11285 3, line_orientation[21]);
11287 new_quads[4]->set_line_orientation(
11288 0, line_orientation[18]);
11289 new_quads[4]->set_line_orientation(
11290 2, line_orientation[0]);
11292 new_quads[5]->set_line_orientation(
11293 1, line_orientation[22]);
11294 new_quads[5]->set_line_orientation(
11295 2, line_orientation[1]);
11297 new_quads[6]->set_line_orientation(
11298 0, line_orientation[19]);
11299 new_quads[6]->set_line_orientation(
11300 3, line_orientation[4]);
11302 new_quads[7]->set_line_orientation(
11303 1, line_orientation[23]);
11304 new_quads[7]->set_line_orientation(
11305 3, line_orientation[5]);
11307 new_quads[8]->set_line_orientation(
11308 0, line_orientation[2]);
11309 new_quads[8]->set_line_orientation(
11310 2, line_orientation[8]);
11312 new_quads[9]->set_line_orientation(
11313 1, line_orientation[6]);
11314 new_quads[9]->set_line_orientation(
11315 2, line_orientation[9]);
11317 new_quads[10]->set_line_orientation(
11318 0, line_orientation[3]);
11319 new_quads[10]->set_line_orientation(
11320 3, line_orientation[12]);
11322 new_quads[11]->set_line_orientation(
11323 1, line_orientation[7]);
11324 new_quads[11]->set_line_orientation(
11325 3, line_orientation[13]);
11366 const int quad_indices[36] = {
11367 new_quads[0]->index(),
11368 new_quads[1]->index(),
11369 new_quads[2]->index(),
11370 new_quads[3]->index(),
11371 new_quads[4]->index(),
11372 new_quads[5]->index(),
11373 new_quads[6]->index(),
11374 new_quads[7]->index(),
11375 new_quads[8]->index(),
11376 new_quads[9]->index(),
11377 new_quads[10]->index(),
11378 new_quads[11]->index(),
11380 hex->face(0)->isotropic_child_index(
11382 0, f_or[0], f_fl[0], f_ro[0])),
11383 hex->face(0)->isotropic_child_index(
11385 1, f_or[0], f_fl[0], f_ro[0])),
11386 hex->face(0)->isotropic_child_index(
11388 2, f_or[0], f_fl[0], f_ro[0])),
11389 hex->face(0)->isotropic_child_index(
11391 3, f_or[0], f_fl[0], f_ro[0])),
11393 hex->face(1)->isotropic_child_index(
11395 0, f_or[1], f_fl[1], f_ro[1])),
11396 hex->face(1)->isotropic_child_index(
11398 1, f_or[1], f_fl[1], f_ro[1])),
11399 hex->face(1)->isotropic_child_index(
11401 2, f_or[1], f_fl[1], f_ro[1])),
11402 hex->face(1)->isotropic_child_index(
11404 3, f_or[1], f_fl[1], f_ro[1])),
11406 hex->face(2)->isotropic_child_index(
11408 0, f_or[2], f_fl[2], f_ro[2])),
11409 hex->face(2)->isotropic_child_index(
11411 1, f_or[2], f_fl[2], f_ro[2])),
11412 hex->face(2)->isotropic_child_index(
11414 2, f_or[2], f_fl[2], f_ro[2])),
11415 hex->face(2)->isotropic_child_index(
11417 3, f_or[2], f_fl[2], f_ro[2])),
11419 hex->face(3)->isotropic_child_index(
11421 0, f_or[3], f_fl[3], f_ro[3])),
11422 hex->face(3)->isotropic_child_index(
11424 1, f_or[3], f_fl[3], f_ro[3])),
11425 hex->face(3)->isotropic_child_index(
11427 2, f_or[3], f_fl[3], f_ro[3])),
11428 hex->face(3)->isotropic_child_index(
11430 3, f_or[3], f_fl[3], f_ro[3])),
11432 hex->face(4)->isotropic_child_index(
11434 0, f_or[4], f_fl[4], f_ro[4])),
11435 hex->face(4)->isotropic_child_index(
11437 1, f_or[4], f_fl[4], f_ro[4])),
11438 hex->face(4)->isotropic_child_index(
11440 2, f_or[4], f_fl[4], f_ro[4])),
11441 hex->face(4)->isotropic_child_index(
11443 3, f_or[4], f_fl[4], f_ro[4])),
11445 hex->face(5)->isotropic_child_index(
11447 0, f_or[5], f_fl[5], f_ro[5])),
11448 hex->face(5)->isotropic_child_index(
11450 1, f_or[5], f_fl[5], f_ro[5])),
11451 hex->face(5)->isotropic_child_index(
11453 2, f_or[5], f_fl[5], f_ro[5])),
11454 hex->face(5)->isotropic_child_index(
11456 3, f_or[5], f_fl[5], f_ro[5]))};
11459 new_hexes[0]->set_bounding_object_indices(
11466 new_hexes[1]->set_bounding_object_indices(
11473 new_hexes[2]->set_bounding_object_indices(
11479 quad_indices[10]});
11480 new_hexes[3]->set_bounding_object_indices(
11486 quad_indices[11]});
11489 new_hexes[4]->set_bounding_object_indices(
11495 quad_indices[32]});
11496 new_hexes[5]->set_bounding_object_indices(
11502 quad_indices[33]});
11503 new_hexes[6]->set_bounding_object_indices(
11509 quad_indices[34]});
11510 new_hexes[7]->set_bounding_object_indices(
11516 quad_indices[35]});
11550 for (
unsigned int s = 0;
11557 const unsigned int current_child =
11566 ref_case, f, f_or[f], f_fl[f], f_ro[f]));
11567 new_hexes[current_child]->set_combined_face_orientation(
11573 if (check_for_distorted_cells &&
11574 has_distorted_children<dim, spacedim>(hex))
11594 return cells_with_distorted_children;