50 namespace TriangulationImplementation
124 template <
int dim,
int spacedim>
126 cell_is_patch_level_1(
131 unsigned int n_active_children = 0;
132 for (
unsigned int i = 0; i < cell->n_children(); ++i)
133 if (cell->child(i)->is_active())
136 return (n_active_children == 0) ||
137 (n_active_children == cell->n_children());
147 template <
int dim,
int spacedim>
149 cell_will_be_coarsened(
155 if (cell->has_children())
157 unsigned int children_to_coarsen = 0;
158 const unsigned int n_children = cell->n_children();
160 for (
unsigned int c = 0; c < n_children; ++c)
161 if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
162 ++children_to_coarsen;
163 if (children_to_coarsen == n_children)
166 for (
unsigned int c = 0; c < n_children; ++c)
167 if (cell->child(c)->is_active())
168 cell->child(c)->clear_coarsen_flag();
204 template <
int dim,
int spacedim>
206 face_will_be_refined_by_neighbor_internal(
208 const unsigned int face_no,
217 cell->neighbor(face_no);
224 if (neighbor->has_children())
229 if (cell_will_be_coarsened(neighbor))
238 expected_face_ref_case = cell->face(face_no)->refinement_case();
250 const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
258 neighbor->face_orientation(neighbor_neighbor),
259 neighbor->face_flip(neighbor_neighbor),
260 neighbor->face_rotation(neighbor_neighbor));
264 neighbor_face = neighbor->face(neighbor_neighbor);
265 const int this_face_index = cell->face_index(face_no);
271 if (neighbor_face->index() == this_face_index)
277 expected_face_ref_case = face_ref_case;
295 for (
unsigned int c = 0; c < neighbor_face->n_children(); ++c)
296 if (neighbor_face->child_index(c) == this_face_index)
305 if ((neighbor_face->refinement_case() | face_ref_case) ==
306 neighbor_face->refinement_case())
325 expected_face_ref_case =
326 ~neighbor_face->refinement_case();
353 template <
int dim,
int spacedim>
355 face_will_be_refined_by_neighbor(
357 const unsigned int face_no)
360 return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
367 template <
int dim,
int spacedim>
369 face_will_be_refined_by_neighbor(
371 const unsigned int face_no,
374 return face_will_be_refined_by_neighbor_internal(cell,
376 expected_face_ref_case);
381 template <
int dim,
int spacedim>
383 satisfies_level1_at_vertex_rule(
386 std::vector<unsigned int> min_adjacent_cell_level(
388 std::vector<unsigned int> max_adjacent_cell_level(
394 min_adjacent_cell_level[cell->vertex_index(v)] =
395 std::min<unsigned int>(
396 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
397 max_adjacent_cell_level[cell->vertex_index(v)] =
398 std::max<unsigned int>(
399 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
402 for (
unsigned int k = 0; k < triangulation.
n_vertices(); ++k)
404 if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
417 template <
int dim,
int spacedim>
418 std::vector<unsigned int>
423 std::vector<unsigned int> line_cell_count(triangulation.
n_raw_lines(),
426 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
427 ++line_cell_count[cell->line_index(
l)];
428 return line_cell_count;
431 return std::vector<unsigned int>();
442 template <
int dim,
int spacedim>
443 std::vector<unsigned int>
448 std::vector<unsigned int> quad_cell_count(triangulation.
n_raw_quads(),
452 ++quad_cell_count[cell->quad_index(q)];
453 return quad_cell_count;
456 return std::vector<unsigned int>();
480 void reorder_compatibility(std::vector<
CellData<2>> &cells,
483 for (
auto &cell : cells)
485 std::swap(cell.vertices[2], cell.vertices[3]);
489 void reorder_compatibility(std::vector<
CellData<3>> &cells,
493 for (
auto &cell : cells)
497 tmp[i] = cell.vertices[i];
505 std::swap(boundary_quad.vertices[2], boundary_quad.vertices[3]);
527 template <
int dim,
int spacedim>
532 if (line->has_children())
533 return line->child(0)->vertex_index(1);
538 template <
int dim,
int spacedim>
543 switch (static_cast<unsigned char>(quad->refinement_case()))
546 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
549 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
552 return quad->child(0)->vertex_index(3);
561 template <
int dim,
int spacedim>
566 switch (static_cast<unsigned char>(hex->refinement_case()))
569 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
572 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
575 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
578 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
581 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
584 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
587 return hex->child(0)->vertex_index(7);
608 template <
class TRIANGULATION>
609 inline typename TRIANGULATION::DistortedCellList
610 collect_distorted_coarse_cells(
const TRIANGULATION &)
612 return typename TRIANGULATION::DistortedCellList();
634 vertices[i] = cell->vertex(i);
640 if (determinants[i] <= 1
e-9 * std::pow(cell->diameter(), 1. * dim))
642 distorted_cells.distorted_cells.push_back(cell);
647 return distorted_cells;
659 has_distorted_children(
664 for (
unsigned int c = 0; c < cell->n_children(); ++c)
668 vertices[i] = cell->child(c)->vertex(i);
674 if (determinants[i] <=
675 1
e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
690 template <
int dim,
int spacedim>
692 has_distorted_children(
699 template <
int dim,
int spacedim>
701 update_periodic_face_map_recursively(
704 unsigned int n_face_1,
705 unsigned int n_face_2,
706 const std::bitset<3> & orientation,
712 std::bitset<3>>> &periodic_face_map)
715 const FaceIterator face_1 = cell_1->face(n_face_1);
716 const FaceIterator face_2 = cell_2->face(n_face_2);
718 const bool face_orientation = orientation[0];
719 const bool face_flip = orientation[1];
720 const bool face_rotation = orientation[2];
722 Assert((dim != 1) || (face_orientation ==
true && face_flip ==
false &&
723 face_rotation ==
false),
725 "(face_orientation, face_flip, face_rotation) " 726 "is invalid for 1D"));
728 Assert((dim != 2) || (face_orientation ==
true && face_rotation ==
false),
730 "(face_orientation, face_flip, face_rotation) " 731 "is invalid for 2D"));
735 Assert(face_1->at_boundary() && face_2->at_boundary(),
736 ExcMessage(
"Periodic faces must be on the boundary"));
742 Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
747 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
749 const CellFace cell_face_1(cell_1, n_face_1);
750 const CellFace cell_face_2(cell_2, n_face_2);
751 const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
752 cell_face_2, orientation);
754 const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
755 periodic_faces(cell_face_1, cell_face_orientation_2);
759 periodic_face_map.insert(periodic_faces);
763 if (cell_1->has_children())
765 if (cell_2->has_children())
767 update_periodic_face_map_recursively<dim, spacedim>(
768 cell_1->child(n_face_1),
769 cell_2->child(n_face_2),
777 update_periodic_face_map_recursively<dim, spacedim>(
778 cell_1->child(n_face_1),
793 static const int lookup_table_2d[2][2] =
800 static const int lookup_table_3d[2][2][2][4] =
819 if (cell_1->has_children())
821 if (cell_2->has_children())
827 Assert(face_1->n_children() ==
833 for (
unsigned int i = 0;
834 i < GeometryInfo<dim>::max_children_per_face;
842 j = lookup_table_2d[face_flip][i];
845 j = lookup_table_3d[face_orientation][face_flip]
853 unsigned int child_cell_1 =
855 cell_1->refinement_case(),
858 cell_1->face_orientation(n_face_1),
859 cell_1->face_flip(n_face_1),
860 cell_1->face_rotation(n_face_1),
861 face_1->refinement_case());
862 unsigned int child_cell_2 =
864 cell_2->refinement_case(),
867 cell_2->face_orientation(n_face_2),
868 cell_2->face_flip(n_face_2),
869 cell_2->face_rotation(n_face_2),
870 face_2->refinement_case());
872 Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
875 Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
881 update_periodic_face_map_recursively<dim, spacedim>(
882 cell_1->child(child_cell_1),
883 cell_2->child(child_cell_2),
892 for (
unsigned int i = 0;
893 i < GeometryInfo<dim>::max_children_per_face;
897 unsigned int child_cell_1 =
899 cell_1->refinement_case(),
902 cell_1->face_orientation(n_face_1),
903 cell_1->face_flip(n_face_1),
904 cell_1->face_rotation(n_face_1),
905 face_1->refinement_case());
908 update_periodic_face_map_recursively<dim, spacedim>(
909 cell_1->child(child_cell_1),
927 namespace TriangulationImplementation
943 <<
"Something went wrong when making cell " << arg1
944 <<
". Read the docs and the source code " 945 <<
"for more information.");
952 <<
"Something went wrong upon construction of cell " 966 <<
" has negative measure. This typically " 967 <<
"indicates some distortion in the cell, or a mistakenly " 968 <<
"swapped pair of vertices in the input to " 969 <<
"Triangulation::create_triangulation().");
981 <<
"Error while creating cell " << arg1
982 <<
": the vertex index " << arg2 <<
" must be between 0 and " 991 <<
"While trying to assign a boundary indicator to a line: " 992 <<
"the line with end vertices " << arg1 <<
" and " << arg2
993 <<
" does not exist.");
1003 <<
"While trying to assign a boundary indicator to a quad: " 1004 <<
"the quad with bounding lines " << arg1 <<
", " << arg2
1005 <<
", " << arg3 <<
", " << arg4 <<
" does not exist.");
1015 <<
"The input data for creating a triangulation contained " 1016 <<
"information about a line with indices " << arg1 <<
" and " << arg2
1017 <<
" that is described to have boundary indicator " 1018 << static_cast<int>(arg3)
1019 <<
". However, this is an internal line not located on the " 1020 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1022 <<
"If this happened at a place where you call " 1023 <<
"Triangulation::create_triangulation() yourself, you need " 1024 <<
"to check the SubCellData object you pass to this function." 1027 <<
"If this happened in a place where you are reading a mesh " 1028 <<
"from a file, then you need to investigate why such a line " 1029 <<
"ended up in the input file. A typical case is a geometry " 1030 <<
"that consisted of multiple parts and for which the mesh " 1031 <<
"generator program assumes that the interface between " 1032 <<
"two parts is a boundary when that isn't supposed to be " 1033 <<
"the case, or where the mesh generator simply assigns " 1034 <<
"'geometry indicators' to lines at the perimeter of " 1035 <<
"a part that are not supposed to be interpreted as " 1036 <<
"'boundary indicators'.");
1048 <<
"The input data for creating a triangulation contained " 1049 <<
"information about a quad with indices " << arg1 <<
", " << arg2
1050 <<
", " << arg3 <<
", and " << arg4
1051 <<
" that is described to have boundary indicator " 1052 << static_cast<int>(arg5)
1053 <<
". However, this is an internal quad not located on the " 1054 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1056 <<
"If this happened at a place where you call " 1057 <<
"Triangulation::create_triangulation() yourself, you need " 1058 <<
"to check the SubCellData object you pass to this function." 1061 <<
"If this happened in a place where you are reading a mesh " 1062 <<
"from a file, then you need to investigate why such a quad " 1063 <<
"ended up in the input file. A typical case is a geometry " 1064 <<
"that consisted of multiple parts and for which the mesh " 1065 <<
"generator program assumes that the interface between " 1066 <<
"two parts is a boundary when that isn't supposed to be " 1067 <<
"the case, or where the mesh generator simply assigns " 1068 <<
"'geometry indicators' to quads at the surface of " 1069 <<
"a part that are not supposed to be interpreted as " 1070 <<
"'boundary indicators'.");
1079 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1080 <<
" and " << arg2 <<
" appears more than once. " 1081 <<
"This is not allowed.");
1091 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1092 <<
" and " << arg2 <<
" appears multiple times with different (valid) " 1093 << arg3 <<
". This is not allowed.");
1105 <<
"In SubCellData the quad info of the quad with line indices " << arg1
1106 <<
", " << arg2 <<
", " << arg3 <<
" and " << arg4
1107 <<
" appears multiple times with different (valid) " << arg5
1108 <<
". This is not allowed.");
1119 const unsigned int new_quads_in_pairs,
1120 const unsigned int new_quads_single)
1126 unsigned int next_free_single = 0;
1127 unsigned int next_free_pair = 0;
1131 unsigned int n_quads = 0;
1132 unsigned int n_unused_pairs = 0;
1133 unsigned int n_unused_singles = 0;
1134 for (
unsigned int i = 0; i < tria_faces.
quads.
used.size(); ++i)
1138 else if (i + 1 < tria_faces.
quads.
used.size())
1143 if (next_free_single == 0)
1144 next_free_single = i;
1149 if (next_free_pair == 0)
1157 Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
1162 const int additional_single_quads = new_quads_single - n_unused_singles;
1164 unsigned int new_size =
1165 tria_faces.
quads.
used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
1166 if (additional_single_quads > 0)
1167 new_size += additional_single_quads;
1177 new_size * GeometryInfo<2>::lines_per_cell -
1206 const unsigned int total_cells,
1207 const unsigned int dimension,
1208 const unsigned int space_dimension)
1260 if (dimension < space_dimension)
1271 tria_level.
parents.reserve((total_cells + 1) / 2);
1273 (total_cells + 1) / 2 -
1277 tria_level.
neighbors.reserve(total_cells * (2 * dimension));
1279 total_cells * (2 * dimension) -
1281 std::make_pair(-1, -1));
1284 if (tria_level.
dim == 3)
1290 total_cells * GeometryInfo<3>::faces_per_cell -
1294 else if (tria_level.
dim == 2)
1300 total_cells * GeometryInfo<2>::faces_per_cell -
1305 if (tria_level.
dim == 2 || tria_level.
dim == 3)
1325 <<
"The containers have sizes " << arg1 <<
" and " << arg2
1326 <<
", which is not as expected.");
1335 const unsigned int true_dimension)
1338 (void)true_dimension;
1365 const unsigned int new_objects_in_pairs,
1366 const unsigned int new_objects_single = 0)
1378 unsigned int n_objects = 0;
1379 unsigned int n_unused_pairs = 0;
1380 unsigned int n_unused_singles = 0;
1381 for (
unsigned int i = 0; i < tria_objects.
used.size(); ++i)
1383 if (tria_objects.
used[i])
1385 else if (i + 1 < tria_objects.
used.size())
1387 if (tria_objects.
used[i + 1])
1404 Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
1405 tria_objects.
used.size(),
1410 const int additional_single_objects =
1411 new_objects_single - n_unused_singles;
1413 unsigned int new_size = tria_objects.
used.size() +
1414 new_objects_in_pairs - 2 * n_unused_pairs;
1415 if (additional_single_objects > 0)
1416 new_size += additional_single_objects;
1419 if (new_size > tria_objects.
n_objects())
1421 unsigned int faces_per_cell = 1;
1422 unsigned int max_children_per_cell = 1;
1442 tria_objects.
cells.reserve(new_size * faces_per_cell);
1443 tria_objects.
cells.insert(tria_objects.
cells.end(),
1448 tria_objects.
used.reserve(new_size);
1449 tria_objects.
used.insert(tria_objects.
used.end(),
1450 new_size - tria_objects.
used.size(),
1459 const unsigned int factor = max_children_per_cell / 2;
1460 tria_objects.
children.reserve(factor * new_size);
1480 tria_objects.
user_data.reserve(new_size);
1481 tria_objects.
user_data.resize(new_size);
1490 if (n_unused_singles == 0)
1498 const unsigned int new_hexes = new_objects_in_pairs;
1500 const unsigned int new_size =
1501 new_hexes + std::count(tria_objects.
used.begin(),
1502 tria_objects.
used.end(),
1506 if (new_size > tria_objects.
n_objects())
1508 unsigned int faces_per_cell = 1;
1519 tria_objects.
cells.reserve(new_size * faces_per_cell);
1520 tria_objects.
cells.insert(tria_objects.
cells.end(),
1525 tria_objects.
used.reserve(new_size);
1526 tria_objects.
used.insert(tria_objects.
used.end(),
1527 new_size - tria_objects.
used.size(),
1536 tria_objects.
children.reserve(4 * new_size);
1555 tria_objects.
user_data.reserve(new_size);
1556 tria_objects.
user_data.resize(new_size);
1580 tria_object.
used.size()));
1621 template <
int dim,
int spacedim>
1628 virtual ~
Policy() =
default;
1643 std::vector<unsigned int> & line_cell_count,
1644 std::vector<unsigned int> &quad_cell_count) = 0;
1651 const bool check_for_distorted_cells) = 0;
1657 prevent_distorted_boundary_cells(
1664 prepare_refinement_dim_dependent(
1680 virtual std::unique_ptr<Policy<dim, spacedim>>
1691 template <
int dim,
int spacedim,
typename T>
1698 T::update_neighbors(tria);
1705 std::vector<unsigned int> & line_cell_count,
1706 std::vector<unsigned int> &quad_cell_count)
override 1708 T::delete_children(tria, cell, line_cell_count, quad_cell_count);
1713 const bool check_for_distorted_cells)
override 1715 return T::execute_refinement(triangulation, check_for_distorted_cells);
1722 T::prevent_distorted_boundary_cells(triangulation);
1729 T::prepare_refinement_dim_dependent(triangulation);
1737 return T::template coarsening_allowed<dim, spacedim>(cell);
1740 std::unique_ptr<Policy<dim, spacedim>>
1743 return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
1857 template <
int dim,
int spacedim>
1861 const unsigned int level_objects,
1864 using line_iterator =
1868 if (level_objects > 0)
1899 line_iterator(triangulation.
end_line()) :
1901 for (; line != endc; ++line)
1904 if (line->has_children() ==
false)
1920 line_iterator line = triangulation.
begin_line(),
1922 for (; line != endc; ++line)
1925 if (line->has_children() ==
false)
1945 template <
int dim,
int spacedim>
1949 const unsigned int level_objects,
1960 &compute_number_cache<dim, spacedim>),
1966 using quad_iterator =
1983 unsigned int n_levels = 0;
1984 if (level_objects > 0)
1988 n_levels =
level + 1;
2001 (
level == n_levels - 1 ?
2002 quad_iterator(triangulation.
end_quad()) :
2004 for (; quad != endc; ++quad)
2007 if (quad->has_children() ==
false)
2023 quad_iterator quad = triangulation.
begin_quad(),
2025 for (; quad != endc; ++quad)
2028 if (quad->has_children() ==
false)
2034 update_lines.
join();
2052 template <
int dim,
int spacedim>
2056 const unsigned int level_objects,
2067 &compute_number_cache<dim, spacedim>),
2073 using hex_iterator =
2091 unsigned int n_levels = 0;
2092 if (level_objects > 0)
2096 n_levels =
level + 1;
2108 endc = (
level == n_levels - 1 ?
2109 hex_iterator(triangulation.
end_hex()) :
2111 for (; hex != endc; ++hex)
2114 if (hex->has_children() ==
false)
2130 hex_iterator hex = triangulation.
begin_hex(),
2131 endc = triangulation.
end_hex();
2132 for (; hex != endc; ++hex)
2135 if (hex->has_children() ==
false)
2141 update_quads_and_lines.
join();
2146 template <
int spacedim>
2151 template <
int dim,
int spacedim>
2206 static const unsigned int left_right_offset[2][6][2] = {
2215 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
2226 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
2227 adjacent_cells(2 * triangulation.
n_raw_faces(), dummy);
2230 for (
auto f : cell->face_indices())
2235 const unsigned int offset =
2236 (cell->direction_flag() ?
2237 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2239 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2241 adjacent_cells[2 * face->index() + offset] = cell;
2251 if (cell->is_active() && face->has_children())
2253 adjacent_cells[2 * face->child(0)->index() + offset] =
2255 adjacent_cells[2 * face->child(1)->index() + offset] =
2281 if (face->has_children() &&
2282 (cell->is_active() ||
2284 cell->refinement_case(), f) ==
2287 for (
unsigned int c = 0; c < face->n_children(); ++c)
2288 adjacent_cells[2 * face->child(c)->index() + offset] =
2290 if (face->child(0)->has_children())
2292 adjacent_cells[2 * face->child(0)->child(0)->index() +
2294 adjacent_cells[2 * face->child(0)->child(1)->index() +
2297 if (face->child(1)->has_children())
2299 adjacent_cells[2 * face->child(1)->child(0)->index() +
2301 adjacent_cells[2 * face->child(1)->child(1)->index() +
2314 for (
auto f : cell->face_indices())
2316 const unsigned int offset =
2317 (cell->direction_flag() ?
2318 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2320 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2322 f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
2330 template <
int dim,
int spacedim>
2344 if (dim == spacedim)
2345 for (
unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2353 const double cell_measure = GridTools::cell_measure<spacedim>(
2368 tria.
faces = std::make_unique<
2376 const auto connectivity = build_connectivity<unsigned int>(cells);
2377 const unsigned int n_cell = cells.size();
2382 auto &lines_0 = tria.
faces->lines;
2385 const auto & crs = connectivity.entity_to_entities(1, 0);
2386 const unsigned int n_lines = crs.ptr.size() - 1;
2389 reserve_space_(lines_0, n_lines);
2392 for (
unsigned int line = 0; line < n_lines; ++line)
2393 for (
unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
2402 auto &quads_0 = tria.
faces->quads;
2403 auto &faces = *tria.
faces;
2406 const auto & crs = connectivity.entity_to_entities(2, 1);
2407 const unsigned int n_quads = crs.ptr.size() - 1;
2410 reserve_space_(quads_0, n_quads);
2411 reserve_space_(faces, 2 , n_quads);
2414 for (
unsigned int q = 0, k = 0; q < n_quads; ++q)
2417 faces.quad_reference_cell[q] = connectivity.entity_types(2)[q];
2420 for (
unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
2428 faces.quads_line_orientations
2430 connectivity.entity_orientations(1)[k];
2437 auto &cells_0 = tria.
levels[0]->cells;
2441 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2442 const auto &nei = connectivity.entity_to_entities(dim, dim);
2446 const bool orientation_needed =
2449 std::any_of(connectivity.entity_orientations(1).begin(),
2450 connectivity.entity_orientations(1).end(),
2451 [](
const auto &i) {
return i == 0; }));
2454 reserve_space_(cells_0, n_cell);
2455 reserve_space_(
level, spacedim, n_cell, orientation_needed);
2458 for (
unsigned int cell = 0; cell < n_cell; ++cell)
2461 cells_0.boundary_or_material_id[cell].material_id =
2462 cells[cell].material_id;
2465 cells_0.manifold_id[cell] = cells[cell].manifold_id;
2468 level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
2471 for (
unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2475 if (nei.col[i] != static_cast<unsigned int>(-1))
2477 j] = {0, nei.col[i]};
2484 if (orientation_needed)
2485 level.face_orientations
2487 connectivity.entity_orientations(dim - 1)[i];
2495 auto &bids_face = dim == 3 ?
2496 tria.
faces->quads.boundary_or_material_id :
2497 tria.
faces->lines.boundary_or_material_id;
2500 std::vector<unsigned int> count(bids_face.size(), 0);
2503 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2506 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
2507 for (
unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
2508 count[crs.col[i]]++;
2511 for (
unsigned int face = 0; face < count.size(); ++face)
2513 if (count[face] != 1)
2517 bids_face[face].boundary_id = 0;
2523 const auto &crs = connectivity.entity_to_entities(2, 1);
2524 for (
unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
2525 tria.
faces->lines.boundary_or_material_id[crs.col[i]]
2531 static const unsigned int t_tba =
static_cast<unsigned int>(-1);
2532 static const unsigned int t_inner =
static_cast<unsigned int>(-2);
2534 std::vector<unsigned int> type(
vertices.size(), t_tba);
2536 const auto &crs = connectivity.entity_to_entities(1, 0);
2538 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
2539 for (
unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2541 if (type[crs.col[i]] != t_inner)
2542 type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
2544 for (
unsigned int face = 0; face < type.size(); ++face)
2549 if (type[face] != t_inner && type[face] != t_tba)
2556 process_subcelldata(connectivity.entity_to_entities(1, 0),
2562 process_subcelldata(connectivity.entity_to_entities(2, 0),
2568 template <
int structdim,
typename T>
2577 if (boundary_objects_in.size() == 0)
2581 auto boundary_objects = boundary_objects_in;
2584 for (
auto &boundary_object : boundary_objects)
2585 std::sort(boundary_object.vertices.begin(),
2586 boundary_object.vertices.end());
2589 std::sort(boundary_objects.begin(),
2590 boundary_objects.end(),
2591 [](
const auto &a,
const auto &
b) {
2592 return a.vertices <
b.vertices;
2595 unsigned int counter = 0;
2597 std::vector<unsigned int> key;
2600 for (
unsigned int o = 0; o < obj.
n_objects(); ++o)
2614 key.assign(crs.
col.data() + crs.
ptr[o],
2615 crs.
col.data() + crs.
ptr[o + 1]);
2616 std::sort(key.begin(), key.end());
2619 const auto subcell_object =
2621 boundary_objects.end(),
2623 [&](
const auto &cell,
const auto &key) {
2624 return cell.vertices < key;
2628 if (subcell_object == boundary_objects.end() ||
2629 subcell_object->vertices != key)
2638 if (subcell_object->boundary_id !=
2656 const unsigned structdim,
2657 const unsigned int size)
2659 const unsigned int dim = faces.
dim;
2661 const unsigned int faces_per_cell =
2666 if (dim == 3 && structdim == 2)
2681 const unsigned int spacedim,
2682 const unsigned int size,
2683 const bool orientation_needed)
2685 const unsigned int dim = level.
dim;
2687 const unsigned int faces_per_cell =
2699 level.
parents.assign((size + 1) / 2, -1);
2704 level.
neighbors.assign(size * faces_per_cell, {-1, -1});
2708 if (orientation_needed)
2722 const unsigned int structdim = obj.
structdim;
2724 const unsigned int max_children_per_cell =
2729 const unsigned int faces_per_cell =
2734 obj.
used.assign(size,
true);
2738 BoundaryOrMaterialId());
2746 obj.
children.assign(max_children_per_cell / 2 * size, -1);
2748 obj.
cells.assign(faces_per_cell * size, -1);
2778 template <
int spacedim>
2782 std::vector<unsigned int> &,
2783 std::vector<unsigned int> &)
2785 const unsigned int dim = 1;
2797 Assert(!cell->child(0)->has_children() &&
2798 !cell->child(1)->has_children(),
2804 if (cell->neighbor(0)->has_children())
2811 neighbor = neighbor->child(1);
2814 Assert(neighbor->neighbor(1) == cell->child(0),
2816 neighbor->set_neighbor(1, cell);
2822 if (neighbor->has_children())
2823 neighbor = neighbor->child(1);
2832 if (cell->neighbor(1)->has_children())
2839 neighbor = neighbor->child(0);
2842 Assert(neighbor->neighbor(0) == cell->child(1),
2844 neighbor->set_neighbor(0, cell);
2850 if (neighbor->has_children())
2851 neighbor = neighbor->child(0);
2861 triangulation.
vertices_used[cell->child(0)->vertex_index(1)] =
false;
2867 for (
unsigned int child = 0; child < cell->n_children(); ++child)
2870 cell->child(child)->clear_user_flag();
2871 cell->child(child)->clear_used_flag();
2876 cell->clear_children();
2877 cell->clear_user_flag();
2882 template <
int spacedim>
2886 std::vector<unsigned int> & line_cell_count,
2887 std::vector<unsigned int> &)
2889 const unsigned int dim = 2;
2897 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2900 lines_to_delete.reserve(4 * 2 + 4);
2905 for (
unsigned int c = 0; c < cell->n_children(); ++c)
2909 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
2910 --line_cell_count[child->line_index(
l)];
2925 .
vertices_used[cell->child(0)->line(1)->vertex_index(1)] =
false;
2927 lines_to_delete.push_back(cell->child(0)->line(1));
2928 lines_to_delete.push_back(cell->child(0)->line(3));
2929 lines_to_delete.push_back(cell->child(3)->line(0));
2930 lines_to_delete.push_back(cell->child(3)->line(2));
2934 unsigned int inner_face_no =
2939 lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2943 for (
unsigned int child = 0; child < cell->n_children(); ++child)
2945 cell->child(child)->clear_user_data();
2946 cell->child(child)->clear_user_flag();
2947 cell->child(child)->clear_used_flag();
2952 cell->clear_children();
2953 cell->clear_refinement_case();
2954 cell->clear_user_flag();
2960 for (
unsigned int line_no = 0;
2961 line_no < GeometryInfo<dim>::lines_per_cell;
2965 cell->line(line_no);
2967 if (line->has_children())
2972 Assert((line_cell_count[line->child_index(0)] == 0 &&
2973 line_cell_count[line->child_index(1)] == 0) ||
2974 (line_cell_count[line->child_index(0)] > 0 &&
2975 line_cell_count[line->child_index(1)] > 0),
2978 if (line_cell_count[line->child_index(0)] == 0)
2980 for (
unsigned int c = 0; c < 2; ++c)
2981 Assert(!line->child(c)->has_children(),
2991 lines_to_delete.push_back(line->child(0));
2992 lines_to_delete.push_back(line->child(1));
2994 line->clear_children();
3006 typename std::vector<
3008 line = lines_to_delete.
begin(),
3009 endline = lines_to_delete.end();
3010 for (; line != endline; ++line)
3012 (*line)->clear_user_data();
3013 (*line)->clear_user_flag();
3014 (*line)->clear_used_flag();
3020 template <
int spacedim>
3024 std::vector<unsigned int> & line_cell_count,
3025 std::vector<unsigned int> & quad_cell_count)
3027 const unsigned int dim = 3;
3039 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3041 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3044 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3045 quads_to_delete.reserve(6 * 4 + 12);
3049 for (
unsigned int c = 0; c < cell->n_children(); ++c)
3053 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
3054 --line_cell_count[child->line_index(
l)];
3056 --quad_cell_count[child->quad_index(f)];
3070 quads_to_delete.push_back(cell->child(0)->face(1));
3073 quads_to_delete.push_back(cell->child(0)->face(3));
3076 quads_to_delete.push_back(cell->child(0)->face(5));
3079 quads_to_delete.push_back(cell->child(0)->face(1));
3080 quads_to_delete.push_back(cell->child(0)->face(3));
3081 quads_to_delete.push_back(cell->child(3)->face(0));
3082 quads_to_delete.push_back(cell->child(3)->face(2));
3084 lines_to_delete.push_back(cell->child(0)->line(11));
3087 quads_to_delete.push_back(cell->child(0)->face(1));
3088 quads_to_delete.push_back(cell->child(0)->face(5));
3089 quads_to_delete.push_back(cell->child(3)->face(0));
3090 quads_to_delete.push_back(cell->child(3)->face(4));
3092 lines_to_delete.push_back(cell->child(0)->line(5));
3095 quads_to_delete.push_back(cell->child(0)->face(3));
3096 quads_to_delete.push_back(cell->child(0)->face(5));
3097 quads_to_delete.push_back(cell->child(3)->face(2));
3098 quads_to_delete.push_back(cell->child(3)->face(4));
3100 lines_to_delete.push_back(cell->child(0)->line(7));
3103 quads_to_delete.push_back(cell->child(0)->face(1));
3104 quads_to_delete.push_back(cell->child(2)->face(1));
3105 quads_to_delete.push_back(cell->child(4)->face(1));
3106 quads_to_delete.push_back(cell->child(6)->face(1));
3108 quads_to_delete.push_back(cell->child(0)->face(3));
3109 quads_to_delete.push_back(cell->child(1)->face(3));
3110 quads_to_delete.push_back(cell->child(4)->face(3));
3111 quads_to_delete.push_back(cell->child(5)->face(3));
3113 quads_to_delete.push_back(cell->child(0)->face(5));
3114 quads_to_delete.push_back(cell->child(1)->face(5));
3115 quads_to_delete.push_back(cell->child(2)->face(5));
3116 quads_to_delete.push_back(cell->child(3)->face(5));
3118 lines_to_delete.push_back(cell->child(0)->line(5));
3119 lines_to_delete.push_back(cell->child(0)->line(7));
3120 lines_to_delete.push_back(cell->child(0)->line(11));
3121 lines_to_delete.push_back(cell->child(7)->line(0));
3122 lines_to_delete.push_back(cell->child(7)->line(2));
3123 lines_to_delete.push_back(cell->child(7)->line(8));
3129 triangulation.
vertices_used[cell->child(0)->vertex_index(7)] =
3141 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3143 cell->child(child)->clear_user_data();
3144 cell->child(child)->clear_user_flag();
3151 cell->child(child)->set_face_orientation(f,
true);
3152 cell->child(child)->set_face_flip(f,
false);
3153 cell->child(child)->set_face_rotation(f,
false);
3156 cell->child(child)->clear_used_flag();
3161 cell->clear_children();
3162 cell->clear_refinement_case();
3163 cell->clear_user_flag();
3176 cell->face(quad_no);
3180 quad->has_children()) ||
3185 switch (quad->refinement_case())
3197 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3198 quad_cell_count[quad->child_index(1)] == 0) ||
3199 (quad_cell_count[quad->child_index(0)] > 0 &&
3200 quad_cell_count[quad->child_index(1)] > 0),
3206 unsigned int deleted_grandchildren = 0;
3207 unsigned int number_of_child_refinements = 0;
3209 for (
unsigned int c = 0; c < 2; ++c)
3210 if (quad->child(c)->has_children())
3212 ++number_of_child_refinements;
3217 (quad_cell_count[quad->child(c)->child_index(0)] ==
3219 quad_cell_count[quad->child(c)->child_index(1)] ==
3221 (quad_cell_count[quad->child(c)->child_index(0)] >
3223 quad_cell_count[quad->child(c)->child_index(1)] >
3226 if (quad_cell_count[quad->child(c)->child_index(0)] ==
3233 Assert(quad->refinement_case() +
3234 quad->child(c)->refinement_case() ==
3242 quads_to_delete.push_back(
3243 quad->child(c)->child(0));
3244 quads_to_delete.push_back(
3245 quad->child(c)->child(1));
3246 if (quad->child(c)->refinement_case() ==
3248 lines_to_delete.push_back(
3249 quad->child(c)->child(0)->line(1));
3251 lines_to_delete.push_back(
3252 quad->child(c)->child(0)->line(3));
3253 quad->child(c)->clear_children();
3254 quad->child(c)->clear_refinement_case();
3255 ++deleted_grandchildren;
3263 if (number_of_child_refinements > 0 &&
3264 deleted_grandchildren == number_of_child_refinements)
3269 middle_line = quad->child(0)->line(1);
3271 middle_line = quad->child(0)->line(3);
3273 lines_to_delete.push_back(middle_line->child(0));
3274 lines_to_delete.push_back(middle_line->child(1));
3277 middle_line)] =
false;
3278 middle_line->clear_children();
3283 if (quad_cell_count[quad->child_index(0)] == 0)
3289 quads_to_delete.push_back(quad->child(0));
3290 quads_to_delete.push_back(quad->child(1));
3292 lines_to_delete.push_back(quad->child(0)->line(1));
3294 lines_to_delete.push_back(quad->child(0)->line(3));
3314 if (quad->child(0)->has_children())
3315 if (quad->refinement_case() ==
3349 quad_iterator switch_1 =
3350 quad->child(0)->child(1),
3352 quad->child(1)->child(0);
3354 Assert(!switch_1->has_children(),
3356 Assert(!switch_2->has_children(),
3359 const int switch_1_index = switch_1->index();
3360 const int switch_2_index = switch_2->index();
3361 for (
unsigned int l = 0;
3362 l < triangulation.
levels.size();
3364 for (
unsigned int h = 0;
3366 triangulation.
levels[
l]->cells.n_objects();
3368 for (
const unsigned int q :
3373 ->cells.get_bounding_object_indices(
3375 if (index == switch_1_index)
3377 ->cells.get_bounding_object_indices(
3378 h)[q] = switch_2_index;
3379 else if (index == switch_2_index)
3381 ->cells.get_bounding_object_indices(
3382 h)[q] = switch_1_index;
3387 const int switch_1_lines[4] = {
3388 static_cast<signed int>(
3389 switch_1->line_index(0)),
3390 static_cast<signed int>(
3391 switch_1->line_index(1)),
3392 static_cast<signed int>(
3393 switch_1->line_index(2)),
3394 static_cast<signed int>(
3395 switch_1->line_index(3))};
3396 const bool switch_1_line_orientations[4] = {
3397 switch_1->line_orientation(0),
3398 switch_1->line_orientation(1),
3399 switch_1->line_orientation(2),
3400 switch_1->line_orientation(3)};
3402 switch_1->boundary_id();
3403 const unsigned int switch_1_user_index =
3404 switch_1->user_index();
3405 const bool switch_1_user_flag =
3406 switch_1->user_flag_set();
3408 switch_1->set_bounding_object_indices(
3409 {switch_2->line_index(0),
3410 switch_2->line_index(1),
3411 switch_2->line_index(2),
3412 switch_2->line_index(3)});
3413 switch_1->set_line_orientation(
3414 0, switch_2->line_orientation(0));
3415 switch_1->set_line_orientation(
3416 1, switch_2->line_orientation(1));
3417 switch_1->set_line_orientation(
3418 2, switch_2->line_orientation(2));
3419 switch_1->set_line_orientation(
3420 3, switch_2->line_orientation(3));
3421 switch_1->set_boundary_id_internal(
3422 switch_2->boundary_id());
3423 switch_1->set_manifold_id(
3424 switch_2->manifold_id());
3425 switch_1->set_user_index(switch_2->user_index());
3426 if (switch_2->user_flag_set())
3427 switch_1->set_user_flag();
3429 switch_1->clear_user_flag();
3431 switch_2->set_bounding_object_indices(
3435 switch_1_lines[3]});
3436 switch_2->set_line_orientation(
3437 0, switch_1_line_orientations[0]);
3438 switch_2->set_line_orientation(
3439 1, switch_1_line_orientations[1]);
3440 switch_2->set_line_orientation(
3441 2, switch_1_line_orientations[2]);
3442 switch_2->set_line_orientation(
3443 3, switch_1_line_orientations[3]);
3444 switch_2->set_boundary_id_internal(
3445 switch_1_boundary_id);
3446 switch_2->set_manifold_id(
3447 switch_1->manifold_id());
3448 switch_2->set_user_index(switch_1_user_index);
3449 if (switch_1_user_flag)
3450 switch_2->set_user_flag();
3452 switch_2->clear_user_flag();
3454 const unsigned int child_0 =
3455 quad->child(0)->child_index(0);
3456 const unsigned int child_2 =
3457 quad->child(1)->child_index(0);
3458 quad->clear_children();
3459 quad->clear_refinement_case();
3460 quad->set_refinement_case(
3462 quad->set_children(0, child_0);
3463 quad->set_children(2, child_2);
3465 quad_cell_count[child_2]);
3480 const unsigned int child_0 =
3481 quad->child(0)->child_index(0);
3482 const unsigned int child_2 =
3483 quad->child(1)->child_index(0);
3484 quad->clear_children();
3485 quad->clear_refinement_case();
3486 quad->set_refinement_case(
3488 quad->set_children(0, child_0);
3489 quad->set_children(2, child_2);
3493 quad->clear_children();
3494 quad->clear_refinement_case();
3505 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3506 quad_cell_count[quad->child_index(1)] == 0 &&
3507 quad_cell_count[quad->child_index(2)] == 0 &&
3508 quad_cell_count[quad->child_index(3)] == 0) ||
3509 (quad_cell_count[quad->child_index(0)] > 0 &&
3510 quad_cell_count[quad->child_index(1)] > 0 &&
3511 quad_cell_count[quad->child_index(2)] > 0 &&
3512 quad_cell_count[quad->child_index(3)] > 0),
3515 if (quad_cell_count[quad->child_index(0)] == 0)
3521 lines_to_delete.push_back(quad->child(0)->line(1));
3522 lines_to_delete.push_back(quad->child(3)->line(0));
3523 lines_to_delete.push_back(quad->child(0)->line(3));
3524 lines_to_delete.push_back(quad->child(3)->line(2));
3526 for (
unsigned int child = 0; child < quad->n_children();
3528 quads_to_delete.push_back(quad->child(child));
3534 quad->clear_children();
3535 quad->clear_refinement_case();
3554 for (
unsigned int line_no = 0;
3555 line_no < GeometryInfo<dim>::lines_per_cell;
3559 cell->line(line_no);
3563 line->has_children()) ||
3568 if (line->has_children())
3573 Assert((line_cell_count[line->child_index(0)] == 0 &&
3574 line_cell_count[line->child_index(1)] == 0) ||
3575 (line_cell_count[line->child_index(0)] > 0 &&
3576 line_cell_count[line->child_index(1)] > 0),
3579 if (line_cell_count[line->child_index(0)] == 0)
3581 for (
unsigned int c = 0; c < 2; ++c)
3582 Assert(!line->child(c)->has_children(),
3592 lines_to_delete.push_back(line->child(0));
3593 lines_to_delete.push_back(line->child(1));
3595 line->clear_children();
3607 typename std::vector<
3609 line = lines_to_delete.begin(),
3610 endline = lines_to_delete.end();
3611 for (; line != endline; ++line)
3613 (*line)->clear_user_data();
3614 (*line)->clear_user_flag();
3615 (*line)->clear_used_flag();
3618 typename std::vector<
3620 quad = quads_to_delete.begin(),
3621 endquad = quads_to_delete.end();
3622 for (; quad != endquad; ++quad)
3624 (*quad)->clear_user_data();
3625 (*quad)->clear_children();
3626 (*quad)->clear_refinement_case();
3627 (*quad)->clear_user_flag();
3628 (*quad)->clear_used_flag();
3650 template <
int spacedim>
3653 unsigned int & next_unused_vertex,
3660 const unsigned int dim = 2;
3663 cell->clear_refine_flag();
3760 int new_vertices[9];
3761 for (
unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3762 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3763 for (
unsigned int line_no = 0; line_no < 4; ++line_no)
3764 if (cell->line(line_no)->has_children())
3765 new_vertices[4 + line_no] =
3766 cell->line(line_no)->child(0)->vertex_index(1);
3775 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
3776 ++next_unused_vertex;
3778 next_unused_vertex < triangulation.
vertices.size(),
3780 "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."));
3783 new_vertices[8] = next_unused_vertex;
3795 if (dim == spacedim)
3798 triangulation.
vertices[next_unused_vertex] = cell->center(
true);
3805 if (cell->user_flag_set())
3808 cell->clear_user_flag();
3809 triangulation.
vertices[next_unused_vertex] =
3810 cell->center(
true,
true);
3820 cell->clear_user_flag();
3825 triangulation.
vertices[next_unused_vertex] =
3826 cell->center(
true,
true);
3833 unsigned int lmin = 8;
3834 unsigned int lmax = 12;
3841 for (
unsigned int l = lmin;
l < lmax; ++
l)
3843 while (next_unused_line->used() ==
true)
3845 new_lines[
l] = next_unused_line;
3849 new_lines[
l]->used() ==
false,
3851 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3866 for (
unsigned int c = 0; c < 2; ++c, ++
l)
3867 new_lines[l] = cell->line(face_no)->child(c);
3870 new_lines[8]->set_bounding_object_indices(
3871 {new_vertices[6], new_vertices[8]});
3872 new_lines[9]->set_bounding_object_indices(
3873 {new_vertices[8], new_vertices[7]});
3874 new_lines[10]->set_bounding_object_indices(
3875 {new_vertices[4], new_vertices[8]});
3876 new_lines[11]->set_bounding_object_indices(
3877 {new_vertices[8], new_vertices[5]});
3886 new_lines[0] = cell->line(0);
3887 new_lines[1] = cell->line(1);
3888 new_lines[2] = cell->line(2)->child(0);
3889 new_lines[3] = cell->line(2)->child(1);
3890 new_lines[4] = cell->line(3)->child(0);
3891 new_lines[5] = cell->line(3)->child(1);
3892 new_lines[6]->set_bounding_object_indices(
3893 {new_vertices[6], new_vertices[7]});
3903 new_lines[0] = cell->line(0)->child(0);
3904 new_lines[1] = cell->line(0)->child(1);
3905 new_lines[2] = cell->line(1)->child(0);
3906 new_lines[3] = cell->line(1)->child(1);
3907 new_lines[4] = cell->line(2);
3908 new_lines[5] = cell->line(3);
3909 new_lines[6]->set_bounding_object_indices(
3910 {new_vertices[4], new_vertices[5]});
3913 for (
unsigned int l = lmin;
l < lmax; ++
l)
3915 new_lines[
l]->set_used_flag();
3916 new_lines[
l]->clear_user_flag();
3918 new_lines[
l]->clear_children();
3920 new_lines[
l]->set_boundary_id_internal(
3922 new_lines[
l]->set_manifold_id(cell->manifold_id());
3929 while (next_unused_cell->used() ==
true)
3933 for (
unsigned int i = 0; i < n_children; ++i)
3936 next_unused_cell->used() ==
false,
3938 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3939 subcells[i] = next_unused_cell;
3941 if (i % 2 == 1 && i < n_children - 1)
3942 while (next_unused_cell->used() ==
true)
3960 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3961 new_lines[8]->index(),
3962 new_lines[4]->index(),
3963 new_lines[10]->index()});
3964 subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3965 new_lines[2]->index(),
3966 new_lines[5]->index(),
3967 new_lines[11]->index()});
3968 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3969 new_lines[9]->index(),
3970 new_lines[10]->index(),
3971 new_lines[6]->index()});
3972 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3973 new_lines[3]->index(),
3974 new_lines[11]->index(),
3975 new_lines[7]->index()});
3991 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3992 new_lines[6]->index(),
3993 new_lines[2]->index(),
3994 new_lines[4]->index()});
3995 subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3996 new_lines[1]->index(),
3997 new_lines[3]->index(),
3998 new_lines[5]->index()});
4015 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4016 new_lines[2]->index(),
4017 new_lines[4]->index(),
4018 new_lines[6]->index()});
4019 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4020 new_lines[3]->index(),
4021 new_lines[6]->index(),
4022 new_lines[5]->index()});
4027 for (
unsigned int i = 0; i < n_children; ++i)
4029 subcells[i]->set_used_flag();
4030 subcells[i]->clear_refine_flag();
4031 subcells[i]->clear_user_flag();
4032 subcells[i]->clear_user_data();
4033 subcells[i]->clear_children();
4036 subcells[i]->set_material_id(cell->material_id());
4037 subcells[i]->set_manifold_id(cell->manifold_id());
4038 subcells[i]->set_subdomain_id(subdomainid);
4041 subcells[i]->set_parent(cell->index());
4049 for (
unsigned int i = 0; i < n_children / 2; ++i)
4050 cell->set_children(2 * i, subcells[2 * i]->index());
4052 cell->set_refinement_case(ref_case);
4060 for (
unsigned int c = 0; c < n_children; ++c)
4061 cell->child(c)->set_direction_flag(cell->direction_flag());
4066 template <
int dim,
int spacedim>
4069 const bool check_for_distorted_cells)
4076 endc = triangulation.
end();
4077 for (; cell != endc; ++cell)
4079 if (cell->refine_flag_set())
4081 triangulation.
levels.push_back(
4093 line->clear_user_flag();
4094 line->clear_user_data();
4097 unsigned int n_single_lines = 0;
4098 unsigned int n_lines_in_pairs = 0;
4099 unsigned int needed_vertices = 0;
4105 unsigned int needed_cells = 0;
4107 for (
const auto &cell :
4109 if (cell->refine_flag_set())
4111 if (cell->reference_cell() ==
4115 needed_vertices += 0;
4116 n_single_lines += 3;
4118 else if (cell->reference_cell() ==
4122 needed_vertices += 1;
4123 n_single_lines += 4;
4130 for (
const auto line_no : cell->face_indices())
4132 auto line = cell->line(line_no);
4133 if (line->has_children() ==
false)
4134 line->set_user_flag();
4139 const unsigned int used_cells =
4140 std::count(triangulation.
levels[
level + 1]->cells.used.begin(),
4141 triangulation.
levels[
level + 1]->cells.used.end(),
4146 used_cells + needed_cells,
4158 if (line->user_flag_set())
4161 n_lines_in_pairs += 2;
4162 needed_vertices += 1;
4167 needed_vertices += std::count(triangulation.
vertices_used.begin(),
4171 if (needed_vertices > triangulation.
vertices.size())
4177 unsigned int next_unused_vertex = 0;
4186 for (; line != endl; ++line)
4187 if (line->user_flag_set())
4193 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
4194 ++next_unused_vertex;
4196 next_unused_vertex < triangulation.
vertices.size(),
4198 "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."));
4201 triangulation.
vertices[next_unused_vertex] = line->center(
true);
4203 bool pair_found =
false;
4205 for (; next_unused_line != endl; ++next_unused_line)
4206 if (!next_unused_line->used() &&
4207 !(++next_unused_line)->used())
4215 line->set_children(0, next_unused_line->index());
4218 children[2] = {next_unused_line, ++next_unused_line};
4221 children[0]->used() ==
false,
4223 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4225 children[1]->used() ==
false,
4227 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4229 children[0]->set_bounding_object_indices(
4230 {line->vertex_index(0), next_unused_vertex});
4231 children[1]->set_bounding_object_indices(
4232 {next_unused_vertex, line->vertex_index(1)});
4234 children[0]->set_used_flag();
4235 children[1]->set_used_flag();
4236 children[0]->clear_children();
4237 children[1]->clear_children();
4240 children[0]->clear_user_flag();
4241 children[1]->clear_user_flag();
4244 children[0]->set_boundary_id_internal(line->boundary_id());
4245 children[1]->set_boundary_id_internal(line->boundary_id());
4247 children[0]->set_manifold_id(line->manifold_id());
4248 children[1]->set_manifold_id(line->manifold_id());
4250 line->clear_user_flag();
4257 cells_with_distorted_children;
4263 unsigned int &next_unused_vertex,
4264 auto & next_unused_line,
4265 auto & next_unused_cell,
4266 const auto & cell) {
4267 const auto ref_case = cell->refine_flag_set();
4268 cell->clear_refine_flag();
4270 unsigned int n_new_vertices = 0;
4274 else if (cell->reference_cell() ==
4280 std::vector<int> new_vertices(n_new_vertices);
4281 for (
unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
4283 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4284 for (
unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
4285 if (cell->line(line_no)->has_children())
4286 new_vertices[cell->n_vertices() + line_no] =
4287 cell->line(line_no)->child(0)->vertex_index(1);
4291 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
4292 ++next_unused_vertex;
4294 next_unused_vertex < triangulation.
vertices.size(),
4296 "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."));
4299 new_vertices[8] = next_unused_vertex;
4301 if (dim == spacedim)
4303 triangulation.
vertices[next_unused_vertex] =
4306 if (cell->user_flag_set())
4308 cell->clear_user_flag();
4309 triangulation.
vertices[next_unused_vertex] =
4310 cell->center(
true,
true);
4315 cell->clear_user_flag();
4317 triangulation.
vertices[next_unused_vertex] =
4318 cell->center(
true,
true);
4322 std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
4325 unsigned int lmin = 0;
4326 unsigned int lmax = 0;
4333 else if (cell->reference_cell() ==
4344 for (
unsigned int l = lmin;
l < lmax; ++
l)
4346 while (next_unused_line->used() ==
true)
4348 new_lines[
l] = next_unused_line;
4352 new_lines[
l]->used() ==
false,
4354 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4362 const auto ref = [&](
const unsigned int face_no,
4363 const unsigned int vertex_no) {
4364 if (cell->line(face_no)->child(0)->vertex_index(0) ==
4365 static_cast<unsigned int>(new_vertices[vertex_no]) ||
4366 cell->line(face_no)->child(0)->vertex_index(1) ==
4367 static_cast<unsigned int>(new_vertices[vertex_no]))
4369 new_lines[2 * face_no + 0] =
4370 cell->line(face_no)->child(0);
4371 new_lines[2 * face_no + 1] =
4372 cell->line(face_no)->child(1);
4376 new_lines[2 * face_no + 0] =
4377 cell->line(face_no)->child(1);
4378 new_lines[2 * face_no + 1] =
4379 cell->line(face_no)->child(0);
4387 new_lines[6]->set_bounding_object_indices(
4388 {new_vertices[3], new_vertices[4]});
4389 new_lines[7]->set_bounding_object_indices(
4390 {new_vertices[4], new_vertices[5]});
4391 new_lines[8]->set_bounding_object_indices(
4392 {new_vertices[5], new_vertices[3]});
4394 else if (cell->reference_cell() ==
4398 for (
const unsigned int face_no : cell->face_indices())
4399 for (
unsigned int c = 0; c < 2; ++c, ++
l)
4400 new_lines[l] = cell->line(face_no)->child(c);
4402 new_lines[8]->set_bounding_object_indices(
4403 {new_vertices[6], new_vertices[8]});
4404 new_lines[9]->set_bounding_object_indices(
4405 {new_vertices[8], new_vertices[7]});
4406 new_lines[10]->set_bounding_object_indices(
4407 {new_vertices[4], new_vertices[8]});
4408 new_lines[11]->set_bounding_object_indices(
4409 {new_vertices[8], new_vertices[5]});
4418 for (
unsigned int l = lmin;
l < lmax; ++
l)
4420 new_lines[
l]->set_used_flag();
4421 new_lines[
l]->clear_user_flag();
4422 new_lines[
l]->clear_user_data();
4423 new_lines[
l]->clear_children();
4425 new_lines[
l]->set_boundary_id_internal(
4427 new_lines[
l]->set_manifold_id(cell->manifold_id());
4432 while (next_unused_cell->used() ==
true)
4435 unsigned int n_children = 0;
4439 else if (cell->reference_cell() ==
4445 for (
unsigned int i = 0; i < n_children; ++i)
4448 next_unused_cell->used() ==
false,
4450 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4451 subcells[i] = next_unused_cell;
4453 if (i % 2 == 1 && i < n_children - 1)
4454 while (next_unused_cell->used() ==
true)
4461 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4462 new_lines[8]->index(),
4463 new_lines[5]->index()});
4464 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4465 new_lines[2]->index(),
4466 new_lines[6]->index()});
4467 subcells[2]->set_bounding_object_indices({new_lines[7]->index(),
4468 new_lines[3]->index(),
4469 new_lines[4]->index()});
4470 subcells[3]->set_bounding_object_indices({new_lines[6]->index(),
4471 new_lines[7]->index(),
4472 new_lines[8]->index()});
4476 const auto ref = [&](
const unsigned int line_no,
4477 const unsigned int vertex_no,
4478 const unsigned int subcell_no,
4479 const unsigned int subcell_line_no) {
4480 if (new_lines[line_no]->vertex_index(1) !=
4481 static_cast<unsigned int>(new_vertices[vertex_no]))
4482 triangulation.
levels[subcells[subcell_no]->level()]
4483 ->face_orientations[subcells[subcell_no]->index() *
4485 subcell_line_no] = 0;
4509 else if ((dim == 2) && (cell->reference_cell() ==
4512 subcells[0]->set_bounding_object_indices(
4513 {new_lines[0]->index(),
4514 new_lines[8]->index(),
4515 new_lines[4]->index(),
4516 new_lines[10]->index()});
4517 subcells[1]->set_bounding_object_indices(
4518 {new_lines[8]->index(),
4519 new_lines[2]->index(),
4520 new_lines[5]->index(),
4521 new_lines[11]->index()});
4522 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
4523 new_lines[9]->index(),
4524 new_lines[10]->index(),
4525 new_lines[6]->index()});
4526 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
4527 new_lines[3]->index(),
4528 new_lines[11]->index(),
4529 new_lines[7]->index()});
4538 for (
unsigned int i = 0; i < n_children; ++i)
4540 subcells[i]->set_used_flag();
4541 subcells[i]->clear_refine_flag();
4542 subcells[i]->clear_user_flag();
4543 subcells[i]->clear_user_data();
4544 subcells[i]->clear_children();
4547 subcells[i]->set_material_id(cell->material_id());
4548 subcells[i]->set_manifold_id(cell->manifold_id());
4549 subcells[i]->set_subdomain_id(subdomainid);
4553 triangulation.
levels[subcells[i]->level()]
4554 ->reference_cell[subcells[i]->index()] = cell->reference_cell();
4557 subcells[i]->set_parent(cell->index());
4560 for (
unsigned int i = 0; i < n_children / 2; ++i)
4561 cell->set_children(2 * i, subcells[2 * i]->index());
4563 cell->set_refinement_case(ref_case);
4566 for (
unsigned int c = 0; c < n_children; ++c)
4567 cell->child(c)->set_direction_flag(cell->direction_flag());
4571 level < static_cast<int>(triangulation.
levels.size()) - 1;
4577 for (
const auto &cell :
4579 if (cell->refine_flag_set())
4581 if (cell->at_boundary())
4582 cell->set_user_flag();
4584 create_children(triangulation,
4590 if (cell->reference_cell() ==
4592 check_for_distorted_cells &&
4593 has_distorted_children<dim, spacedim>(cell))
4597 triangulation.
signals.post_refinement_on_cell(cell);
4601 return cells_with_distorted_children;
4610 template <
int spacedim>
4615 const unsigned int dim = 1;
4623 endc = triangulation.
end();
4624 for (; cell != endc; ++cell)
4626 if (cell->refine_flag_set())
4628 triangulation.
levels.push_back(
4640 unsigned int needed_vertices = 0;
4645 unsigned int flagged_cells = 0;
4647 for (
const auto &acell :
4649 if (acell->refine_flag_set())
4654 const unsigned int used_cells =
4655 std::count(triangulation.
levels[
level + 1]->cells.used.begin(),
4656 triangulation.
levels[
level + 1]->cells.used.end(),
4674 needed_vertices += flagged_cells;
4679 needed_vertices += std::count(triangulation.
vertices_used.begin(),
4685 if (needed_vertices > triangulation.
vertices.size())
4696 unsigned int next_unused_vertex = 0;
4703 for (
const auto &cell :
4705 if (cell->refine_flag_set())
4708 cell->clear_refine_flag();
4714 ++next_unused_vertex;
4716 next_unused_vertex < triangulation.
vertices.size(),
4718 "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."));
4723 triangulation.
vertices[next_unused_vertex] =
4733 while (next_unused_cell->used() ==
true)
4735 first_child = next_unused_cell;
4736 first_child->set_used_flag();
4737 first_child->clear_user_data();
4740 next_unused_cell->used() ==
false,
4742 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4743 second_child = next_unused_cell;
4744 second_child->set_used_flag();
4745 second_child->clear_user_data();
4750 cell->set_children(0, first_child->index());
4751 first_child->clear_children();
4752 first_child->set_bounding_object_indices(
4753 {cell->vertex_index(0), next_unused_vertex});
4754 first_child->set_material_id(cell->material_id());
4755 first_child->set_manifold_id(cell->manifold_id());
4756 first_child->set_subdomain_id(subdomainid);
4757 first_child->set_direction_flag(cell->direction_flag());
4759 first_child->set_parent(cell->index());
4763 first_child->face(1)->set_manifold_id(cell->manifold_id());
4768 first_child->set_neighbor(1, second_child);
4770 first_child->set_neighbor(0, cell->neighbor(0));
4771 else if (cell->neighbor(0)->is_active())
4777 Assert(cell->neighbor(0)->level() <= cell->level(),
4779 first_child->set_neighbor(0, cell->neighbor(0));
4785 const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4786 first_child->set_neighbor(0,
4787 cell->neighbor(0)->child(nbnb));
4792 left_neighbor = cell->neighbor(0);
4793 while (left_neighbor->has_children())
4795 left_neighbor = left_neighbor->child(nbnb);
4796 left_neighbor->set_neighbor(nbnb, first_child);
4801 second_child->clear_children();
4802 second_child->set_bounding_object_indices(
4803 {next_unused_vertex, cell->vertex_index(1)});
4804 second_child->set_neighbor(0, first_child);
4805 second_child->set_material_id(cell->material_id());
4806 second_child->set_manifold_id(cell->manifold_id());
4807 second_child->set_subdomain_id(subdomainid);
4808 second_child->set_direction_flag(cell->direction_flag());
4811 second_child->set_neighbor(1, cell->neighbor(1));
4812 else if (cell->neighbor(1)->is_active())
4814 Assert(cell->neighbor(1)->level() <= cell->level(),
4816 second_child->set_neighbor(1, cell->neighbor(1));
4821 const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4822 second_child->set_neighbor(
4823 1, cell->neighbor(1)->child(nbnb));
4826 right_neighbor = cell->neighbor(1);
4827 while (right_neighbor->has_children())
4829 right_neighbor = right_neighbor->child(nbnb);
4830 right_neighbor->set_neighbor(nbnb, second_child);
4834 triangulation.
signals.post_refinement_on_cell(cell);
4850 template <
int spacedim>
4853 const bool check_for_distorted_cells)
4855 const unsigned int dim = 2;
4859 bool flag_isotropic_mesh =
true;
4861 cell = triangulation.
begin(),
4862 endc = triangulation.
end();
4863 for (; cell != endc; ++cell)
4868 flag_isotropic_mesh =
false;
4872 if (flag_isotropic_mesh)
4873 return execute_refinement_isotropic(triangulation,
4874 check_for_distorted_cells);
4883 endc = triangulation.
end();
4884 for (; cell != endc; ++cell)
4886 if (cell->refine_flag_set())
4888 triangulation.
levels.push_back(
4902 line->clear_user_flag();
4903 line->clear_user_data();
4908 unsigned int n_single_lines = 0;
4912 unsigned int n_lines_in_pairs = 0;
4918 unsigned int needed_vertices = 0;
4923 unsigned int needed_cells = 0;
4925 for (
const auto &cell :
4927 if (cell->refine_flag_set())
4938 n_single_lines += 4;
4950 n_single_lines += 1;
4958 for (
const unsigned int line_no :
4962 cell->refine_flag_set(), line_no) ==
4966 line = cell->line(line_no);
4967 if (line->has_children() ==
false)
4968 line->set_user_flag();
4975 const unsigned int used_cells =
4976 std::count(triangulation.
levels[
level + 1]->cells.used.begin(),
4977 triangulation.
levels[
level + 1]->cells.used.end(),
4985 used_cells + needed_cells,
5001 if (line->user_flag_set())
5004 n_lines_in_pairs += 2;
5005 needed_vertices += 1;
5017 needed_vertices += std::count(triangulation.
vertices_used.begin(),
5023 if (needed_vertices > triangulation.
vertices.size())
5034 unsigned int next_unused_vertex = 0;
5046 for (; line != endl; ++line)
5047 if (line->user_flag_set())
5053 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5054 ++next_unused_vertex;
5056 next_unused_vertex < triangulation.
vertices.size(),
5058 "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."));
5061 triangulation.
vertices[next_unused_vertex] = line->center(
true);
5066 bool pair_found =
false;
5068 for (; next_unused_line != endl; ++next_unused_line)
5069 if (!next_unused_line->used() &&
5070 !(++next_unused_line)->used())
5083 line->set_children(0, next_unused_line->index());
5087 children[2] = {next_unused_line, ++next_unused_line};
5091 children[0]->used() ==
false,
5093 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5095 children[1]->used() ==
false,
5097 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5099 children[0]->set_bounding_object_indices(
5100 {line->vertex_index(0), next_unused_vertex});
5101 children[1]->set_bounding_object_indices(
5102 {next_unused_vertex, line->vertex_index(1)});
5104 children[0]->set_used_flag();
5105 children[1]->set_used_flag();
5106 children[0]->clear_children();
5107 children[1]->clear_children();
5110 children[0]->clear_user_flag();
5111 children[1]->clear_user_flag();
5114 children[0]->set_boundary_id_internal(line->boundary_id());
5115 children[1]->set_boundary_id_internal(line->boundary_id());
5117 children[0]->set_manifold_id(line->manifold_id());
5118 children[1]->set_manifold_id(line->manifold_id());
5122 line->clear_user_flag();
5134 cells_with_distorted_children;
5142 level < static_cast<int>(triangulation.
levels.size()) - 1;
5148 for (
const auto &cell :
5150 if (cell->refine_flag_set())
5158 if (cell->at_boundary())
5159 cell->set_user_flag();
5163 create_children(triangulation,
5169 if (check_for_distorted_cells &&
5170 has_distorted_children<dim, spacedim>(cell))
5174 triangulation.
signals.post_refinement_on_cell(cell);
5178 return cells_with_distorted_children;
5186 template <
int spacedim>
5189 const bool check_for_distorted_cells)
5191 const unsigned int dim = 3;
5205 endc = triangulation.
end();
5206 for (; cell != endc; ++cell)
5208 if (cell->refine_flag_set())
5210 triangulation.
levels.push_back(
5220 triangulation.
faces->quads.clear_user_data();
5226 line->clear_user_flag();
5231 quad->clear_user_flag();
5254 unsigned int needed_vertices = 0;
5255 unsigned int needed_lines_single = 0;
5256 unsigned int needed_quads_single = 0;
5257 unsigned int needed_lines_pair = 0;
5258 unsigned int needed_quads_pair = 0;
5263 unsigned int new_cells = 0;
5265 for (
const auto &acell :
5267 if (acell->refine_flag_set())
5277 ++needed_quads_single;
5285 ++needed_lines_single;
5286 needed_quads_single += 4;
5293 needed_lines_single += 6;
5294 needed_quads_single += 12;
5308 for (
const unsigned int face :
5312 aface = acell->face(face);
5319 acell->face_orientation(face),
5320 acell->face_flip(face),
5321 acell->face_rotation(face));
5326 if (face_ref_case ==
5329 if (aface->number_of_children() < 4)
5332 aface->set_user_flag();
5334 else if (aface->refinement_case() != face_ref_case)
5345 Assert(aface->refinement_case() ==
5347 dim - 1>::isotropic_refinement ||
5348 aface->refinement_case() ==
5351 aface->set_user_index(face_ref_case);
5357 for (
unsigned int line = 0;
5358 line < GeometryInfo<dim>::lines_per_cell;
5362 !acell->line(line)->has_children())
5363 acell->line(line)->set_user_flag();
5369 const unsigned int used_cells =
5370 std::count(triangulation.
levels[
level + 1]->cells.used.begin(),
5371 triangulation.
levels[
level + 1]->cells.used.end(),
5379 used_cells + new_cells,
5393 if (quad->user_flag_set())
5399 needed_quads_pair += 4;
5400 needed_lines_pair += 4;
5401 needed_vertices += 1;
5403 if (quad->user_index())
5407 needed_quads_pair += 2;
5408 needed_lines_single += 1;
5421 if (quad->has_children())
5423 Assert(quad->refinement_case() ==
5426 if ((face_refinement_cases[quad->user_index()] ==
5428 (quad->child(0)->line_index(1) + 1 !=
5429 quad->child(2)->line_index(1))) ||
5430 (face_refinement_cases[quad->user_index()] ==
5432 (quad->child(0)->line_index(3) + 1 !=
5433 quad->child(1)->line_index(3))))
5434 needed_lines_pair += 2;
5443 if (line->user_flag_set())
5445 needed_lines_pair += 2;
5446 needed_vertices += 1;
5452 needed_lines_single);
5456 needed_quads_single);
5459 needed_quads_single);
5463 needed_vertices += std::count(triangulation.
vertices_used.begin(),
5469 if (needed_vertices > triangulation.
vertices.size())
5488 if (!cell->refine_flag_set())
5489 for (
unsigned int line = 0;
5490 line < GeometryInfo<dim>::lines_per_cell;
5492 if (cell->line(line)->has_children())
5493 for (
unsigned int c = 0; c < 2; ++c)
5494 Assert(cell->line(line)->child(c)->user_flag_set() ==
false,
5506 unsigned int next_unused_vertex = 0;
5517 for (; line != endl; ++line)
5518 if (line->user_flag_set())
5524 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5525 ++next_unused_vertex;
5527 next_unused_vertex < triangulation.
vertices.size(),
5529 "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."));
5532 triangulation.
vertices[next_unused_vertex] = line->center(
true);
5538 triangulation.
faces->lines.template next_free_pair_object<1>(
5546 line->set_children(0, next_unused_line->index());
5550 children[2] = {next_unused_line, ++next_unused_line};
5555 children[0]->used() ==
false,
5557 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5559 children[1]->used() ==
false,
5561 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5563 children[0]->set_bounding_object_indices(
5564 {line->vertex_index(0), next_unused_vertex});
5565 children[1]->set_bounding_object_indices(
5566 {next_unused_vertex, line->vertex_index(1)});
5568 children[0]->set_used_flag();
5569 children[1]->set_used_flag();
5570 children[0]->clear_children();
5571 children[1]->clear_children();
5574 children[0]->clear_user_flag();
5575 children[1]->clear_user_flag();
5577 children[0]->set_boundary_id_internal(line->boundary_id());
5578 children[1]->set_boundary_id_internal(line->boundary_id());
5580 children[0]->set_manifold_id(line->manifold_id());
5581 children[1]->set_manifold_id(line->manifold_id());
5586 line->clear_user_flag();
5629 for (; quad != endq; ++quad)
5631 if (quad->user_index())
5634 face_refinement_cases[quad->user_index()];
5643 if (aniso_quad_ref_case == quad->refinement_case())
5646 Assert(quad->refinement_case() ==
5648 quad->refinement_case() ==
5653 Assert(quad->user_index() ==
5655 quad->user_index() ==
5664 triangulation.
faces->lines
5667 new_line->used() ==
false,
5669 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5684 unsigned int vertex_indices[2];
5688 quad->line(2)->child(0)->vertex_index(1);
5690 quad->line(3)->child(0)->vertex_index(1);
5695 quad->line(0)->child(0)->vertex_index(1);
5697 quad->line(1)->child(0)->vertex_index(1);
5700 new_line->set_bounding_object_indices(
5701 {vertex_indices[0], vertex_indices[1]});
5702 new_line->set_used_flag();
5703 new_line->clear_user_flag();
5705 new_line->clear_children();
5706 new_line->set_boundary_id_internal(quad->boundary_id());
5707 new_line->set_manifold_id(quad->manifold_id());
5715 const unsigned int index[2][2] = {
5725 triangulation.
faces->quads
5727 new_quads[0] = next_unused_quad;
5729 new_quads[0]->used() ==
false,
5731 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5734 new_quads[1] = next_unused_quad;
5736 new_quads[1]->used() ==
false,
5738 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5743 new_quads[0]->set_bounding_object_indices(
5744 {
static_cast<int>(quad->line_index(0)),
5747 ->child(index[0][quad->line_orientation(2)])
5750 ->child(index[0][quad->line_orientation(3)])
5752 new_quads[1]->set_bounding_object_indices(
5754 static_cast<int>(quad->line_index(1)),
5756 ->child(index[1][quad->line_orientation(2)])
5759 ->child(index[1][quad->line_orientation(3)])
5764 new_quads[0]->set_bounding_object_indices(
5766 ->child(index[0][quad->line_orientation(0)])
5769 ->child(index[0][quad->line_orientation(1)])
5771 static_cast<int>(quad->line_index(2)),
5772 new_line->index()});
5773 new_quads[1]->set_bounding_object_indices(
5775 ->child(index[1][quad->line_orientation(0)])
5778 ->child(index[1][quad->line_orientation(1)])
5781 static_cast<int>(quad->line_index(3))});
5784 for (
const auto &new_quad : new_quads)
5786 new_quad->set_used_flag();
5787 new_quad->clear_user_flag();
5788 new_quad->clear_user_data();
5789 new_quad->clear_children();
5790 new_quad->set_boundary_id_internal(quad->boundary_id());
5791 new_quad->set_manifold_id(quad->manifold_id());
5795 for (
unsigned int j = 0;
5796 j < GeometryInfo<dim>::lines_per_face;
5798 new_quad->set_line_orientation(j,
true);
5804 new_quads[0]->set_line_orientation(
5805 0, quad->line_orientation(0));
5806 new_quads[0]->set_line_orientation(
5807 2, quad->line_orientation(2));
5808 new_quads[1]->set_line_orientation(
5809 1, quad->line_orientation(1));
5810 new_quads[1]->set_line_orientation(
5811 3, quad->line_orientation(3));
5814 new_quads[0]->set_line_orientation(
5815 3, quad->line_orientation(3));
5816 new_quads[1]->set_line_orientation(
5817 2, quad->line_orientation(2));
5821 new_quads[0]->set_line_orientation(
5822 1, quad->line_orientation(1));
5823 new_quads[1]->set_line_orientation(
5824 0, quad->line_orientation(0));
5830 if (quad->refinement_case() ==
5853 if (aniso_quad_ref_case ==
5856 old_child[0] = quad->child(0)->line(1);
5857 old_child[1] = quad->child(2)->line(1);
5861 Assert(aniso_quad_ref_case ==
5865 old_child[0] = quad->child(0)->line(3);
5866 old_child[1] = quad->child(1)->line(3);
5869 if (old_child[0]->index() + 1 != old_child[1]->index())
5875 spacedim>::raw_line_iterator
5878 new_child[0] = new_child[1] =
5879 triangulation.
faces->lines
5880 .template next_free_pair_object<1>(
5884 new_child[0]->set_used_flag();
5885 new_child[1]->set_used_flag();
5887 const int old_index_0 = old_child[0]->index(),
5888 old_index_1 = old_child[1]->index(),
5889 new_index_0 = new_child[0]->index(),
5890 new_index_1 = new_child[1]->index();
5894 for (
unsigned int q = 0;
5895 q < triangulation.
faces->quads.n_objects();
5897 for (
unsigned int l = 0;
5898 l < GeometryInfo<dim>::lines_per_face;
5901 const int this_index =
5902 triangulation.
faces->quads
5903 .get_bounding_object_indices(q)[
l];
5904 if (this_index == old_index_0)
5905 triangulation.
faces->quads
5906 .get_bounding_object_indices(q)[
l] =
5908 else if (this_index == old_index_1)
5909 triangulation.
faces->quads
5910 .get_bounding_object_indices(q)[
l] =
5915 for (
unsigned int i = 0; i < 2; ++i)
5917 Assert(!old_child[i]->has_children(),
5920 new_child[i]->set_bounding_object_indices(
5921 {old_child[i]->vertex_index(0),
5922 old_child[i]->vertex_index(1)});
5923 new_child[i]->set_boundary_id_internal(
5925 new_child[i]->set_manifold_id(
5927 new_child[i]->set_user_index(
5928 old_child[i]->user_index());
5929 if (old_child[i]->user_flag_set())
5930 new_child[i]->set_user_flag();
5932 new_child[i]->clear_user_flag();
5934 new_child[i]->clear_children();
5936 old_child[i]->clear_user_flag();
5937 old_child[i]->clear_user_index();
5938 old_child[i]->clear_used_flag();
5944 if (aniso_quad_ref_case ==
5947 new_line->set_children(
5948 0, quad->child(0)->line_index(1));
5949 Assert(new_line->child(1) ==
5950 quad->child(2)->line(1),
5987 quad_iterator switch_1 = quad->child(1),
5988 switch_2 = quad->child(2);
5989 const int switch_1_index = switch_1->index();
5990 const int switch_2_index = switch_2->index();
5991 for (
unsigned int l = 0;
5992 l < triangulation.
levels.size();
5994 for (
unsigned int h = 0;
5996 triangulation.
levels[
l]->cells.n_objects();
5998 for (
const unsigned int q :
6001 const int face_index =
6003 ->cells.get_bounding_object_indices(
6005 if (face_index == switch_1_index)
6007 ->cells.get_bounding_object_indices(
6008 h)[q] = switch_2_index;
6009 else if (face_index == switch_2_index)
6011 ->cells.get_bounding_object_indices(
6012 h)[q] = switch_1_index;
6016 const unsigned int switch_1_lines[4] = {
6017 switch_1->line_index(0),
6018 switch_1->line_index(1),
6019 switch_1->line_index(2),
6020 switch_1->line_index(3)};
6021 const bool switch_1_line_orientations[4] = {
6022 switch_1->line_orientation(0),
6023 switch_1->line_orientation(1),
6024 switch_1->line_orientation(2),
6025 switch_1->line_orientation(3)};
6027 switch_1->boundary_id();
6028 const unsigned int switch_1_user_index =
6029 switch_1->user_index();
6030 const bool switch_1_user_flag =
6031 switch_1->user_flag_set();
6033 switch_1_refinement_case =
6034 switch_1->refinement_case();
6035 const int switch_1_first_child_pair =
6036 (switch_1_refinement_case ?
6037 switch_1->child_index(0) :
6039 const int switch_1_second_child_pair =
6040 (switch_1_refinement_case ==
6042 switch_1->child_index(2) :
6045 switch_1->set_bounding_object_indices(
6046 {switch_2->line_index(0),
6047 switch_2->line_index(1),
6048 switch_2->line_index(2),
6049 switch_2->line_index(3)});
6050 switch_1->set_line_orientation(
6051 0, switch_2->line_orientation(0));
6052 switch_1->set_line_orientation(
6053 1, switch_2->line_orientation(1));
6054 switch_1->set_line_orientation(
6055 2, switch_2->line_orientation(2));
6056 switch_1->set_line_orientation(
6057 3, switch_2->line_orientation(3));
6058 switch_1->set_boundary_id_internal(
6059 switch_2->boundary_id());
6060 switch_1->set_manifold_id(switch_2->manifold_id());
6061 switch_1->set_user_index(switch_2->user_index());
6062 if (switch_2->user_flag_set())
6063 switch_1->set_user_flag();
6065 switch_1->clear_user_flag();
6066 switch_1->clear_refinement_case();
6067 switch_1->set_refinement_case(
6068 switch_2->refinement_case());
6069 switch_1->clear_children();
6070 if (switch_2->refinement_case())
6071 switch_1->set_children(0,
6072 switch_2->child_index(0));
6073 if (switch_2->refinement_case() ==
6075 switch_1->set_children(2,
6076 switch_2->child_index(2));
6078 switch_2->set_bounding_object_indices(
6082 switch_1_lines[3]});
6083 switch_2->set_line_orientation(
6084 0, switch_1_line_orientations[0]);
6085 switch_2->set_line_orientation(
6086 1, switch_1_line_orientations[1]);
6087 switch_2->set_line_orientation(
6088 2, switch_1_line_orientations[2]);
6089 switch_2->set_line_orientation(
6090 3, switch_1_line_orientations[3]);
6091 switch_2->set_boundary_id_internal(
6092 switch_1_boundary_id);
6093 switch_2->set_manifold_id(switch_1->manifold_id());
6094 switch_2->set_user_index(switch_1_user_index);
6095 if (switch_1_user_flag)
6096 switch_2->set_user_flag();
6098 switch_2->clear_user_flag();
6099 switch_2->clear_refinement_case();
6100 switch_2->set_refinement_case(
6101 switch_1_refinement_case);
6102 switch_2->clear_children();
6103 switch_2->set_children(0,
6104 switch_1_first_child_pair);
6105 switch_2->set_children(2,
6106 switch_1_second_child_pair);
6108 new_quads[0]->set_refinement_case(
6110 new_quads[0]->set_children(0, quad->child_index(0));
6111 new_quads[1]->set_refinement_case(
6113 new_quads[1]->set_children(0, quad->child_index(2));
6117 new_quads[0]->set_refinement_case(
6119 new_quads[0]->set_children(0, quad->child_index(0));
6120 new_quads[1]->set_refinement_case(
6122 new_quads[1]->set_children(0, quad->child_index(2));
6123 new_line->set_children(
6124 0, quad->child(0)->line_index(3));
6125 Assert(new_line->child(1) ==
6126 quad->child(1)->line(3),
6129 quad->clear_children();
6133 quad->set_children(0, new_quads[0]->index());
6135 quad->set_refinement_case(aniso_quad_ref_case);
6142 if (quad->user_flag_set())
6154 ++next_unused_vertex;
6156 next_unused_vertex < triangulation.
vertices.size(),
6158 "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."));
6166 quad->refinement_case();
6172 quad->child(0)->set_user_index(
6174 quad->child(1)->set_user_index(
6180 middle_line = quad->child(0)->line(1);
6182 middle_line = quad->child(0)->line(3);
6189 if (!middle_line->has_children())
6196 triangulation.
vertices[next_unused_vertex] =
6197 middle_line->center(
true);
6204 triangulation.
faces->lines
6205 .template next_free_pair_object<1>(
6210 middle_line->set_children(
6211 0, next_unused_line->index());
6215 raw_line_iterator children[2] = {
6216 next_unused_line, ++next_unused_line};
6222 children[0]->used() ==
false,
6224 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6226 children[1]->used() ==
false,
6228 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6230 children[0]->set_bounding_object_indices(
6231 {middle_line->vertex_index(0),
6232 next_unused_vertex});
6233 children[1]->set_bounding_object_indices(
6234 {next_unused_vertex,
6235 middle_line->vertex_index(1)});
6237 children[0]->set_used_flag();
6238 children[1]->set_used_flag();
6239 children[0]->clear_children();
6240 children[1]->clear_children();
6243 children[0]->clear_user_flag();
6244 children[1]->clear_user_flag();
6246 children[0]->set_boundary_id_internal(
6247 middle_line->boundary_id());
6248 children[1]->set_boundary_id_internal(
6249 middle_line->boundary_id());
6251 children[0]->set_manifold_id(
6252 middle_line->manifold_id());
6253 children[1]->set_manifold_id(
6254 middle_line->manifold_id());
6260 quad->clear_user_flag();
6289 triangulation.
vertices[next_unused_vertex] =
6290 quad->center(
true,
true);
6299 for (
unsigned int i = 0; i < 4; ++i)
6309 triangulation.
faces->lines
6312 new_lines[i] = next_unused_line;
6316 new_lines[i]->used() ==
false,
6318 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6338 const unsigned int vertex_indices[5] = {
6339 quad->line(0)->child(0)->vertex_index(1),
6340 quad->line(1)->child(0)->vertex_index(1),
6341 quad->line(2)->child(0)->vertex_index(1),
6342 quad->line(3)->child(0)->vertex_index(1),
6343 next_unused_vertex};
6345 new_lines[0]->set_bounding_object_indices(
6346 {vertex_indices[2], vertex_indices[4]});
6347 new_lines[1]->set_bounding_object_indices(
6348 {vertex_indices[4], vertex_indices[3]});
6349 new_lines[2]->set_bounding_object_indices(
6350 {vertex_indices[0], vertex_indices[4]});
6351 new_lines[3]->set_bounding_object_indices(
6352 {vertex_indices[4], vertex_indices[1]});
6354 for (
const auto &new_line : new_lines)
6356 new_line->set_used_flag();
6357 new_line->clear_user_flag();
6359 new_line->clear_children();
6360 new_line->set_boundary_id_internal(quad->boundary_id());
6361 new_line->set_manifold_id(quad->manifold_id());
6380 const unsigned int index[2][2] = {
6384 const int line_indices[12] = {
6386 ->child(index[0][quad->line_orientation(0)])
6389 ->child(index[1][quad->line_orientation(0)])
6392 ->child(index[0][quad->line_orientation(1)])
6395 ->child(index[1][quad->line_orientation(1)])
6398 ->child(index[0][quad->line_orientation(2)])
6401 ->child(index[1][quad->line_orientation(2)])
6404 ->child(index[0][quad->line_orientation(3)])
6407 ->child(index[1][quad->line_orientation(3)])
6409 new_lines[0]->index(),
6410 new_lines[1]->index(),
6411 new_lines[2]->index(),
6412 new_lines[3]->index()};
6421 triangulation.
faces->quads
6424 new_quads[0] = next_unused_quad;
6426 new_quads[0]->used() ==
false,
6428 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6431 new_quads[1] = next_unused_quad;
6433 new_quads[1]->used() ==
false,
6435 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6438 triangulation.
faces->quads
6440 new_quads[2] = next_unused_quad;
6442 new_quads[2]->used() ==
false,
6444 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6447 new_quads[3] = next_unused_quad;
6449 new_quads[3]->used() ==
false,
6451 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6454 quad->set_children(0, new_quads[0]->index());
6455 quad->set_children(2, new_quads[2]->index());
6458 new_quads[0]->set_bounding_object_indices(
6463 new_quads[1]->set_bounding_object_indices(
6468 new_quads[2]->set_bounding_object_indices(
6473 new_quads[3]->set_bounding_object_indices(
6478 for (
const auto &new_quad : new_quads)
6480 new_quad->set_used_flag();
6481 new_quad->clear_user_flag();
6483 new_quad->clear_children();
6484 new_quad->set_boundary_id_internal(quad->boundary_id());
6485 new_quad->set_manifold_id(quad->manifold_id());
6489 for (
unsigned int j = 0;
6490 j < GeometryInfo<dim>::lines_per_face;
6492 new_quad->set_line_orientation(j,
true);
6498 new_quads[0]->set_line_orientation(
6499 0, quad->line_orientation(0));
6500 new_quads[0]->set_line_orientation(
6501 2, quad->line_orientation(2));
6502 new_quads[1]->set_line_orientation(
6503 1, quad->line_orientation(1));
6504 new_quads[1]->set_line_orientation(
6505 2, quad->line_orientation(2));
6506 new_quads[2]->set_line_orientation(
6507 0, quad->line_orientation(0));
6508 new_quads[2]->set_line_orientation(
6509 3, quad->line_orientation(3));
6510 new_quads[3]->set_line_orientation(
6511 1, quad->line_orientation(1));
6512 new_quads[3]->set_line_orientation(
6513 3, quad->line_orientation(3));
6517 quad->clear_user_flag();
6528 cells_with_distorted_children;
6542 for (; hex != endh; ++hex)
6543 if (hex->refine_flag_set())
6551 hex->clear_refine_flag();
6552 hex->set_refinement_case(ref_case);
6563 unsigned int n_new_lines = 0;
6564 unsigned int n_new_quads = 0;
6565 unsigned int n_new_hexes = 0;
6596 new_lines(n_new_lines);
6597 for (
unsigned int i = 0; i < n_new_lines; ++i)
6600 triangulation.
faces->lines
6604 new_lines[i]->used() ==
false,