6391 const bool check_for_distorted_cells)
6393 constexpr int dim = 3;
6395 using raw_line_iterator =
6397 using raw_quad_iterator =
6409 triangulation.
levels.size() - 1))
6410 if (cell->refine_flag_set())
6415 triangulation.
levels.push_back(
6418 triangulation.
strides.max_children_per_cell,
6419 triangulation.
strides.max_faces_per_cell,
6420 triangulation.
strides.max_vertices_per_cell));
6426 triangulation.
faces->quads.clear_user_data();
6427 triangulation.
faces->lines.clear_user_flags();
6428 triangulation.
faces->quads.clear_user_flags();
6443 unsigned int needed_vertices = 0;
6444 unsigned int needed_lines_single = 0;
6445 unsigned int needed_faces_single = 0;
6446 unsigned int needed_lines_pair = 0;
6447 unsigned int needed_faces_pair = 0;
6448 bool contains_pyramid =
false;
6449 for (
int level_no = triangulation.
levels.size() - 2; level_no >= 0;
6452 std::size_t needed_cells = 0;
6454 for (
const auto &cell :
6456 if (cell->refine_flag_set())
6459 Assert(cell->refine_flag_set() ==
6464 const auto cell_reference_cell = cell->reference_cell();
6468 needed_cells += cell_reference_cell.n_isotropic_children();
6470 switch (cell_reference_cell)
6473 needed_lines_single += 1;
6474 needed_faces_single += 8;
6482 needed_lines_single += 4;
6483 needed_faces_single += 13;
6484 contains_pyramid =
true;
6497 needed_lines_single += 3;
6498 needed_faces_single += 10;
6504 needed_lines_single += 6;
6505 needed_faces_single += 12;
6518 for (
const auto face : cell->face_indices())
6519 if (cell->face(face)->n_children() == 0)
6520 cell->face(face)->set_user_flag();
6523 Assert(cell->face(face)->n_children() ==
6524 cell_reference_cell.face_reference_cell(face)
6525 .n_isotropic_children(),
6529 for (
const auto line : cell->line_indices())
6530 if (cell->line(line)->has_children() ==
false)
6531 cell->line(line)->set_user_flag();
6533 Assert(cell->line(line)->n_children() == 2,
6537 auto &next_level = *triangulation.
levels[level_no + 1];
6538 const auto used_cells = std::count(next_level.cells.used.begin(),
6539 next_level.cells.used.end(),
6541 if (used_cells + needed_cells > next_level.size())
6542 next_level.allocate_end(
6543 (used_cells + needed_cells) - next_level.size(),
6553 next_level.cells.allocate_end(needed_cells, 0);
6562 if (contains_pyramid)
6580 if (quad->user_flag_set() ==
false)
6585 needed_faces_pair += 4;
6586 needed_lines_pair += 4;
6587 needed_vertices += 1;
6591 needed_faces_pair += 4;
6592 needed_lines_single += 3;
6605 if (line->user_flag_set() ==
false)
6608 needed_lines_pair += 2;
6609 needed_vertices += 1;
6612 triangulation.
faces->allocate_end(needed_faces_pair,
6613 needed_faces_single);
6614 triangulation.
faces->lines.allocate_end(needed_lines_pair,
6615 needed_lines_single);
6618 needed_vertices += std::count(triangulation.
vertices_used.begin(),
6622 if (needed_vertices > triangulation.
vertices.size())
6642 if (!cell->refine_flag_set())
6643 for (
unsigned int line_n = 0; line_n < cell->n_lines();
6645 if (cell->line(line_n)->has_children())
6646 for (
unsigned int c = 0; c < 2; ++c)
6647 Assert(cell->line(line_n)->child(c)->user_flag_set() ==
6652 unsigned int current_vertex = 0;
6656 auto get_next_unused_vertex = [](
const unsigned int current_vertex,
6657 std::vector<bool> &vertices_used) {
6658 unsigned int next_vertex = current_vertex;
6659 while (next_vertex < vertices_used.size() &&
6660 vertices_used[next_vertex] ==
true)
6663 vertices_used[next_vertex] =
true;
6673 raw_line_iterator next_unused_line = triangulation.
begin_raw_line();
6675 for (; line != endl; ++line)
6677 if (line->user_flag_set() ==
false)
6681 triangulation.
faces->lines.template next_free_pair_object<1>(
6689 line->set_children(0, next_unused_line->index());
6691 const std::array<raw_line_iterator, 2> children{
6692 {next_unused_line, ++next_unused_line}};
6698 get_next_unused_vertex(current_vertex,
6700 triangulation.
vertices[current_vertex] = line->center(
true);
6702 children[0]->set_bounding_object_indices(
6703 {line->vertex_index(0), current_vertex});
6704 children[1]->set_bounding_object_indices(
6705 {current_vertex, line->vertex_index(1)});
6707 const auto manifold_id = line->manifold_id();
6708 const auto boundary_id = line->boundary_id();
6709 for (
const auto &child : children)
6711 child->set_used_flag();
6712 child->clear_children();
6714 child->clear_user_flag();
6715 child->set_boundary_id_internal(boundary_id);
6716 child->set_manifold_id(manifold_id);
6719 line->clear_user_flag();
6729 for (; face != endf; ++face)
6731 if (face->user_flag_set() ==
false)
6734 const auto reference_face_type = face->reference_cell();
6735 const auto n_lines_per_face = reference_face_type.n_lines();
6739 std::array<raw_line_iterator, 4> new_lines;
6740 switch (reference_face_type)
6743 for (
unsigned int l = 0; l < 3; ++l)
6745 triangulation.
faces->lines
6746 .template next_free_single_object<1>(triangulation);
6750 for (
unsigned int l = 0; l < 2; ++l)
6752 auto next_unused_line =
6753 triangulation.
faces->lines
6754 .template next_free_pair_object<1>(triangulation);
6755 new_lines[2 * l] = next_unused_line;
6756 new_lines[2 * l + 1] = ++next_unused_line;
6766 for (
const unsigned int line : face->line_indices())
6775 std::array<raw_quad_iterator, 4> new_faces;
6776 for (
unsigned int f = 0; f < 2; ++f)
6778 auto next_unused_quad =
6779 triangulation.
faces->quads
6780 .template next_free_pair_object<2>(triangulation);
6782 new_faces[2 * f] = next_unused_quad;
6783 new_faces[2 * f + 1] = ++next_unused_quad;
6785 face->set_children(2 * f, new_faces[2 * f]->
index());
6792 for (
const auto &quad : new_faces)
6807 for (
const auto i : face->vertex_indices())
6810 for (
unsigned int l = 0; l < n_lines_per_face; ++l)
6816 get_next_unused_vertex(current_vertex,
6820 triangulation.
vertices[current_vertex] =
6821 face->center(
true,
true);
6827 std::array<raw_line_iterator, 12> lines;
6828 unsigned int n_lines = 0;
6829 for (
unsigned int l = 0; l < n_lines_per_face; ++l)
6830 for (
unsigned int c = 0; c < 2; ++c)
6831 lines[n_lines++] = face->line(l)->child(
6832 child_line_index(c, face->line_orientation(l)));
6834 for (
unsigned int l = 0; l < n_lines_per_face; ++l)
6835 lines[n_lines++] = new_lines[l];
6837 std::array<int, 12> line_indices;
6838 for (
unsigned int i = 0; i < n_lines; ++i)
6839 line_indices[i] = lines[i]->
index();
6846 constexpr ::ndarray<unsigned int, 12, 2> line_vertices_quad{
6860 constexpr ::ndarray<unsigned int, 4, 4> quad_lines_quad{
6875 constexpr ::ndarray<unsigned int, 12, 2> line_vertices_tri{
6889 constexpr ::ndarray<unsigned int, 4, 4> tri_lines_tri{
6899 constexpr ::ndarray<unsigned int, 4, 4, 2>
6900 tri_line_vertices_tri{
6901 {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
6902 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
6903 {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
6904 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
6907 const auto &line_vertices =
6909 line_vertices_quad :
6911 const auto &face_lines =
6919 for (
unsigned int i = 0, j = 2 * n_lines_per_face;
6920 i < n_lines_per_face;
6923 auto &new_line = new_lines[i];
6924 new_line->set_bounding_object_indices(
6927 new_line->set_used_flag();
6928 new_line->clear_user_flag();
6929 new_line->clear_user_data();
6930 new_line->clear_children();
6931 new_line->set_boundary_id_internal(face->boundary_id());
6932 new_line->set_manifold_id(face->manifold_id());
6936 for (
unsigned int i = 0; i < new_faces.size(); ++i)
6939 auto &new_face = new_faces[i];
6943 triangulation.
faces->set_quad_type(new_face->index(),
6944 reference_face_type);
6946 switch (reference_face_type)
6949 new_face->set_bounding_object_indices(
6950 {line_indices[face_lines[i][0]],
6951 line_indices[face_lines[i][1]],
6952 line_indices[face_lines[i][2]]});
6956 new_face->set_bounding_object_indices(
6957 {line_indices[face_lines[i][0]],
6958 line_indices[face_lines[i][1]],
6959 line_indices[face_lines[i][2]],
6960 line_indices[face_lines[i][3]]});
6967 new_face->set_used_flag();
6968 new_face->clear_user_flag();
6969 new_face->clear_user_data();
6970 new_face->clear_children();
6971 new_face->set_boundary_id_internal(face->boundary_id());
6972 new_face->set_manifold_id(face->manifold_id());
6974 [[maybe_unused]] std::set<unsigned int> s;
6982 for (
const auto f : new_face->line_indices())
6984 const std::array<unsigned int, 2> vertices_0 = {
6985 {lines[face_lines[i][f]]->vertex_index(0),
6986 lines[face_lines[i][f]]->vertex_index(1)}};
6988 const std::array<unsigned int, 2> vertices_1 = {
6992 const auto orientation =
7000 for (
const auto i : vertices_0)
7002 for (
const auto i : vertices_1)
7006 new_face->set_line_orientation(f, orientation);
7020 constexpr ::ndarray<unsigned int, 4, 2>
7021 quad_child_line_vertices{{{{0, 2}},
7026 for (
unsigned int i = 0; i < 4; ++i)
7027 for (
unsigned int j = 0; j < 2; ++j)
7028 new_faces[quad_child_line_vertices[i][j]]
7029 ->set_line_orientation(i, face->line_orientation(i));
7032 face->clear_user_flag();
7038 cells_with_distorted_children;
7048 cell->level() >=
static_cast<int>(
level),
7052 for (; cell != triangulation.
end() &&
7053 cell->level() ==
static_cast<int>(
level);
7056 if (cell->refine_flag_set() ==
7063 cell->clear_refine_flag();
7064 cell->set_refinement_case(ref_case);
7066 const auto cell_reference_cell = cell->reference_cell();
7068 unsigned int n_new_lines = 0;
7069 unsigned int n_new_faces = 0;
7070 const unsigned int n_new_cells =
7071 cell_reference_cell.n_isotropic_children();
7075 switch (cell_reference_cell)
7103 constexpr unsigned int max_n_new_lines = 6;
7104 constexpr unsigned int max_n_new_faces = 13;
7105 constexpr unsigned int max_n_new_children = 10;
7106 constexpr unsigned int max_relevant_vertices = 27;
7107 constexpr unsigned int max_relevant_lines = 30;
7108 constexpr unsigned int max_face_indices = 36;
7110 std::array<raw_line_iterator, max_n_new_lines> new_lines;
7111 for (
unsigned int i = 0; i < n_new_lines; ++i)
7114 triangulation.
faces->lines
7115 .template next_free_single_object<1>(triangulation);
7118 new_lines[i]->set_used_flag();
7119 new_lines[i]->clear_user_flag();
7120 new_lines[i]->clear_user_data();
7121 new_lines[i]->clear_children();
7122 new_lines[i]->set_boundary_id_internal(
7124 new_lines[i]->set_manifold_id(cell->manifold_id());
7127 std::array<raw_quad_iterator, max_n_new_faces> new_faces;
7128 for (
unsigned int i = 0; i < n_new_faces; ++i)
7131 triangulation.
faces->quads
7132 .template next_free_single_object<2>(triangulation);
7134 auto &new_face = new_faces[i];
7137 new_face->set_used_flag();
7138 new_face->clear_user_flag();
7139 new_face->clear_user_data();
7140 new_face->clear_children();
7141 new_face->set_boundary_id_internal(
7143 new_face->set_manifold_id(cell->manifold_id());
7148 for (
unsigned int j = 0;
7149 j < triangulation.
faces->lines_per_quad;
7151 new_face->set_line_orientation(
7204 std::array<unsigned int, max_relevant_vertices>
7213 const unsigned int n_vertices =
7217 for (
unsigned int i = 0; i < n_vertices; ++i)
7224 const unsigned int n_lines =
std::min(cell->n_lines(), 12u);
7226 const std::array<unsigned int, 12> line_indices =
7227 TriaAccessorImplementation::Implementation::
7228 get_line_indices_of_cell(*cell);
7230 for (
unsigned int l = 0; l < n_lines; ++l)
7232 raw_line_iterator line(&triangulation,
7240 for (
auto quad_face_index :
7241 cell_reference_cell.face_indices_by_type(
7254 get_next_unused_vertex(current_vertex,
7258 triangulation.
vertices[current_vertex] =
7259 cell->center(
true,
true);
7264 unsigned int chosen_line_tetrahedron = 0;
7267 switch (cell_reference_cell)
7275 constexpr ::ndarray<unsigned int, 3, 2>
7276 new_line_vertices = {{{{6, 8}},
7282 std::uint8_t refinement_choice =
7283 cell->refine_choice();
7284 if (refinement_choice ==
7288 const auto &vertices =
7290 double min_distance =
7291 std::numeric_limits<double>::infinity();
7292 for (
unsigned int i = 0;
7293 i < new_line_vertices.size();
7296 const double current_distance =
7298 [new_line_vertices[i][0]]]
7301 [new_line_vertices[i][1]]]);
7302 if (current_distance < min_distance)
7304 chosen_line_tetrahedron = i;
7305 min_distance = current_distance;
7309 else if (refinement_choice ==
7312 chosen_line_tetrahedron = 0;
7313 else if (refinement_choice ==
7316 chosen_line_tetrahedron = 1;
7317 else if (refinement_choice ==
7320 chosen_line_tetrahedron = 2;
7324 cell->set_refinement_case(
7327 new_lines[0]->set_bounding_object_indices(
7329 [new_line_vertices[chosen_line_tetrahedron][0]],
7331 [chosen_line_tetrahedron][1]]});
7340 static constexpr ::ndarray<unsigned int, 4, 2>
7341 new_line_vertices = {{{{13, 9}},
7346 for (
unsigned int i = 0; i < n_new_lines; ++i)
7347 new_lines[i]->set_bounding_object_indices(
7357 constexpr ::ndarray<unsigned int, 3, 2>
7358 new_line_vertices = {{{{15, 16}},
7362 for (
unsigned int i = 0; i < n_new_lines; ++i)
7363 new_lines[i]->set_bounding_object_indices(
7371 constexpr ::ndarray<unsigned int, 6, 2>
7372 new_line_vertices = {{{{22, 26}},
7378 for (
unsigned int i = 0; i < n_new_lines; ++i)
7379 new_lines[i]->set_bounding_object_indices(
7392 auto standard_to_real_quad_child_index =
7393 [](
const unsigned int child,
7404 switch (face_ref_cell_type)
7462 std::array<raw_line_iterator, max_relevant_lines>
7464 std::array<int, max_relevant_lines> relevant_line_indices;
7465 unsigned int relevant_lines_counter = 0;
7467 for (
auto f : cell_reference_cell.face_indices_by_type(
7487 constexpr ::ndarray<unsigned int, 4, 2>
7488 quad_relevant_children_and_lines = {{{{0, 1}},
7492 for (
unsigned int l = 0; l < 4;
7493 ++l, ++relevant_lines_counter)
7497 const unsigned int relevant_child =
7498 standard_to_real_quad_child_index(
7499 quad_relevant_children_and_lines[l][0],
7500 cell->combined_face_orientation(f),
7506 const unsigned int relevant_child_line =
7508 .standard_to_real_face_line(
7509 quad_relevant_children_and_lines[l][1],
7512 cell->combined_face_orientation(f));
7516 relevant_lines[relevant_lines_counter] =
7518 ->isotropic_child(relevant_child)
7519 ->line(relevant_child_line);
7524 for (
auto f : cell_reference_cell.face_indices_by_type(
7527 for (
unsigned int l = 0; l < 3;
7528 ++l, ++relevant_lines_counter)
7533 const std::array<std::array<unsigned int, 3>, 6>
7544 tri_line_perm = {{{{0, 1, 2}},
7551 const auto combined_orientation =
7552 cell->combined_face_orientation(f);
7553 relevant_lines[relevant_lines_counter] =
7556 ->line(tri_line_perm[combined_orientation][l]);
7561 for (
unsigned int i = 0; i < n_new_lines;
7562 ++i, ++relevant_lines_counter)
7563 relevant_lines[relevant_lines_counter] = new_lines[i];
7567 for (
unsigned int i = 0; i < relevant_lines_counter; ++i)
7568 relevant_line_indices[i] = relevant_lines[i]->
index();
7577 const auto &new_face_lines =
7578 cell_reference_cell.new_isotropic_child_face_lines(
7579 chosen_line_tetrahedron);
7587 const auto &new_face_lines_vert =
7589 .new_isotropic_child_face_line_vertices(
7590 chosen_line_tetrahedron);
7596 constexpr ::ndarray<unsigned int, 4, 2>
7597 representative_lines{{{{0, 2}},
7602 for (
unsigned int q = 0; q < n_new_faces; ++q)
7604 auto &new_face = new_faces[q];
7606 if (new_face_lines[q][3] == X)
7610 triangulation.
faces->set_quad_type(
7613 new_face->set_bounding_object_indices(
7614 {relevant_line_indices[new_face_lines[q][0]],
7615 relevant_line_indices[new_face_lines[q][1]],
7616 relevant_line_indices[new_face_lines[q][2]]});
7620 triangulation.
faces->set_quad_type(
7623 new_face->set_bounding_object_indices(
7624 {relevant_line_indices[new_face_lines[q][0]],
7625 relevant_line_indices[new_face_lines[q][1]],
7626 relevant_line_indices[new_face_lines[q][2]],
7627 relevant_line_indices[new_face_lines[q][3]]});
7634 const unsigned int n_compute_lines =
7637 new_face->n_lines();
7638 for (
unsigned int line = 0; line < n_compute_lines;
7641 const unsigned int l =
7642 (cell_reference_cell ==
7644 representative_lines[q % 4][0] :
7647 const std::array<unsigned int, 2> vertices_0 = {
7648 {relevant_lines[new_face_lines[q][l]]
7650 relevant_lines[new_face_lines[q][l]]
7651 ->vertex_index(1)}};
7653 const std::array<unsigned int, 2> vertices_1 = {
7657 const auto orientation =
7662 new_face->set_line_orientation(l, orientation);
7667 if (cell_reference_cell ==
7669 new_faces[representative_lines[q % 4][1] + q -
7671 ->set_line_orientation(l, orientation);
7693 std::array<int, max_face_indices> face_indices;
7695 unsigned int face_indices_counter = 0;
7698 for (; face_indices_counter < n_new_faces;
7699 ++face_indices_counter)
7700 face_indices[face_indices_counter] =
7701 new_faces[face_indices_counter]->
index();
7704 for (
auto f_q : cell_reference_cell.face_indices_by_type(
7707 for (
unsigned int c = 0; c < 4;
7708 ++c, ++face_indices_counter)
7710 const auto combined_orientation =
7711 cell->combined_face_orientation(f_q);
7714 face_indices[face_indices_counter] =
7715 cell->face(f_q)->isotropic_child_index(
7716 standard_to_real_quad_child_index(
7718 combined_orientation,
7724 for (
auto f_t : cell_reference_cell.face_indices_by_type(
7727 for (
unsigned int c = 0; c < 4;
7728 ++c, ++face_indices_counter)
7730 const auto combined_orientation =
7731 cell->combined_face_orientation(f_t);
7734 face_indices[face_indices_counter] =
7735 cell->face(f_t)->child_index(
7736 standard_to_real_quad_child_index(
7738 combined_orientation,
7747 const auto cell_faces =
7748 cell_reference_cell.new_isotropic_child_cell_faces(
7749 chosen_line_tetrahedron);
7755 const auto new_cells_vertices =
7756 cell_reference_cell.new_isotropic_child_cell_vertices(
7757 chosen_line_tetrahedron);
7762 const ::types::subdomain_id cell_subdomain_id =
7763 cell->subdomain_id();
7770 for (
unsigned int c = 0; c < n_new_cells; ++c)
7778 ->cells.next_free_hex(triangulation,
level + 1);
7782 new_cells[c] = next_unused_cell;
7784 auto &new_cell = new_cells[c];
7789 new_cell->set_used_flag();
7790 new_cell->clear_user_flag();
7791 new_cell->clear_user_data();
7792 new_cell->clear_children();
7793 new_cell->set_material_id(cell->material_id());
7794 new_cell->set_manifold_id(cell->manifold_id());
7795 new_cell->set_subdomain_id(cell_subdomain_id);
7801 new_cell->set_parent(cell->index());
7802 cell->set_children(c, new_cell->index());
7815 triangulation.
levels[new_cell->level()]
7816 ->reference_cell[new_cell->index()] =
7817 child_reference_cell;
7823 switch (child_reference_cell)
7826 new_cell->set_bounding_object_indices(
7827 {face_indices[cell_faces[c][0]],
7828 face_indices[cell_faces[c][1]],
7829 face_indices[cell_faces[c][2]],
7830 face_indices[cell_faces[c][3]]});
7836 new_cell->set_bounding_object_indices(
7837 {face_indices[cell_faces[c][0]],
7838 face_indices[cell_faces[c][1]],
7839 face_indices[cell_faces[c][2]],
7840 face_indices[cell_faces[c][3]],
7841 face_indices[cell_faces[c][4]]});
7845 new_cell->set_bounding_object_indices(
7846 {face_indices[cell_faces[c][0]],
7847 face_indices[cell_faces[c][1]],
7848 face_indices[cell_faces[c][2]],
7849 face_indices[cell_faces[c][3]],
7850 face_indices[cell_faces[c][4]],
7851 face_indices[cell_faces[c][5]]});
7870 for (
const auto f : new_cell->face_indices())
7871 new_cell->set_combined_face_orientation(
7876 for (
const auto f : new_cell->face_indices())
7878 const auto &face = new_cell->face(f);
7881 auto new_cell_vertices = new_cells_vertices[c];
7884 std::array<unsigned int, 4> vertices_0,
7889 for (
unsigned int face_vertex_no :
7890 face->vertex_indices())
7895 const auto cell_vertex_no =
7896 child_reference_cell
7901 default_geometric_orientation);
7915 [new_cell_vertices[cell_vertex_no]];
7919 for (
const auto i : face->vertex_indices())
7920 vertices_1[i] = face->vertex_index(i);
7925 new_cell->set_combined_face_orientation(
7927 face->reference_cell()
7928 .get_combined_orientation(
7930 vertices_0.cbegin() +
7931 face->n_vertices()),
7933 vertices_1.cbegin() +
7934 face->n_vertices())));
7943 constexpr ::ndarray<unsigned int, 6, 4>
7944 face_to_child_indices_hex{{{{0, 2, 4, 6}},
7951 for (
const auto f : cell->face_indices())
7953 const auto combined_orientation =
7954 cell->combined_face_orientation(f);
7955 for (
unsigned int c = 0; c < 4; ++c)
7956 new_cells[face_to_child_indices_hex[f][c]]
7957 ->set_combined_face_orientation(
7958 f, combined_orientation);
7964 if (check_for_distorted_cells &&
7965 has_distorted_children<dim, spacedim>(cell))
7972 triangulation.
faces->quads.clear_user_data();
7974 return cells_with_distorted_children;
7984 const bool check_for_distorted_cells)
7986 const unsigned int dim = 3;
7989 bool flag_isotropic_mesh =
true;
7991 cell = triangulation.
begin(),
7992 endc = triangulation.
end();
7993 for (; cell != endc; ++cell)
8003 flag_isotropic_mesh =
false;
8007 if (flag_isotropic_mesh)
8008 return execute_refinement_isotropic(triangulation,
8009 check_for_distorted_cells);
8021 triangulation.
levels.size() - 1))
8022 if (cell->refine_flag_set())
8027 triangulation.
levels.push_back(
8030 triangulation.
strides.max_children_per_cell,
8031 triangulation.
strides.max_faces_per_cell,
8032 triangulation.
strides.max_vertices_per_cell));
8039 triangulation.
faces->quads.clear_user_data();
8045 line->clear_user_flag();
8050 quad->clear_user_flag();
8073 unsigned int needed_vertices = 0;
8074 unsigned int needed_lines_single = 0;
8075 unsigned int needed_quads_single = 0;
8076 unsigned int needed_lines_pair = 0;
8077 unsigned int needed_quads_pair = 0;
8078 for (
int level_no = triangulation.
levels.size() - 2; level_no >= 0;
8083 std::size_t needed_cells = 0;
8085 for (
const auto &acell :
8087 if (acell->refine_flag_set())
8097 ++needed_quads_single;
8105 ++needed_lines_single;
8106 needed_quads_single += 4;
8113 needed_lines_single += 6;
8114 needed_quads_single += 12;
8128 for (
const unsigned int face :
8132 aface = acell->face(face);
8139 acell->face_orientation(face),
8140 acell->face_flip(face),
8141 acell->face_rotation(face));
8146 if (face_ref_case ==
8149 if (aface->n_active_descendants() < 4)
8152 aface->set_user_flag();
8154 else if (aface->refinement_case() != face_ref_case)
8165 Assert(aface->refinement_case() ==
8168 aface->refinement_case() ==
8171 aface->set_user_index(face_ref_case);
8177 for (
unsigned int line = 0;
8178 line < GeometryInfo<dim>::lines_per_cell;
8182 !acell->line(line)->has_children())
8183 acell->line(line)->set_user_flag();
8188 auto &next_level = *triangulation.
levels[level_no + 1];
8189 const auto used_cells = std::count(next_level.cells.used.begin(),
8190 next_level.cells.used.end(),
8192 if (used_cells + needed_cells > next_level.size())
8193 next_level.allocate_end(
8194 (used_cells + needed_cells) - next_level.size(),
8204 next_level.cells.allocate_end(needed_cells, 0);
8213 if (quad->user_flag_set())
8219 needed_quads_pair += 4;
8220 needed_lines_pair += 4;
8221 needed_vertices += 1;
8223 if (quad->user_index())
8227 needed_quads_pair += 2;
8228 needed_lines_single += 1;
8241 if (quad->has_children())
8243 Assert(quad->refinement_case() ==
8246 if ((face_refinement_cases[quad->user_index()] ==
8248 (quad->child(0)->line_index(1) + 1 !=
8249 quad->child(2)->line_index(1))) ||
8250 (face_refinement_cases[quad->user_index()] ==
8252 (quad->child(0)->line_index(3) + 1 !=
8253 quad->child(1)->line_index(3))))
8254 needed_lines_pair += 2;
8263 if (line->user_flag_set())
8265 needed_lines_pair += 2;
8266 needed_vertices += 1;
8270 triangulation.
faces->lines.allocate_end(needed_lines_pair,
8271 needed_lines_single);
8274 triangulation.
faces->allocate_end(needed_quads_pair,
8275 needed_quads_single);
8278 needed_vertices += std::count(triangulation.
vertices_used.begin(),
8284 if (needed_vertices > triangulation.
vertices.size())
8304 if (!cell->refine_flag_set())
8305 for (
unsigned int line = 0;
8306 line < GeometryInfo<dim>::lines_per_cell;
8308 if (cell->line(line)->has_children())
8309 for (
unsigned int c = 0; c < 2; ++c)
8310 Assert(cell->line(line)->child(c)->user_flag_set() ==
8323 unsigned int next_unused_vertex = 0;
8334 for (; line != endl; ++line)
8335 if (line->user_flag_set())
8341 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
8342 ++next_unused_vertex;
8344 next_unused_vertex < triangulation.
vertices.size(),
8346 "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."));
8349 triangulation.
vertices[next_unused_vertex] = line->center(
true);
8355 triangulation.
faces->lines.template next_free_pair_object<1>(
8363 line->set_children(0, next_unused_line->index());
8367 children[2] = {next_unused_line, ++next_unused_line};
8374 children[0]->set_bounding_object_indices(
8375 {line->vertex_index(0), next_unused_vertex});
8376 children[1]->set_bounding_object_indices(
8377 {next_unused_vertex, line->vertex_index(1)});
8379 children[0]->set_used_flag();
8380 children[1]->set_used_flag();
8381 children[0]->clear_children();
8382 children[1]->clear_children();
8385 children[0]->clear_user_flag();
8386 children[1]->clear_user_flag();
8388 children[0]->set_boundary_id_internal(line->boundary_id());
8389 children[1]->set_boundary_id_internal(line->boundary_id());
8391 children[0]->set_manifold_id(line->manifold_id());
8392 children[1]->set_manifold_id(line->manifold_id());
8397 line->clear_user_flag();
8427 for (
unsigned int loop = 0; loop < 2; ++loop)
8440 for (; quad != endq; ++quad)
8442 if (quad->user_index())
8445 face_refinement_cases[quad->user_index()];
8454 if (aniso_quad_ref_case == quad->refinement_case())
8457 Assert(quad->refinement_case() ==
8459 quad->refinement_case() ==
8464 Assert(quad->user_index() ==
8466 quad->user_index() ==
8475 triangulation.
faces->lines
8476 .template next_free_single_object<1>(triangulation);
8496 quad->line(2)->child(0)->vertex_index(1);
8498 quad->line(3)->child(0)->vertex_index(1);
8503 quad->line(0)->child(0)->vertex_index(1);
8505 quad->line(1)->child(0)->vertex_index(1);
8508 new_line->set_bounding_object_indices(
8510 new_line->set_used_flag();
8511 new_line->clear_user_flag();
8513 new_line->clear_children();
8514 new_line->set_boundary_id_internal(quad->boundary_id());
8515 new_line->set_manifold_id(quad->manifold_id());
8523 triangulation.
faces->quads
8524 .template next_free_pair_object<2>(triangulation);
8525 new_quads[0] = next_unused_quad;
8529 new_quads[1] = next_unused_quad;
8534 new_quads[0]->set_bounding_object_indices(
8535 {
static_cast<int>(quad->line_index(0)),
8539 child_line_index(0, quad->line_orientation(2)))
8543 child_line_index(0, quad->line_orientation(3)))
8545 new_quads[1]->set_bounding_object_indices(
8547 static_cast<int>(quad->line_index(1)),
8550 child_line_index(1, quad->line_orientation(2)))
8554 child_line_index(1, quad->line_orientation(3)))
8559 new_quads[0]->set_bounding_object_indices(
8562 child_line_index(0, quad->line_orientation(0)))
8566 child_line_index(0, quad->line_orientation(1)))
8568 static_cast<int>(quad->line_index(2)),
8569 new_line->index()});
8570 new_quads[1]->set_bounding_object_indices(
8573 child_line_index(1, quad->line_orientation(0)))
8577 child_line_index(1, quad->line_orientation(1)))
8580 static_cast<int>(quad->line_index(3))});
8583 for (
const auto &new_quad : new_quads)
8585 new_quad->set_used_flag();
8586 new_quad->clear_user_flag();
8588 new_quad->clear_children();
8589 new_quad->set_boundary_id_internal(quad->boundary_id());
8590 new_quad->set_manifold_id(quad->manifold_id());
8594 for (
unsigned int j = 0;
8595 j < GeometryInfo<dim>::lines_per_face;
8597 new_quad->set_line_orientation(
8604 new_quads[0]->set_line_orientation(
8605 0, quad->line_orientation(0));
8606 new_quads[0]->set_line_orientation(
8607 2, quad->line_orientation(2));
8608 new_quads[1]->set_line_orientation(
8609 1, quad->line_orientation(1));
8610 new_quads[1]->set_line_orientation(
8611 3, quad->line_orientation(3));
8614 new_quads[0]->set_line_orientation(
8615 3, quad->line_orientation(3));
8616 new_quads[1]->set_line_orientation(
8617 2, quad->line_orientation(2));
8621 new_quads[0]->set_line_orientation(
8622 1, quad->line_orientation(1));
8623 new_quads[1]->set_line_orientation(
8624 0, quad->line_orientation(0));
8630 if (quad->refinement_case() ==
8653 if (aniso_quad_ref_case ==
8656 old_child[0] = quad->child(0)->line(1);
8657 old_child[1] = quad->child(2)->line(1);
8661 Assert(aniso_quad_ref_case ==
8665 old_child[0] = quad->child(0)->line(3);
8666 old_child[1] = quad->child(1)->line(3);
8669 if (old_child[0]->
index() + 1 != old_child[1]->
index())
8675 spacedim>::raw_line_iterator
8678 new_child[0] = new_child[1] =
8679 triangulation.
faces->lines
8680 .template next_free_pair_object<1>(
8684 new_child[0]->set_used_flag();
8685 new_child[1]->set_used_flag();
8687 const int old_index_0 = old_child[0]->index(),
8688 old_index_1 = old_child[1]->index(),
8689 new_index_0 = new_child[0]->index(),
8690 new_index_1 = new_child[1]->index();
8694 for (
unsigned int q = 0;
8695 q < triangulation.
faces->quads.n_objects();
8697 for (
unsigned int l = 0;
8698 l < GeometryInfo<dim>::lines_per_face;
8701 const int this_index =
8702 triangulation.
faces->quads
8703 .get_bounding_object_indices(q)[l];
8704 if (this_index == old_index_0)
8705 triangulation.
faces->quads
8706 .get_bounding_object_indices(q)[l] =
8708 else if (this_index == old_index_1)
8709 triangulation.
faces->quads
8710 .get_bounding_object_indices(q)[l] =
8715 for (
unsigned int i = 0; i < 2; ++i)
8717 Assert(!old_child[i]->has_children(),
8720 new_child[i]->set_bounding_object_indices(
8721 {old_child[i]->vertex_index(0),
8722 old_child[i]->vertex_index(1)});
8723 new_child[i]->set_boundary_id_internal(
8724 old_child[i]->boundary_id());
8725 new_child[i]->set_manifold_id(
8726 old_child[i]->manifold_id());
8727 new_child[i]->set_user_index(
8728 old_child[i]->user_index());
8729 if (old_child[i]->user_flag_set())
8730 new_child[i]->set_user_flag();
8732 new_child[i]->clear_user_flag();
8734 new_child[i]->clear_children();
8736 old_child[i]->clear_user_flag();
8737 old_child[i]->clear_user_index();
8738 old_child[i]->clear_used_flag();
8744 if (aniso_quad_ref_case ==
8747 new_line->set_children(
8748 0, quad->child(0)->line_index(1));
8749 Assert(new_line->child(1) ==
8750 quad->child(2)->line(1),
8787 quad_iterator switch_1 = quad->child(1),
8788 switch_2 = quad->child(2);
8789 const int switch_1_index = switch_1->index();
8790 const int switch_2_index = switch_2->index();
8791 for (
unsigned int l = 0;
8792 l < triangulation.
levels.size();
8794 for (
unsigned int h = 0;
8796 triangulation.
levels[l]->cells.n_objects();
8798 for (
const unsigned int q :
8801 const int face_index =
8803 ->cells.get_bounding_object_indices(
8805 if (face_index == switch_1_index)
8807 ->cells.get_bounding_object_indices(
8808 h)[q] = switch_2_index;
8809 else if (face_index == switch_2_index)
8811 ->cells.get_bounding_object_indices(
8812 h)[q] = switch_1_index;
8816 int switch_1_lines[4], switch_2_lines[4];
8817 for (
int i = 0; i < 4; ++i)
8819 switch_1_lines[i] = switch_1->line(i)->index();
8820 switch_2_lines[i] = switch_2->line(i)->index();
8823 switch_1_line_orientations[4] = {
8824 switch_1->line_orientation(0),
8825 switch_1->line_orientation(1),
8826 switch_1->line_orientation(2),
8827 switch_1->line_orientation(3)};
8829 switch_1->boundary_id();
8830 const unsigned int switch_1_user_index =
8831 switch_1->user_index();
8832 const bool switch_1_user_flag =
8833 switch_1->user_flag_set();
8835 switch_1_refinement_case =
8836 switch_1->refinement_case();
8837 const int switch_1_first_child_pair =
8838 (switch_1_refinement_case ?
8839 switch_1->child_index(0) :
8841 const int switch_1_second_child_pair =
8842 (switch_1_refinement_case ==
8844 switch_1->child_index(2) :
8847 switch_1->set_bounding_object_indices(
8851 switch_2_lines[3]});
8852 switch_1->set_line_orientation(
8853 0, switch_2->line_orientation(0));
8854 switch_1->set_line_orientation(
8855 1, switch_2->line_orientation(1));
8856 switch_1->set_line_orientation(
8857 2, switch_2->line_orientation(2));
8858 switch_1->set_line_orientation(
8859 3, switch_2->line_orientation(3));
8860 switch_1->set_boundary_id_internal(
8861 switch_2->boundary_id());
8862 switch_1->set_manifold_id(switch_2->manifold_id());
8863 switch_1->set_user_index(switch_2->user_index());
8864 if (switch_2->user_flag_set())
8865 switch_1->set_user_flag();
8867 switch_1->clear_user_flag();
8868 switch_1->clear_refinement_case();
8869 switch_1->set_refinement_case(
8870 switch_2->refinement_case());
8871 switch_1->clear_children();
8872 if (switch_2->refinement_case())
8873 switch_1->set_children(0,
8874 switch_2->child_index(0));
8875 if (switch_2->refinement_case() ==
8877 switch_1->set_children(2,
8878 switch_2->child_index(2));
8880 switch_2->set_bounding_object_indices(
8884 switch_1_lines[3]});
8885 switch_2->set_line_orientation(
8886 0, switch_1_line_orientations[0]);
8887 switch_2->set_line_orientation(
8888 1, switch_1_line_orientations[1]);
8889 switch_2->set_line_orientation(
8890 2, switch_1_line_orientations[2]);
8891 switch_2->set_line_orientation(
8892 3, switch_1_line_orientations[3]);
8893 switch_2->set_boundary_id_internal(
8894 switch_1_boundary_id);
8895 switch_2->set_manifold_id(switch_1->manifold_id());
8896 switch_2->set_user_index(switch_1_user_index);
8897 if (switch_1_user_flag)
8898 switch_2->set_user_flag();
8900 switch_2->clear_user_flag();
8901 switch_2->clear_refinement_case();
8902 switch_2->set_refinement_case(
8903 switch_1_refinement_case);
8904 switch_2->clear_children();
8905 switch_2->set_children(0,
8906 switch_1_first_child_pair);
8907 switch_2->set_children(2,
8908 switch_1_second_child_pair);
8910 new_quads[0]->set_refinement_case(
8912 new_quads[0]->set_children(0, quad->child_index(0));
8913 new_quads[1]->set_refinement_case(
8915 new_quads[1]->set_children(0, quad->child_index(2));
8919 new_quads[0]->set_refinement_case(
8921 new_quads[0]->set_children(0, quad->child_index(0));
8922 new_quads[1]->set_refinement_case(
8924 new_quads[1]->set_children(0, quad->child_index(2));
8925 new_line->set_children(
8926 0, quad->child(0)->line_index(3));
8927 Assert(new_line->child(1) ==
8928 quad->child(1)->line(3),
8931 quad->clear_children();
8935 quad->set_children(0, new_quads[0]->
index());
8937 quad->set_refinement_case(aniso_quad_ref_case);
8944 if (quad->user_flag_set())
8956 ++next_unused_vertex;
8958 next_unused_vertex < triangulation.
vertices.size(),
8960 "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."));
8968 quad->refinement_case();
8974 quad->child(0)->set_user_index(
8976 quad->child(1)->set_user_index(
8982 middle_line = quad->child(0)->line(1);
8984 middle_line = quad->child(0)->line(3);
8991 if (!middle_line->has_children())
8998 triangulation.
vertices[next_unused_vertex] =
8999 middle_line->center(
true);
9006 triangulation.
faces->lines
9007 .template next_free_pair_object<1>(
9012 middle_line->set_children(
9013 0, next_unused_line->index());
9017 raw_line_iterator children[2] = {
9018 next_unused_line, ++next_unused_line};
9026 children[0]->set_bounding_object_indices(
9027 {middle_line->vertex_index(0),
9028 next_unused_vertex});
9029 children[1]->set_bounding_object_indices(
9030 {next_unused_vertex,
9031 middle_line->vertex_index(1)});
9033 children[0]->set_used_flag();
9034 children[1]->set_used_flag();
9035 children[0]->clear_children();
9036 children[1]->clear_children();
9039 children[0]->clear_user_flag();
9040 children[1]->clear_user_flag();
9042 children[0]->set_boundary_id_internal(
9043 middle_line->boundary_id());
9044 children[1]->set_boundary_id_internal(
9045 middle_line->boundary_id());
9047 children[0]->set_manifold_id(
9048 middle_line->manifold_id());
9049 children[1]->set_manifold_id(
9050 middle_line->manifold_id());
9056 quad->clear_user_flag();
9085 triangulation.
vertices[next_unused_vertex] =
9086 quad->center(
true,
true);
9095 for (
unsigned int i = 0; i < 4; ++i)
9105 triangulation.
faces->lines
9106 .template next_free_pair_object<1>(triangulation);
9108 new_lines[i] = next_unused_line;
9132 quad->line(0)->child(0)->vertex_index(1),
9133 quad->line(1)->child(0)->vertex_index(1),
9134 quad->line(2)->child(0)->vertex_index(1),
9135 quad->line(3)->child(0)->vertex_index(1),
9136 next_unused_vertex};
9138 new_lines[0]->set_bounding_object_indices(
9140 new_lines[1]->set_bounding_object_indices(
9142 new_lines[2]->set_bounding_object_indices(
9144 new_lines[3]->set_bounding_object_indices(
9147 for (
const auto &new_line : new_lines)
9149 new_line->set_used_flag();
9150 new_line->clear_user_flag();
9152 new_line->clear_children();
9153 new_line->set_boundary_id_internal(quad->boundary_id());
9154 new_line->set_manifold_id(quad->manifold_id());
9167 const int line_indices[12] = {
9169 ->child(child_line_index(0, quad->line_orientation(0)))
9172 ->child(child_line_index(1, quad->line_orientation(0)))
9175 ->child(child_line_index(0, quad->line_orientation(1)))
9178 ->child(child_line_index(1, quad->line_orientation(1)))
9181 ->child(child_line_index(0, quad->line_orientation(2)))
9184 ->child(child_line_index(1, quad->line_orientation(2)))
9187 ->child(child_line_index(0, quad->line_orientation(3)))
9190 ->child(child_line_index(1, quad->line_orientation(3)))
9192 new_lines[0]->index(),
9193 new_lines[1]->index(),
9194 new_lines[2]->index(),
9195 new_lines[3]->index()};
9204 triangulation.
faces->quads
9205 .template next_free_pair_object<2>(triangulation);
9207 new_quads[0] = next_unused_quad;
9211 new_quads[1] = next_unused_quad;
9215 triangulation.
faces->quads
9216 .template next_free_pair_object<2>(triangulation);
9217 new_quads[2] = next_unused_quad;
9221 new_quads[3] = next_unused_quad;
9225 quad->set_children(0, new_quads[0]->
index());
9226 quad->set_children(2, new_quads[2]->
index());
9229 new_quads[0]->set_bounding_object_indices(
9234 new_quads[1]->set_bounding_object_indices(
9239 new_quads[2]->set_bounding_object_indices(
9244 new_quads[3]->set_bounding_object_indices(
9249 for (
const auto &new_quad : new_quads)
9251 new_quad->set_used_flag();
9252 new_quad->clear_user_flag();
9254 new_quad->clear_children();
9255 new_quad->set_boundary_id_internal(quad->boundary_id());
9256 new_quad->set_manifold_id(quad->manifold_id());
9260 for (
unsigned int j = 0;
9261 j < GeometryInfo<dim>::lines_per_face;
9263 new_quad->set_line_orientation(
9270 new_quads[0]->set_line_orientation(
9271 0, quad->line_orientation(0));
9272 new_quads[0]->set_line_orientation(
9273 2, quad->line_orientation(2));
9274 new_quads[1]->set_line_orientation(
9275 1, quad->line_orientation(1));
9276 new_quads[1]->set_line_orientation(
9277 2, quad->line_orientation(2));
9278 new_quads[2]->set_line_orientation(
9279 0, quad->line_orientation(0));
9280 new_quads[2]->set_line_orientation(
9281 3, quad->line_orientation(3));
9282 new_quads[3]->set_line_orientation(
9283 1, quad->line_orientation(1));
9284 new_quads[3]->set_line_orientation(
9285 3, quad->line_orientation(3));
9289 quad->clear_user_flag();
9300 cells_with_distorted_children;
9314 for (; hex != endh; ++hex)
9315 if (hex->refine_flag_set())
9323 hex->clear_refine_flag();
9324 hex->set_refinement_case(ref_case);
9335 unsigned int n_new_lines = 0;
9336 unsigned int n_new_quads = 0;
9337 unsigned int n_new_hexes = 0;
9368 new_lines(n_new_lines);
9369 for (
unsigned int i = 0; i < n_new_lines; ++i)
9372 triangulation.
faces->lines
9373 .template next_free_single_object<1>(triangulation);
9376 new_lines[i]->set_used_flag();
9377 new_lines[i]->clear_user_flag();
9378 new_lines[i]->clear_user_data();
9379 new_lines[i]->clear_children();
9381 new_lines[i]->set_boundary_id_internal(
9385 new_lines[i]->set_manifold_id(hex->manifold_id());
9392 new_quads(n_new_quads);
9393 for (
unsigned int i = 0; i < n_new_quads; ++i)
9396 triangulation.
faces->quads
9397 .template next_free_single_object<2>(triangulation);
9400 new_quads[i]->set_used_flag();
9401 new_quads[i]->clear_user_flag();
9402 new_quads[i]->clear_user_data();
9403 new_quads[i]->clear_children();
9405 new_quads[i]->set_boundary_id_internal(
9409 new_quads[i]->set_manifold_id(hex->manifold_id());
9412 for (
unsigned int j = 0;
9413 j < GeometryInfo<dim>::lines_per_face;
9415 new_quads[i]->set_line_orientation(
9425 new_hexes(n_new_hexes);
9426 for (
unsigned int i = 0; i < n_new_hexes; ++i)
9430 triangulation.
levels[
level + 1]->cells.next_free_hex(
9431 triangulation,
level + 1);
9435 new_hexes[i] = next_unused_hex;
9438 new_hexes[i]->set_used_flag();
9439 new_hexes[i]->clear_user_flag();
9440 new_hexes[i]->clear_user_data();
9441 new_hexes[i]->clear_children();
9444 new_hexes[i]->set_material_id(hex->material_id());
9445 new_hexes[i]->set_manifold_id(hex->manifold_id());
9446 new_hexes[i]->set_subdomain_id(subdomainid);
9453 new_hexes[i]->set_parent(hex->index());
9466 for (
const unsigned int f :
9468 new_hexes[i]->set_combined_face_orientation(
9472 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
9473 hex->set_children(2 * i, new_hexes[2 * i]->index());
9480 const bool f_or[6] = {hex->face_orientation(0),
9481 hex->face_orientation(1),
9482 hex->face_orientation(2),
9483 hex->face_orientation(3),
9484 hex->face_orientation(4),
9485 hex->face_orientation(5)};
9488 const bool f_fl[6] = {hex->face_flip(0),
9496 const bool f_ro[6] = {hex->face_rotation(0),
9497 hex->face_rotation(1),
9498 hex->face_rotation(2),
9499 hex->face_rotation(3),
9500 hex->face_rotation(4),
9501 hex->face_rotation(5)};
9505 hex->combined_face_orientation(0),
9506 hex->combined_face_orientation(1),
9507 hex->combined_face_orientation(2),
9508 hex->combined_face_orientation(3),
9509 hex->combined_face_orientation(4),
9510 hex->combined_face_orientation(5)};
9520 const unsigned int child_at_origin[2][2][2] = {
9607 raw_line_iterator lines[4] = {
9608 hex->face(2)->child(0)->line(
9609 (hex->face(2)->refinement_case() ==
9613 hex->face(3)->child(0)->line(
9614 (hex->face(3)->refinement_case() ==
9618 hex->face(4)->child(0)->line(
9619 (hex->face(4)->refinement_case() ==
9623 hex->face(5)->child(0)->line(
9624 (hex->face(5)->refinement_case() ==
9630 int line_indices[4];
9631 for (
unsigned int i = 0; i < 4; ++i)
9632 line_indices[i] = lines[i]->
index();
9645 const unsigned int middle_vertices[2] = {
9646 hex->line(2)->child(0)->vertex_index(1),
9647 hex->line(7)->child(0)->vertex_index(1)};
9649 for (
unsigned int i = 0; i < 4; ++i)
9650 if (lines[i]->vertex_index(i % 2) ==
9651 middle_vertices[i % 2])
9652 line_orientation[i] =
9657 Assert(lines[i]->vertex_index((i + 1) % 2) ==
9658 middle_vertices[i % 2],
9660 line_orientation[i] =
9666 new_quads[0]->set_bounding_object_indices(
9672 new_quads[0]->set_line_orientation(
9673 0, line_orientation[0]);
9674 new_quads[0]->set_line_orientation(
9675 1, line_orientation[1]);
9676 new_quads[0]->set_line_orientation(
9677 2, line_orientation[2]);
9678 new_quads[0]->set_line_orientation(
9679 3, line_orientation[3]);
9712 const int quad_indices[11] = {
9713 new_quads[0]->index(),
9715 hex->face(0)->index(),
9717 hex->face(1)->index(),
9719 hex->face(2)->child_index(
9720 child_at_origin[hex->face(2)->refinement_case() -
9721 1][f_fl[2]][f_ro[2]]),
9722 hex->face(2)->child_index(
9724 child_at_origin[hex->face(2)->refinement_case() -
9725 1][f_fl[2]][f_ro[2]]),
9727 hex->face(3)->child_index(
9728 child_at_origin[hex->face(3)->refinement_case() -
9729 1][f_fl[3]][f_ro[3]]),
9730 hex->face(3)->child_index(
9732 child_at_origin[hex->face(3)->refinement_case() -
9733 1][f_fl[3]][f_ro[3]]),
9735 hex->face(4)->child_index(
9736 child_at_origin[hex->face(4)->refinement_case() -
9737 1][f_fl[4]][f_ro[4]]),
9738 hex->face(4)->child_index(
9740 child_at_origin[hex->face(4)->refinement_case() -
9741 1][f_fl[4]][f_ro[4]]),
9743 hex->face(5)->child_index(
9744 child_at_origin[hex->face(5)->refinement_case() -
9745 1][f_fl[5]][f_ro[5]]),
9746 hex->face(5)->child_index(
9748 child_at_origin[hex->face(5)->refinement_case() -
9749 1][f_fl[5]][f_ro[5]])
9753 new_hexes[0]->set_bounding_object_indices(
9760 new_hexes[1]->set_bounding_object_indices(
9835 raw_line_iterator lines[4] = {
9836 hex->face(0)->child(0)->line(
9837 (hex->face(0)->refinement_case() ==
9841 hex->face(1)->child(0)->line(
9842 (hex->face(1)->refinement_case() ==
9846 hex->face(4)->child(0)->line(
9847 (hex->face(4)->refinement_case() ==
9851 hex->face(5)->child(0)->line(
9852 (hex->face(5)->refinement_case() ==
9858 int line_indices[4];
9859 for (
unsigned int i = 0; i < 4; ++i)
9860 line_indices[i] = lines[i]->
index();
9873 const unsigned int middle_vertices[2] = {
9874 hex->line(0)->child(0)->vertex_index(1),
9875 hex->line(5)->child(0)->vertex_index(1)};
9877 for (
unsigned int i = 0; i < 4; ++i)
9878 if (lines[i]->vertex_index(i % 2) ==
9879 middle_vertices[i % 2])
9880 line_orientation[i] =
9885 Assert(lines[i]->vertex_index((i + 1) % 2) ==
9886 middle_vertices[i % 2],
9888 line_orientation[i] =
9894 new_quads[0]->set_bounding_object_indices(
9900 new_quads[0]->set_line_orientation(
9901 0, line_orientation[2]);
9902 new_quads[0]->set_line_orientation(
9903 1, line_orientation[3]);
9904 new_quads[0]->set_line_orientation(
9905 2, line_orientation[0]);
9906 new_quads[0]->set_line_orientation(
9907 3, line_orientation[1]);
9940 const int quad_indices[11] = {
9941 new_quads[0]->index(),
9943 hex->face(0)->child_index(
9944 child_at_origin[hex->face(0)->refinement_case() -
9945 1][f_fl[0]][f_ro[0]]),
9946 hex->face(0)->child_index(
9948 child_at_origin[hex->face(0)->refinement_case() -
9949 1][f_fl[0]][f_ro[0]]),
9951 hex->face(1)->child_index(
9952 child_at_origin[hex->face(1)->refinement_case() -
9953 1][f_fl[1]][f_ro[1]]),
9954 hex->face(1)->child_index(
9956 child_at_origin[hex->face(1)->refinement_case() -
9957 1][f_fl[1]][f_ro[1]]),
9959 hex->face(2)->index(),
9961 hex->face(3)->index(),
9963 hex->face(4)->child_index(
9964 child_at_origin[hex->face(4)->refinement_case() -
9965 1][f_fl[4]][f_ro[4]]),
9966 hex->face(4)->child_index(
9968 child_at_origin[hex->face(4)->refinement_case() -
9969 1][f_fl[4]][f_ro[4]]),
9971 hex->face(5)->child_index(
9972 child_at_origin[hex->face(5)->refinement_case() -
9973 1][f_fl[5]][f_ro[5]]),
9974 hex->face(5)->child_index(
9976 child_at_origin[hex->face(5)->refinement_case() -
9977 1][f_fl[5]][f_ro[5]])
9981 new_hexes[0]->set_bounding_object_indices(
9988 new_hexes[1]->set_bounding_object_indices(
10065 raw_line_iterator lines[4] = {
10066 hex->face(0)->child(0)->line(
10067 (hex->face(0)->refinement_case() ==
10071 hex->face(1)->child(0)->line(
10072 (hex->face(1)->refinement_case() ==
10076 hex->face(2)->child(0)->line(
10077 (hex->face(2)->refinement_case() ==
10081 hex->face(3)->child(0)->line(
10082 (hex->face(3)->refinement_case() ==
10088 int line_indices[4];
10089 for (
unsigned int i = 0; i < 4; ++i)
10090 line_indices[i] = lines[i]->
index();
10103 const unsigned int middle_vertices[2] = {
10104 middle_vertex_index<dim, spacedim>(hex->line(8)),
10105 middle_vertex_index<dim, spacedim>(hex->line(11))};
10107 for (
unsigned int i = 0; i < 4; ++i)
10108 if (lines[i]->vertex_index(i % 2) ==
10109 middle_vertices[i % 2])
10110 line_orientation[i] =
10115 Assert(lines[i]->vertex_index((i + 1) % 2) ==
10116 middle_vertices[i % 2],
10118 line_orientation[i] =
10124 new_quads[0]->set_bounding_object_indices(
10130 new_quads[0]->set_line_orientation(
10131 0, line_orientation[0]);
10132 new_quads[0]->set_line_orientation(
10133 1, line_orientation[1]);
10134 new_quads[0]->set_line_orientation(
10135 2, line_orientation[2]);
10136 new_quads[0]->set_line_orientation(
10137 3, line_orientation[3]);
10171 const int quad_indices[11] = {
10172 new_quads[0]->index(),
10174 hex->face(0)->child_index(
10175 child_at_origin[hex->face(0)->refinement_case() -
10176 1][f_fl[0]][f_ro[0]]),
10177 hex->face(0)->child_index(
10179 child_at_origin[hex->face(0)->refinement_case() -
10180 1][f_fl[0]][f_ro[0]]),
10182 hex->face(1)->child_index(
10183 child_at_origin[hex->face(1)->refinement_case() -
10184 1][f_fl[1]][f_ro[1]]),
10185 hex->face(1)->child_index(
10187 child_at_origin[hex->face(1)->refinement_case() -
10188 1][f_fl[1]][f_ro[1]]),
10190 hex->face(2)->child_index(
10191 child_at_origin[hex->face(2)->refinement_case() -
10192 1][f_fl[2]][f_ro[2]]),
10193 hex->face(2)->child_index(
10195 child_at_origin[hex->face(2)->refinement_case() -
10196 1][f_fl[2]][f_ro[2]]),
10198 hex->face(3)->child_index(
10199 child_at_origin[hex->face(3)->refinement_case() -
10200 1][f_fl[3]][f_ro[3]]),
10201 hex->face(3)->child_index(
10203 child_at_origin[hex->face(3)->refinement_case() -
10204 1][f_fl[3]][f_ro[3]]),
10206 hex->face(4)->index(),
10208 hex->face(5)->index()
10211 new_hexes[0]->set_bounding_object_indices(
10218 new_hexes[1]->set_bounding_object_indices(
10224 quad_indices[10]});
10250 new_lines[0]->set_bounding_object_indices(
10251 {middle_vertex_index<dim, spacedim>(hex->face(4)),
10252 middle_vertex_index<dim, spacedim>(hex->face(5))});
10320 spacedim>::raw_line_iterator lines[13] = {
10321 hex->face(0)->child(0)->line(
10322 (hex->face(0)->refinement_case() ==
10326 hex->face(1)->child(0)->line(
10327 (hex->face(1)->refinement_case() ==
10331 hex->face(2)->child(0)->line(
10332 (hex->face(2)->refinement_case() ==
10336 hex->face(3)->child(0)->line(
10337 (hex->face(3)->refinement_case() ==
10345 0, f_or[4], f_fl[4], f_ro[4]))
10348 1, f_or[4], f_fl[4], f_ro[4])),
10352 3, f_or[4], f_fl[4], f_ro[4]))
10355 0, f_or[4], f_fl[4], f_ro[4])),
10359 0, f_or[4], f_fl[4], f_ro[4]))
10362 3, f_or[4], f_fl[4], f_ro[4])),
10366 3, f_or[4], f_fl[4], f_ro[4]))
10369 2, f_or[4], f_fl[4], f_ro[4])),
10374 0, f_or[5], f_fl[5], f_ro[5]))
10377 1, f_or[5], f_fl[5], f_ro[5])),
10381 3, f_or[5], f_fl[5], f_ro[5]))
10384 0, f_or[5], f_fl[5], f_ro[5])),
10388 0, f_or[5], f_fl[5], f_ro[5]))
10391 3, f_or[5], f_fl[5], f_ro[5])),
10395 3, f_or[5], f_fl[5], f_ro[5]))
10398 2, f_or[5], f_fl[5], f_ro[5])),
10403 int line_indices[13];
10404 for (
unsigned int i = 0; i < 13; ++i)
10405 line_indices[i] = lines[i]->
index();
10416 const unsigned int middle_vertices[4] = {
10417 hex->line(0)->child(0)->vertex_index(1),
10418 hex->line(1)->child(0)->vertex_index(1),
10419 hex->line(2)->child(0)->vertex_index(1),
10420 hex->line(3)->child(0)->vertex_index(1),
10426 for (
unsigned int i = 0; i < 4; ++i)
10427 if (lines[i]->vertex_index(0) == middle_vertices[i])
10428 line_orientation[i] =
10433 Assert(lines[i]->vertex_index(1) ==
10434 middle_vertices[i],
10436 line_orientation[i] =
10446 for (
unsigned int i = 4; i < 12; ++i)
10447 if (lines[i]->vertex_index((i + 1) % 2) ==
10448 middle_vertex_index<dim, spacedim>(
10449 hex->face(3 + i / 4)))
10450 line_orientation[i] =
10455 Assert(lines[i]->vertex_index(i % 2) ==
10456 (middle_vertex_index<dim, spacedim>(
10457 hex->face(3 + i / 4))),
10459 line_orientation[i] =
10465 line_orientation[12] =
10491 new_quads[0]->set_bounding_object_indices(
10496 new_quads[1]->set_bounding_object_indices(
10501 new_quads[2]->set_bounding_object_indices(
10505 line_indices[12]});
10506 new_quads[3]->set_bounding_object_indices(
10512 new_quads[0]->set_line_orientation(
10513 0, line_orientation[2]);
10514 new_quads[0]->set_line_orientation(
10515 2, line_orientation[4]);
10516 new_quads[0]->set_line_orientation(
10517 3, line_orientation[8]);
10519 new_quads[1]->set_line_orientation(
10520 1, line_orientation[3]);
10521 new_quads[1]->set_line_orientation(
10522 2, line_orientation[5]);
10523 new_quads[1]->set_line_orientation(
10524 3, line_orientation[9]);
10526 new_quads[2]->set_line_orientation(
10527 0, line_orientation[6]);
10528 new_quads[2]->set_line_orientation(
10529 1, line_orientation[10]);
10530 new_quads[2]->set_line_orientation(
10531 2, line_orientation[0]);
10533 new_quads[3]->set_line_orientation(
10534 0, line_orientation[7]);
10535 new_quads[3]->set_line_orientation(
10536 1, line_orientation[11]);
10537 new_quads[3]->set_line_orientation(
10538 3, line_orientation[1]);
10571 const int quad_indices[20] = {
10572 new_quads[0]->index(),
10573 new_quads[1]->index(),
10574 new_quads[2]->index(),
10575 new_quads[3]->index(),
10577 hex->face(0)->child_index(
10578 child_at_origin[hex->face(0)->refinement_case() -
10579 1][f_fl[0]][f_ro[0]]),
10580 hex->face(0)->child_index(
10582 child_at_origin[hex->face(0)->refinement_case() -
10583 1][f_fl[0]][f_ro[0]]),
10585 hex->face(1)->child_index(
10586 child_at_origin[hex->face(1)->refinement_case() -
10587 1][f_fl[1]][f_ro[1]]),
10588 hex->face(1)->child_index(
10590 child_at_origin[hex->face(1)->refinement_case() -
10591 1][f_fl[1]][f_ro[1]]),
10593 hex->face(2)->child_index(
10594 child_at_origin[hex->face(2)->refinement_case() -
10595 1][f_fl[2]][f_ro[2]]),
10596 hex->face(2)->child_index(
10598 child_at_origin[hex->face(2)->refinement_case() -
10599 1][f_fl[2]][f_ro[2]]),
10601 hex->face(3)->child_index(
10602 child_at_origin[hex->face(3)->refinement_case() -
10603 1][f_fl[3]][f_ro[3]]),
10604 hex->face(3)->child_index(
10606 child_at_origin[hex->face(3)->refinement_case() -
10607 1][f_fl[3]][f_ro[3]]),
10609 hex->face(4)->isotropic_child_index(
10611 0, f_or[4], f_fl[4], f_ro[4])),
10612 hex->face(4)->isotropic_child_index(
10614 1, f_or[4], f_fl[4], f_ro[4])),
10615 hex->face(4)->isotropic_child_index(
10617 2, f_or[4], f_fl[4], f_ro[4])),
10618 hex->face(4)->isotropic_child_index(
10620 3, f_or[4], f_fl[4], f_ro[4])),
10622 hex->face(5)->isotropic_child_index(
10624 0, f_or[5], f_fl[5], f_ro[5])),
10625 hex->face(5)->isotropic_child_index(
10627 1, f_or[5], f_fl[5], f_ro[5])),
10628 hex->face(5)->isotropic_child_index(
10630 2, f_or[5], f_fl[5], f_ro[5])),
10631 hex->face(5)->isotropic_child_index(
10633 3, f_or[5], f_fl[5], f_ro[5]))};
10635 new_hexes[0]->set_bounding_object_indices(
10641 quad_indices[16]});
10642 new_hexes[1]->set_bounding_object_indices(
10648 quad_indices[17]});
10649 new_hexes[2]->set_bounding_object_indices(
10655 quad_indices[18]});
10656 new_hexes[3]->set_bounding_object_indices(
10662 quad_indices[19]});
10688 new_lines[0]->set_bounding_object_indices(
10689 {middle_vertex_index<dim, spacedim>(hex->face(2)),
10690 middle_vertex_index<dim, spacedim>(hex->face(3))});
10758 spacedim>::raw_line_iterator lines[13] = {
10759 hex->face(0)->child(0)->line(
10760 (hex->face(0)->refinement_case() ==
10764 hex->face(1)->child(0)->line(
10765 (hex->face(1)->refinement_case() ==
10769 hex->face(4)->child(0)->line(
10770 (hex->face(4)->refinement_case() ==
10774 hex->face(5)->child(0)->line(
10775 (hex->face(5)->refinement_case() ==
10783 0, f_or[2], f_fl[2], f_ro[2]))
10786 3, f_or[2], f_fl[2], f_ro[2])),
10790 3, f_or[2], f_fl[2], f_ro[2]))
10793 2, f_or[2], f_fl[2], f_ro[2])),
10797 0, f_or[2], f_fl[2], f_ro[2]))
10800 1, f_or[2], f_fl[2], f_ro[2])),
10804 3, f_or[2], f_fl[2], f_ro[2]))
10807 0, f_or[2], f_fl[2], f_ro[2])),
10812 0, f_or[3], f_fl[3], f_ro[3]))
10815 3, f_or[3], f_fl[3], f_ro[3])),
10819 3, f_or[3], f_fl[3], f_ro[3]))
10822 2, f_or[3], f_fl[3], f_ro[3])),
10826 0, f_or[3], f_fl[3], f_ro[3]))
10829 1, f_or[3], f_fl[3], f_ro[3])),
10833 3, f_or[3], f_fl[3], f_ro[3]))
10836 0, f_or[3], f_fl[3], f_ro[3])),
10841 int line_indices[13];
10842 for (
unsigned int i = 0; i < 13; ++i)
10843 line_indices[i] = lines[i]->
index();
10854 const unsigned int middle_vertices[4] = {
10855 hex->line(8)->child(0)->vertex_index(1),
10856 hex->line(9)->child(0)->vertex_index(1),
10857 hex->line(2)->child(0)->vertex_index(1),
10858 hex->line(6)->child(0)->vertex_index(1),
10863 for (
unsigned int i = 0; i < 4; ++i)
10864 if (lines[i]->vertex_index(0) == middle_vertices[i])
10865 line_orientation[i] =
10870 Assert(lines[i]->vertex_index(1) ==
10871 middle_vertices[i],
10873 line_orientation[i] =
10883 for (
unsigned int i = 4; i < 12; ++i)
10884 if (lines[i]->vertex_index((i + 1) % 2) ==
10885 middle_vertex_index<dim, spacedim>(
10886 hex->face(1 + i / 4)))
10887 line_orientation[i] =
10893 Assert(lines[i]->vertex_index(i % 2) ==
10894 (middle_vertex_index<dim, spacedim>(
10895 hex->face(1 + i / 4))),
10897 line_orientation[i] =
10903 line_orientation[12] =
10930 new_quads[0]->set_bounding_object_indices(
10934 line_indices[10]});
10935 new_quads[1]->set_bounding_object_indices(
10939 line_indices[11]});
10940 new_quads[2]->set_bounding_object_indices(
10944 line_indices[12]});
10945 new_quads[3]->set_bounding_object_indices(
10951 new_quads[0]->set_line_orientation(
10952 0, line_orientation[0]);
10953 new_quads[0]->set_line_orientation(
10954 2, line_orientation[6]);
10955 new_quads[0]->set_line_orientation(
10956 3, line_orientation[10]);
10958 new_quads[1]->set_line_orientation(
10959 1, line_orientation[1]);
10960 new_quads[1]->set_line_orientation(
10961 2, line_orientation[7]);
10962 new_quads[1]->set_line_orientation(
10963 3, line_orientation[11]);
10965 new_quads[2]->set_line_orientation(
10966 0, line_orientation[4]);
10967 new_quads[2]->set_line_orientation(
10968 1, line_orientation[8]);
10969 new_quads[2]->set_line_orientation(
10970 2, line_orientation[2]);
10972 new_quads[3]->set_line_orientation(
10973 0, line_orientation[5]);
10974 new_quads[3]->set_line_orientation(
10975 1, line_orientation[9]);
10976 new_quads[3]->set_line_orientation(
10977 3, line_orientation[3]);
11011 const int quad_indices[20] = {
11012 new_quads[0]->index(),
11013 new_quads[1]->index(),
11014 new_quads[2]->index(),
11015 new_quads[3]->index(),
11017 hex->face(0)->child_index(
11018 child_at_origin[hex->face(0)->refinement_case() -
11019 1][f_fl[0]][f_ro[0]]),
11020 hex->face(0)->child_index(
11022 child_at_origin[hex->face(0)->refinement_case() -
11023 1][f_fl[0]][f_ro[0]]),
11025 hex->face(1)->child_index(
11026 child_at_origin[hex->face(1)->refinement_case() -
11027 1][f_fl[1]][f_ro[1]]),
11028 hex->face(1)->child_index(
11030 child_at_origin[hex->face(1)->refinement_case() -
11031 1][f_fl[1]][f_ro[1]]),
11033 hex->face(2)->isotropic_child_index(
11035 0, f_or[2], f_fl[2], f_ro[2])),
11036 hex->face(2)->isotropic_child_index(
11038 1, f_or[2], f_fl[2], f_ro[2])),
11039 hex->face(2)->isotropic_child_index(
11041 2, f_or[2], f_fl[2], f_ro[2])),
11042 hex->face(2)->isotropic_child_index(
11044 3, f_or[2], f_fl[2], f_ro[2])),
11046 hex->face(3)->isotropic_child_index(
11048 0, f_or[3], f_fl[3], f_ro[3])),
11049 hex->face(3)->isotropic_child_index(
11051 1, f_or[3], f_fl[3], f_ro[3])),
11052 hex->face(3)->isotropic_child_index(
11054 2, f_or[3], f_fl[3], f_ro[3])),
11055 hex->face(3)->isotropic_child_index(
11057 3, f_or[3], f_fl[3], f_ro[3])),
11059 hex->face(4)->child_index(
11060 child_at_origin[hex->face(4)->refinement_case() -
11061 1][f_fl[4]][f_ro[4]]),
11062 hex->face(4)->child_index(
11064 child_at_origin[hex->face(4)->refinement_case() -
11065 1][f_fl[4]][f_ro[4]]),
11067 hex->face(5)->child_index(
11068 child_at_origin[hex->face(5)->refinement_case() -
11069 1][f_fl[5]][f_ro[5]]),
11070 hex->face(5)->child_index(
11072 child_at_origin[hex->face(5)->refinement_case() -
11073 1][f_fl[5]][f_ro[5]])};
11084 new_hexes[0]->set_bounding_object_indices(
11091 new_hexes[1]->set_bounding_object_indices(
11097 quad_indices[18]});
11098 new_hexes[2]->set_bounding_object_indices(
11105 new_hexes[3]->set_bounding_object_indices(
11111 quad_indices[19]});
11138 new_lines[0]->set_bounding_object_indices(
11140 {middle_vertex_index<dim, spacedim>(hex->face(0)),
11141 middle_vertex_index<dim, spacedim>(hex->face(1))});
11210 spacedim>::raw_line_iterator lines[13] = {
11211 hex->face(2)->child(0)->line(
11212 (hex->face(2)->refinement_case() ==
11216 hex->face(3)->child(0)->line(
11217 (hex->face(3)->refinement_case() ==
11221 hex->face(4)->child(0)->line(
11222 (hex->face(4)->refinement_case() ==
11226 hex->face(5)->child(0)->line(
11227 (hex->face(5)->refinement_case() ==
11235 0, f_or[0], f_fl[0], f_ro[0]))
11238 1, f_or[0], f_fl[0], f_ro[0])),
11242 3, f_or[0], f_fl[0], f_ro[0]))
11245 0, f_or[0], f_fl[0], f_ro[0])),
11249 0, f_or[0], f_fl[0], f_ro[0]))
11252 3, f_or[0], f_fl[0], f_ro[0])),
11256 3, f_or[0], f_fl[0], f_ro[0]))
11259 2, f_or[0], f_fl[0], f_ro[0])),
11264 0, f_or[1], f_fl[1], f_ro[1]))
11267 1, f_or[1], f_fl[1], f_ro[1])),
11271 3, f_or[1], f_fl[1], f_ro[1]))
11274 0, f_or[1], f_fl[1], f_ro[1])),
11278 0, f_or[1], f_fl[1], f_ro[1]))
11281 3, f_or[1], f_fl[1], f_ro[1])),
11285 3, f_or[1], f_fl[1], f_ro[1]))
11288 2, f_or[1], f_fl[1], f_ro[1])),
11293 int line_indices[13];
11294 for (
unsigned int i = 0; i < 13; ++i)
11295 line_indices[i] = lines[i]->
index();
11306 const unsigned int middle_vertices[4] = {
11307 hex->line(8)->child(0)->vertex_index(1),
11308 hex->line(10)->child(0)->vertex_index(1),
11309 hex->line(0)->child(0)->vertex_index(1),
11310 hex->line(4)->child(0)->vertex_index(1),
11315 for (
unsigned int i = 0; i < 4; ++i)
11316 if (lines[i]->vertex_index(0) == middle_vertices[i])
11317 line_orientation[i] =
11322 Assert(lines[i]->vertex_index(1) ==
11323 middle_vertices[i],
11325 line_orientation[i] =
11335 for (
unsigned int i = 4; i < 12; ++i)
11336 if (lines[i]->vertex_index((i + 1) % 2) ==
11337 middle_vertex_index<dim, spacedim>(
11338 hex->face(i / 4 - 1)))
11339 line_orientation[i] =
11344 Assert(lines[i]->vertex_index(i % 2) ==
11345 (middle_vertex_index<dim, spacedim>(
11346 hex->face(i / 4 - 1))),
11348 line_orientation[i] =
11353 line_orientation[12] =
11374 new_quads[0]->set_bounding_object_indices(
11378 line_indices[12]});
11379 new_quads[1]->set_bounding_object_indices(
11384 new_quads[2]->set_bounding_object_indices(
11389 new_quads[3]->set_bounding_object_indices(
11395 new_quads[0]->set_line_orientation(
11396 0, line_orientation[6]);
11397 new_quads[0]->set_line_orientation(
11398 1, line_orientation[10]);
11399 new_quads[0]->set_line_orientation(
11400 2, line_orientation[0]);
11402 new_quads[1]->set_line_orientation(
11403 0, line_orientation[7]);
11404 new_quads[1]->set_line_orientation(
11405 1, line_orientation[11]);
11406 new_quads[1]->set_line_orientation(
11407 3, line_orientation[1]);
11409 new_quads[2]->set_line_orientation(
11410 0, line_orientation[2]);
11411 new_quads[2]->set_line_orientation(
11412 2, line_orientation[4]);
11413 new_quads[2]->set_line_orientation(
11414 3, line_orientation[8]);
11416 new_quads[3]->set_line_orientation(
11417 1, line_orientation[3]);
11418 new_quads[3]->set_line_orientation(
11419 2, line_orientation[5]);
11420 new_quads[3]->set_line_orientation(
11421 3, line_orientation[9]);
11455 const int quad_indices[20] = {
11456 new_quads[0]->index(),
11457 new_quads[1]->index(),
11458 new_quads[2]->index(),
11459 new_quads[3]->index(),
11461 hex->face(0)->isotropic_child_index(
11463 0, f_or[0], f_fl[0], f_ro[0])),
11464 hex->face(0)->isotropic_child_index(
11466 1, f_or[0], f_fl[0], f_ro[0])),
11467 hex->face(0)->isotropic_child_index(
11469 2, f_or[0], f_fl[0], f_ro[0])),
11470 hex->face(0)->isotropic_child_index(
11472 3, f_or[0], f_fl[0], f_ro[0])),
11474 hex->face(1)->isotropic_child_index(
11476 0, f_or[1], f_fl[1], f_ro[1])),
11477 hex->face(1)->isotropic_child_index(
11479 1, f_or[1], f_fl[1], f_ro[1])),
11480 hex->face(1)->isotropic_child_index(
11482 2, f_or[1], f_fl[1], f_ro[1])),
11483 hex->face(1)->isotropic_child_index(
11485 3, f_or[1], f_fl[1], f_ro[1])),
11487 hex->face(2)->child_index(
11488 child_at_origin[hex->face(2)->refinement_case() -
11489 1][f_fl[2]][f_ro[2]]),
11490 hex->face(2)->child_index(
11492 child_at_origin[hex->face(2)->refinement_case() -
11493 1][f_fl[2]][f_ro[2]]),
11495 hex->face(3)->child_index(
11496 child_at_origin[hex->face(3)->refinement_case() -
11497 1][f_fl[3]][f_ro[3]]),
11498 hex->face(3)->child_index(
11500 child_at_origin[hex->face(3)->refinement_case() -
11501 1][f_fl[3]][f_ro[3]]),
11503 hex->face(4)->child_index(
11504 child_at_origin[hex->face(4)->refinement_case() -
11505 1][f_fl[4]][f_ro[4]]),
11506 hex->face(4)->child_index(
11508 child_at_origin[hex->face(4)->refinement_case() -
11509 1][f_fl[4]][f_ro[4]]),
11511 hex->face(5)->child_index(
11512 child_at_origin[hex->face(5)->refinement_case() -
11513 1][f_fl[5]][f_ro[5]]),
11514 hex->face(5)->child_index(
11516 child_at_origin[hex->face(5)->refinement_case() -
11517 1][f_fl[5]][f_ro[5]])};
11519 new_hexes[0]->set_bounding_object_indices(
11526 new_hexes[1]->set_bounding_object_indices(
11533 new_hexes[2]->set_bounding_object_indices(
11539 quad_indices[18]});
11540 new_hexes[3]->set_bounding_object_indices(
11546 quad_indices[19]});
11578 ++next_unused_vertex;
11580 next_unused_vertex < triangulation.
vertices.size(),
11582 "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."));
11592 triangulation.
vertices[next_unused_vertex] =
11593 hex->center(
true,
true);
11615 middle_vertex_index<dim, spacedim>(hex->face(0)),
11616 middle_vertex_index<dim, spacedim>(hex->face(1)),
11617 middle_vertex_index<dim, spacedim>(hex->face(2)),
11618 middle_vertex_index<dim, spacedim>(hex->face(3)),
11619 middle_vertex_index<dim, spacedim>(hex->face(4)),
11620 middle_vertex_index<dim, spacedim>(hex->face(5)),
11621 next_unused_vertex};
11623 new_lines[0]->set_bounding_object_indices(
11625 new_lines[1]->set_bounding_object_indices(
11627 new_lines[2]->set_bounding_object_indices(
11629 new_lines[3]->set_bounding_object_indices(
11631 new_lines[4]->set_bounding_object_indices(
11633 new_lines[5]->set_bounding_object_indices(
11703 spacedim>::raw_line_iterator lines[30] = {
11707 0, f_or[0], f_fl[0], f_ro[0]))
11710 1, f_or[0], f_fl[0], f_ro[0])),
11714 3, f_or[0], f_fl[0], f_ro[0]))
11717 0, f_or[0], f_fl[0], f_ro[0])),
11721 0, f_or[0], f_fl[0], f_ro[0]))
11724 3, f_or[0], f_fl[0], f_ro[0])),
11728 3, f_or[0], f_fl[0], f_ro[0]))
11731 2, f_or[0], f_fl[0], f_ro[0])),
11736 0, f_or[1], f_fl[1], f_ro[1]))
11739 1, f_or[1], f_fl[1], f_ro[1])),
11743 3, f_or[1], f_fl[1], f_ro[1]))
11746 0, f_or[1], f_fl[1], f_ro[1])),
11750 0, f_or[1], f_fl[1], f_ro[1]))
11753 3, f_or[1], f_fl[1], f_ro[1])),
11757 3, f_or[1], f_fl[1], f_ro[1]))
11760 2, f_or[1], f_fl[1], f_ro[1])),
11765 0, f_or[2], f_fl[2], f_ro[2]))
11768 1, f_or[2], f_fl[2], f_ro[2])),
11772 3, f_or[2], f_fl[2], f_ro[2]))
11775 0, f_or[2], f_fl[2], f_ro[2])),
11779 0, f_or[2], f_fl[2], f_ro[2]))
11782 3, f_or[2], f_fl[2], f_ro[2])),
11786 3, f_or[2], f_fl[2], f_ro[2]))
11789 2, f_or[2], f_fl[2], f_ro[2])),
11794 0, f_or[3], f_fl[3], f_ro[3]))
11797 1, f_or[3], f_fl[3], f_ro[3])),
11801 3, f_or[3], f_fl[3], f_ro[3]))
11804 0, f_or[3], f_fl[3], f_ro[3])),
11808 0, f_or[3], f_fl[3], f_ro[3]))
11811 3, f_or[3], f_fl[3], f_ro[3])),
11815 3, f_or[3], f_fl[3], f_ro[3]))
11818 2, f_or[3], f_fl[3], f_ro[3])),
11823 0, f_or[4], f_fl[4], f_ro[4]))
11826 1, f_or[4], f_fl[4], f_ro[4])),
11830 3, f_or[4], f_fl[4], f_ro[4]))
11833 0, f_or[4], f_fl[4], f_ro[4])),
11837 0, f_or[4], f_fl[4], f_ro[4]))
11840 3, f_or[4], f_fl[4], f_ro[4])),
11844 3, f_or[4], f_fl[4], f_ro[4]))
11847 2, f_or[4], f_fl[4], f_ro[4])),
11852 0, f_or[5], f_fl[5], f_ro[5]))
11855 1, f_or[5], f_fl[5], f_ro[5])),
11859 3, f_or[5], f_fl[5], f_ro[5]))
11862 0, f_or[5], f_fl[5], f_ro[5])),
11866 0, f_or[5], f_fl[5], f_ro[5]))
11869 3, f_or[5], f_fl[5], f_ro[5])),
11873 3, f_or[5], f_fl[5], f_ro[5]))
11876 2, f_or[5], f_fl[5], f_ro[5])),
11886 int line_indices[30];
11887 for (
unsigned int i = 0; i < 30; ++i)
11888 line_indices[i] = lines[i]->
index();
11903 for (
unsigned int i = 0; i < 24; ++i)
11904 if (lines[i]->vertex_index((i + 1) % 2) ==
11906 line_orientation[i] =
11912 Assert(lines[i]->vertex_index(i % 2) ==
11915 line_orientation[i] =
11921 for (
unsigned int i = 24; i < 30; ++i)
11922 line_orientation[i] =
11955 new_quads[0]->set_bounding_object_indices(
11959 line_indices[24]});
11960 new_quads[1]->set_bounding_object_indices(
11964 line_indices[25]});
11965 new_quads[2]->set_bounding_object_indices(
11969 line_indices[20]});
11970 new_quads[3]->set_bounding_object_indices(
11974 line_indices[21]});
11975 new_quads[4]->set_bounding_object_indices(
11979 line_indices[28]});
11980 new_quads[5]->set_bounding_object_indices(
11984 line_indices[29]});
11985 new_quads[6]->set_bounding_object_indices(
11990 new_quads[7]->set_bounding_object_indices(
11995 new_quads[8]->set_bounding_object_indices(
11999 line_indices[26]});
12000 new_quads[9]->set_bounding_object_indices(
12004 line_indices[27]});
12005 new_quads[10]->set_bounding_object_indices(
12009 line_indices[12]});
12010 new_quads[11]->set_bounding_object_indices(
12014 line_indices[13]});
12019 new_quads[0]->set_line_orientation(
12020 0, line_orientation[10]);
12021 new_quads[0]->set_line_orientation(
12022 2, line_orientation[16]);
12024 new_quads[1]->set_line_orientation(
12025 1, line_orientation[14]);
12026 new_quads[1]->set_line_orientation(
12027 2, line_orientation[17]);
12029 new_quads[2]->set_line_orientation(
12030 0, line_orientation[11]);
12031 new_quads[2]->set_line_orientation(
12032 3, line_orientation[20]);
12034 new_quads[3]->set_line_orientation(
12035 1, line_orientation[15]);
12036 new_quads[3]->set_line_orientation(
12037 3, line_orientation[21]);
12039 new_quads[4]->set_line_orientation(
12040 0, line_orientation[18]);
12041 new_quads[4]->set_line_orientation(
12042 2, line_orientation[0]);
12044 new_quads[5]->set_line_orientation(
12045 1, line_orientation[22]);
12046 new_quads[5]->set_line_orientation(
12047 2, line_orientation[1]);
12049 new_quads[6]->set_line_orientation(
12050 0, line_orientation[19]);
12051 new_quads[6]->set_line_orientation(
12052 3, line_orientation[4]);
12054 new_quads[7]->set_line_orientation(
12055 1, line_orientation[23]);
12056 new_quads[7]->set_line_orientation(
12057 3, line_orientation[5]);
12059 new_quads[8]->set_line_orientation(
12060 0, line_orientation[2]);
12061 new_quads[8]->set_line_orientation(
12062 2, line_orientation[8]);
12064 new_quads[9]->set_line_orientation(
12065 1, line_orientation[6]);
12066 new_quads[9]->set_line_orientation(
12067 2, line_orientation[9]);
12069 new_quads[10]->set_line_orientation(
12070 0, line_orientation[3]);
12071 new_quads[10]->set_line_orientation(
12072 3, line_orientation[12]);
12074 new_quads[11]->set_line_orientation(
12075 1, line_orientation[7]);
12076 new_quads[11]->set_line_orientation(
12077 3, line_orientation[13]);
12118 const int quad_indices[36] = {
12119 new_quads[0]->index(),
12120 new_quads[1]->index(),
12121 new_quads[2]->index(),
12122 new_quads[3]->index(),
12123 new_quads[4]->index(),
12124 new_quads[5]->index(),
12125 new_quads[6]->index(),
12126 new_quads[7]->index(),
12127 new_quads[8]->index(),
12128 new_quads[9]->index(),
12129 new_quads[10]->index(),
12130 new_quads[11]->index(),
12132 hex->face(0)->isotropic_child_index(
12134 0, f_or[0], f_fl[0], f_ro[0])),
12135 hex->face(0)->isotropic_child_index(
12137 1, f_or[0], f_fl[0], f_ro[0])),
12138 hex->face(0)->isotropic_child_index(
12140 2, f_or[0], f_fl[0], f_ro[0])),
12141 hex->face(0)->isotropic_child_index(
12143 3, f_or[0], f_fl[0], f_ro[0])),
12145 hex->face(1)->isotropic_child_index(
12147 0, f_or[1], f_fl[1], f_ro[1])),
12148 hex->face(1)->isotropic_child_index(
12150 1, f_or[1], f_fl[1], f_ro[1])),
12151 hex->face(1)->isotropic_child_index(
12153 2, f_or[1], f_fl[1], f_ro[1])),
12154 hex->face(1)->isotropic_child_index(
12156 3, f_or[1], f_fl[1], f_ro[1])),
12158 hex->face(2)->isotropic_child_index(
12160 0, f_or[2], f_fl[2], f_ro[2])),
12161 hex->face(2)->isotropic_child_index(
12163 1, f_or[2], f_fl[2], f_ro[2])),
12164 hex->face(2)->isotropic_child_index(
12166 2, f_or[2], f_fl[2], f_ro[2])),
12167 hex->face(2)->isotropic_child_index(
12169 3, f_or[2], f_fl[2], f_ro[2])),
12171 hex->face(3)->isotropic_child_index(
12173 0, f_or[3], f_fl[3], f_ro[3])),
12174 hex->face(3)->isotropic_child_index(
12176 1, f_or[3], f_fl[3], f_ro[3])),
12177 hex->face(3)->isotropic_child_index(
12179 2, f_or[3], f_fl[3], f_ro[3])),
12180 hex->face(3)->isotropic_child_index(
12182 3, f_or[3], f_fl[3], f_ro[3])),
12184 hex->face(4)->isotropic_child_index(
12186 0, f_or[4], f_fl[4], f_ro[4])),
12187 hex->face(4)->isotropic_child_index(
12189 1, f_or[4], f_fl[4], f_ro[4])),
12190 hex->face(4)->isotropic_child_index(
12192 2, f_or[4], f_fl[4], f_ro[4])),
12193 hex->face(4)->isotropic_child_index(
12195 3, f_or[4], f_fl[4], f_ro[4])),
12197 hex->face(5)->isotropic_child_index(
12199 0, f_or[5], f_fl[5], f_ro[5])),
12200 hex->face(5)->isotropic_child_index(
12202 1, f_or[5], f_fl[5], f_ro[5])),
12203 hex->face(5)->isotropic_child_index(
12205 2, f_or[5], f_fl[5], f_ro[5])),
12206 hex->face(5)->isotropic_child_index(
12208 3, f_or[5], f_fl[5], f_ro[5]))};
12211 new_hexes[0]->set_bounding_object_indices(
12218 new_hexes[1]->set_bounding_object_indices(
12225 new_hexes[2]->set_bounding_object_indices(
12231 quad_indices[10]});
12232 new_hexes[3]->set_bounding_object_indices(
12238 quad_indices[11]});
12241 new_hexes[4]->set_bounding_object_indices(
12247 quad_indices[32]});
12248 new_hexes[5]->set_bounding_object_indices(
12254 quad_indices[33]});
12255 new_hexes[6]->set_bounding_object_indices(
12261 quad_indices[34]});
12262 new_hexes[7]->set_bounding_object_indices(
12268 quad_indices[35]});
12302 for (
unsigned int s = 0;
12309 const unsigned int current_child =
12318 ref_case, f, f_or[f], f_fl[f], f_ro[f]));
12319 new_hexes[current_child]->set_combined_face_orientation(
12325 if (check_for_distorted_cells &&
12326 has_distorted_children<dim, spacedim>(hex))
12343 triangulation.
faces->quads.clear_user_data();
12346 return cells_with_distorted_children;