47 template <
int structdim>
59 template <
int structdim>
97 namespace TriangulationImplementation
171 template <
int dim,
int spacedim>
173 cell_is_patch_level_1(
178 unsigned int n_active_children = 0;
179 for (
unsigned int i = 0; i < cell->n_children(); ++i)
180 if (cell->child(i)->is_active())
183 return (n_active_children == 0) ||
184 (n_active_children == cell->n_children());
194 template <
int dim,
int spacedim>
196 cell_will_be_coarsened(
202 if (cell->has_children())
204 unsigned int children_to_coarsen = 0;
205 const unsigned int n_children = cell->n_children();
207 for (
unsigned int c = 0; c < n_children; ++c)
208 if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
209 ++children_to_coarsen;
210 if (children_to_coarsen == n_children)
213 for (
unsigned int c = 0; c < n_children; ++c)
214 if (cell->child(c)->is_active())
215 cell->child(c)->clear_coarsen_flag();
251 template <
int dim,
int spacedim>
253 face_will_be_refined_by_neighbor_internal(
255 const unsigned int face_no,
264 cell->neighbor(face_no);
271 if (neighbor->has_children())
276 if (cell_will_be_coarsened(neighbor))
285 expected_face_ref_case = cell->face(face_no)->refinement_case();
297 const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
305 neighbor->face_orientation(neighbor_neighbor),
306 neighbor->face_flip(neighbor_neighbor),
307 neighbor->face_rotation(neighbor_neighbor));
311 neighbor_face = neighbor->face(neighbor_neighbor);
312 const int this_face_index = cell->face_index(face_no);
318 if (neighbor_face->index() == this_face_index)
324 expected_face_ref_case = face_ref_case;
342 for (
unsigned int c = 0; c < neighbor_face->n_children(); ++c)
343 if (neighbor_face->child_index(c) == this_face_index)
352 if ((neighbor_face->refinement_case() | face_ref_case) ==
353 neighbor_face->refinement_case())
372 expected_face_ref_case =
373 ~neighbor_face->refinement_case();
400 template <
int dim,
int spacedim>
402 face_will_be_refined_by_neighbor(
404 const unsigned int face_no)
407 return face_will_be_refined_by_neighbor_internal(cell, face_no,
dummy);
414 template <
int dim,
int spacedim>
416 face_will_be_refined_by_neighbor(
418 const unsigned int face_no,
421 return face_will_be_refined_by_neighbor_internal(cell,
423 expected_face_ref_case);
428 template <
int dim,
int spacedim>
430 satisfies_level1_at_vertex_rule(
433 std::vector<unsigned int> min_adjacent_cell_level(
435 std::vector<unsigned int> max_adjacent_cell_level(
438 for (
const auto &cell :
triangulation.active_cell_iterators())
441 min_adjacent_cell_level[cell->vertex_index(v)] =
442 std::min<unsigned int>(
443 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
444 max_adjacent_cell_level[cell->vertex_index(v)] =
445 std::max<unsigned int>(
446 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
449 for (
unsigned int k = 0; k <
triangulation.n_vertices(); ++k)
451 if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
464 template <
int dim,
int spacedim>
465 std::vector<unsigned int>
470 std::vector<unsigned int> line_cell_count(
triangulation.n_raw_lines(),
475 for (; cell != endc; ++cell)
476 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
477 ++line_cell_count[cell->line_index(
l)];
478 return line_cell_count;
481 return std::vector<unsigned int>();
492 template <
int dim,
int spacedim>
493 std::vector<unsigned int>
498 std::vector<unsigned int> quad_cell_count(
triangulation.n_raw_quads(),
503 for (; cell != endc; ++cell)
505 ++quad_cell_count[cell->quad_index(q)];
506 return quad_cell_count;
509 return std::vector<unsigned int>();
533 void reorder_compatibility(std::vector<
CellData<2>> &cells,
536 for (
auto &cell : cells)
537 std::swap(cell.vertices[2], cell.vertices[3]);
541 void reorder_compatibility(std::vector<
CellData<3>> &cells,
545 for (
auto &cell : cells)
548 tmp[i] = cell.vertices[i];
554 std::vector<CellData<2>>::iterator boundary_quad =
556 std::vector<CellData<2>>::iterator end_quad =
558 for (
unsigned int quad_no = 0; boundary_quad != end_quad;
559 ++boundary_quad, ++quad_no)
560 std::swap(boundary_quad->vertices[2], boundary_quad->vertices[3]);
582 template <
int dim,
int spacedim>
587 if (line->has_children())
588 return line->child(0)->vertex_index(1);
593 template <
int dim,
int spacedim>
598 switch (
static_cast<unsigned char>(quad->refinement_case()))
601 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
604 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
607 return quad->child(0)->vertex_index(3);
616 template <
int dim,
int spacedim>
621 switch (
static_cast<unsigned char>(hex->refinement_case()))
624 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
627 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
630 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
633 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
636 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
639 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
642 return hex->child(0)->vertex_index(7);
663 template <
class TRIANGULATION>
664 inline typename TRIANGULATION::DistortedCellList
665 collect_distorted_coarse_cells(
const TRIANGULATION &)
667 return typename TRIANGULATION::DistortedCellList();
698 if (determinants[i] <= 1
e-9 * std::pow(cell->diameter(), 1. * dim))
700 distorted_cells.distorted_cells.push_back(cell);
705 return distorted_cells;
717 has_distorted_children(
719 std::integral_constant<int, dim>,
720 std::integral_constant<int, dim>)
724 for (
unsigned int c = 0; c < cell->n_children(); ++c)
728 vertices[i] = cell->child(c)->vertex(i);
734 if (determinants[i] <=
735 1
e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
750 template <
int dim,
int spacedim>
752 has_distorted_children(
754 std::integral_constant<int, dim>,
755 std::integral_constant<int, spacedim>)
766 template <
int spacedim>
771 template <
int dim,
int spacedim>
826 static const unsigned int left_right_offset[2][6][2] = {
835 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
846 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
853 for (; cell != endc; ++cell)
859 const unsigned int offset =
860 (cell->direction_flag() ?
861 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
862 1 - left_right_offset[dim - 2][f][cell->face_orientation(f)]);
864 adjacent_cells[2 * face->index() + offset] = cell;
874 if (cell->is_active() && face->has_children())
876 adjacent_cells[2 * face->child(0)->index() + offset] = cell;
877 adjacent_cells[2 * face->child(1)->index() + offset] = cell;
902 if (face->has_children() &&
903 (cell->is_active() ||
905 cell->refinement_case(), f) ==
908 for (
unsigned int c = 0; c < face->n_children(); ++c)
909 adjacent_cells[2 * face->child(c)->index() + offset] = cell;
910 if (face->child(0)->has_children())
912 adjacent_cells[2 * face->child(0)->child(0)->index() +
914 adjacent_cells[2 * face->child(0)->child(1)->index() +
917 if (face->child(1)->has_children())
919 adjacent_cells[2 * face->child(1)->child(0)->index() +
921 adjacent_cells[2 * face->child(1)->child(1)->index() +
936 const unsigned int offset =
937 (cell->direction_flag() ?
938 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
939 1 - left_right_offset[dim - 2][f][cell->face_orientation(f)]);
941 f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
946 template <
int dim,
int spacedim>
948 update_periodic_face_map_recursively(
951 unsigned int n_face_1,
952 unsigned int n_face_2,
953 const std::bitset<3> & orientation,
959 std::bitset<3>>> &periodic_face_map)
962 const FaceIterator face_1 = cell_1->face(n_face_1);
963 const FaceIterator face_2 = cell_2->face(n_face_2);
965 const bool face_orientation = orientation[0];
966 const bool face_flip = orientation[1];
967 const bool face_rotation = orientation[2];
969 Assert((dim != 1) || (face_orientation ==
true && face_flip ==
false &&
970 face_rotation ==
false),
972 "(face_orientation, face_flip, face_rotation) "
973 "is invalid for 1D"));
975 Assert((dim != 2) || (face_orientation ==
true && face_rotation ==
false),
977 "(face_orientation, face_flip, face_rotation) "
978 "is invalid for 2D"));
982 Assert(face_1->at_boundary() && face_2->at_boundary(),
983 ExcMessage(
"Periodic faces must be on the boundary"));
989 Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
994 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
996 const CellFace cell_face_1(cell_1, n_face_1);
997 const CellFace cell_face_2(cell_2, n_face_2);
998 const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
999 cell_face_2, orientation);
1001 const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
1002 periodic_faces(cell_face_1, cell_face_orientation_2);
1006 periodic_face_map.insert(periodic_faces);
1010 if (cell_1->has_children())
1012 if (cell_2->has_children())
1014 update_periodic_face_map_recursively<dim, spacedim>(
1015 cell_1->child(n_face_1),
1016 cell_2->child(n_face_2),
1024 update_periodic_face_map_recursively<dim, spacedim>(
1025 cell_1->child(n_face_1),
1040 static const int lookup_table_2d[2][2] =
1047 static const int lookup_table_3d[2][2][2][4] =
1066 if (cell_1->has_children())
1068 if (cell_2->has_children())
1074 Assert(face_1->n_children() ==
1080 for (
unsigned int i = 0;
1081 i < GeometryInfo<dim>::max_children_per_face;
1089 j = lookup_table_2d[face_flip][i];
1092 j = lookup_table_3d[face_orientation][face_flip]
1100 unsigned int child_cell_1 =
1102 cell_1->refinement_case(),
1105 cell_1->face_orientation(n_face_1),
1106 cell_1->face_flip(n_face_1),
1107 cell_1->face_rotation(n_face_1),
1108 face_1->refinement_case());
1109 unsigned int child_cell_2 =
1111 cell_2->refinement_case(),
1114 cell_2->face_orientation(n_face_2),
1115 cell_2->face_flip(n_face_2),
1116 cell_2->face_rotation(n_face_2),
1117 face_2->refinement_case());
1119 Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
1122 Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
1128 update_periodic_face_map_recursively<dim, spacedim>(
1129 cell_1->child(child_cell_1),
1130 cell_2->child(child_cell_2),
1139 for (
unsigned int i = 0;
1140 i < GeometryInfo<dim>::max_children_per_face;
1144 unsigned int child_cell_1 =
1146 cell_1->refinement_case(),
1149 cell_1->face_orientation(n_face_1),
1150 cell_1->face_flip(n_face_1),
1151 cell_1->face_rotation(n_face_1),
1152 face_1->refinement_case());
1155 update_periodic_face_map_recursively<dim, spacedim>(
1156 cell_1->child(child_cell_1),
1174 namespace TriangulationImplementation
1190 <<
"Something went wrong when making cell " << arg1
1191 <<
". Read the docs and the source code "
1192 <<
"for more information.");
1199 <<
"Something went wrong upon construction of cell "
1213 <<
" has negative measure. This typically "
1214 <<
"indicates some distortion in the cell, or a mistakenly "
1215 <<
"swapped pair of vertices in the input to "
1216 <<
"Triangulation::create_triangulation().");
1228 <<
"Error while creating cell " << arg1
1229 <<
": the vertex index " << arg2 <<
" must be between 0 and "
1238 <<
"While trying to assign a boundary indicator to a line: "
1239 <<
"the line with end vertices " << arg1 <<
" and " << arg2
1240 <<
" does not exist.");
1250 <<
"While trying to assign a boundary indicator to a quad: "
1251 <<
"the quad with bounding lines " << arg1 <<
", " << arg2
1252 <<
", " << arg3 <<
", " << arg4 <<
" does not exist.");
1262 <<
"The input data for creating a triangulation contained "
1263 <<
"information about a line with indices " << arg1 <<
" and " << arg2
1264 <<
" that is described to have boundary indicator "
1265 <<
static_cast<int>(arg3)
1266 <<
". However, this is an internal line not located on the "
1267 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1269 <<
"If this happened at a place where you call "
1270 <<
"Triangulation::create_triangulation() yourself, you need "
1271 <<
"to check the SubCellData object you pass to this function."
1274 <<
"If this happened in a place where you are reading a mesh "
1275 <<
"from a file, then you need to investigate why such a line "
1276 <<
"ended up in the input file. A typical case is a geometry "
1277 <<
"that consisted of multiple parts and for which the mesh "
1278 <<
"generator program assumes that the interface between "
1279 <<
"two parts is a boundary when that isn't supposed to be "
1280 <<
"the case, or where the mesh generator simply assigns "
1281 <<
"'geometry indicators' to lines at the perimeter of "
1282 <<
"a part that are not supposed to be interpreted as "
1283 <<
"'boundary indicators'.");
1295 <<
"The input data for creating a triangulation contained "
1296 <<
"information about a quad with indices " << arg1 <<
", " << arg2
1297 <<
", " << arg3 <<
", and " << arg4
1298 <<
" that is described to have boundary indicator "
1299 <<
static_cast<int>(arg5)
1300 <<
". However, this is an internal quad not located on the "
1301 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1303 <<
"If this happened at a place where you call "
1304 <<
"Triangulation::create_triangulation() yourself, you need "
1305 <<
"to check the SubCellData object you pass to this function."
1308 <<
"If this happened in a place where you are reading a mesh "
1309 <<
"from a file, then you need to investigate why such a quad "
1310 <<
"ended up in the input file. A typical case is a geometry "
1311 <<
"that consisted of multiple parts and for which the mesh "
1312 <<
"generator program assumes that the interface between "
1313 <<
"two parts is a boundary when that isn't supposed to be "
1314 <<
"the case, or where the mesh generator simply assigns "
1315 <<
"'geometry indicators' to quads at the surface of "
1316 <<
"a part that are not supposed to be interpreted as "
1317 <<
"'boundary indicators'.");
1326 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1327 <<
" and " << arg2 <<
" appears more than once. "
1328 <<
"This is not allowed.");
1338 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1339 <<
" and " << arg2 <<
" appears multiple times with different (valid) "
1340 << arg3 <<
". This is not allowed.");
1352 <<
"In SubCellData the quad info of the quad with line indices " << arg1
1353 <<
", " << arg2 <<
", " << arg3 <<
" and " << arg4
1354 <<
" appears multiple times with different (valid) " << arg5
1355 <<
". This is not allowed.");
1465 template <
int dim,
int spacedim>
1469 const unsigned int level_objects,
1472 using line_iterator =
1476 if (level_objects > 0)
1509 for (; line != endc; ++line)
1512 if (line->has_children() ==
false)
1530 for (; line != endc; ++line)
1533 if (line->has_children() ==
false)
1553 template <
int dim,
int spacedim>
1557 const unsigned int level_objects,
1568 &compute_number_cache<dim, spacedim>),
1574 using quad_iterator =
1591 unsigned int n_levels = 0;
1592 if (level_objects > 0)
1596 n_levels =
level + 1;
1609 (
level == n_levels - 1 ?
1612 for (; quad != endc; ++quad)
1615 if (quad->has_children() ==
false)
1633 for (; quad != endc; ++quad)
1636 if (quad->has_children() ==
false)
1642 update_lines.
join();
1660 template <
int dim,
int spacedim>
1664 const unsigned int level_objects,
1675 &compute_number_cache<dim, spacedim>),
1681 using hex_iterator =
1699 unsigned int n_levels = 0;
1700 if (level_objects > 0)
1704 n_levels =
level + 1;
1716 endc = (
level == n_levels - 1 ?
1719 for (; hex != endc; ++hex)
1722 if (hex->has_children() ==
false)
1740 for (; hex != endc; ++hex)
1743 if (hex->has_children() ==
false)
1749 update_quads_and_lines.
join();
1760 template <
int spacedim>
1776 const unsigned int dim = 1;
1780 triangulation.vertices_used = std::vector<bool>(v.size(),
true);
1788 if (dim == spacedim)
1790 for (
unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
1800 cells[cell_no].vertices);
1812 std::vector<std::vector<int>> lines_at_vertex(v.size());
1816 std_cxx14::make_unique<
1818 triangulation.levels[0]->reserve_space(cells.size(), dim, spacedim);
1819 triangulation.levels[0]->cells.reserve_space(0, cells.size());
1824 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
1826 while (next_free_line->used())
1829 next_free_line->set(
1832 next_free_line->set_used_flag();
1833 next_free_line->set_material_id(cells[cell].
material_id);
1834 next_free_line->set_manifold_id(cells[cell].
manifold_id);
1836 next_free_line->set_subdomain_id(0);
1840 lines_at_vertex[cells[cell].vertices[0]].push_back(cell);
1841 lines_at_vertex[cells[cell].vertices[1]].push_back(cell);
1847 unsigned int boundary_nodes = 0;
1848 for (
const auto &line : lines_at_vertex)
1849 switch (line.size())
1862 "You have a vertex in your triangulation "
1863 "at which more than two cells come together. "
1864 "(For one dimensional triangulation, cells are "
1867 "This is not currently supported because the "
1868 "Triangulation class makes the assumption that "
1869 "every cell has zero or one neighbors behind "
1870 "each face (here, behind each vertex), but in your "
1871 "situation there would be more than one."
1873 "Support for this is not currently implemented. "
1874 "If you need to work with triangulations where "
1875 "more than two cells come together at a vertex, "
1876 "duplicate the vertices once per cell (i.e., put "
1877 "multiple vertices at the same physical location, "
1878 "but using different vertex indices for each) "
1879 "and then ensure continuity of the solution by "
1880 "explicitly creating constraints that the degrees "
1881 "of freedom at these vertices have the same "
1882 "value, using the AffineConstraints class."));
1900 if (lines_at_vertex[line->vertex_index(vertex)][0] == line->index())
1901 if (lines_at_vertex[line->vertex_index(vertex)].size() == 2)
1906 lines_at_vertex[line->vertex_index(vertex)][1]);
1907 line->set_neighbor(vertex, neighbor);
1922 lines_at_vertex[line->vertex_index(vertex)][0]);
1923 line->set_neighbor(vertex, neighbor);
1932 for (
const auto &cell :
triangulation.active_cell_iterators())
1938 if (cell->at_boundary(f))
1940 [cell->face(f)->vertex_index()] = f;
1952 template <
int spacedim>
1962 const unsigned int dim = 2;
1966 triangulation.vertices_used = std::vector<bool>(v.size(),
true);
1974 if (dim == spacedim)
1976 for (
unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
1983 cells[cell_no].vertices);
2004 std::map<std::pair<int, int>,
2007 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
2009 for (
const auto vertex : cells[cell].
vertices)
2020 std::pair<int, int> line_vertices(
2043 line_vertices.second, line_vertices.first)) ==
2059 std::vector<unsigned short int> vertex_touch_count(v.size(), 0);
2061 std::pair<int, int>,
2063 for (i = needed_lines.begin(); i != needed_lines.end(); ++i)
2067 ++vertex_touch_count[i->first.first];
2068 ++vertex_touch_count[i->first.second];
2075 AssertThrow(*(std::min_element(vertex_touch_count.begin(),
2076 vertex_touch_count.end())) >= 2,
2078 "During creation of a triangulation, a part of the "
2079 "algorithm encountered a vertex that is part of only "
2080 "a single adjacent line. However, in 2d, every vertex "
2081 "needs to be at least part of two lines."));
2086 std_cxx14::make_unique<
2090 triangulation.levels[0]->reserve_space(cells.size(), dim, spacedim);
2091 triangulation.faces->lines.reserve_space(0, needed_lines.size());
2092 triangulation.levels[0]->cells.reserve_space(0, cells.size());
2099 std::pair<int, int>,
2101 for (i = needed_lines.begin(); line !=
triangulation.end_line();
2105 i->first.first, i->first.second));
2106 line->set_used_flag();
2107 line->clear_user_flag();
2118 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>>
2125 for (
unsigned int c = 0; c < cells.size(); ++c, ++cell)
2129 for (
unsigned int line = 0;
2130 line < GeometryInfo<dim>::lines_per_cell;
2132 lines[line] = needed_lines[std::make_pair(
2142 lines[3]->index()));
2144 cell->set_used_flag();
2147 cell->clear_user_data();
2148 cell->set_subdomain_id(0);
2153 for (
const auto &line : lines)
2154 adjacent_cells[line->index()].push_back(cell);
2164 const unsigned int n_adj_cells =
2165 adjacent_cells[line->index()].size();
2173 AssertThrow((n_adj_cells >= 1) && (n_adj_cells <= 2),
2179 (n_adj_cells >= 1) && (n_adj_cells <= 2),
2180 ExcMessage(
"You have a line in your triangulation at which "
2181 "more than two cells come together."
2183 "This is not currently supported because the "
2184 "Triangulation class makes the assumption that "
2185 "every cell has zero or one neighbors behind each "
2186 "face (here, behind each line), but in your "
2187 "situation there would be more than one."
2189 "Support for this is not currently implemented. "
2190 "If you need to work with triangulations where "
2191 "more than two cells come together at a line, "
2192 "duplicate the vertices once per cell (i.e., put "
2193 "multiple vertices at the same physical location, "
2194 "but using different vertex indices for each) "
2195 "and then ensure continuity of the solution by "
2196 "explicitly creating constraints that the degrees "
2197 "of freedom at these lines have the same "
2198 "value, using the AffineConstraints class."));
2203 line->set_boundary_id_internal(
2212 std::pair<int, int> line_vertices(
2213 std::make_pair(subcell_line.vertices[0],
2214 subcell_line.vertices[1]));
2215 if (needed_lines.find(line_vertices) != needed_lines.end())
2217 line = needed_lines[line_vertices];
2221 std::swap(line_vertices.first, line_vertices.second);
2222 if (needed_lines.find(line_vertices) != needed_lines.end())
2223 line = needed_lines[line_vertices];
2228 line_vertices.second));
2233 line->boundary_id() !=
2236 line_vertices.second));
2241 line->manifold_id() == subcell_line.manifold_id,
2243 line_vertices.second,
2245 line->set_manifold_id(subcell_line.manifold_id);
2254 line->vertex_index(1),
2255 subcell_line.boundary_id));
2256 line->set_boundary_id_internal(subcell_line.boundary_id);
2266 for (
unsigned int side = 0; side < 4; ++side)
2267 if (adjacent_cells[cell->line(side)->index()][0] == cell)
2271 if (adjacent_cells[cell->line(side)->index()].size() == 2)
2275 side, adjacent_cells[cell->line(side)->index()][1]);
2281 cell->set_neighbor(side,
2282 adjacent_cells[cell->line(side)->index()][0]);
2330 template <
int spacedim>
2340 const unsigned int dim = 3;
2344 triangulation.vertices_used = std::vector<bool>(v.size(),
true);
2350 for (
unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2357 cells[cell_no].vertices);
2382 typename std::map<std::pair<int, int>,
2385 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
2389 for (
const auto vertex : cells[cell].
vertices)
2395 for (
unsigned int line = 0;
2396 line < GeometryInfo<dim>::lines_per_cell;
2405 std::pair<int, int> line_vertices(
2414 if ((needed_lines.find(std::make_pair(line_vertices.second,
2415 line_vertices.first)) ==
2416 needed_lines.end()))
2434 std::vector<unsigned short int> vertex_touch_count(v.size(), 0);
2436 std::pair<int, int>,
2438 for (i = needed_lines.begin(); i != needed_lines.end(); ++i)
2442 ++vertex_touch_count[i->first.first];
2443 ++vertex_touch_count[i->first.second];
2451 *(std::min_element(vertex_touch_count.begin(),
2452 vertex_touch_count.end())) >= 3,
2454 "During creation of a triangulation, a part of the "
2455 "algorithm encountered a vertex that is part of only "
2456 "one or two adjacent lines. However, in 3d, every vertex "
2457 "needs to be at least part of three lines."));
2466 std_cxx14::make_unique<
2470 triangulation.levels[0]->reserve_space(cells.size(), dim, spacedim);
2471 triangulation.faces->lines.reserve_space(0, needed_lines.size());
2478 std::pair<int, int>,
2480 for (i = needed_lines.begin(); line !=
triangulation.end_line();
2484 i->first.first, i->first.second));
2485 line->set_used_flag();
2486 line->clear_user_flag();
2510 std::map<internal::TriangulationImplementation::TriaObject<2>,
2511 std::pair<typename Triangulation<dim, spacedim>::quad_iterator,
2512 std::array<bool, GeometryInfo<dim>::lines_per_face>>,
2515 for (
const auto &cell : cells)
2543 std::array<bool, GeometryInfo<dim>::lines_per_face> orientation;
2545 for (
unsigned int line = 0;
2546 line < GeometryInfo<dim>::lines_per_cell;
2549 line_list[line] = std::pair<int, int>(
2554 inverse_line_list[line] = std::pair<int, int>(
2571 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_face;
2573 if (needed_lines.find(
2575 face,
l)]) == needed_lines.end())
2579 dim>::face_to_cell_lines(face,
l)]]
2581 orientation[
l] =
true;
2587 dim>::face_to_cell_lines(face,
l)]]
2589 orientation[
l] =
false;
2638 test_quad_1(quad.
face(2),
2644 test_quad_2(quad.
face(0),
2650 test_quad_3(quad.
face(3),
2656 test_quad_4(quad.
face(1),
2662 test_quad_5(quad.
face(2),
2668 test_quad_6(quad.
face(1),
2674 test_quad_7(quad.
face(3),
2680 if (needed_quads.find(test_quad_1) == needed_quads.end() &&
2681 needed_quads.find(test_quad_2) == needed_quads.end() &&
2682 needed_quads.find(test_quad_3) == needed_quads.end() &&
2683 needed_quads.find(test_quad_4) == needed_quads.end() &&
2684 needed_quads.find(test_quad_5) == needed_quads.end() &&
2685 needed_quads.find(test_quad_6) == needed_quads.end() &&
2686 needed_quads.find(test_quad_7) == needed_quads.end())
2687 needed_quads[quad] =
2698 triangulation.faces->quads.reserve_space(0, needed_quads.size());
2705 std::pair<typename Triangulation<dim, spacedim>::quad_iterator,
2706 std::array<bool, GeometryInfo<dim>::lines_per_face>>,
2708 for (q = needed_quads.begin(); quad !=
triangulation.end_quad();
2711 quad->set(q->first);
2712 quad->set_used_flag();
2713 quad->clear_user_flag();
2716 quad->set_line_orientation(0, q->second.second[0]);
2717 quad->set_line_orientation(1, q->second.second[1]);
2718 quad->set_line_orientation(2, q->second.second[2]);
2719 quad->set_line_orientation(3, q->second.second[3]);
2724 q->second.first = quad;
2736 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>>
2743 for (
unsigned int c = 0; c < cells.size(); ++c, ++cell)
2758 unsigned int face_line_list[4];
2759 for (
unsigned int line = 0;
2760 line < GeometryInfo<dim>::lines_per_cell;
2763 line_list[line] = std::make_pair(
2768 inverse_line_list[line] = std::pair<int, int>(
2787 for (
unsigned int l = 0;
2788 l < GeometryInfo<dim>::lines_per_face;
2791 dim>::face_to_cell_lines(face,
l)]) ==
2795 dim>::face_to_cell_lines(face,
l)]]
2800 dim>::face_to_cell_lines(face,
l)]]
2809 if (needed_quads.find(quad) != needed_quads.end())
2820 face_iterator[face] = needed_quads[quad].first;
2821 face_orientation[face] =
true;
2822 face_flip[face] =
false;
2823 face_rotation[face] =
false;
2851 test_quad_4(quad.
face(1),
2869 test_quad_7(quad.
face(3),
2875 if (needed_quads.find(test_quad_1) != needed_quads.end())
2877 face_iterator[face] = needed_quads[test_quad_1].first;
2878 face_orientation[face] =
false;
2879 face_flip[face] =
false;
2880 face_rotation[face] =
false;
2882 else if (needed_quads.find(test_quad_2) !=
2885 face_iterator[face] = needed_quads[test_quad_2].first;
2886 face_orientation[face] =
false;
2887 face_flip[face] =
false;
2888 face_rotation[face] =
true;
2890 else if (needed_quads.find(test_quad_3) !=
2893 face_iterator[face] = needed_quads[test_quad_3].first;
2894 face_orientation[face] =
false;
2895 face_flip[face] =
true;
2896 face_rotation[face] =
false;
2898 else if (needed_quads.find(test_quad_4) !=
2901 face_iterator[face] = needed_quads[test_quad_4].first;
2902 face_orientation[face] =
false;
2903 face_flip[face] =
true;
2904 face_rotation[face] =
true;
2906 else if (needed_quads.find(test_quad_5) !=
2909 face_iterator[face] = needed_quads[test_quad_5].first;
2910 face_orientation[face] =
true;
2911 face_flip[face] =
false;
2912 face_rotation[face] =
true;
2914 else if (needed_quads.find(test_quad_6) !=
2917 face_iterator[face] = needed_quads[test_quad_6].first;
2918 face_orientation[face] =
true;
2919 face_flip[face] =
true;
2920 face_rotation[face] =
false;
2922 else if (needed_quads.find(test_quad_7) !=
2925 face_iterator[face] = needed_quads[test_quad_7].first;
2926 face_orientation[face] =
true;
2927 face_flip[face] =
true;
2928 face_rotation[face] =
true;
2943 face_iterator[0]->index(),
2944 face_iterator[1]->index(),
2945 face_iterator[2]->index(),
2946 face_iterator[3]->index(),
2947 face_iterator[4]->index(),
2948 face_iterator[5]->index()));
2950 cell->set_used_flag();
2953 cell->clear_user_flag();
2954 cell->clear_user_data();
2955 cell->set_subdomain_id(0);
2961 cell->set_face_orientation(quad, face_orientation[quad]);
2962 cell->set_face_flip(quad, face_flip[quad]);
2963 cell->set_face_rotation(quad, face_rotation[quad]);
2970 for (
const auto &quad : face_iterator)
2971 adjacent_cells[quad->index()].push_back(cell);
2988 std::multimap<unsigned int, std::pair<unsigned int, unsigned int>>
2991 for (
unsigned int line = 0;
2992 line < GeometryInfo<dim>::lines_per_face;
2994 cell_to_face_lines.insert(
2995 std::pair<
unsigned int,
2996 std::pair<unsigned int, unsigned int>>(
2998 std::pair<unsigned int, unsigned int>(face, line)));
2999 std::multimap<
unsigned int,
3000 std::pair<unsigned int, unsigned int>>::
3001 const_iterator map_iter = cell_to_face_lines.begin();
3003 for (; map_iter != cell_to_face_lines.end(); ++map_iter)
3005 const unsigned int cell_line = map_iter->first;
3006 const unsigned int face1 = map_iter->second.first;
3007 const unsigned int line1 = map_iter->second.second;
3009 Assert(map_iter != cell_to_face_lines.end(),
3011 Assert(map_iter->first == cell_line,
3013 const unsigned int face2 = map_iter->second.first;
3014 const unsigned int line2 = map_iter->second.second;
3021 Assert(face_iterator[face1]->line(
3024 face_orientation[face1],
3026 face_rotation[face1])) ==
3027 face_iterator[face2]->line(
3030 face_orientation[face2],
3032 face_rotation[face2])),
3048 const unsigned int n_adj_cells =
3049 adjacent_cells[quad->index()].size();
3052 AssertThrow((n_adj_cells >= 1) && (n_adj_cells <= 2),
3057 quad->set_boundary_id_internal(
3102 if (quad->at_boundary())
3104 for (
unsigned int l = 0;
l < 4; ++
l)
3108 line->set_boundary_id_internal(0);
3120 std::pair<int, int> line_vertices(
3121 std::make_pair(subcell_line.vertices[0],
3122 subcell_line.vertices[1]));
3123 if (needed_lines.find(line_vertices) != needed_lines.end())
3126 line = needed_lines[line_vertices];
3132 std::swap(line_vertices.first, line_vertices.second);
3133 if (needed_lines.find(line_vertices) != needed_lines.end())
3134 line = needed_lines[line_vertices];
3139 line_vertices.second));
3142 if (line->at_boundary())
3147 line->boundary_id() == subcell_line.boundary_id,
3149 line_vertices.second,
3155 if (subcell_line.boundary_id !=
3157 line->set_boundary_id(subcell_line.boundary_id);
3161 line->manifold_id() == subcell_line.manifold_id,
3163 line_vertices.second,
3165 line->set_manifold_id(subcell_line.manifold_id);
3182 for (
unsigned int i = 0; i < 4; ++i)
3184 std::pair<int, int> line_vertices(
3194 if (needed_lines.find(line_vertices) != needed_lines.end())
3195 line[i] = needed_lines[line_vertices];
3200 std::swap(line_vertices.first, line_vertices.second);
3201 if (needed_lines.find(line_vertices) != needed_lines.end())
3202 line[i] = needed_lines[line_vertices];
3208 line_vertices.second));
3247 static const unsigned int lex2cclock[4] = {3, 1, 0, 2};
3254 line_counterclock[4];
3255 for (
unsigned int i = 0; i < 4; ++i)
3256 line_counterclock[lex2cclock[i]] = line[i];
3257 unsigned int n_rotations = 0;
3258 bool not_found_quad_1;
3259 while ((not_found_quad_1 = (needed_quads.find(quad_compare_1) ==
3260 needed_quads.end())) &&
3261 (needed_quads.find(quad_compare_2) == needed_quads.end()) &&
3267 line_counterclock + 1,
3268 line_counterclock + 4);
3272 for (
unsigned int i = 0; i < 4; ++i)
3275 i, line_counterclock[lex2cclock[i]]->index());
3277 (i + 2) % 4, line_counterclock[lex2cclock[i]]->index());
3289 if (not_found_quad_1)
3290 quad = needed_quads[quad_compare_2].first;
3292 quad = needed_quads[quad_compare_1].first;
3296 if (quad->at_boundary())
3301 quad->boundary_id() == subcell_quad.boundary_id,
3311 if (subcell_quad.boundary_id !=
3313 quad->set_boundary_id(subcell_quad.boundary_id);
3317 AssertThrow(quad->manifold_id() == subcell_quad.manifold_id,
3324 quad->set_manifold_id(subcell_quad.manifold_id);
3334 for (
unsigned int face = 0; face < 6; ++face)
3335 if (adjacent_cells[cell->quad(face)->index()][0] == cell)
3339 if (adjacent_cells[cell->quad(face)->index()].size() == 2)
3343 face, adjacent_cells[cell->quad(face)->index()][1]);
3349 cell->set_neighbor(face,
3350 adjacent_cells[cell->quad(face)->index()][0]);
3369 template <
int spacedim>
3373 std::vector<unsigned int> &,
3374 std::vector<unsigned int> &)
3376 const unsigned int dim = 1;
3388 Assert(!cell->child(0)->has_children() &&
3389 !cell->child(1)->has_children(),
3395 if (cell->neighbor(0)->has_children())
3402 neighbor = neighbor->child(1);
3405 Assert(neighbor->neighbor(1) == cell->child(0),
3407 neighbor->set_neighbor(1, cell);
3413 if (neighbor->has_children())
3414 neighbor = neighbor->child(1);
3423 if (cell->neighbor(1)->has_children())
3430 neighbor = neighbor->child(0);
3433 Assert(neighbor->neighbor(0) == cell->child(1),
3435 neighbor->set_neighbor(0, cell);
3441 if (neighbor->has_children())
3442 neighbor = neighbor->child(0);
3452 triangulation.vertices_used[cell->child(0)->vertex_index(1)] =
false;
3458 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3461 cell->child(child)->clear_user_flag();
3462 cell->child(child)->clear_used_flag();
3467 cell->clear_children();
3468 cell->clear_user_flag();
3473 template <
int spacedim>
3477 std::vector<unsigned int> & line_cell_count,
3478 std::vector<unsigned int> &)
3480 const unsigned int dim = 2;
3488 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3491 lines_to_delete.reserve(4 * 2 + 4);
3496 for (
unsigned int c = 0; c < cell->n_children(); ++c)
3500 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
3501 --line_cell_count[child->line_index(
l)];
3516 .vertices_used[cell->child(0)->line(1)->vertex_index(1)] =
false;
3518 lines_to_delete.push_back(cell->child(0)->line(1));
3519 lines_to_delete.push_back(cell->child(0)->line(3));
3520 lines_to_delete.push_back(cell->child(3)->line(0));
3521 lines_to_delete.push_back(cell->child(3)->line(2));
3525 unsigned int inner_face_no =
3530 lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
3534 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3536 cell->child(child)->clear_user_data();
3537 cell->child(child)->clear_user_flag();
3538 cell->child(child)->clear_used_flag();
3543 cell->clear_children();
3544 cell->clear_refinement_case();
3545 cell->clear_user_flag();
3551 for (
unsigned int line_no = 0;
3552 line_no < GeometryInfo<dim>::lines_per_cell;
3556 cell->line(line_no);
3558 if (line->has_children())
3563 Assert((line_cell_count[line->child_index(0)] == 0 &&
3564 line_cell_count[line->child_index(1)] == 0) ||
3565 (line_cell_count[line->child_index(0)] > 0 &&
3566 line_cell_count[line->child_index(1)] > 0),
3569 if (line_cell_count[line->child_index(0)] == 0)
3571 for (
unsigned int c = 0; c < 2; ++c)
3572 Assert(!line->child(c)->has_children(),
3580 .vertices_used[line->child(0)->vertex_index(1)] =
false;
3582 lines_to_delete.push_back(line->child(0));
3583 lines_to_delete.push_back(line->child(1));
3585 line->clear_children();
3597 typename std::vector<
3599 line = lines_to_delete.begin(),
3600 endline = lines_to_delete.end();
3601 for (; line != endline; ++line)
3603 (*line)->clear_user_data();
3604 (*line)->clear_user_flag();
3605 (*line)->clear_used_flag();
3611 template <
int spacedim>
3615 std::vector<unsigned int> & line_cell_count,
3616 std::vector<unsigned int> & quad_cell_count)
3618 const unsigned int dim = 3;
3630 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3632 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3635 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3636 quads_to_delete.reserve(6 * 4 + 12);
3640 for (
unsigned int c = 0; c < cell->n_children(); ++c)
3644 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
3645 --line_cell_count[child->line_index(
l)];
3647 --quad_cell_count[child->quad_index(f)];
3661 quads_to_delete.push_back(cell->child(0)->face(1));
3664 quads_to_delete.push_back(cell->child(0)->face(3));
3667 quads_to_delete.push_back(cell->child(0)->face(5));
3670 quads_to_delete.push_back(cell->child(0)->face(1));
3671 quads_to_delete.push_back(cell->child(0)->face(3));
3672 quads_to_delete.push_back(cell->child(3)->face(0));
3673 quads_to_delete.push_back(cell->child(3)->face(2));
3675 lines_to_delete.push_back(cell->child(0)->line(11));
3678 quads_to_delete.push_back(cell->child(0)->face(1));
3679 quads_to_delete.push_back(cell->child(0)->face(5));
3680 quads_to_delete.push_back(cell->child(3)->face(0));
3681 quads_to_delete.push_back(cell->child(3)->face(4));
3683 lines_to_delete.push_back(cell->child(0)->line(5));
3686 quads_to_delete.push_back(cell->child(0)->face(3));
3687 quads_to_delete.push_back(cell->child(0)->face(5));
3688 quads_to_delete.push_back(cell->child(3)->face(2));
3689 quads_to_delete.push_back(cell->child(3)->face(4));
3691 lines_to_delete.push_back(cell->child(0)->line(7));
3694 quads_to_delete.push_back(cell->child(0)->face(1));
3695 quads_to_delete.push_back(cell->child(2)->face(1));
3696 quads_to_delete.push_back(cell->child(4)->face(1));
3697 quads_to_delete.push_back(cell->child(6)->face(1));
3699 quads_to_delete.push_back(cell->child(0)->face(3));
3700 quads_to_delete.push_back(cell->child(1)->face(3));
3701 quads_to_delete.push_back(cell->child(4)->face(3));
3702 quads_to_delete.push_back(cell->child(5)->face(3));
3704 quads_to_delete.push_back(cell->child(0)->face(5));
3705 quads_to_delete.push_back(cell->child(1)->face(5));
3706 quads_to_delete.push_back(cell->child(2)->face(5));
3707 quads_to_delete.push_back(cell->child(3)->face(5));
3709 lines_to_delete.push_back(cell->child(0)->line(5));
3710 lines_to_delete.push_back(cell->child(0)->line(7));
3711 lines_to_delete.push_back(cell->child(0)->line(11));
3712 lines_to_delete.push_back(cell->child(7)->line(0));
3713 lines_to_delete.push_back(cell->child(7)->line(2));
3714 lines_to_delete.push_back(cell->child(7)->line(8));
3720 triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3732 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3734 cell->child(child)->clear_user_data();
3735 cell->child(child)->clear_user_flag();
3742 cell->child(child)->set_face_orientation(f,
true);
3743 cell->child(child)->set_face_flip(f,
false);
3744 cell->child(child)->set_face_rotation(f,
false);
3747 cell->child(child)->clear_used_flag();
3752 cell->clear_children();
3753 cell->clear_refinement_case();
3754 cell->clear_user_flag();
3767 cell->face(quad_no);
3771 quad->has_children()) ||
3776 switch (quad->refinement_case())
3788 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3789 quad_cell_count[quad->child_index(1)] == 0) ||
3790 (quad_cell_count[quad->child_index(0)] > 0 &&
3791 quad_cell_count[quad->child_index(1)] > 0),
3797 unsigned int deleted_grandchildren = 0;
3798 unsigned int number_of_child_refinements = 0;
3800 for (
unsigned int c = 0; c < 2; ++c)
3801 if (quad->child(c)->has_children())
3803 ++number_of_child_refinements;
3808 (quad_cell_count[quad->child(c)->child_index(0)] ==
3810 quad_cell_count[quad->child(c)->child_index(1)] ==
3812 (quad_cell_count[quad->child(c)->child_index(0)] >
3814 quad_cell_count[quad->child(c)->child_index(1)] >
3817 if (quad_cell_count[quad->child(c)->child_index(0)] ==
3824 Assert(quad->refinement_case() +
3825 quad->child(c)->refinement_case() ==
3833 quads_to_delete.push_back(
3834 quad->child(c)->child(0));
3835 quads_to_delete.push_back(
3836 quad->child(c)->child(1));
3837 if (quad->child(c)->refinement_case() ==
3839 lines_to_delete.push_back(
3840 quad->child(c)->child(0)->line(1));
3842 lines_to_delete.push_back(
3843 quad->child(c)->child(0)->line(3));
3844 quad->child(c)->clear_children();
3845 quad->child(c)->clear_refinement_case();
3846 ++deleted_grandchildren;
3854 if (number_of_child_refinements > 0 &&
3855 deleted_grandchildren == number_of_child_refinements)
3860 middle_line = quad->child(0)->line(1);
3862 middle_line = quad->child(0)->line(3);
3864 lines_to_delete.push_back(middle_line->child(0));
3865 lines_to_delete.push_back(middle_line->child(1));
3867 .vertices_used[middle_vertex_index<dim, spacedim>(
3868 middle_line)] =
false;
3869 middle_line->clear_children();
3874 if (quad_cell_count[quad->child_index(0)] == 0)
3880 quads_to_delete.push_back(quad->child(0));
3881 quads_to_delete.push_back(quad->child(1));
3883 lines_to_delete.push_back(quad->child(0)->line(1));
3885 lines_to_delete.push_back(quad->child(0)->line(3));
3905 if (quad->child(0)->has_children())
3906 if (quad->refinement_case() ==
3940 quad_iterator switch_1 =
3941 quad->child(0)->child(1),
3943 quad->child(1)->child(0);
3945 Assert(!switch_1->has_children(),
3947 Assert(!switch_2->has_children(),
3950 const int switch_1_index = switch_1->index();
3951 const int switch_2_index = switch_2->index();
3952 for (
unsigned int l = 0;
3955 for (
unsigned int h = 0;
3957 ->cells.cells.size();
3959 for (
const unsigned int q :
3965 if (index == switch_1_index)
3968 .set_face(q, switch_2_index);
3969 else if (index == switch_2_index)
3972 .set_face(q, switch_1_index);
3977 const int switch_1_lines[4] = {
3978 static_cast<signed int>(
3979 switch_1->line_index(0)),
3980 static_cast<signed int>(
3981 switch_1->line_index(1)),
3982 static_cast<signed int>(
3983 switch_1->line_index(2)),
3984 static_cast<signed int>(
3985 switch_1->line_index(3))};
3986 const bool switch_1_line_orientations[4] = {
3987 switch_1->line_orientation(0),
3988 switch_1->line_orientation(1),
3989 switch_1->line_orientation(2),
3990 switch_1->line_orientation(3)};
3992 switch_1->boundary_id();
3993 const unsigned int switch_1_user_index =
3994 switch_1->user_index();
3995 const bool switch_1_user_flag =
3996 switch_1->user_flag_set();
4001 switch_2->line_index(1),
4002 switch_2->line_index(2),
4003 switch_2->line_index(3)));
4004 switch_1->set_line_orientation(
4005 0, switch_2->line_orientation(0));
4006 switch_1->set_line_orientation(
4007 1, switch_2->line_orientation(1));
4008 switch_1->set_line_orientation(
4009 2, switch_2->line_orientation(2));
4010 switch_1->set_line_orientation(
4011 3, switch_2->line_orientation(3));
4012 switch_1->set_boundary_id_internal(
4013 switch_2->boundary_id());
4014 switch_1->set_manifold_id(
4015 switch_2->manifold_id());
4016 switch_1->set_user_index(switch_2->user_index());
4017 if (switch_2->user_flag_set())
4018 switch_1->set_user_flag();
4020 switch_1->clear_user_flag();
4027 switch_1_lines[3]));
4028 switch_2->set_line_orientation(
4029 0, switch_1_line_orientations[0]);
4030 switch_2->set_line_orientation(
4031 1, switch_1_line_orientations[1]);
4032 switch_2->set_line_orientation(
4033 2, switch_1_line_orientations[2]);
4034 switch_2->set_line_orientation(
4035 3, switch_1_line_orientations[3]);
4036 switch_2->set_boundary_id_internal(
4037 switch_1_boundary_id);
4038 switch_2->set_manifold_id(
4039 switch_1->manifold_id());
4040 switch_2->set_user_index(switch_1_user_index);
4041 if (switch_1_user_flag)
4042 switch_2->set_user_flag();
4044 switch_2->clear_user_flag();
4046 const unsigned int child_0 =
4047 quad->child(0)->child_index(0);
4048 const unsigned int child_2 =
4049 quad->child(1)->child_index(0);
4050 quad->clear_children();
4051 quad->clear_refinement_case();
4052 quad->set_refinement_case(
4054 quad->set_children(0, child_0);
4055 quad->set_children(2, child_2);
4057 quad_cell_count[child_2]);
4072 const unsigned int child_0 =
4073 quad->child(0)->child_index(0);
4074 const unsigned int child_2 =
4075 quad->child(1)->child_index(0);
4076 quad->clear_children();
4077 quad->clear_refinement_case();
4078 quad->set_refinement_case(
4080 quad->set_children(0, child_0);
4081 quad->set_children(2, child_2);
4085 quad->clear_children();
4086 quad->clear_refinement_case();
4097 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
4098 quad_cell_count[quad->child_index(1)] == 0 &&
4099 quad_cell_count[quad->child_index(2)] == 0 &&
4100 quad_cell_count[quad->child_index(3)] == 0) ||
4101 (quad_cell_count[quad->child_index(0)] > 0 &&
4102 quad_cell_count[quad->child_index(1)] > 0 &&
4103 quad_cell_count[quad->child_index(2)] > 0 &&
4104 quad_cell_count[quad->child_index(3)] > 0),
4107 if (quad_cell_count[quad->child_index(0)] == 0)
4113 lines_to_delete.push_back(quad->child(0)->line(1));
4114 lines_to_delete.push_back(quad->child(3)->line(0));
4115 lines_to_delete.push_back(quad->child(0)->line(3));
4116 lines_to_delete.push_back(quad->child(3)->line(2));
4118 for (
unsigned int child = 0; child < quad->n_children();
4120 quads_to_delete.push_back(quad->child(child));
4123 .vertices_used[quad->child(0)->vertex_index(3)] =
4126 quad->clear_children();
4127 quad->clear_refinement_case();
4146 for (
unsigned int line_no = 0;
4147 line_no < GeometryInfo<dim>::lines_per_cell;
4151 cell->line(line_no);
4155 line->has_children()) ||
4160 if (line->has_children())
4165 Assert((line_cell_count[line->child_index(0)] == 0 &&
4166 line_cell_count[line->child_index(1)] == 0) ||
4167 (line_cell_count[line->child_index(0)] > 0 &&
4168 line_cell_count[line->child_index(1)] > 0),
4171 if (line_cell_count[line->child_index(0)] == 0)
4173 for (
unsigned int c = 0; c < 2; ++c)
4174 Assert(!line->child(c)->has_children(),
4182 .vertices_used[line->child(0)->vertex_index(1)] =
false;
4184 lines_to_delete.push_back(line->child(0));
4185 lines_to_delete.push_back(line->child(1));
4187 line->clear_children();
4199 typename std::vector<
4201 line = lines_to_delete.begin(),
4202 endline = lines_to_delete.end();
4203 for (; line != endline; ++line)
4205 (*line)->clear_user_data();
4206 (*line)->clear_user_flag();
4207 (*line)->clear_used_flag();
4210 typename std::vector<
4212 quad = quads_to_delete.begin(),
4213 endquad = quads_to_delete.end();
4214 for (; quad != endquad; ++quad)
4216 (*quad)->clear_user_data();
4217 (*quad)->clear_children();
4218 (*quad)->clear_refinement_case();
4219 (*quad)->clear_user_flag();
4220 (*quad)->clear_used_flag();
4242 template <
int spacedim>
4245 unsigned int & next_unused_vertex,
4252 const unsigned int dim = 2;
4255 cell->clear_refine_flag();
4352 int new_vertices[9];
4353 for (
unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
4354 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4355 for (
unsigned int line_no = 0; line_no < 4; ++line_no)
4356 if (cell->line(line_no)->has_children())
4357 new_vertices[4 + line_no] =
4358 cell->line(line_no)->child(0)->vertex_index(1);
4367 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
4368 ++next_unused_vertex;
4372 "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."));
4375 new_vertices[8] = next_unused_vertex;
4387 if (dim == spacedim)
4390 triangulation.vertices[next_unused_vertex] = cell->center(
true);
4397 if (cell->user_flag_set())
4400 cell->clear_user_flag();
4402 cell->center(
true,
true);
4412 cell->clear_user_flag();
4418 cell->center(
true,
true);
4425 unsigned int lmin = 8;
4426 unsigned int lmax = 12;
4433 for (
unsigned int l = lmin;
l < lmax; ++
l)
4435 while (next_unused_line->used() ==
true)
4437 new_lines[
l] = next_unused_line;
4441 new_lines[
l]->used() ==
false,
4443 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4458 for (
unsigned int c = 0; c < 2; ++c, ++
l)
4459 new_lines[
l] = cell->line(face_no)->child(c);
4464 new_vertices[6], new_vertices[8]));
4467 new_vertices[8], new_vertices[7]));
4470 new_vertices[4], new_vertices[8]));
4473 new_vertices[8], new_vertices[5]));
4482 new_lines[0] = cell->line(0);
4483 new_lines[1] = cell->line(1);
4484 new_lines[2] = cell->line(2)->child(0);
4485 new_lines[3] = cell->line(2)->child(1);
4486 new_lines[4] = cell->line(3)->child(0);
4487 new_lines[5] = cell->line(3)->child(1);
4490 new_vertices[6], new_vertices[7]));
4500 new_lines[0] = cell->line(0)->child(0);
4501 new_lines[1] = cell->line(0)->child(1);
4502 new_lines[2] = cell->line(1)->child(0);
4503 new_lines[3] = cell->line(1)->child(1);
4504 new_lines[4] = cell->line(2);
4505 new_lines[5] = cell->line(3);
4508 new_vertices[4], new_vertices[5]));
4511 for (
unsigned int l = lmin;
l < lmax; ++
l)
4513 new_lines[
l]->set_used_flag();
4514 new_lines[
l]->clear_user_flag();
4516 new_lines[
l]->clear_children();
4518 new_lines[
l]->set_boundary_id_internal(
4520 new_lines[
l]->set_manifold_id(cell->manifold_id());
4527 while (next_unused_cell->used() ==
true)
4531 for (
unsigned int i = 0; i < n_children; ++i)
4534 next_unused_cell->used() ==
false,
4536 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4537 subcells[i] = next_unused_cell;
4539 if (i % 2 == 1 && i < n_children - 1)
4540 while (next_unused_cell->used() ==
true)
4560 new_lines[0]->index(),
4561 new_lines[8]->index(),
4562 new_lines[4]->index(),
4563 new_lines[10]->index()));
4566 new_lines[8]->index(),
4567 new_lines[2]->index(),
4568 new_lines[5]->index(),
4569 new_lines[11]->index()));
4572 new_lines[1]->index(),
4573 new_lines[9]->index(),
4574 new_lines[10]->index(),
4575 new_lines[6]->index()));
4578 new_lines[9]->index(),
4579 new_lines[3]->index(),
4580 new_lines[11]->index(),
4581 new_lines[7]->index()));
4599 new_lines[0]->index(),
4600 new_lines[6]->index(),
4601 new_lines[2]->index(),
4602 new_lines[4]->index()));
4605 new_lines[6]->index(),
4606 new_lines[1]->index(),
4607 new_lines[3]->index(),
4608 new_lines[5]->index()));
4627 new_lines[0]->index(),
4628 new_lines[2]->index(),
4629 new_lines[4]->index(),
4630 new_lines[6]->index()));
4633 new_lines[1]->index(),
4634 new_lines[3]->index(),
4635 new_lines[6]->index(),
4636 new_lines[5]->index()));
4641 for (
unsigned int i = 0; i < n_children; ++i)
4643 subcells[i]->set_used_flag();
4644 subcells[i]->clear_refine_flag();
4645 subcells[i]->clear_user_flag();
4646 subcells[i]->clear_user_data();
4647 subcells[i]->clear_children();
4650 subcells[i]->set_material_id(cell->material_id());
4651 subcells[i]->set_manifold_id(cell->manifold_id());
4652 subcells[i]->set_subdomain_id(subdomainid);
4655 subcells[i]->set_parent(cell->index());
4663 for (
unsigned int i = 0; i < n_children / 2; ++i)
4664 cell->set_children(2 * i, subcells[2 * i]->index());
4666 cell->set_refinement_case(ref_case);
4674 for (
unsigned int c = 0; c < n_children; ++c)
4675 cell->child(c)->set_direction_flag(cell->direction_flag());
4684 template <
int spacedim>
4689 const unsigned int dim = 1;
4698 for (; cell != endc; ++cell)
4700 if (cell->refine_flag_set())
4703 std_cxx14::make_unique<
4714 unsigned int needed_vertices = 0;
4719 unsigned int flagged_cells = 0;
4723 for (; acell != aendc; ++acell)
4724 if (acell->refine_flag_set())
4729 const unsigned int used_cells =
4747 needed_vertices += flagged_cells;
4752 needed_vertices += std::count(
triangulation.vertices_used.begin(),
4769 unsigned int next_unused_vertex = 0;
4780 for (; (cell != endc) && (cell->level() ==
level); ++cell)
4781 if (cell->refine_flag_set())
4784 cell->clear_refine_flag();
4790 ++next_unused_vertex;
4794 "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."));
4809 while (next_unused_cell->used() ==
true)
4811 first_child = next_unused_cell;
4812 first_child->set_used_flag();
4813 first_child->clear_user_data();
4816 next_unused_cell->used() ==
false,
4818 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4819 second_child = next_unused_cell;
4820 second_child->set_used_flag();
4821 second_child->clear_user_data();
4826 cell->set_children(0, first_child->index());
4827 first_child->clear_children();
4830 cell->vertex_index(0), next_unused_vertex));
4831 first_child->set_material_id(cell->material_id());
4832 first_child->set_manifold_id(cell->manifold_id());
4833 first_child->set_subdomain_id(subdomainid);
4834 first_child->set_direction_flag(cell->direction_flag());
4836 first_child->set_parent(cell->index());
4840 first_child->face(1)->set_manifold_id(cell->manifold_id());
4845 first_child->set_neighbor(1, second_child);
4847 first_child->set_neighbor(0, cell->neighbor(0));
4848 else if (cell->neighbor(0)->is_active())
4854 Assert(cell->neighbor(0)->level() <= cell->level(),
4856 first_child->set_neighbor(0, cell->neighbor(0));
4862 const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4863 first_child->set_neighbor(0,
4864 cell->neighbor(0)->child(nbnb));
4869 left_neighbor = cell->neighbor(0);
4870 while (left_neighbor->has_children())
4872 left_neighbor = left_neighbor->child(nbnb);
4873 left_neighbor->set_neighbor(nbnb, first_child);
4878 second_child->clear_children();
4881 next_unused_vertex, cell->vertex_index(1)));
4882 second_child->set_neighbor(0, first_child);
4883 second_child->set_material_id(cell->material_id());
4884 second_child->set_manifold_id(cell->manifold_id());
4885 second_child->set_subdomain_id(subdomainid);
4886 second_child->set_direction_flag(cell->direction_flag());
4889 second_child->set_neighbor(1, cell->neighbor(1));
4890 else if (cell->neighbor(1)->is_active())
4892 Assert(cell->neighbor(1)->level() <= cell->level(),
4894 second_child->set_neighbor(1, cell->neighbor(1));
4899 const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4900 second_child->set_neighbor(
4901 1, cell->neighbor(1)->child(nbnb));
4904 right_neighbor = cell->neighbor(1);
4905 while (right_neighbor->has_children())
4907 right_neighbor = right_neighbor->child(nbnb);
4908 right_neighbor->set_neighbor(nbnb, second_child);
4928 template <
int spacedim>
4931 const bool check_for_distorted_cells)
4933 const unsigned int dim = 2;
4942 for (; cell != endc; ++cell)
4944 if (cell->refine_flag_set())
4947 std_cxx14::make_unique<
4960 line->clear_user_flag();
4966 unsigned int n_single_lines = 0;
4970 unsigned int n_lines_in_pairs = 0;
4976 unsigned int needed_vertices = 0;
4981 unsigned int needed_cells = 0;
4986 for (; cell != endc; ++cell)
4987 if (cell->refine_flag_set())
4998 n_single_lines += 4;
5010 n_single_lines += 1;
5018 for (
const unsigned int line_no :
5022 cell->refine_flag_set(), line_no) ==
5026 line = cell->line(line_no);
5027 if (line->has_children() ==
false)
5028 line->set_user_flag();
5035 const unsigned int used_cells =
5045 used_cells + needed_cells, 2, spacedim);
5058 if (line->user_flag_set())
5061 n_lines_in_pairs += 2;
5062 needed_vertices += 1;
5071 triangulation.faces->lines.reserve_space(n_lines_in_pairs, 0);
5074 needed_vertices += std::count(
triangulation.vertices_used.begin(),
5091 unsigned int next_unused_vertex = 0;
5103 for (; line != endl; ++line)
5104 if (line->user_flag_set())
5110 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
5111 ++next_unused_vertex;
5115 "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."));
5118 triangulation.vertices[next_unused_vertex] = line->center(
true);
5123 bool pair_found =
false;
5125 for (; next_unused_line != endl; ++next_unused_line)
5126 if (!next_unused_line->used() &&
5127 !(++next_unused_line)->used())
5140 line->set_children(0, next_unused_line->index());
5144 children[2] = {next_unused_line, ++next_unused_line};
5148 children[0]->used() ==
false,
5150 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5152 children[1]->used() ==
false,
5154 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5158 line->vertex_index(0), next_unused_vertex));
5161 next_unused_vertex, line->vertex_index(1)));
5163 children[0]->set_used_flag();
5164 children[1]->set_used_flag();
5165 children[0]->clear_children();
5166 children[1]->clear_children();
5169 children[0]->clear_user_flag();
5170 children[1]->clear_user_flag();
5173 children[0]->set_boundary_id_internal(line->boundary_id());
5174 children[1]->set_boundary_id_internal(line->boundary_id());
5176 children[0]->set_manifold_id(line->manifold_id());
5177 children[1]->set_manifold_id(line->manifold_id());
5181 line->clear_user_flag();
5190 triangulation.faces->lines.reserve_space(0, n_single_lines);
5193 cells_with_distorted_children;
5213 for (; cell != endc; ++cell)
5214 if (cell->refine_flag_set())
5222 if (cell->at_boundary())
5223 cell->set_user_flag();
5233 if ((check_for_distorted_cells ==
true) &&
5234 has_distorted_children(
5236 std::integral_constant<int, dim>(),
5237 std::integral_constant<int, spacedim>()))
5245 return cells_with_distorted_children;
5253 template <
int spacedim>
5256 const bool check_for_distorted_cells)
5258 const unsigned int dim = 3;
5273 for (; cell != endc; ++cell)
5275 if (cell->refine_flag_set())
5278 std_cxx14::make_unique<
5293 line->clear_user_flag();
5298 quad->clear_user_flag();
5321 unsigned int needed_vertices = 0;
5322 unsigned int needed_lines_single = 0;
5323 unsigned int needed_quads_single = 0;
5324 unsigned int needed_lines_pair = 0;
5325 unsigned int needed_quads_pair = 0;
5330 unsigned int new_cells = 0;
5335 for (; acell != aendc; ++acell)
5336 if (acell->refine_flag_set())
5346 ++needed_quads_single;
5354 ++needed_lines_single;
5355 needed_quads_single += 4;
5362 needed_lines_single += 6;
5363 needed_quads_single += 12;
5377 for (
const unsigned int face :
5381 aface = acell->face(face);
5388 acell->face_orientation(face),
5389 acell->face_flip(face),
5390 acell->face_rotation(face));
5395 if (face_ref_case ==
5398 if (aface->number_of_children() < 4)
5401 aface->set_user_flag();
5403 else if (aface->refinement_case() != face_ref_case)
5414 Assert(aface->refinement_case() ==
5416 dim - 1>::isotropic_refinement ||
5417 aface->refinement_case() ==
5420 aface->set_user_index(face_ref_case);
5426 for (
unsigned int line = 0;
5427 line < GeometryInfo<dim>::lines_per_cell;
5431 !acell->line(line)->has_children())
5432 acell->line(line)->set_user_flag();
5438 const unsigned int used_cells =
5448 used_cells + new_cells, 3, spacedim);
5460 if (quad->user_flag_set())
5466 needed_quads_pair += 4;
5467 needed_lines_pair += 4;
5468 needed_vertices += 1;
5470 if (quad->user_index())
5474 needed_quads_pair += 2;
5475 needed_lines_single += 1;
5488 if (quad->has_children())
5490 Assert(quad->refinement_case() ==
5493 if ((face_refinement_cases[quad->user_index()] ==
5495 (quad->child(0)->line_index(1) + 1 !=
5496 quad->child(2)->line_index(1))) ||
5497 (face_refinement_cases[quad->user_index()] ==
5499 (quad->child(0)->line_index(3) + 1 !=
5500 quad->child(1)->line_index(3))))
5501 needed_lines_pair += 2;
5510 if (line->user_flag_set())
5512 needed_lines_pair += 2;
5513 needed_vertices += 1;
5518 needed_lines_single);
5521 needed_quads_single);
5525 needed_vertices += std::count(
triangulation.vertices_used.begin(),
5549 for (
const auto &cell :
triangulation.active_cell_iterators())
5550 if (!cell->refine_flag_set())
5551 for (
unsigned int line = 0;
5552 line < GeometryInfo<dim>::lines_per_cell;
5554 if (cell->line(line)->has_children())
5555 for (
unsigned int c = 0; c < 2; ++c)
5556 Assert(cell->line(line)->child(c)->user_flag_set() ==
false,
5568 unsigned int next_unused_vertex = 0;
5579 for (; line != endl; ++line)
5580 if (line->user_flag_set())
5586 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
5587 ++next_unused_vertex;
5591 "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."));
5594 triangulation.vertices[next_unused_vertex] = line->center(
true);
5608 line->set_children(0, next_unused_line->index());
5612 children[2] = {next_unused_line, ++next_unused_line};
5617 children[0]->used() ==
false,
5619 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5621 children[1]->used() ==
false,
5623 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5627 line->vertex_index(0), next_unused_vertex));
5630 next_unused_vertex, line->vertex_index(1)));
5632 children[0]->set_used_flag();
5633 children[1]->set_used_flag();
5634 children[0]->clear_children();
5635 children[1]->clear_children();
5638 children[0]->clear_user_flag();
5639 children[1]->clear_user_flag();
5641 children[0]->set_boundary_id_internal(line->boundary_id());
5642 children[1]->set_boundary_id_internal(line->boundary_id());
5644 children[0]->set_manifold_id(line->manifold_id());
5645 children[1]->set_manifold_id(line->manifold_id());
5650 line->clear_user_flag();
5693 for (; quad != endq; ++quad)
5695 if (quad->user_index())
5698 face_refinement_cases[quad->user_index()];
5707 if (aniso_quad_ref_case == quad->refinement_case())
5710 Assert(quad->refinement_case() ==
5712 quad->refinement_case() ==
5717 Assert(quad->user_index() ==
5719 quad->user_index() ==
5731 new_line->used() ==
false,
5733 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5752 quad->line(2)->child(0)->vertex_index(1);
5754 quad->line(3)->child(0)->vertex_index(1);
5759 quad->line(0)->child(0)->vertex_index(1);
5761 quad->line(1)->child(0)->vertex_index(1);
5767 new_line->set_used_flag();
5768 new_line->clear_user_flag();
5770 new_line->clear_children();
5771 new_line->set_boundary_id_internal(quad->boundary_id());
5772 new_line->set_manifold_id(quad->manifold_id());
5780 const unsigned int index[2][2] = {
5792 new_quads[0] = next_unused_quad;
5794 new_quads[0]->used() ==
false,
5796 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5799 new_quads[1] = next_unused_quad;
5801 new_quads[1]->used() ==
false,
5803 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5810 quad->line_index(0),
5813 ->child(index[0][quad->line_orientation(2)])
5816 ->child(index[0][quad->line_orientation(3)])
5821 quad->line_index(1),
5823 ->child(index[1][quad->line_orientation(2)])
5826 ->child(index[1][quad->line_orientation(3)])
5834 ->child(index[0][quad->line_orientation(0)])
5837 ->child(index[0][quad->line_orientation(1)])
5839 quad->line_index(2),
5840 new_line->index()));
5844 ->child(index[1][quad->line_orientation(0)])
5847 ->child(index[1][quad->line_orientation(1)])
5850 quad->line_index(3)));
5853 for (
const auto &new_quad : new_quads)
5855 new_quad->set_used_flag();
5856 new_quad->clear_user_flag();
5858 new_quad->clear_children();
5859 new_quad->set_boundary_id_internal(quad->boundary_id());
5860 new_quad->set_manifold_id(quad->manifold_id());
5864 for (
unsigned int j = 0;
5865 j < GeometryInfo<dim>::lines_per_face;
5867 new_quad->set_line_orientation(j,
true);
5873 new_quads[0]->set_line_orientation(
5874 0, quad->line_orientation(0));
5875 new_quads[0]->set_line_orientation(
5876 2, quad->line_orientation(2));
5877 new_quads[1]->set_line_orientation(
5878 1, quad->line_orientation(1));
5879 new_quads[1]->set_line_orientation(
5880 3, quad->line_orientation(3));
5883 new_quads[0]->set_line_orientation(
5884 3, quad->line_orientation(3));
5885 new_quads[1]->set_line_orientation(
5886 2, quad->line_orientation(2));
5890 new_quads[0]->set_line_orientation(
5891 1, quad->line_orientation(1));
5892 new_quads[1]->set_line_orientation(
5893 0, quad->line_orientation(0));
5899 if (quad->refinement_case() ==
5922 if (aniso_quad_ref_case ==
5925 old_child[0] = quad->child(0)->line(1);
5926 old_child[1] = quad->child(2)->line(1);
5930 Assert(aniso_quad_ref_case ==
5934 old_child[0] = quad->child(0)->line(3);
5935 old_child[1] = quad->child(1)->line(3);
5938 if (old_child[0]->index() + 1 != old_child[1]->index())
5944 spacedim>::raw_line_iterator
5947 new_child[0] = new_child[1] =
5952 new_child[0]->set_used_flag();
5953 new_child[1]->set_used_flag();
5955 const int old_index_0 = old_child[0]->index(),
5956 old_index_1 = old_child[1]->index(),
5957 new_index_0 = new_child[0]->index(),
5958 new_index_1 = new_child[1]->index();
5962 for (
unsigned int q = 0;
5965 for (
unsigned int l = 0;
5966 l < GeometryInfo<dim>::lines_per_face;
5969 const int this_index =
5971 if (this_index == old_index_0)
5973 .set_face(
l, new_index_0);
5974 else if (this_index == old_index_1)
5976 .set_face(
l, new_index_1);
5980 for (
unsigned int i = 0; i < 2; ++i)
5982 Assert(!old_child[i]->has_children(),
5988 old_child[i]->vertex_index(
5990 new_child[i]->set_boundary_id_internal(
5992 new_child[i]->set_manifold_id(
5994 new_child[i]->set_user_index(
5995 old_child[i]->user_index());
5996 if (old_child[i]->user_flag_set())
5997 new_child[i]->set_user_flag();
5999 new_child[i]->clear_user_flag();
6001 new_child[i]->clear_children();
6003 old_child[i]->clear_user_flag();
6004 old_child[i]->clear_user_index();
6005 old_child[i]->clear_used_flag();
6011 if (aniso_quad_ref_case ==
6014 new_line->set_children(
6015 0, quad->child(0)->line_index(1));
6016 Assert(new_line->child(1) ==
6017 quad->child(2)->line(1),
6054 quad_iterator switch_1 = quad->child(1),
6055 switch_2 = quad->child(2);
6056 const int switch_1_index = switch_1->index();
6057 const int switch_2_index = switch_2->index();
6058 for (
unsigned int l = 0;
6061 for (
unsigned int h = 0;
6065 for (
const unsigned int q :
6068 const int face_index =
6072 if (face_index == switch_1_index)
6075 .set_face(q, switch_2_index);
6076 else if (face_index == switch_2_index)
6079 .set_face(q, switch_1_index);
6083 const unsigned int switch_1_lines[4] = {
6084 switch_1->line_index(0),
6085 switch_1->line_index(1),
6086 switch_1->line_index(2),
6087 switch_1->line_index(3)};
6088 const bool switch_1_line_orientations[4] = {
6089 switch_1->line_orientation(0),
6090 switch_1->line_orientation(1),
6091 switch_1->line_orientation(2),
6092 switch_1->line_orientation(3)};
6094 switch_1->boundary_id();
6095 const unsigned int switch_1_user_index =
6096 switch_1->user_index();
6097 const bool switch_1_user_flag =
6098 switch_1->user_flag_set();
6100 switch_1_refinement_case =
6101 switch_1->refinement_case();
6102 const int switch_1_first_child_pair =
6103 (switch_1_refinement_case ?
6104 switch_1->child_index(0) :
6106 const int switch_1_second_child_pair =
6107 (switch_1_refinement_case ==
6109 switch_1->child_index(2) :
6114 2>(switch_2->line_index(0),
6115 switch_2->line_index(1),
6116 switch_2->line_index(2),
6117 switch_2->line_index(3)));
6118 switch_1->set_line_orientation(
6119 0, switch_2->line_orientation(0));
6120 switch_1->set_line_orientation(
6121 1, switch_2->line_orientation(1));
6122 switch_1->set_line_orientation(
6123 2, switch_2->line_orientation(2));
6124 switch_1->set_line_orientation(
6125 3, switch_2->line_orientation(3));
6126 switch_1->set_boundary_id_internal(
6127 switch_2->boundary_id());
6128 switch_1->set_manifold_id(switch_2->manifold_id());
6129 switch_1->set_user_index(switch_2->user_index());
6130 if (switch_2->user_flag_set())
6131 switch_1->set_user_flag();
6133 switch_1->clear_user_flag();
6134 switch_1->clear_refinement_case();
6135 switch_1->set_refinement_case(
6136 switch_2->refinement_case());
6137 switch_1->clear_children();
6138 if (switch_2->refinement_case())
6139 switch_1->set_children(0,
6140 switch_2->child_index(0));
6141 if (switch_2->refinement_case() ==
6143 switch_1->set_children(2,
6144 switch_2->child_index(2));
6148 2>(switch_1_lines[0],
6151 switch_1_lines[3]));
6152 switch_2->set_line_orientation(
6153 0, switch_1_line_orientations[0]);
6154 switch_2->set_line_orientation(
6155 1, switch_1_line_orientations[1]);
6156 switch_2->set_line_orientation(
6157 2, switch_1_line_orientations[2]);
6158 switch_2->set_line_orientation(
6159 3, switch_1_line_orientations[3]);
6160 switch_2->set_boundary_id_internal(
6161 switch_1_boundary_id);
6162 switch_2->set_manifold_id(switch_1->manifold_id());
6163 switch_2->set_user_index(switch_1_user_index);
6164 if (switch_1_user_flag)
6165 switch_2->set_user_flag();
6167 switch_2->clear_user_flag();
6168 switch_2->clear_refinement_case();
6169 switch_2->set_refinement_case(
6170 switch_1_refinement_case);
6171 switch_2->clear_children();
6172 switch_2->set_children(0,
6173 switch_1_first_child_pair);
6174 switch_2->set_children(2,
6175 switch_1_second_child_pair);
6177 new_quads[0]->set_refinement_case(
6179 new_quads[0]->set_children(0, quad->child_index(0));
6180 new_quads[1]->set_refinement_case(
6182 new_quads[1]->set_children(0, quad->child_index(2));
6186 new_quads[0]->set_refinement_case(
6188 new_quads[0]->set_children(0, quad->child_index(0));
6189 new_quads[1]->set_refinement_case(
6191 new_quads[1]->set_children(0, quad->child_index(2));
6192 new_line->set_children(
6193 0, quad->child(0)->line_index(3));
6194 Assert(new_line->child(1) ==
6195 quad->child(1)->line(3),
6198 quad->clear_children();
6202 quad->set_children(0, new_quads[0]->index());
6204 quad->set_refinement_case(aniso_quad_ref_case);
6211 if (quad->user_flag_set())
6223 ++next_unused_vertex;
6227 "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."));
6235 quad->refinement_case();
6241 quad->child(0)->set_user_index(
6243 quad->child(1)->set_user_index(
6249 middle_line = quad->child(0)->line(1);
6251 middle_line = quad->child(0)->line(3);
6258 if (!middle_line->has_children())
6266 middle_line->center(
true);
6278 middle_line->set_children(
6279 0, next_unused_line->index());
6283 raw_line_iterator children[2] = {
6284 next_unused_line, ++next_unused_line};
6290 children[0]->used() ==
false,
6292 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6294 children[1]->used() ==
false,
6296 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6300 1>(middle_line->vertex_index(0),
6301 next_unused_vertex));
6304 1>(next_unused_vertex,
6305 middle_line->vertex_index(1)));
6307 children[0]->set_used_flag();
6308 children[1]->set_used_flag();
6309 children[0]->clear_children();
6310 children[1]->clear_children();
6313 children[0]->clear_user_flag();
6314 children[1]->clear_user_flag();
6316 children[0]->set_boundary_id_internal(
6317 middle_line->boundary_id());
6318 children[1]->set_boundary_id_internal(
6319 middle_line->boundary_id());
6321 children[0]->set_manifold_id(
6322 middle_line->manifold_id());
6323 children[1]->set_manifold_id(
6324 middle_line->manifold_id());
6330 quad->clear_user_flag();
6360 quad->center(
true,
true);
6369 for (
unsigned int i = 0; i < 4; ++i)
6382 new_lines[i] = next_unused_line;
6386 new_lines[i]->used() ==
false,
6388 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6409 quad->line(0)->child(0)->vertex_index(1),
6410 quad->line(1)->child(0)->vertex_index(1),
6411 quad->line(2)->child(0)->vertex_index(1),
6412 quad->line(3)->child(0)->vertex_index(1),
6413 next_unused_vertex};
6428 for (
const auto &new_line : new_lines)
6430 new_line->set_used_flag();
6431 new_line->clear_user_flag();
6433 new_line->clear_children();
6434 new_line->set_boundary_id_internal(quad->boundary_id());
6435 new_line->set_manifold_id(quad->manifold_id());
6454 const unsigned int index[2][2] = {
6458 const int line_indices[12] = {
6460 ->child(index[0][quad->line_orientation(0)])
6463 ->child(index[1][quad->line_orientation(0)])
6466 ->child(index[0][quad->line_orientation(1)])
6469 ->child(index[1][quad->line_orientation(1)])
6472 ->child(index[0][quad->line_orientation(2)])
6475 ->child(index[1][quad->line_orientation(2)])
6478 ->child(index[0][quad->line_orientation(3)])
6481 ->child(index[1][quad->line_orientation(3)])
6483 new_lines[0]->index(),
6484 new_lines[1]->index(),
6485 new_lines[2]->index(),
6486 new_lines[3]->index()};
6498 new_quads[0] = next_unused_quad;
6500 new_quads[0]->used() ==
false,
6502 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6505 new_quads[1] = next_unused_quad;
6507 new_quads[1]->used() ==
false,
6509 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6514 new_quads[2] = next_unused_quad;
6516 new_quads[2]->used() ==
false,
6518 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6521 new_quads[3] = next_unused_quad;
6523 new_quads[3]->used() ==
false,
6525 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6528 quad->set_children(0, new_quads[0]->index());
6529 quad->set_children(2, new_quads[2]->index());
6563 for (
const auto &new_quad : new_quads)
6565 new_quad->set_used_flag();
6566 new_quad->clear_user_flag();
6568 new_quad->clear_children();
6569 new_quad->set_boundary_id_internal(quad->boundary_id());
6570 new_quad->set_manifold_id(quad->manifold_id());
6574 for (
unsigned int j = 0;
6575 j < GeometryInfo<dim>::lines_per_face;
6577 new_quad->set_line_orientation(j,
true);
6583 new_quads[0]->set_line_orientation(
6584 0, quad->line_orientation(0));
6585 new_quads[0]->set_line_orientation(
6586 2, quad->line_orientation(2));
6587 new_quads[1]->set_line_orientation(
6588 1, quad->line_orientation(1));
6589 new_quads[1]->set_line_orientation(
6590 2, quad->line_orientation(2));
6591 new_quads[2]->set_line_orientation(
6592 0, quad->line_orientation(0));
6593 new_quads[2]->set_line_orientation(
6594 3, quad->line_orientation(3));
6595 new_quads[3]->set_line_orientation(
6596 1, quad->line_orientation(1));
6597 new_quads[3]->set_line_orientation(
6598 3, quad->line_orientation(3));