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));
6602 quad->clear_user_flag();
6613 cells_with_distorted_children;
6627 for (; hex != endh; ++hex)
6628 if (hex->refine_flag_set())
6636 hex->clear_refine_flag();
6637 hex->set_refinement_case(ref_case);
6648 unsigned int n_new_lines = 0;
6649 unsigned int n_new_quads = 0;
6650 unsigned int n_new_hexes = 0;
6681 new_lines(n_new_lines);
6682 for (
unsigned int i = 0; i < n_new_lines; ++i)
6689 new_lines[i]->used() ==
false,
6691 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6692 new_lines[i]->set_used_flag();
6693 new_lines[i]->clear_user_flag();
6694 new_lines[i]->clear_user_data();
6695 new_lines[i]->clear_children();
6697 new_lines[i]->set_boundary_id_internal(
6701 new_lines[i]->set_manifold_id(hex->manifold_id());
6708 new_quads(n_new_quads);
6709 for (
unsigned int i = 0; i < n_new_quads; ++i)
6716 new_quads[i]->used() ==
false,
6718 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6719 new_quads[i]->set_used_flag();
6720 new_quads[i]->clear_user_flag();
6721 new_quads[i]->clear_user_data();
6722 new_quads[i]->clear_children();
6724 new_quads[i]->set_boundary_id_internal(
6728 new_quads[i]->set_manifold_id(hex->manifold_id());
6731 for (
unsigned int j = 0;
6732 j < GeometryInfo<dim>::lines_per_face;
6734 new_quads[i]->set_line_orientation(j,
true);
6743 new_hexes(n_new_hexes);
6744 for (
unsigned int i = 0; i < n_new_hexes; ++i)
6753 new_hexes[i] = next_unused_hex;
6756 new_hexes[i]->used() ==
false,
6758 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6759 new_hexes[i]->set_used_flag();
6760 new_hexes[i]->clear_user_flag();
6761 new_hexes[i]->clear_user_data();
6762 new_hexes[i]->clear_children();
6765 new_hexes[i]->set_material_id(hex->material_id());
6766 new_hexes[i]->set_manifold_id(hex->manifold_id());
6767 new_hexes[i]->set_subdomain_id(subdomainid);
6770 new_hexes[i]->set_parent(hex->index());
6782 for (
const unsigned int f :
6785 new_hexes[i]->set_face_orientation(f,
true);
6786 new_hexes[i]->set_face_flip(f,
false);
6787 new_hexes[i]->set_face_rotation(f,
false);
6791 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
6792 hex->set_children(2 * i, new_hexes[2 * i]->index());
6799 const bool f_or[6] = {hex->face_orientation(0),
6800 hex->face_orientation(1),
6801 hex->face_orientation(2),
6802 hex->face_orientation(3),
6803 hex->face_orientation(4),
6804 hex->face_orientation(5)};
6807 const bool f_fl[6] = {hex->face_flip(0),
6815 const bool f_ro[6] = {hex->face_rotation(0),
6816 hex->face_rotation(1),
6817 hex->face_rotation(2),
6818 hex->face_rotation(3),
6819 hex->face_rotation(4),
6820 hex->face_rotation(5)};
6830 const unsigned int child_at_origin[2][2][2] = {
6917 raw_line_iterator lines[4] = {
6918 hex->face(2)->child(0)->line(
6919 (hex->face(2)->refinement_case() ==
6923 hex->face(3)->child(0)->line(
6924 (hex->face(3)->refinement_case() ==
6928 hex->face(4)->child(0)->line(
6929 (hex->face(4)->refinement_case() ==
6933 hex->face(5)->child(0)->line(
6934 (hex->face(5)->refinement_case() ==
6940 unsigned int line_indices[4];
6941 for (
unsigned int i = 0; i < 4; ++i)
6942 line_indices[i] = lines[i]->index();
6949 bool line_orientation[4];
6955 const unsigned int middle_vertices[2] = {
6956 hex->line(2)->child(0)->vertex_index(1),
6957 hex->line(7)->child(0)->vertex_index(1)};
6959 for (
unsigned int i = 0; i < 4; ++i)
6960 if (lines[i]->vertex_index(i % 2) ==
6961 middle_vertices[i % 2])
6962 line_orientation[i] =
true;
6967 Assert(lines[i]->vertex_index((i + 1) % 2) ==
6968 middle_vertices[i % 2],
6970 line_orientation[i] =
false;
6982 new_quads[0]->set_line_orientation(
6983 0, line_orientation[0]);
6984 new_quads[0]->set_line_orientation(
6985 1, line_orientation[1]);
6986 new_quads[0]->set_line_orientation(
6987 2, line_orientation[2]);
6988 new_quads[0]->set_line_orientation(
6989 3, line_orientation[3]);
7022 const int quad_indices[11] = {
7023 new_quads[0]->index(),
7025 hex->face(0)->index(),
7027 hex->face(1)->index(),
7029 hex->face(2)->child_index(
7030 child_at_origin[hex->face(2)->refinement_case() -
7031 1][f_fl[2]][f_ro[2]]),
7032 hex->face(2)->child_index(
7034 child_at_origin[hex->face(2)->refinement_case() -
7035 1][f_fl[2]][f_ro[2]]),
7037 hex->face(3)->child_index(
7038 child_at_origin[hex->face(3)->refinement_case() -
7039 1][f_fl[3]][f_ro[3]]),
7040 hex->face(3)->child_index(
7042 child_at_origin[hex->face(3)->refinement_case() -
7043 1][f_fl[3]][f_ro[3]]),
7045 hex->face(4)->child_index(
7046 child_at_origin[hex->face(4)->refinement_case() -
7047 1][f_fl[4]][f_ro[4]]),
7048 hex->face(4)->child_index(
7050 child_at_origin[hex->face(4)->refinement_case() -
7051 1][f_fl[4]][f_ro[4]]),
7053 hex->face(5)->child_index(
7054 child_at_origin[hex->face(5)->refinement_case() -
7055 1][f_fl[5]][f_ro[5]]),
7056 hex->face(5)->child_index(
7058 child_at_origin[hex->face(5)->refinement_case() -
7059 1][f_fl[5]][f_ro[5]])
7147 raw_line_iterator lines[4] = {
7148 hex->face(0)->child(0)->line(
7149 (hex->face(0)->refinement_case() ==
7153 hex->face(1)->child(0)->line(
7154 (hex->face(1)->refinement_case() ==
7158 hex->face(4)->child(0)->line(
7159 (hex->face(4)->refinement_case() ==
7163 hex->face(5)->child(0)->line(
7164 (hex->face(5)->refinement_case() ==
7170 unsigned int line_indices[4];
7171 for (
unsigned int i = 0; i < 4; ++i)
7172 line_indices[i] = lines[i]->index();
7179 bool line_orientation[4];
7185 const unsigned int middle_vertices[2] = {
7186 hex->line(0)->child(0)->vertex_index(1),
7187 hex->line(5)->child(0)->vertex_index(1)};
7189 for (
unsigned int i = 0; i < 4; ++i)
7190 if (lines[i]->vertex_index(i % 2) ==
7191 middle_vertices[i % 2])
7192 line_orientation[i] =
true;
7196 Assert(lines[i]->vertex_index((i + 1) % 2) ==
7197 middle_vertices[i % 2],
7199 line_orientation[i] =
false;
7211 new_quads[0]->set_line_orientation(
7212 0, line_orientation[2]);
7213 new_quads[0]->set_line_orientation(
7214 1, line_orientation[3]);
7215 new_quads[0]->set_line_orientation(
7216 2, line_orientation[0]);
7217 new_quads[0]->set_line_orientation(
7218 3, line_orientation[1]);
7251 const int quad_indices[11] = {
7252 new_quads[0]->index(),
7254 hex->face(0)->child_index(
7255 child_at_origin[hex->face(0)->refinement_case() -
7256 1][f_fl[0]][f_ro[0]]),
7257 hex->face(0)->child_index(
7259 child_at_origin[hex->face(0)->refinement_case() -
7260 1][f_fl[0]][f_ro[0]]),
7262 hex->face(1)->child_index(
7263 child_at_origin[hex->face(1)->refinement_case() -
7264 1][f_fl[1]][f_ro[1]]),
7265 hex->face(1)->child_index(
7267 child_at_origin[hex->face(1)->refinement_case() -
7268 1][f_fl[1]][f_ro[1]]),
7270 hex->face(2)->index(),
7272 hex->face(3)->index(),
7274 hex->face(4)->child_index(
7275 child_at_origin[hex->face(4)->refinement_case() -
7276 1][f_fl[4]][f_ro[4]]),
7277 hex->face(4)->child_index(
7279 child_at_origin[hex->face(4)->refinement_case() -
7280 1][f_fl[4]][f_ro[4]]),
7282 hex->face(5)->child_index(
7283 child_at_origin[hex->face(5)->refinement_case() -
7284 1][f_fl[5]][f_ro[5]]),
7285 hex->face(5)->child_index(
7287 child_at_origin[hex->face(5)->refinement_case() -
7288 1][f_fl[5]][f_ro[5]])
7378 raw_line_iterator lines[4] = {
7379 hex->face(0)->child(0)->line(
7380 (hex->face(0)->refinement_case() ==
7384 hex->face(1)->child(0)->line(
7385 (hex->face(1)->refinement_case() ==
7389 hex->face(2)->child(0)->line(
7390 (hex->face(2)->refinement_case() ==
7394 hex->face(3)->child(0)->line(
7395 (hex->face(3)->refinement_case() ==
7401 unsigned int line_indices[4];
7402 for (
unsigned int i = 0; i < 4; ++i)
7403 line_indices[i] = lines[i]->index();
7410 bool line_orientation[4];
7416 const unsigned int middle_vertices[2] = {
7417 middle_vertex_index<dim, spacedim>(hex->line(8)),
7418 middle_vertex_index<dim, spacedim>(hex->line(11))};
7420 for (
unsigned int i = 0; i < 4; ++i)
7421 if (lines[i]->vertex_index(i % 2) ==
7422 middle_vertices[i % 2])
7423 line_orientation[i] =
true;
7427 Assert(lines[i]->vertex_index((i + 1) % 2) ==
7428 middle_vertices[i % 2],
7430 line_orientation[i] =
false;
7442 new_quads[0]->set_line_orientation(
7443 0, line_orientation[0]);
7444 new_quads[0]->set_line_orientation(
7445 1, line_orientation[1]);
7446 new_quads[0]->set_line_orientation(
7447 2, line_orientation[2]);
7448 new_quads[0]->set_line_orientation(
7449 3, line_orientation[3]);
7483 const int quad_indices[11] = {
7484 new_quads[0]->index(),
7486 hex->face(0)->child_index(
7487 child_at_origin[hex->face(0)->refinement_case() -
7488 1][f_fl[0]][f_ro[0]]),
7489 hex->face(0)->child_index(
7491 child_at_origin[hex->face(0)->refinement_case() -
7492 1][f_fl[0]][f_ro[0]]),
7494 hex->face(1)->child_index(
7495 child_at_origin[hex->face(1)->refinement_case() -
7496 1][f_fl[1]][f_ro[1]]),
7497 hex->face(1)->child_index(
7499 child_at_origin[hex->face(1)->refinement_case() -
7500 1][f_fl[1]][f_ro[1]]),
7502 hex->face(2)->child_index(
7503 child_at_origin[hex->face(2)->refinement_case() -
7504 1][f_fl[2]][f_ro[2]]),
7505 hex->face(2)->child_index(
7507 child_at_origin[hex->face(2)->refinement_case() -
7508 1][f_fl[2]][f_ro[2]]),
7510 hex->face(3)->child_index(
7511 child_at_origin[hex->face(3)->refinement_case() -
7512 1][f_fl[3]][f_ro[3]]),
7513 hex->face(3)->child_index(
7515 child_at_origin[hex->face(3)->refinement_case() -
7516 1][f_fl[3]][f_ro[3]]),
7518 hex->face(4)->index(),
7520 hex->face(5)->index()
7566 1>(middle_vertex_index<dim, spacedim>(
7568 middle_vertex_index<dim, spacedim>(
7637 spacedim>::raw_line_iterator lines[13] = {
7638 hex->face(0)->child(0)->line(
7639 (hex->face(0)->refinement_case() ==
7643 hex->face(1)->child(0)->line(
7644 (hex->face(1)->refinement_case() ==
7648 hex->face(2)->child(0)->line(
7649 (hex->face(2)->refinement_case() ==
7653 hex->face(3)->child(0)->line(
7654 (hex->face(3)->refinement_case() ==
7662 0, f_or[4], f_fl[4], f_ro[4]))
7665 1, f_or[4], f_fl[4], f_ro[4])),
7669 3, f_or[4], f_fl[4], f_ro[4]))
7672 0, f_or[4], f_fl[4], f_ro[4])),
7676 0, f_or[4], f_fl[4], f_ro[4]))
7679 3, f_or[4], f_fl[4], f_ro[4])),
7683 3, f_or[4], f_fl[4], f_ro[4]))
7686 2, f_or[4], f_fl[4], f_ro[4])),
7691 0, f_or[5], f_fl[5], f_ro[5]))
7694 1, f_or[5], f_fl[5], f_ro[5])),
7698 3, f_or[5], f_fl[5], f_ro[5]))
7701 0, f_or[5], f_fl[5], f_ro[5])),
7705 0, f_or[5], f_fl[5], f_ro[5]))
7708 3, f_or[5], f_fl[5], f_ro[5])),
7712 3, f_or[5], f_fl[5], f_ro[5]))
7715 2, f_or[5], f_fl[5], f_ro[5])),
7720 unsigned int line_indices[13];
7721 for (
unsigned int i = 0; i < 13; ++i)
7722 line_indices[i] = lines[i]->index();
7729 bool line_orientation[13];
7733 const unsigned int middle_vertices[4] = {
7734 hex->line(0)->child(0)->vertex_index(1),
7735 hex->line(1)->child(0)->vertex_index(1),
7736 hex->line(2)->child(0)->vertex_index(1),
7737 hex->line(3)->child(0)->vertex_index(1),
7743 for (
unsigned int i = 0; i < 4; ++i)
7744 if (lines[i]->vertex_index(0) == middle_vertices[i])
7745 line_orientation[i] =
true;
7749 Assert(lines[i]->vertex_index(1) ==
7752 line_orientation[i] =
false;
7761 for (
unsigned int i = 4; i < 12; ++i)
7762 if (lines[i]->vertex_index((i + 1) % 2) ==
7763 middle_vertex_index<dim, spacedim>(
7764 hex->face(3 + i / 4)))
7765 line_orientation[i] =
true;
7770 Assert(lines[i]->vertex_index(i % 2) ==
7771 (middle_vertex_index<dim, spacedim>(
7772 hex->face(3 + i / 4))),
7774 line_orientation[i] =
false;
7779 line_orientation[12] =
true;
7812 2>(line_indices[12],
7829 new_quads[0]->set_line_orientation(
7830 0, line_orientation[2]);
7831 new_quads[0]->set_line_orientation(
7832 2, line_orientation[4]);
7833 new_quads[0]->set_line_orientation(
7834 3, line_orientation[8]);
7836 new_quads[1]->set_line_orientation(
7837 1, line_orientation[3]);
7838 new_quads[1]->set_line_orientation(
7839 2, line_orientation[5]);
7840 new_quads[1]->set_line_orientation(
7841 3, line_orientation[9]);
7843 new_quads[2]->set_line_orientation(
7844 0, line_orientation[6]);
7845 new_quads[2]->set_line_orientation(
7846 1, line_orientation[10]);
7847 new_quads[2]->set_line_orientation(
7848 2, line_orientation[0]);
7850 new_quads[3]->set_line_orientation(
7851 0, line_orientation[7]);
7852 new_quads[3]->set_line_orientation(
7853 1, line_orientation[11]);
7854 new_quads[3]->set_line_orientation(
7855 3, line_orientation[1]);
7888 const int quad_indices[20] = {
7889 new_quads[0]->index(),
7890 new_quads[1]->index(),
7891 new_quads[2]->index(),
7892 new_quads[3]->index(),
7894 hex->face(0)->child_index(
7895 child_at_origin[hex->face(0)->refinement_case() -
7896 1][f_fl[0]][f_ro[0]]),
7897 hex->face(0)->child_index(
7899 child_at_origin[hex->face(0)->refinement_case() -
7900 1][f_fl[0]][f_ro[0]]),
7902 hex->face(1)->child_index(
7903 child_at_origin[hex->face(1)->refinement_case() -
7904 1][f_fl[1]][f_ro[1]]),
7905 hex->face(1)->child_index(
7907 child_at_origin[hex->face(1)->refinement_case() -
7908 1][f_fl[1]][f_ro[1]]),
7910 hex->face(2)->child_index(
7911 child_at_origin[hex->face(2)->refinement_case() -
7912 1][f_fl[2]][f_ro[2]]),
7913 hex->face(2)->child_index(
7915 child_at_origin[hex->face(2)->refinement_case() -
7916 1][f_fl[2]][f_ro[2]]),
7918 hex->face(3)->child_index(
7919 child_at_origin[hex->face(3)->refinement_case() -
7920 1][f_fl[3]][f_ro[3]]),
7921 hex->face(3)->child_index(
7923 child_at_origin[hex->face(3)->refinement_case() -
7924 1][f_fl[3]][f_ro[3]]),
7926 hex->face(4)->isotropic_child_index(
7928 0, f_or[4], f_fl[4], f_ro[4])),
7929 hex->face(4)->isotropic_child_index(
7931 1, f_or[4], f_fl[4], f_ro[4])),
7932 hex->face(4)->isotropic_child_index(
7934 2, f_or[4], f_fl[4], f_ro[4])),
7935 hex->face(4)->isotropic_child_index(
7937 3, f_or[4], f_fl[4], f_ro[4])),
7939 hex->face(5)->isotropic_child_index(
7941 0, f_or[5], f_fl[5], f_ro[5])),
7942 hex->face(5)->isotropic_child_index(
7944 1, f_or[5], f_fl[5], f_ro[5])),
7945 hex->face(5)->isotropic_child_index(
7947 2, f_or[5], f_fl[5], f_ro[5])),
7948 hex->face(5)->isotropic_child_index(
7950 3, f_or[5], f_fl[5], f_ro[5]))};
8011 1>(middle_vertex_index<dim, spacedim>(
8013 middle_vertex_index<dim, spacedim>(
8082 spacedim>::raw_line_iterator lines[13] = {
8083 hex->face(0)->child(0)->line(
8084 (hex->face(0)->refinement_case() ==
8088 hex->face(1)->child(0)->line(
8089 (hex->face(1)->refinement_case() ==
8093 hex->face(4)->child(0)->line(
8094 (hex->face(4)->refinement_case() ==
8098 hex->face(5)->child(0)->line(
8099 (hex->face(5)->refinement_case() ==
8107 0, f_or[2], f_fl[2], f_ro[2]))
8110 3, f_or[2], f_fl[2], f_ro[2])),
8114 3, f_or[2], f_fl[2], f_ro[2]))
8117 2, f_or[2], f_fl[2], f_ro[2])),
8121 0, f_or[2], f_fl[2], f_ro[2]))
8124 1, f_or[2], f_fl[2], f_ro[2])),
8128 3, f_or[2], f_fl[2], f_ro[2]))
8131 0, f_or[2], f_fl[2], f_ro[2])),
8136 0, f_or[3], f_fl[3], f_ro[3]))
8139 3, f_or[3], f_fl[3], f_ro[3])),
8143 3, f_or[3], f_fl[3], f_ro[3]))
8146 2, f_or[3], f_fl[3], f_ro[3])),
8150 0, f_or[3], f_fl[3], f_ro[3]))
8153 1, f_or[3], f_fl[3], f_ro[3])),
8157 3, f_or[3], f_fl[3], f_ro[3]))
8160 0, f_or[3], f_fl[3], f_ro[3])),
8165 unsigned int line_indices[13];
8166 for (
unsigned int i = 0; i < 13; ++i)
8167 line_indices[i] = lines[i]->index();
8174 bool line_orientation[13];
8178 const unsigned int middle_vertices[4] = {
8179 hex->line(8)->child(0)->vertex_index(1),
8180 hex->line(9)->child(0)->vertex_index(1),
8181 hex->line(2)->child(0)->vertex_index(1),
8182 hex->line(6)->child(0)->vertex_index(1),
8187 for (
unsigned int i = 0; i < 4; ++i)
8188 if (lines[i]->vertex_index(0) == middle_vertices[i])
8189 line_orientation[i] =
true;
8193 Assert(lines[i]->vertex_index(1) ==
8196 line_orientation[i] =
false;
8205 for (
unsigned int i = 4; i < 12; ++i)
8206 if (lines[i]->vertex_index((i + 1) % 2) ==
8207 middle_vertex_index<dim, spacedim>(
8208 hex->face(1 + i / 4)))
8209 line_orientation[i] =
true;
8214 Assert(lines[i]->vertex_index(i % 2) ==
8215 (middle_vertex_index<dim, spacedim>(
8216 hex->face(1 + i / 4))),
8218 line_orientation[i] =
false;
8223 line_orientation[12] =
true;
8257 2>(line_indices[12],
8274 new_quads[0]->set_line_orientation(
8275 0, line_orientation[0]);
8276 new_quads[0]->set_line_orientation(
8277 2, line_orientation[6]);
8278 new_quads[0]->set_line_orientation(
8279 3, line_orientation[10]);
8281 new_quads[1]->set_line_orientation(
8282 1, line_orientation[1]);
8283 new_quads[1]->set_line_orientation(
8284 2, line_orientation[7]);
8285 new_quads[1]->set_line_orientation(
8286 3, line_orientation[11]);
8288 new_quads[2]->set_line_orientation(
8289 0, line_orientation[4]);
8290 new_quads[2]->set_line_orientation(
8291 1, line_orientation[8]);
8292 new_quads[2]->set_line_orientation(
8293 2, line_orientation[2]);
8295 new_quads[3]->set_line_orientation(
8296 0, line_orientation[5]);
8297 new_quads[3]->set_line_orientation(
8298 1, line_orientation[9]);
8299 new_quads[3]->set_line_orientation(
8300 3, line_orientation[3]);
8334 const int quad_indices[20] = {
8335 new_quads[0]->index(),
8336 new_quads[1]->index(),
8337 new_quads[2]->index(),
8338 new_quads[3]->index(),
8340 hex->face(0)->child_index(
8341 child_at_origin[hex->face(0)->refinement_case() -
8342 1][f_fl[0]][f_ro[0]]),
8343 hex->face(0)->child_index(
8345 child_at_origin[hex->face(0)->refinement_case() -
8346 1][f_fl[0]][f_ro[0]]),
8348 hex->face(1)->child_index(
8349 child_at_origin[hex->face(1)->refinement_case() -
8350 1][f_fl[1]][f_ro[1]]),
8351 hex->face(1)->child_index(
8353 child_at_origin[hex->face(1)->refinement_case() -
8354 1][f_fl[1]][f_ro[1]]),
8356 hex->face(2)->isotropic_child_index(
8358 0, f_or[2], f_fl[2], f_ro[2])),
8359 hex->face(2)->isotropic_child_index(
8361 1, f_or[2], f_fl[2], f_ro[2])),
8362 hex->face(2)->isotropic_child_index(
8364 2, f_or[2], f_fl[2], f_ro[2])),
8365 hex->face(2)->isotropic_child_index(
8367 3, f_or[2], f_fl[2], f_ro[2])),
8369 hex->face(3)->isotropic_child_index(
8371 0, f_or[3], f_fl[3], f_ro[3])),
8372 hex->face(3)->isotropic_child_index(
8374 1, f_or[3], f_fl[3], f_ro[3])),
8375 hex->face(3)->isotropic_child_index(
8377 2, f_or[3], f_fl[3], f_ro[3])),
8378 hex->face(3)->isotropic_child_index(
8380 3, f_or[3], f_fl[3], f_ro[3])),
8382 hex->face(4)->child_index(
8383 child_at_origin[hex->face(4)->refinement_case() -
8384 1][f_fl[4]][f_ro[4]]),
8385 hex->face(4)->child_index(
8387 child_at_origin[hex->face(4)->refinement_case() -
8388 1][f_fl[4]][f_ro[4]]),
8390 hex->face(5)->child_index(
8391 child_at_origin[hex->face(5)->refinement_case() -
8392 1][f_fl[5]][f_ro[5]]),
8393 hex->face(5)->child_index(
8395 child_at_origin[hex->face(5)->refinement_case() -
8396 1][f_fl[5]][f_ro[5]])};
8467 1>(middle_vertex_index<dim, spacedim>(
8469 middle_vertex_index<dim, spacedim>(
8539 spacedim>::raw_line_iterator lines[13] = {
8540 hex->face(2)->child(0)->line(
8541 (hex->face(2)->refinement_case() ==
8545 hex->face(3)->child(0)->line(
8546 (hex->face(3)->refinement_case() ==
8550 hex->face(4)->child(0)->line(
8551 (hex->face(4)->refinement_case() ==
8555 hex->face(5)->child(0)->line(
8556 (hex->face(5)->refinement_case() ==
8564 0, f_or[0], f_fl[0], f_ro[0]))
8567 1, f_or[0], f_fl[0], f_ro[0])),
8571 3, f_or[0], f_fl[0], f_ro[0]))
8574 0, f_or[0], f_fl[0], f_ro[0])),
8578 0, f_or[0], f_fl[0], f_ro[0]))
8581 3, f_or[0], f_fl[0], f_ro[0])),
8585 3, f_or[0], f_fl[0], f_ro[0]))
8588 2, f_or[0], f_fl[0], f_ro[0])),
8593 0, f_or[1], f_fl[1], f_ro[1]))
8596 1, f_or[1], f_fl[1], f_ro[1])),
8600 3, f_or[1], f_fl[1], f_ro[1]))
8603 0, f_or[1], f_fl[1], f_ro[1])),
8607 0, f_or[1], f_fl[1], f_ro[1]))
8610 3, f_or[1], f_fl[1], f_ro[1])),
8614 3, f_or[1], f_fl[1], f_ro[1]))
8617 2, f_or[1], f_fl[1], f_ro[1])),
8622 unsigned int line_indices[13];
8624 for (
unsigned int i = 0; i < 13; ++i)
8625 line_indices[i] = lines[i]->index();
8632 bool line_orientation[13];
8636 const unsigned int middle_vertices[4] = {
8637 hex->line(8)->child(0)->vertex_index(1),
8638 hex->line(10)->child(0)->vertex_index(1),
8639 hex->line(0)->child(0)->vertex_index(1),
8640 hex->line(4)->child(0)->vertex_index(1),
8645 for (
unsigned int i = 0; i < 4; ++i)
8646 if (lines[i]->vertex_index(0) == middle_vertices[i])
8647 line_orientation[i] =
true;
8651 Assert(lines[i]->vertex_index(1) ==
8654 line_orientation[i] =
false;
8663 for (
unsigned int i = 4; i < 12; ++i)
8664 if (lines[i]->vertex_index((i + 1) % 2) ==
8665 middle_vertex_index<dim, spacedim>(
8666 hex->face(i / 4 - 1)))
8667 line_orientation[i] =
true;
8672 Assert(lines[i]->vertex_index(i % 2) ==
8673 (middle_vertex_index<dim, spacedim>(
8674 hex->face(i / 4 - 1))),
8676 line_orientation[i] =
false;
8681 line_orientation[12] =
true;
8721 2>(line_indices[12],
8726 new_quads[0]->set_line_orientation(
8727 0, line_orientation[6]);
8728 new_quads[0]->set_line_orientation(
8729 1, line_orientation[10]);
8730 new_quads[0]->set_line_orientation(
8731 2, line_orientation[0]);
8733 new_quads[1]->set_line_orientation(
8734 0, line_orientation[7]);
8735 new_quads[1]->set_line_orientation(
8736 1, line_orientation[11]);
8737 new_quads[1]->set_line_orientation(
8738 3, line_orientation[1]);
8740 new_quads[2]->set_line_orientation(
8741 0, line_orientation[2]);
8742 new_quads[2]->set_line_orientation(
8743 2, line_orientation[4]);
8744 new_quads[2]->set_line_orientation(
8745 3, line_orientation[8]);
8747 new_quads[3]->set_line_orientation(
8748 1, line_orientation[3]);
8749 new_quads[3]->set_line_orientation(
8750 2, line_orientation[5]);
8751 new_quads[3]->set_line_orientation(
8752 3, line_orientation[9]);
8786 const int quad_indices[20] = {
8787 new_quads[0]->index(),
8788 new_quads[1]->index(),
8789 new_quads[2]->index(),
8790 new_quads[3]->index(),
8792 hex->face(0)->isotropic_child_index(
8794 0, f_or[0], f_fl[0], f_ro[0])),
8795 hex->face(0)->isotropic_child_index(
8797 1, f_or[0], f_fl[0], f_ro[0])),
8798 hex->face(0)->isotropic_child_index(
8800 2, f_or[0], f_fl[0], f_ro[0])),
8801 hex->face(0)->isotropic_child_index(
8803 3, f_or[0], f_fl[0], f_ro[0])),
8805 hex->face(1)->isotropic_child_index(
8807 0, f_or[1], f_fl[1], f_ro[1])),
8808 hex->face(1)->isotropic_child_index(
8810 1, f_or[1], f_fl[1], f_ro[1])),
8811 hex->face(1)->isotropic_child_index(
8813 2, f_or[1], f_fl[1], f_ro[1])),
8814 hex->face(1)->isotropic_child_index(
8816 3, f_or[1], f_fl[1], f_ro[1])),
8818 hex->face(2)->child_index(
8819 child_at_origin[hex->face(2)->refinement_case() -
8820 1][f_fl[2]][f_ro[2]]),
8821 hex->face(2)->child_index(
8823 child_at_origin[hex->face(2)->refinement_case() -
8824 1][f_fl[2]][f_ro[2]]),
8826 hex->face(3)->child_index(
8827 child_at_origin[hex->face(3)->refinement_case() -
8828 1][f_fl[3]][f_ro[3]]),
8829 hex->face(3)->child_index(
8831 child_at_origin[hex->face(3)->refinement_case() -
8832 1][f_fl[3]][f_ro[3]]),
8834 hex->face(4)->child_index(
8835 child_at_origin[hex->face(4)->refinement_case() -
8836 1][f_fl[4]][f_ro[4]]),
8837 hex->face(4)->child_index(
8839 child_at_origin[hex->face(4)->refinement_case() -
8840 1][f_fl[4]][f_ro[4]]),
8842 hex->face(5)->child_index(
8843 child_at_origin[hex->face(5)->refinement_case() -
8844 1][f_fl[5]][f_ro[5]]),
8845 hex->face(5)->child_index(
8847 child_at_origin[hex->face(5)->refinement_case() -
8848 1][f_fl[5]][f_ro[5]])};
8913 ++next_unused_vertex;
8917 "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."));
8928 hex->center(
true,
true);
8950 middle_vertex_index<dim, spacedim>(hex->face(0)),
8951 middle_vertex_index<dim, spacedim>(hex->face(1)),
8952 middle_vertex_index<dim, spacedim>(hex->face(2)),
8953 middle_vertex_index<dim, spacedim>(hex->face(3)),
8954 middle_vertex_index<dim, spacedim>(hex->face(4)),
8955 middle_vertex_index<dim, spacedim>(hex->face(5)),
8956 next_unused_vertex};
9044 spacedim>::raw_line_iterator lines[30] = {
9048 0, f_or[0], f_fl[0], f_ro[0]))
9051 1, f_or[0], f_fl[0], f_ro[0])),
9055 3, f_or[0], f_fl[0], f_ro[0]))
9058 0, f_or[0], f_fl[0], f_ro[0])),
9062 0, f_or[0], f_fl[0], f_ro[0]))
9065 3, f_or[0], f_fl[0], f_ro[0])),
9069 3, f_or[0], f_fl[0], f_ro[0]))
9072 2, f_or[0], f_fl[0], f_ro[0])),
9077 0, f_or[1], f_fl[1], f_ro[1]))
9080 1, f_or[1], f_fl[1], f_ro[1])),
9084 3, f_or[1], f_fl[1], f_ro[1]))
9087 0, f_or[1], f_fl[1], f_ro[1])),
9091 0, f_or[1], f_fl[1], f_ro[1]))
9094 3, f_or[1], f_fl[1], f_ro[1])),
9098 3, f_or[1], f_fl[1], f_ro[1]))
9101 2, f_or[1], f_fl[1], f_ro[1])),
9106 0, f_or[2], f_fl[2], f_ro[2]))
9109 1, f_or[2], f_fl[2], f_ro[2])),
9113 3, f_or[2], f_fl[2], f_ro[2]))
9116 0, f_or[2], f_fl[2], f_ro[2])),
9120 0, f_or[2], f_fl[2], f_ro[2]))
9123 3, f_or[2], f_fl[2], f_ro[2])),
9127 3, f_or[2], f_fl[2], f_ro[2]))
9130 2, f_or[2], f_fl[2], f_ro[2])),
9135 0, f_or[3], f_fl[3], f_ro[3]))
9138 1, f_or[3], f_fl[3], f_ro[3])),
9142 3, f_or[3], f_fl[3], f_ro[3]))
9145 0, f_or[3], f_fl[3], f_ro[3])),
9149 0, f_or[3], f_fl[3], f_ro[3]))
9152 3, f_or[3], f_fl[3], f_ro[3])),
9156 3, f_or[3], f_fl[3], f_ro[3]))
9159 2, f_or[3], f_fl[3], f_ro[3])),
9164 0, f_or[4], f_fl[4], f_ro[4]))
9167 1, f_or[4], f_fl[4], f_ro[4])),
9171 3, f_or[4], f_fl[4], f_ro[4]))
9174 0, f_or[4], f_fl[4], f_ro[4])),
9178 0, f_or[4], f_fl[4], f_ro[4]))
9181 3, f_or[4], f_fl[4], f_ro[4])),
9185 3, f_or[4], f_fl[4], f_ro[4]))
9188 2, f_or[4], f_fl[4], f_ro[4])),
9193 0, f_or[5], f_fl[5], f_ro[5]))
9196 1, f_or[5], f_fl[5], f_ro[5])),
9200 3, f_or[5], f_fl[5], f_ro[5]))
9203 0, f_or[5], f_fl[5], f_ro[5])),
9207 0, f_or[5], f_fl[5], f_ro[5]))
9210 3, f_or[5], f_fl[5], f_ro[5])),
9214 3, f_or[5], f_fl[5], f_ro[5]))
9217 2, f_or[5], f_fl[5], f_ro[5])),
9227 unsigned int line_indices[30];
9228 for (
unsigned int i = 0; i < 30; ++i)
9229 line_indices[i] = lines[i]->index();
9236 bool line_orientation[30];
9244 for (
unsigned int i = 0; i < 24; ++i)
9245 if (lines[i]->vertex_index((i + 1) % 2) ==
9247 line_orientation[i] =
true;
9252 Assert(lines[i]->vertex_index(i % 2) ==
9255 line_orientation[i] =
false;
9260 for (
unsigned int i = 24; i < 30; ++i)
9261 line_orientation[i] =
true;
9295 2>(line_indices[10],
9301 2>(line_indices[28],
9307 2>(line_indices[11],
9313 2>(line_indices[29],
9319 2>(line_indices[18],
9325 2>(line_indices[26],
9331 2>(line_indices[19],
9337 2>(line_indices[27],
9349 2>(line_indices[24],
9361 2>(line_indices[25],
9369 new_quads[0]->set_line_orientation(
9370 0, line_orientation[10]);
9371 new_quads[0]->set_line_orientation(
9372 2, line_orientation[16]);
9374 new_quads[1]->set_line_orientation(
9375 1, line_orientation[14]);
9376 new_quads[1]->set_line_orientation(
9377 2, line_orientation[17]);
9379 new_quads[2]->set_line_orientation(
9380 0, line_orientation[11]);
9381 new_quads[2]->set_line_orientation(
9382 3, line_orientation[20]);
9384 new_quads[3]->set_line_orientation(
9385 1, line_orientation[15]);
9386 new_quads[3]->set_line_orientation(
9387 3, line_orientation[21]);
9389 new_quads[4]->set_line_orientation(
9390 0, line_orientation[18]);
9391 new_quads[4]->set_line_orientation(
9392 2, line_orientation[0]);
9394 new_quads[5]->set_line_orientation(
9395 1, line_orientation[22]);
9396 new_quads[5]->set_line_orientation(
9397 2, line_orientation[1]);
9399 new_quads[6]->set_line_orientation(
9400 0, line_orientation[19]);
9401 new_quads[6]->set_line_orientation(
9402 3, line_orientation[4]);
9404 new_quads[7]->set_line_orientation(
9405 1, line_orientation[23]);
9406 new_quads[7]->set_line_orientation(
9407 3, line_orientation[5]);
9409 new_quads[8]->set_line_orientation(
9410 0, line_orientation[2]);
9411 new_quads[8]->set_line_orientation(
9412 2, line_orientation[8]);
9414 new_quads[9]->set_line_orientation(
9415 1, line_orientation[6]);
9416 new_quads[9]->set_line_orientation(
9417 2, line_orientation[9]);
9419 new_quads[10]->set_line_orientation(
9420 0, line_orientation[3]);
9421 new_quads[10]->set_line_orientation(
9422 3, line_orientation[12]);
9424 new_quads[11]->set_line_orientation(
9425 1, line_orientation[7]);
9426 new_quads[11]->set_line_orientation(
9427 3, line_orientation[13]);
9468 const int quad_indices[36] = {
9469 new_quads[0]->index(),
9470 new_quads[1]->index(),
9471 new_quads[2]->index(),
9472 new_quads[3]->index(),
9473 new_quads[4]->index(),
9474 new_quads[5]->index(),
9475 new_quads[6]->index(),
9476 new_quads[7]->index(),
9477 new_quads[8]->index(),
9478 new_quads[9]->index(),
9479 new_quads[10]->index(),
9480 new_quads[11]->index(),
9482 hex->face(0)->isotropic_child_index(
9484 0, f_or[0], f_fl[0], f_ro[0])),
9485 hex->face(0)->isotropic_child_index(
9487 1, f_or[0], f_fl[0], f_ro[0])),
9488 hex->face(0)->isotropic_child_index(
9490 2, f_or[0], f_fl[0], f_ro[0])),
9491 hex->face(0)->isotropic_child_index(
9493 3, f_or[0], f_fl[0], f_ro[0])),
9495 hex->face(1)->isotropic_child_index(
9497 0, f_or[1], f_fl[1], f_ro[1])),
9498 hex->face(1)->isotropic_child_index(
9500 1, f_or[1], f_fl[1], f_ro[1])),
9501 hex->face(1)->isotropic_child_index(
9503 2, f_or[1], f_fl[1], f_ro[1])),
9504 hex->face(1)->isotropic_child_index(
9506 3, f_or[1], f_fl[1], f_ro[1])),
9508 hex->face(2)->isotropic_child_index(
9510 0, f_or[2], f_fl[2], f_ro[2])),
9511 hex->face(2)->isotropic_child_index(
9513 1, f_or[2], f_fl[2], f_ro[2])),
9514 hex->face(2)->isotropic_child_index(
9516 2, f_or[2], f_fl[2], f_ro[2])),
9517 hex->face(2)->isotropic_child_index(
9519 3, f_or[2], f_fl[2], f_ro[2])),
9521 hex->face(3)->isotropic_child_index(
9523 0, f_or[3], f_fl[3], f_ro[3])),
9524 hex->face(3)->isotropic_child_index(
9526 1, f_or[3], f_fl[3], f_ro[3])),
9527 hex->face(3)->isotropic_child_index(
9529 2, f_or[3], f_fl[3], f_ro[3])),
9530 hex->face(3)->isotropic_child_index(
9532 3, f_or[3], f_fl[3], f_ro[3])),
9534 hex->face(4)->isotropic_child_index(
9536 0, f_or[4], f_fl[4], f_ro[4])),
9537 hex->face(4)->isotropic_child_index(
9539 1, f_or[4], f_fl[4], f_ro[4])),
9540 hex->face(4)->isotropic_child_index(
9542 2, f_or[4], f_fl[4], f_ro[4])),
9543 hex->face(4)->isotropic_child_index(
9545 3, f_or[4], f_fl[4], f_ro[4])),
9547 hex->face(5)->isotropic_child_index(
9549 0, f_or[5], f_fl[5], f_ro[5])),
9550 hex->face(5)->isotropic_child_index(
9552 1, f_or[5], f_fl[5], f_ro[5])),
9553 hex->face(5)->isotropic_child_index(
9555 2, f_or[5], f_fl[5], f_ro[5])),
9556 hex->face(5)->isotropic_child_index(
9558 3, f_or[5], f_fl[5], f_ro[5]))};
9563 3>(quad_indices[12],
9579 3>(quad_indices[13],
9597 3>(quad_indices[14],
9613 3>(quad_indices[15],
9660 for (
unsigned int s = 0;
9667 const unsigned int current_child =
9676 ref_case, f, f_or[f], f_fl[f], f_ro[f]));
9677 new_hexes[current_child]->set_face_orientation(f,
9679 new_hexes[current_child]->set_face_flip(f, f_fl[f]);
9680 new_hexes[current_child]->set_face_rotation(f, f_ro[f]);
9685 if ((check_for_distorted_cells ==
true) &&
9686 has_distorted_children(
9688 std::integral_constant<int, dim>(),
9689 std::integral_constant<int, spacedim>()))
9709 return cells_with_distorted_children;
9725 template <
int spacedim>
9730 template <
int dim,
int spacedim>
9744 if (cell->at_boundary() && cell->refine_flag_set() &&
9745 cell->refine_flag_set() !=
9752 for (
const unsigned int face_no :
9754 if (cell->face(face_no)->at_boundary())
9765 spacedim>::face_iterator
9766 face = cell->face(face_no);
9774 .transform_real_to_unit_cell(cell, new_bound);
9785 std::fabs(new_unit[face_no / 2] - face_no % 2);
9790 const double allowed = 0.25;
9793 cell->flag_for_face_refinement(face_no);
9812 template <
int dim,
int spacedim>
9817 ExcMessage(
"Wrong function called -- there should "
9818 "be a specialization."));
9822 template <
int spacedim>
9826 const unsigned int dim = 3;
9838 bool mesh_changed =
false;
9842 mesh_changed =
false;
9856 if (cell->refine_flag_set())
9858 for (
unsigned int line = 0;
9859 line < GeometryInfo<dim>::lines_per_cell;
9862 cell->refine_flag_set(), line) ==
9866 cell->line(line)->set_user_flag();
9868 else if (cell->has_children() &&
9869 !cell->child(0)->coarsen_flag_set())
9871 for (
unsigned int line = 0;
9872 line < GeometryInfo<dim>::lines_per_cell;
9875 cell->refinement_case(), line) ==
9879 cell->line(line)->set_user_flag();
9881 else if (cell->has_children() &&
9882 cell->child(0)->coarsen_flag_set())
9883 cell->set_user_flag();
9894 for (
unsigned int line = 0;
9895 line < GeometryInfo<dim>::lines_per_cell;
9898 if (cell->line(line)->has_children())
9907 bool offending_line_found =
false;
9909 for (
unsigned int c = 0; c < 2; ++c)
9911 Assert(cell->line(line)->child(c)->has_children() ==
9915 if (cell->line(line)->child(c)->user_flag_set() &&
9917 cell->refine_flag_set(), line) ==
9921 cell->clear_coarsen_flag();
9928 allow_anisotropic_smoothing)
9929 cell->flag_for_line_refinement(line);
9931 cell->set_refine_flag();
9933 for (
unsigned int l = 0;
9934 l < GeometryInfo<dim>::lines_per_cell;
9937 cell->refine_flag_set(), line) ==
9940 cell->line(
l)->set_user_flag();
9943 offending_line_found =
true;
9949 for (
unsigned int l = 0;
9950 l < GeometryInfo<dim>::lines_per_cell;
9952 if (!cell->line(
l)->has_children() &&
9954 cell->refine_flag_set(),
l) !=
9956 cell->line(
l)->set_user_flag();
9962 if (offending_line_found)
9964 mesh_changed =
true;
9986 if (cell->user_flag_set())
9987 for (
unsigned int line = 0;
9988 line < GeometryInfo<dim>::lines_per_cell;
9990 if (cell->line(line)->has_children() &&
9991 (cell->line(line)->child(0)->user_flag_set() ||
9992 cell->line(line)->child(1)->user_flag_set()))
9994 for (
unsigned int c = 0; c < cell->n_children(); ++c)
9995 cell->child(c)->clear_coarsen_flag();
9996 cell->clear_user_flag();
9997 for (
unsigned int l = 0;
9998 l < GeometryInfo<dim>::lines_per_cell;
10001 cell->refinement_case(),
l) ==
10005 cell->line(
l)->set_user_flag();
10006 mesh_changed =
true;
10011 while (mesh_changed ==
true);
10022 template <
int dim,
int spacedim>
10042 const unsigned int n_subfaces =
10045 if (n_subfaces == 0 || cell->at_boundary(n))
10047 for (
unsigned int c = 0; c < n_subfaces; ++c)
10050 child = cell->child(
10054 child_neighbor = child->neighbor(n);
10055 if (!child->neighbor_is_coarser(n))
10068 if ((child_neighbor->has_children() &&
10069 !child_neighbor->user_flag_set()) ||
10078 child_neighbor->refine_flag_set())
10088 template <
int dim,
int spacedim>
10093 return flat_manifold;
10100 template <
int dim,
int spacedim>
10105 template <
int dim,
int spacedim>
10108 const bool check_for_distorted_cells)
10109 : smooth_grid(smooth_grid)
10110 , anisotropic_refinement(false)
10111 , check_for_distorted_cells(check_for_distorted_cells)
10116 std_cxx14::make_unique<std::map<unsigned int, types::boundary_id>>();
10118 std_cxx14::make_unique<std::map<unsigned int, types::manifold_id>>();
10130 template <
int dim,
int spacedim>
10134 , smooth_grid(tria.smooth_grid)
10135 , periodic_face_pairs_level_0(std::move(tria.periodic_face_pairs_level_0))
10136 , periodic_face_map(std::move(tria.periodic_face_map))
10137 , levels(std::move(tria.levels))
10138 , faces(std::move(tria.faces))
10139 ,
vertices(std::move(tria.vertices))
10140 , vertices_used(std::move(tria.vertices_used))
10141 , manifold(std::move(tria.manifold))
10142 , anisotropic_refinement(tria.anisotropic_refinement)
10143 , check_for_distorted_cells(tria.check_for_distorted_cells)
10144 , number_cache(std::move(tria.number_cache))
10145 , vertex_to_boundary_id_map_1d(std::move(tria.vertex_to_boundary_id_map_1d))
10146 , vertex_to_manifold_id_map_1d(std::move(tria.vertex_to_manifold_id_map_1d))
10152 template <
int dim,
int spacedim>
10159 smooth_grid = tria.smooth_grid;
10160 periodic_face_pairs_level_0 = std::move(tria.periodic_face_pairs_level_0);
10161 periodic_face_map = std::move(tria.periodic_face_map);
10162 levels = std::move(tria.levels);
10163 faces = std::move(tria.faces);
10164 vertices = std::move(tria.vertices);
10165 vertices_used = std::move(tria.vertices_used);
10166 manifold = std::move(tria.manifold);
10167 anisotropic_refinement = tria.anisotropic_refinement;
10168 number_cache = tria.number_cache;
10169 vertex_to_boundary_id_map_1d = std::move(tria.vertex_to_boundary_id_map_1d);
10170 vertex_to_manifold_id_map_1d = std::move(tria.vertex_to_manifold_id_map_1d);
10179 template <
int dim,
int spacedim>
10196 AssertNothrow((dim == 1) || (vertex_to_boundary_id_map_1d ==
nullptr),
10201 AssertNothrow((dim == 1) || (vertex_to_manifold_id_map_1d ==
nullptr),
10207 template <
int dim,
int spacedim>
10215 clear_despite_subscriptions();
10216 periodic_face_pairs_level_0.clear();
10217 periodic_face_map.clear();
10222 template <
int dim,
int spacedim>
10228 ExcTriangulationNotEmpty(
vertices.size(), levels.size()));
10229 smooth_grid = mesh_smoothing;
10234 template <
int dim,
int spacedim>
10238 return smooth_grid;
10243 template <
int dim,
int spacedim>
10251 manifold[m_number] = manifold_object.
clone();
10256 template <
int dim,
int spacedim>
10263 manifold.erase(m_number);
10267 template <
int dim,
int spacedim>
10275 template <
int dim,
int spacedim>
10283 "Error: set_all_manifold_ids() can not be called on an empty Triangulation."));
10286 cell = this->begin_active(),
10287 endc = this->
end();
10289 for (; cell != endc; ++cell)
10290 cell->set_all_manifold_ids(m_number);
10294 template <
int dim,
int spacedim>
10302 "Error: set_all_manifold_ids_on_boundary() can not be called on an empty Triangulation."));
10305 cell = this->begin_active(),
10306 endc = this->
end();
10308 for (; cell != endc; ++cell)
10310 if (cell->face(f)->at_boundary())
10311 cell->face(f)->set_all_manifold_ids(m_number);
10315 template <
int dim,
int spacedim>
10324 "Error: set_all_manifold_ids_on_boundary() can not be called on an empty Triangulation."));
10326 bool boundary_found =
false;
10328 cell = this->begin_active(),
10329 endc = this->
end();
10331 for (; cell != endc; ++cell)
10335 if (cell->face(f)->at_boundary() &&
10336 cell->face(f)->boundary_id() == b_id)
10338 boundary_found =
true;
10339 cell->face(f)->set_manifold_id(m_number);
10344 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++
e)
10345 if (cell->line(
e)->at_boundary() &&
10346 cell->line(
e)->boundary_id() == b_id)
10348 boundary_found =
true;
10349 cell->line(
e)->set_manifold_id(m_number);
10353 (void)boundary_found;
10354 Assert(boundary_found, ExcBoundaryIdNotFound(b_id));
10359 template <
int dim,
int spacedim>
10366 const auto it = manifold.find(m_number);
10368 if (it != manifold.end())
10371 return *(it->second);
10376 return internal::TriangulationImplementation::
10377 get_default_flat_manifold<dim, spacedim>();
10382 template <
int dim,
int spacedim>
10383 std::vector<types::boundary_id>
10390 std::vector<types::boundary_id> boundary_ids;
10391 for (std::map<unsigned int, types::boundary_id>::const_iterator p =
10392 vertex_to_boundary_id_map_1d->begin();
10393 p != vertex_to_boundary_id_map_1d->end();
10395 boundary_ids.push_back(p->second);
10397 return boundary_ids;
10401 std::set<types::boundary_id> b_ids;
10402 for (
auto cell : active_cell_iterators())
10403 if (cell->is_locally_owned())
10405 if (cell->at_boundary(face))
10406 b_ids.insert(cell->face(face)->boundary_id());
10407 std::vector<types::boundary_id> boundary_ids(b_ids.begin(), b_ids.end());
10408 return boundary_ids;
10414 template <
int dim,
int spacedim>
10415 std::vector<types::manifold_id>
10418 std::set<types::manifold_id> m_ids;
10419 for (
auto cell : active_cell_iterators())
10420 if (cell->is_locally_owned())
10422 m_ids.insert(cell->manifold_id());
10425 if (cell->at_boundary(face))
10426 m_ids.insert(cell->face(face)->manifold_id());
10428 std::vector<types::manifold_id> manifold_indicators(m_ids.begin(),
10430 return manifold_indicators;
10436 template <
int dim,
int spacedim>
10441 Assert((
vertices.size() == 0) && (levels.size() == 0) && (faces ==
nullptr),
10442 ExcTriangulationNotEmpty(
vertices.size(), levels.size()));
10444 (dim == 1 || other_tria.
faces !=
nullptr),
10446 "When calling Triangulation::copy_triangulation(), "
10447 "the target triangulation must be empty but the source "
10448 "triangulation (the argument to this function) must contain "
10449 "something. Here, it seems like the source does not "
10450 "contain anything at all."));
10460 faces = std_cxx14::make_unique<
10463 auto bdry_iterator = other_tria.
manifold.begin();
10464 for (; bdry_iterator != other_tria.
manifold.end(); ++bdry_iterator)
10465 manifold[bdry_iterator->first] = bdry_iterator->second->clone();
10468 levels.reserve(other_tria.
levels.size());
10470 levels.push_back(std_cxx14::make_unique<
10478 vertex_to_boundary_id_map_1d =
10479 std_cxx14::make_unique<std::map<unsigned int, types::boundary_id>>(
10482 vertex_to_manifold_id_map_1d =
10483 std_cxx14::make_unique<std::map<unsigned int, types::manifold_id>>(
10488 other_tria.
signals.copy(*
this);
10499 template <
int dim,
int spacedim>
10506 std::vector<CellData<dim>> reordered_cells(cells);
10510 reorder_compatibility(reordered_cells, reordered_subcelldata);
10519 template <
int dim,
int spacedim>
10526 Assert((
vertices.size() == 0) && (levels.size() == 0) && (faces ==
nullptr),
10527 ExcTriangulationNotEmpty(
vertices.size(), levels.size()));
10542 clear_despite_subscriptions();
10549 *
this, levels.size(), number_cache);
10550 reset_active_cell_indices();
10555 if (check_for_distorted_cells ==
true)
10557 DistortedCellList distorted_cells = collect_distorted_coarse_cells(*
this);
10561 AssertThrow(distorted_cells.distorted_cells.size() == 0, distorted_cells);
10594 if (dim < spacedim)
10602 bool values[][2] = {{
false,
true}, {
true,
false}};
10605 correct(i, j) = (values[i][j]);
10610 bool values[][4] = {{
false,
true,
true,
false},
10611 {
true,
false,
false,
true},
10612 {
true,
false,
false,
true},
10613 {
false,
true,
true,
false}};
10616 correct(i, j) = (values[i][j]);
10624 std::list<active_cell_iterator> this_round, next_round;
10627 this_round.push_back(begin_active());
10628 begin_active()->set_direction_flag(
true);
10629 begin_active()->set_user_flag();
10631 while (this_round.size() > 0)
10633 for (
typename std::list<active_cell_iterator>::iterator cell =
10634 this_round.begin();
10635 cell != this_round.end();
10640 if (!((*cell)->face(i)->at_boundary()))
10642 neighbor = (*cell)->neighbor(i);
10644 unsigned int cf = (*cell)->face_index(i);
10645 unsigned int j = 0;
10646 while (neighbor->face_index(j) != cf)
10654 if (neighbor->user_flag_set())
10658 Assert(!(correct(i, j) ^
10659 (neighbor->direction_flag() ==
10660 (*cell)->direction_flag())),
10661 ExcNonOrientableTriangulation());
10665 next_round.push_back(neighbor);
10666 neighbor->set_user_flag();
10667 if ((correct(i, j) ^ (neighbor->direction_flag() ==
10668 (*cell)->direction_flag())))
10669 neighbor->set_direction_flag(
10670 !neighbor->direction_flag());
10680 if (next_round.size() == 0)
10681 for (
const auto &cell : this->active_cell_iterators())
10682 if (cell->user_flag_set() ==
false)
10684 next_round.push_back(cell);
10685 cell->set_direction_flag(
true);
10686 cell->set_user_flag();
10690 this_round = next_round;
10691 next_round.clear();
10701 template <
int dim,
int spacedim>
10712 auto cell_infos = construction_data.
cell_infos;
10715 for (
auto &cell_info : cell_infos)
10721 const CellId a_id(a.id);
10722 const CellId b_id(b.id);
10724 const auto a_coarse_cell_index =
10725 this->coarse_cell_id_to_coarse_cell_index(a_id.get_coarse_cell_id());
10726 const auto b_coarse_cell_index =
10727 this->coarse_cell_id_to_coarse_cell_index(b_id.get_coarse_cell_id());
10734 if (a_coarse_cell_index != b_coarse_cell_index)
10735 return a_coarse_cell_index < b_coarse_cell_index;
10737 return a_id < b_id;
10747 auto cell_info = cell_infos[
level].begin();
10748 for (; cell_info != cell_infos[
level].end(); ++cell_info)
10750 while (cell_info->id != cell->id().template to_binary<dim>())
10753 for (
unsigned int quad = 0;
10754 quad < GeometryInfo<dim>::quads_per_cell;
10756 cell->quad(quad)->set_manifold_id(
10757 cell_info->manifold_quad_ids[quad]);
10760 for (
unsigned int line = 0;
10761 line < GeometryInfo<dim>::lines_per_cell;
10763 cell->line(line)->set_manifold_id(
10764 cell_info->manifold_line_ids[line]);
10766 cell->set_manifold_id(cell_info->manifold_id);
10771 if (
level + 1 != cell_infos.size())
10776 auto fine_cell_info = cell_infos[
level + 1].begin();
10779 for (; fine_cell_info != cell_infos[
level + 1].end();
10784 !coarse_cell->id().is_parent_of(
CellId(fine_cell_info->id)))
10788 coarse_cell->set_refine_flag();
10793 spacedim>::execute_coarsening_and_refinement();
10801 auto cell_info = cell_infos[
level].begin();
10802 for (; cell_info != cell_infos[
level].end(); ++cell_info)
10805 while (cell_info->id != cell->id().template to_binary<dim>())
10809 for (
auto pair : cell_info->boundary_ids)
10811 Assert(cell->at_boundary(pair.first),
10812 ExcMessage(
"Cell face is not on the boundary!"));
10813 cell->face(pair.first)->set_boundary_id(pair.second);
10820 template <
int dim,
int spacedim>
10825 ExcMessage(
"Only works for dim == spacedim-1"));
10826 for (
const auto &cell : this->active_cell_iterators())
10827 cell->set_direction_flag(!cell->direction_flag());
10832 template <
int dim,
int spacedim>
10837 ExcMessage(
"Error: An empty Triangulation can not be refined."));
10840 for (; cell != endc; ++cell)
10842 cell->clear_coarsen_flag();
10843 cell->set_refine_flag();
10849 template <
int dim,
int spacedim>
10853 for (
unsigned int i = 0; i < times; ++i)
10855 set_all_refine_flags();
10856 execute_coarsening_and_refinement();
10866 template <
int dim,
int spacedim>
10871 std::vector<bool>::iterator i = v.begin();
10873 for (; cell != endc; ++cell)
10874 for (
unsigned int j = 0; j < dim; ++j, ++i)
10875 if (cell->refine_flag_set() & (1 << j))
10883 template <
int dim,
int spacedim>
10887 std::vector<bool> v;
10888 save_refine_flags(v);
10897 template <
int dim,
int spacedim>
10901 std::vector<bool> v;
10903 load_refine_flags(v);
10908 template <
int dim,
int spacedim>
10915 std::vector<bool>::const_iterator i = v.begin();
10916 for (; cell != endc; ++cell)
10918 unsigned int ref_case = 0;
10920 for (
unsigned int j = 0; j < dim; ++j, ++i)
10922 ref_case += 1 << j;
10924 ExcGridReadError());
10928 cell->clear_refine_flag();
10936 template <
int dim,
int spacedim>
10941 std::vector<bool>::iterator i = v.begin();
10943 for (; cell != endc; ++cell, ++i)
10944 *i = cell->coarsen_flag_set();
10951 template <
int dim,
int spacedim>
10955 std::vector<bool> v;
10956 save_coarsen_flags(v);
10965 template <
int dim,
int spacedim>
10969 std::vector<bool> v;
10974 load_coarsen_flags(v);
10979 template <
int dim,
int spacedim>
10986 std::vector<bool>::const_iterator i = v.begin();
10987 for (; cell != endc; ++cell, ++i)
10989 cell->set_coarsen_flag();
10991 cell->clear_coarsen_flag();
10997 template <
int dim,
int spacedim>
11001 return anisotropic_refinement;
11020 levels[
level]->cells.clear_user_data();
11034 faces->
lines.clear_user_data();
11041 faces->
lines.clear_user_data();
11047 template <
int dim,
int spacedim>
11052 ::clear_user_data(levels);
11053 ::clear_user_data(faces.get());
11060 void clear_user_flags_line(
11066 for (
const auto &
level : levels)
11067 level->cells.clear_user_flags();
11072 clear_user_flags_line(
11077 faces->lines.clear_user_flags();
11082 template <
int dim,
int spacedim>
11086 ::clear_user_flags_line(levels, faces.get());
11093 void clear_user_flags_quad(
11101 void clear_user_flags_quad(
11107 for (
const auto &
level : levels)
11108 level->cells.clear_user_flags();
11113 clear_user_flags_quad(
11118 faces->quads.clear_user_flags();
11123 template <
int dim,
int spacedim>
11127 ::clear_user_flags_quad(levels, faces.get());
11134 void clear_user_flags_hex(
11143 void clear_user_flags_hex(
11151 void clear_user_flags_hex(
11157 for (
const auto &
level : levels)
11158 level->cells.clear_user_flags();
11163 template <
int dim,
int spacedim>
11167 ::clear_user_flags_hex(levels, faces.get());
11172 template <
int dim,
int spacedim>
11176 clear_user_flags_line();
11177 clear_user_flags_quad();
11178 clear_user_flags_hex();
11183 template <
int dim,
int spacedim>
11187 save_user_flags_line(out);
11190 save_user_flags_quad(out);
11193 save_user_flags_hex(out);
11201 template <
int dim,
int spacedim>
11209 std::vector<bool> tmp;
11211 save_user_flags_line(tmp);
11212 v.insert(v.end(), tmp.begin(), tmp.end());
11216 save_user_flags_quad(tmp);
11217 v.insert(v.end(), tmp.begin(), tmp.end());
11222 save_user_flags_hex(tmp);
11223 v.insert(v.end(), tmp.begin(), tmp.end());
11232 template <
int dim,
int spacedim>
11236 load_user_flags_line(in);
11239 load_user_flags_quad(in);
11242 load_user_flags_hex(in);
11250 template <
int dim,
int spacedim>
11255 std::vector<bool> tmp;
11259 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
11261 load_user_flags_line(tmp);
11266 tmp.insert(tmp.end(),
11267 v.begin() + n_lines(),
11268 v.begin() + n_lines() + n_quads());
11269 load_user_flags_quad(tmp);
11275 tmp.insert(tmp.end(),
11276 v.begin() + n_lines() + n_quads(),
11277 v.begin() + n_lines() + n_quads() + n_hexs());
11278 load_user_flags_hex(tmp);
11287 template <
int dim,
int spacedim>
11291 v.resize(n_lines(),
false);
11292 std::vector<bool>::iterator i = v.begin();
11294 for (; line != endl; ++line, ++i)
11295 *i = line->user_flag_set();
11302 template <
int dim,
int spacedim>
11306 std::vector<bool> v;
11307 save_user_flags_line(v);
11316 template <
int dim,
int spacedim>
11320 std::vector<bool> v;
11325 load_user_flags_line(v);
11330 template <
int dim,
int spacedim>
11334 Assert(v.size() == n_lines(), ExcGridReadError());
11337 std::vector<bool>::const_iterator i = v.begin();
11338 for (; line != endl; ++line, ++i)
11340 line->set_user_flag();
11342 line->clear_user_flag();
11350 template <
typename Iterator>
11354 return i->user_flag_set();
11359 template <
int structdim,
int dim,
int spacedim>
11369 template <
typename Iterator>
11373 i->set_user_flag();
11378 template <
int structdim,
int dim,
int spacedim>
11387 template <
typename Iterator>
11389 clear_user_flag(
const Iterator &i)
11391 i->clear_user_flag();
11396 template <
int structdim,
int dim,
int spacedim>
11406 template <
int dim,
int spacedim>
11410 v.resize(n_quads(),
false);
11414 std::vector<bool>::iterator i = v.begin();
11416 for (; quad != endq; ++quad, ++i)
11417 *i = get_user_flag(quad);
11425 template <
int dim,
int spacedim>
11429 std::vector<bool> v;
11430 save_user_flags_quad(v);
11439 template <
int dim,
int spacedim>
11443 std::vector<bool> v;
11448 load_user_flags_quad(v);
11453 template <
int dim,
int spacedim>
11457 Assert(v.size() == n_quads(), ExcGridReadError());
11462 std::vector<bool>::const_iterator i = v.begin();
11463 for (; quad != endq; ++quad, ++i)
11465 set_user_flag(quad);
11467 clear_user_flag(quad);
11475 template <
int dim,
int spacedim>
11479 v.resize(n_hexs(),
false);
11483 std::vector<bool>::iterator i = v.begin();
11485 for (; hex != endh; ++hex, ++i)
11486 *i = get_user_flag(hex);
11494 template <
int dim,
int spacedim>
11498 std::vector<bool> v;
11499 save_user_flags_hex(v);
11508 template <
int dim,
int spacedim>
11512 std::vector<bool> v;
11517 load_user_flags_hex(v);
11522 template <
int dim,
int spacedim>
11526 Assert(v.size() == n_hexs(), ExcGridReadError());
11531 std::vector<bool>::const_iterator i = v.begin();
11532 for (; hex != endh; ++hex, ++i)
11534 set_user_flag(hex);
11536 clear_user_flag(hex);
11544 template <
int dim,
int spacedim>
11547 std::vector<unsigned int> &v)
const
11553 std::vector<unsigned int> tmp;
11555 save_user_indices_line(tmp);
11556 v.insert(v.end(), tmp.begin(), tmp.end());
11560 save_user_indices_quad(tmp);
11561 v.insert(v.end(), tmp.begin(), tmp.end());
11566 save_user_indices_hex(tmp);
11567 v.insert(v.end(), tmp.begin(), tmp.end());
11576 template <
int dim,
int spacedim>
11579 const std::vector<unsigned int> &v)
11582 std::vector<unsigned int> tmp;
11586 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
11588 load_user_indices_line(tmp);
11593 tmp.insert(tmp.end(),
11594 v.begin() + n_lines(),
11595 v.begin() + n_lines() + n_quads());
11596 load_user_indices_quad(tmp);
11602 tmp.insert(tmp.end(),
11603 v.begin() + n_lines() + n_quads(),
11604 v.begin() + n_lines() + n_quads() + n_hexs());
11605 load_user_indices_hex(tmp);
11616 template <
typename Iterator>
11620 return i->user_index();
11625 template <
int structdim,
int dim,
int spacedim>
11636 template <
typename Iterator>
11638 set_user_index(
const Iterator &i,
const unsigned int x)
11640 i->set_user_index(x);
11645 template <
int structdim,
int dim,
int spacedim>
11649 const unsigned int)
11656 template <
int dim,
int spacedim>
11659 std::vector<unsigned int> &v)
const
11661 v.resize(n_lines(), 0);
11662 std::vector<unsigned int>::iterator i = v.begin();
11664 for (; line != endl; ++line, ++i)
11665 *i = line->user_index();
11670 template <
int dim,
int spacedim>
11673 const std::vector<unsigned int> &v)
11675 Assert(v.size() == n_lines(), ExcGridReadError());
11678 std::vector<unsigned int>::const_iterator i = v.begin();
11679 for (; line != endl; ++line, ++i)
11680 line->set_user_index(*i);
11684 template <
int dim,
int spacedim>
11687 std::vector<unsigned int> &v)
const
11689 v.resize(n_quads(), 0);
11693 std::vector<unsigned int>::iterator i = v.begin();
11695 for (; quad != endq; ++quad, ++i)
11696 *i = get_user_index(quad);
11702 template <
int dim,
int spacedim>
11705 const std::vector<unsigned int> &v)
11707 Assert(v.size() == n_quads(), ExcGridReadError());
11712 std::vector<unsigned int>::const_iterator i = v.begin();
11713 for (; quad != endq; ++quad, ++i)
11714 set_user_index(quad, *i);
11719 template <
int dim,
int spacedim>
11722 std::vector<unsigned int> &v)
const
11724 v.resize(n_hexs(), 0);
11728 std::vector<unsigned int>::iterator i = v.begin();
11730 for (; hex != endh; ++hex, ++i)
11731 *i = get_user_index(hex);
11737 template <
int dim,
int spacedim>
11740 const std::vector<unsigned int> &v)
11742 Assert(v.size() == n_hexs(), ExcGridReadError());
11747 std::vector<unsigned int>::const_iterator i = v.begin();
11748 for (; hex != endh; ++hex, ++i)
11749 set_user_index(hex, *i);
11760 template <
typename Iterator>
11762 get_user_pointer(
const Iterator &i)
11764 return i->user_pointer();
11769 template <
int structdim,
int dim,
int spacedim>
11780 template <
typename Iterator>
11782 set_user_pointer(
const Iterator &i,
void *x)
11784 i->set_user_pointer(x);
11789 template <
int structdim,
int dim,
int spacedim>
11800 template <
int dim,
int spacedim>
11808 std::vector<void *> tmp;
11810 save_user_pointers_line(tmp);
11811 v.insert(v.end(), tmp.begin(), tmp.end());
11815 save_user_pointers_quad(tmp);
11816 v.insert(v.end(), tmp.begin(), tmp.end());
11821 save_user_pointers_hex(tmp);
11822 v.insert(v.end(), tmp.begin(), tmp.end());
11831 template <
int dim,
int spacedim>
11836 std::vector<void *> tmp;
11840 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
11842 load_user_pointers_line(tmp);
11847 tmp.insert(tmp.end(),
11848 v.begin() + n_lines(),
11849 v.begin() + n_lines() + n_quads());
11850 load_user_pointers_quad(tmp);
11856 tmp.insert(tmp.end(),
11857 v.begin() + n_lines() + n_quads(),
11858 v.begin() + n_lines() + n_quads() + n_hexs());
11859 load_user_pointers_hex(tmp);
11868 template <
int dim,
int spacedim>
11871 std::vector<void *> &v)
const
11873 v.resize(n_lines(),
nullptr);
11874 std::vector<void *>::iterator i = v.begin();
11876 for (; line != endl; ++line, ++i)
11877 *i = line->user_pointer();
11882 template <
int dim,
int spacedim>
11885 const std::vector<void *> &v)
11887 Assert(v.size() == n_lines(), ExcGridReadError());
11890 std::vector<void *>::const_iterator i = v.begin();
11891 for (; line != endl; ++line, ++i)
11892 line->set_user_pointer(*i);
11897 template <
int dim,
int spacedim>
11900 std::vector<void *> &v)
const
11902 v.resize(n_quads(),
nullptr);
11906 std::vector<void *>::iterator i = v.begin();
11908 for (; quad != endq; ++quad, ++i)
11909 *i = get_user_pointer(quad);
11915 template <
int dim,
int spacedim>
11918 const std::vector<void *> &v)
11920 Assert(v.size() == n_quads(), ExcGridReadError());
11925 std::vector<void *>::const_iterator i = v.begin();
11926 for (; quad != endq; ++quad, ++i)
11927 set_user_pointer(quad, *i);
11932 template <
int dim,
int spacedim>
11935 std::vector<void *> &v)
const
11937 v.resize(n_hexs(),
nullptr);
11941 std::vector<void *>::iterator i = v.begin();
11943 for (; hex != endh; ++hex, ++i)
11944 *i = get_user_pointer(hex);
11950 template <
int dim,
int spacedim>
11953 const std::vector<void *> &v)
11955 Assert(v.size() == n_hexs(), ExcGridReadError());
11960 std::vector<void *>::const_iterator i = v.begin();
11961 for (; hex != endh; ++hex, ++i)
11962 set_user_pointer(hex, *i);
11971 template <
int dim,
int spacedim>
11978 return begin_raw_line(
level);
11980 return begin_raw_quad(
level);
11982 return begin_raw_hex(
level);
11991 template <
int dim,
int spacedim>
11998 return begin_line(
level);
12000 return begin_quad(
level);
12002 return begin_hex(
level);
12011 template <
int dim,
int spacedim>
12018 return begin_active_line(
level);
12020 return begin_active_quad(
level);
12022 return begin_active_hex(
level);
12031 template <
int dim,
int spacedim>
12035 const unsigned int level = levels.size() - 1;
12036 if (levels[
level]->cells.cells.size() == 0)
12043 levels[
level]->cells.cells.size() - 1);
12046 if (ri->used() ==
true)
12049 if (ri->used() ==
true)
12056 template <
int dim,
int spacedim>
12066 if (cell->is_active() ==
true)
12069 if (cell->is_active() ==
true)
12077 template <
int dim,
int spacedim>
12088 template <
int dim,
int spacedim>
12100 if (
level >= levels.size())
12103 ExcInvalidLevel(
level, n_global_levels()));
12110 if (
level < levels.size() - 1)
12111 return begin_raw(
level + 1);
12117 template <
int dim,
int spacedim>
12129 if (
level >= levels.size())
12132 ExcInvalidLevel(
level, n_global_levels()));
12139 if (
level < levels.size() - 1)
12146 template <
int dim,
int spacedim>
12158 if (
level >= levels.size())
12161 ExcInvalidLevel(
level, n_global_levels()));
12169 begin_active(
level + 1));
12174 template <
int dim,
int spacedim>
12183 template <
int dim,
int spacedim>
12194 template <
int dim,
int spacedim>
12197 const unsigned int level)
const
12205 template <
int dim,
int spacedim>
12208 const unsigned int level)
const
12219 template <
int dim,
int spacedim>
12229 return begin_line();
12231 return begin_quad();
12240 template <
int dim,
int spacedim>
12250 return begin_active_line();
12252 return begin_active_quad();
12261 template <
int dim,
int spacedim>
12282 template <
int dim,
int spacedim>
12288 begin_active_face(), end_face());
12294 template <
int dim,
int spacedim>
12303 while (i->used() ==
false)
12311 template <
int dim,
int spacedim>
12316 return begin_vertex();
12321 template <
int dim,
int spacedim>
12336 template <
int dim,
int spacedim>
12348 if (
level >= levels.size())
12351 ExcInvalidLevel(
level, n_global_levels()));
12362 if (
level >= levels.size() || levels[
level]->cells.cells.size() == 0)
12376 template <
int dim,
int spacedim>
12384 while (ri->used() ==
false)
12392 template <
int dim,
int spacedim>
12400 while (i->has_children())
12408 template <
int dim,
int spacedim>
12422 template <
int dim,
int spacedim>
12434 if (
level >= levels.size())
12437 ExcInvalidLevel(
level, n_global_levels()));
12452 if (
level >= levels.size() || levels[
level]->cells.cells.size() == 0)
12476 template <
int dim,
int spacedim>
12484 while (ri->used() ==
false)
12492 template <
int dim,
int spacedim>
12500 while (i->has_children())
12508 template <
int dim,
int spacedim>
12521 template <
int dim,
int spacedim>
12533 if (
level >= levels.size())
12536 ExcInvalidLevel(
level, n_global_levels()));
12552 if (
level >= levels.size() || levels[
level]->cells.cells.size() == 0)
12567 template <
int dim,
int spacedim>
12575 while (ri->used() ==
false)
12583 template <
int dim,
int spacedim>
12591 while (i->has_children())
12599 template <
int dim,
int spacedim>
12615 namespace TriangulationImplementation
12617 inline unsigned int
12624 inline unsigned int
12632 inline unsigned int
12639 inline unsigned int
12647 inline unsigned int
12654 inline unsigned int
12665 template <
int dim,
int spacedim>
12673 template <
int dim,
int spacedim>
12680 template <
int dim,
int spacedim>
12689 template <
int dim,
int spacedim>
12696 return n_used_vertices();
12708 template <
int dim,
int spacedim>
12715 return n_vertices();
12717 return n_raw_lines();
12719 return n_raw_quads();
12727 template <
int dim,
int spacedim>
12734 return n_used_vertices();
12736 return n_active_lines();
12738 return n_active_quads();
12746 template <
int dim,
int spacedim>
12753 return n_raw_lines(
level);
12755 return n_raw_quads(
level);
12757 return n_raw_hexs(
level);
12766 template <
int dim,
int spacedim>
12773 return n_lines(
level);
12775 return n_quads(
level);
12777 return n_hexs(
level);
12786 template <
int dim,
int spacedim>
12793 return n_active_lines(
level);
12795 return n_active_quads(
level);
12797 return n_active_hexs(
level);
12805 template <
int dim,
int spacedim>
12809 for (
unsigned int lvl = 0; lvl < n_global_levels() - 1; lvl++)
12817 template <
int dim,
int spacedim>
12821 return number_cache.n_lines;
12831 return levels[
level]->cells.cells.size();
12850 return levels[
level]->cells.cells.size();
12868 return levels[
level]->cells.cells.size();
12881 template <
int dim,
int spacedim>
12890 template <
int dim,
int spacedim>
12894 return faces->
lines.cells.size();
12898 template <
int dim,
int spacedim>
12904 return number_cache.n_lines_level[
level];
12908 template <
int dim,
int spacedim>
12912 return number_cache.n_active_lines;
12916 template <
int dim,
int spacedim>
12923 return number_cache.n_active_lines_level[
level];
13073 template <
int dim,
int spacedim>
13077 return number_cache.n_quads;
13081 template <
int dim,
int spacedim>
13087 return number_cache.n_quads_level[
level];
13097 return levels[
level]->cells.cells.size();
13107 return levels[
level]->cells.cells.size();
13121 template <
int dim,
int spacedim>
13140 template <
int dim,
int spacedim>
13144 return number_cache.n_active_quads;
13148 template <
int dim,
int spacedim>
13155 return number_cache.n_active_quads_level[
level];
13159 template <
int dim,
int spacedim>
13168 template <
int dim,
int spacedim>
13177 template <
int dim,
int spacedim>
13185 template <
int dim,
int spacedim>
13194 template <
int dim,
int spacedim>
13206 return number_cache.n_hexes;
13217 return number_cache.n_hexes_level[
level];
13227 return levels[
level]->cells.cells.size();
13235 return number_cache.n_active_hexes;
13246 return number_cache.n_active_hexes_level[
level];
13251 template <
int dim,
int spacedim>
13255 return std::count(vertices_used.begin(), vertices_used.end(),
true);
13260 template <
int dim,
int spacedim>
13261 const std::vector<bool> &
13264 return vertices_used;
13294 template <
int dim,
int spacedim>
13302 unsigned int max_vertex_index = 0;
13303 for (; cell != endc; ++cell)
13305 if (cell->vertex_index(vertex) > max_vertex_index)
13306 max_vertex_index = cell->vertex_index(vertex);
13312 std::vector<unsigned short int> usage_count(max_vertex_index + 1, 0);
13316 for (cell =
begin(); cell != endc; ++cell)
13318 ++usage_count[cell->vertex_index(vertex)];
13321 static_cast<unsigned int>(
13322 *std::max_element(usage_count.begin(), usage_count.end())));
13327 template <
int dim,
int spacedim>
13336 template <
int dim,
int spacedim>
13345 template <
int dim,
int spacedim>
13354 template <
int dim,
int spacedim>
13358 &periodicity_vector)
13360 periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(),
13361 periodicity_vector.begin(),
13362 periodicity_vector.end());
13365 update_periodic_face_map();
13370 template <
int dim,
int spacedim>
13371 const typename std::map<
13372 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
unsigned int>,
13373 std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
13378 return periodic_face_map;
13383 template <
int dim,
int spacedim>
13387 prepare_coarsening_and_refinement();
13392 if (smooth_grid & limit_level_difference_at_vertices)
13396 signals.pre_refinement();
13398 execute_coarsening();
13400 const DistortedCellList cells_with_distorted_children = execute_refinement();
13405 if (smooth_grid & limit_level_difference_at_vertices)
13410 update_neighbors(*
this);
13411 reset_active_cell_indices();
13414 signals.post_refinement();
13416 AssertThrow(cells_with_distorted_children.distorted_cells.size() == 0,
13417 cells_with_distorted_children);
13419 update_periodic_face_map();
13424 template <
int dim,
int spacedim>
13428 unsigned int active_cell_index = 0;
13430 if ((cell->used() ==
false) || cell->has_children())
13434 cell->set_active_cell_index(active_cell_index);
13435 ++active_cell_index;
13442 template <
int dim,
int spacedim>
13447 periodic_face_map.clear();
13449 typename std::vector<
13451 for (it = periodic_face_pairs_level_0.begin();
13452 it != periodic_face_pairs_level_0.end();
13455 update_periodic_face_map_recursively<dim, spacedim>(it->cell[0],
13460 periodic_face_map);
13463 std::bitset<3> inverted_orientation;
13465 bool orientation, flip, rotation;
13466 orientation = it->orientation[0];
13467 rotation = it->orientation[2];
13468 flip = orientation ? rotation ^ it->orientation[1] : it->orientation[1];
13469 inverted_orientation[0] = orientation;
13470 inverted_orientation[1] = flip;
13471 inverted_orientation[2] = rotation;
13473 update_periodic_face_map_recursively<dim, spacedim>(it->cell[1],
13477 inverted_orientation,
13478 periodic_face_map);
13482 typename std::map<std::pair<cell_iterator, unsigned int>,
13483 std::pair<std::pair<cell_iterator, unsigned int>,
13484 std::bitset<3>>>::const_iterator it_test;
13485 for (it_test = periodic_face_map.begin(); it_test != periodic_face_map.end();
13489 it_test->first.first;
13491 it_test->second.first.first;
13492 if (cell_1->level() == cell_2->level())
13496 Assert(periodic_face_map[it_test->second.first].first ==
13505 template <
int dim,
int spacedim>
13513 vertices_used.clear();
13521 template <
int dim,
int spacedim>
13525 const DistortedCellList cells_with_distorted_children =
13527 *
this, check_for_distorted_cells);
13533 *
this, levels.size(), number_cache);
13537 levels[
level]->cells.monitor_memory(dim);
13544 while (cell != endc)
13548 return cells_with_distorted_children;
13553 template <
int dim,
int spacedim>
13561 std::vector<unsigned int> line_cell_count =
13562 count_cells_bounded_by_line(*
this);
13563 std::vector<unsigned int> quad_cell_count =
13564 count_cells_bounded_by_quad(*
this);
13576 clear_user_flags();
13579 for (; cell != endc; ++cell)
13580 if (!cell->is_active())
13581 if (cell->child(0)->coarsen_flag_set())
13583 cell->set_user_flag();
13584 for (
unsigned int child = 0; child < cell->n_children(); ++child)
13586 Assert(cell->child(child)->coarsen_flag_set(),
13588 cell->child(child)->clear_coarsen_flag();
13603 if (levels.size() >= 2)
13604 for (cell = last(); cell != endc; --cell)
13605 if (cell->level() <=
static_cast<int>(levels.size() - 2) &&
13606 cell->user_flag_set())
13609 signals.pre_coarsening_on_cell(cell);
13617 *
this, levels.size(), number_cache);
13621 for (cell =
begin(); cell != endc; ++cell)
13628 template <
int dim,
int spacedim>
13646 save_coarsen_flags(previous_coarsen_flags);
13648 std::vector<int> vertex_level(
vertices.size(), 0);
13650 bool continue_iterating =
true;
13654 if (smooth_grid & limit_level_difference_at_vertices)
13656 Assert(!anisotropic_refinement,
13657 ExcMessage(
"In case of anisotropic refinement the "
13658 "limit_level_difference_at_vertices flag for "
13659 "mesh smoothing must not be set!"));
13663 std::fill(vertex_level.begin(), vertex_level.end(), 0);
13665 for (; cell != endc; ++cell)
13667 if (cell->refine_flag_set())
13668 for (
const unsigned int vertex :
13670 vertex_level[cell->vertex_index(vertex)] =
13671 std::max(vertex_level[cell->vertex_index(vertex)],
13672 cell->level() + 1);
13673 else if (!cell->coarsen_flag_set())
13674 for (
const unsigned int vertex :
13676 vertex_level[cell->vertex_index(vertex)] =
13677 std::max(vertex_level[cell->vertex_index(vertex)],
13688 for (
const unsigned int vertex :
13690 vertex_level[cell->vertex_index(vertex)] =
13691 std::max(vertex_level[cell->vertex_index(vertex)],
13692 cell->level() - 1);
13706 for (cell = last_active(); cell != endc; --cell)
13707 if (cell->refine_flag_set() ==
false)
13709 for (
const unsigned int vertex :
13711 if (vertex_level[cell->vertex_index(vertex)] >=
13715 cell->clear_coarsen_flag();
13720 if (vertex_level[cell->vertex_index(vertex)] >
13723 cell->set_refine_flag();
13725 for (
const unsigned int v :
13727 vertex_level[cell->vertex_index(v)] =
13728 std::max(vertex_level[cell->vertex_index(v)],
13729 cell->level() + 1);
13751 clear_user_flags();
13755 for (; acell != end_ac; ++acell)
13756 acell->clear_coarsen_flag();
13759 for (; cell != endc; ++cell)
13762 if (cell->is_active())
13765 const unsigned int n_children = cell->n_children();
13766 unsigned int flagged_children = 0;
13767 for (
unsigned int child = 0; child < n_children; ++child)
13768 if (cell->child(child)->is_active() &&
13769 cell->child(child)->coarsen_flag_set())
13771 ++flagged_children;
13773 cell->child(child)->clear_coarsen_flag();
13778 if (flagged_children == n_children)
13779 cell->set_user_flag();
13785 for (cell =
begin(); cell != endc; ++cell)
13806 for (cell = last(); cell != endc; --cell)
13807 if (cell->user_flag_set())
13811 template coarsening_allowed<dim, spacedim>(cell))
13812 for (
unsigned int c = 0; c < cell->n_children(); ++c)
13814 Assert(cell->child(c)->refine_flag_set() ==
false,
13817 cell->child(c)->set_coarsen_flag();
13822 clear_user_flags();
13828 save_coarsen_flags(current_coarsen_flags);
13830 continue_iterating = (current_coarsen_flags != previous_coarsen_flags);
13831 previous_coarsen_flags = current_coarsen_flags;
13833 while (continue_iterating ==
true);
13844 std::vector<bool> flags_before;
13845 save_coarsen_flags(flags_before);
13848 fix_coarsen_flags();
13850 std::vector<bool> flags_after;
13851 save_coarsen_flags(flags_after);
13853 return (flags_before != flags_after);
13863 std::vector<bool> flags_before;
13864 save_coarsen_flags(flags_before);
13867 fix_coarsen_flags();
13869 std::vector<bool> flags_after;
13870 save_coarsen_flags(flags_after);
13872 return (flags_before != flags_after);
13882 std::vector<bool> flags_before;
13883 save_coarsen_flags(flags_before);
13886 fix_coarsen_flags();
13888 std::vector<bool> flags_after;
13889 save_coarsen_flags(flags_after);
13891 return (flags_before != flags_after);
13902 template <
int dim,
int spacedim>
13904 possibly_do_not_produce_unrefined_islands(
13909 unsigned int n_neighbors = 0;
13912 unsigned int count = 0;
13920 if (face_will_be_refined_by_neighbor(cell, n))
13927 if (count == n_neighbors ||
13928 (count >= n_neighbors - 1 &&
13931 for (
unsigned int c = 0; c < cell->n_children(); ++c)
13932 cell->child(c)->clear_coarsen_flag();
13935 if (!cell->at_boundary(face) &&
13936 (!cell->neighbor(face)->is_active()) &&
13937 (cell_will_be_coarsened(cell->neighbor(face))))
13938 possibly_do_not_produce_unrefined_islands<dim, spacedim>(
13939 cell->neighbor(face));
13954 template <
int dim,
int spacedim>
13956 possibly_refine_unrefined_island(
13958 const bool allow_anisotropic_smoothing)
13973 if (allow_anisotropic_smoothing ==
false)
13976 unsigned int refined_neighbors = 0, unrefined_neighbors = 0;
13978 if (!cell->at_boundary(face))
13980 if (face_will_be_refined_by_neighbor(cell, face))
13981 ++refined_neighbors;
13983 ++unrefined_neighbors;
13986 if (unrefined_neighbors < refined_neighbors)
13988 cell->clear_coarsen_flag();
13989 cell->set_refine_flag();
13994 if (unrefined_neighbors > 0)
13996 if (!cell->at_boundary(face) &&
13997 (face_will_be_refined_by_neighbor(cell, face) ==
false) &&
13998 (cell->neighbor(face)->has_children() ==
false) &&
13999 (cell->neighbor(face)->refine_flag_set() ==
false))
14000 possibly_refine_unrefined_island<dim, spacedim>(
14001 cell->neighbor(face), allow_anisotropic_smoothing);
14013 for (
unsigned int face_pair = 0;
14014 face_pair < GeometryInfo<dim>::faces_per_cell / 2;
14023 for (
unsigned int face_index = 0; face_index < 2; ++face_index)
14025 unsigned int face = 2 * face_pair + face_index;
14032 face_will_be_refined_by_neighbor<dim, spacedim>(
14033 cell, face, expected_face_ref_case);
14047 directional_cell_refinement_case =
14048 (directional_cell_refinement_case &
14051 expected_face_ref_case,
14053 cell->face_orientation(face),
14054 cell->face_flip(face),
14055 cell->face_rotation(face)));
14061 Assert(directional_cell_refinement_case <
14064 smoothing_cell_refinement_case =
14065 smoothing_cell_refinement_case | directional_cell_refinement_case;
14070 if (smoothing_cell_refinement_case)
14072 cell->clear_coarsen_flag();
14073 cell->set_refine_flag(cell->refine_flag_set() |
14074 smoothing_cell_refinement_case);
14081 template <
int dim,
int spacedim>
14087 std::vector<bool> flags_before[2];
14088 save_coarsen_flags(flags_before[0]);
14089 save_refine_flags(flags_before[1]);
14107 std::vector<bool> flags_before_loop[2] = {flags_before[0], flags_before[1]};
14184 if (((smooth_grid & coarsest_level_1) || (smooth_grid & patch_level_1)) &&
14187 for (
const auto &cell : active_cell_iterators_on_level(1))
14188 cell->clear_coarsen_flag();
14191 bool mesh_changed_in_this_loop =
false;
14200 if (smooth_grid & do_not_produce_unrefined_islands &&
14201 !(smooth_grid & patch_level_1))
14203 for (
const auto &cell : cell_iterators())
14207 if (!cell->is_active() && cell_will_be_coarsened(cell))
14208 possibly_do_not_produce_unrefined_islands<dim, spacedim>(cell);
14232 if (smooth_grid & (eliminate_refined_inner_islands |
14233 eliminate_refined_boundary_islands) &&
14234 !(smooth_grid & patch_level_1))
14236 for (
const auto &cell : cell_iterators())
14237 if (!cell->is_active() ||
14238 (cell->is_active() && cell->refine_flag_set() &&
14239 cell->is_locally_owned()))
14246 bool all_children_active =
true;
14247 if (!cell->is_active())
14248 for (
unsigned int c = 0; c < cell->n_children(); ++c)
14249 if (!cell->child(c)->is_active() ||
14250 cell->child(c)->is_ghost() ||
14251 cell->child(c)->is_artificial())
14253 all_children_active =
false;
14257 if (all_children_active)
14264 unsigned int unrefined_neighbors = 0, total_neighbors = 0;
14266 for (
const unsigned int n :
14274 if (!face_will_be_refined_by_neighbor(cell, n))
14275 ++unrefined_neighbors;
14290 if ((unrefined_neighbors == total_neighbors) &&
14291 (((unrefined_neighbors ==
14293 (smooth_grid & eliminate_refined_inner_islands)) ||
14294 ((unrefined_neighbors <
14297 eliminate_refined_boundary_islands))) &&
14298 (total_neighbors != 0))
14300 if (!cell->is_active())
14301 for (
unsigned int c = 0; c < cell->n_children(); ++c)
14303 cell->child(c)->clear_refine_flag();
14304 cell->child(c)->set_coarsen_flag();
14307 cell->clear_refine_flag();
14322 if (smooth_grid & limit_level_difference_at_vertices)
14324 Assert(!anisotropic_refinement,
14325 ExcMessage(
"In case of anisotropic refinement the "
14326 "limit_level_difference_at_vertices flag for "
14327 "mesh smoothing must not be set!"));
14331 std::vector<int> vertex_level(
vertices.size(), 0);
14332 for (
const auto &cell : active_cell_iterators())
14334 if (cell->refine_flag_set())
14335 for (
const unsigned int vertex :
14337 vertex_level[cell->vertex_index(vertex)] =
14338 std::max(vertex_level[cell->vertex_index(vertex)],
14339 cell->level() + 1);
14340 else if (!cell->coarsen_flag_set())
14341 for (
const unsigned int vertex :
14343 vertex_level[cell->vertex_index(vertex)] =
14344 std::max(vertex_level[cell->vertex_index(vertex)],
14353 for (
const unsigned int vertex :
14355 vertex_level[cell->vertex_index(vertex)] =
14356 std::max(vertex_level[cell->vertex_index(vertex)],
14357 cell->level() - 1);
14372 if (cell->refine_flag_set() ==
false)
14374 for (
const unsigned int vertex :
14376 if (vertex_level[cell->vertex_index(vertex)] >=
14380 cell->clear_coarsen_flag();
14385 if (vertex_level[cell->vertex_index(vertex)] >
14388 cell->set_refine_flag();
14390 for (
const unsigned int v :
14392 vertex_level[cell->vertex_index(v)] =
14393 std::max(vertex_level[cell->vertex_index(v)],
14394 cell->level() + 1);
14415 if (smooth_grid & eliminate_unrefined_islands)
14420 if (cell->refine_flag_set() !=
14422 possibly_refine_unrefined_island<dim, spacedim>(
14423 cell, (smooth_grid & allow_anisotropic_smoothing) != 0);
14451 if (smooth_grid & patch_level_1)
14461 for (
const auto &cell : cell_iterators())
14462 if (!cell->is_active())
14468 if (cell->child(0)->has_children() ==
true)
14475 for (
unsigned int i = 0; i < cell->n_children(); ++i)
14476 combined_ref_case =
14477 combined_ref_case | cell->child(i)->refine_flag_set();
14479 for (
unsigned int i = 0; i < cell->n_children(); ++i)
14483 child->clear_coarsen_flag();
14484 child->set_refine_flag(combined_ref_case);
14497 for (
const auto &cell : cell_iterators())
14503 if (cell->is_active() || cell->child(0)->is_active())
14510 const unsigned int n_children = cell->n_children();
14511 bool has_active_grandchildren =
false;
14513 for (
unsigned int i = 0; i < n_children; ++i)
14514 if (cell->child(i)->child(0)->is_active())
14516 has_active_grandchildren =
true;
14520 if (has_active_grandchildren ==
false)
14526 unsigned int n_grandchildren = 0;
14529 unsigned int n_coarsen_flags = 0;
14534 for (
unsigned int c = 0; c < n_children; ++c)
14541 const unsigned int nn_children = child->n_children();
14542 n_grandchildren += nn_children;
14547 if (child->child(0)->is_active())
14548 for (
unsigned int cc = 0; cc < nn_children; ++cc)
14549 if (child->child(cc)->coarsen_flag_set())
14564 if ((n_coarsen_flags != n_grandchildren) && (n_coarsen_flags > 0))
14565 for (
unsigned int c = 0; c < n_children; ++c)
14568 if (child->child(0)->is_active())
14569 for (
unsigned int cc = 0; cc < child->n_children(); ++cc)
14570 child->child(cc)->clear_coarsen_flag();
14598 bool changed =
true;
14604 for (; cell != endc; --cell)
14605 if (cell->refine_flag_set())
14613 const bool has_periodic_neighbor =
14614 cell->has_periodic_neighbor(i);
14615 const bool has_neighbor_or_periodic_neighbor =
14616 !cell->at_boundary(i) || has_periodic_neighbor;
14617 if (has_neighbor_or_periodic_neighbor &&
14619 cell->refine_flag_set(), i) !=
14632 if (cell->neighbor_or_periodic_neighbor(i)->is_active())
14634 if ((!has_periodic_neighbor &&
14635 cell->neighbor_is_coarser(i)) ||
14636 (has_periodic_neighbor &&
14637 cell->periodic_neighbor_is_coarser(i)))
14639 if (cell->neighbor_or_periodic_neighbor(i)
14640 ->coarsen_flag_set())
14641 cell->neighbor_or_periodic_neighbor(i)
14642 ->clear_coarsen_flag();
14659 allow_anisotropic_smoothing)
14661 has_periodic_neighbor ?
14662 cell->periodic_neighbor(i)
14663 ->flag_for_face_refinement(
14665 ->periodic_neighbor_of_coarser_periodic_neighbor(
14670 ->flag_for_face_refinement(
14672 ->neighbor_of_coarser_neighbor(
14679 ->neighbor_or_periodic_neighbor(
14681 ->refine_flag_set())
14683 cell->neighbor_or_periodic_neighbor(i)
14684 ->set_refine_flag();
14816 std::pair<unsigned int, unsigned int>
14818 has_periodic_neighbor ?
14820 ->periodic_neighbor_of_coarser_periodic_neighbor(
14822 cell->neighbor_of_coarser_neighbor(i);
14823 unsigned int refined_along_x = 0,
14824 refined_along_y = 0,
14825 to_be_refined_along_x = 0,
14826 to_be_refined_along_y = 0;
14828 const int this_face_index =
14829 cell->face_index(i);
14838 const auto parent_face = [&]() {
14839 if (has_periodic_neighbor)
14841 const auto neighbor =
14842 cell->periodic_neighbor(i);
14843 const auto parent_face_no =
14845 ->periodic_neighbor_of_periodic_neighbor(
14848 neighbor->periodic_neighbor(
14850 return parent->face(parent_face_no);
14853 return cell->neighbor(i)->face(
14857 if ((this_face_index ==
14858 parent_face->child_index(0)) ||
14859 (this_face_index ==
14860 parent_face->child_index(1)))
14869 parent_face->refinement_case();
14887 cell->refine_flag_set(),
14889 cell->face_orientation(i),
14890 cell->face_flip(i),
14891 cell->face_rotation(i));
14894 ++to_be_refined_along_x;
14897 ++to_be_refined_along_y;
14902 allow_anisotropic_smoothing) ||
14903 cell->neighbor_or_periodic_neighbor(i)
14904 ->refine_flag_set())
14906 if (refined_along_x +
14907 to_be_refined_along_x >
14911 ->neighbor_or_periodic_neighbor(i)
14912 ->flag_for_face_refinement(
14916 if (refined_along_y +
14917 to_be_refined_along_y >
14921 ->neighbor_or_periodic_neighbor(i)
14922 ->flag_for_face_refinement(
14930 ->neighbor_or_periodic_neighbor(i)
14931 ->refine_flag_set() !=
14933 dim>::isotropic_refinement)
14935 cell->neighbor_or_periodic_neighbor(i)
14936 ->set_refine_flag();
14942 cell->neighbor_or_periodic_neighbor(i);
14945 nb->refine_flag_set(),
14947 nb->face_orientation(nb_indices.first),
14948 nb->face_flip(nb_indices.first),
14949 nb->face_rotation(nb_indices.first));
14951 !(refined_along_x ||
14952 to_be_refined_along_x))
14953 changed |= cell->flag_for_face_refinement(
14957 !(refined_along_y ||
14958 to_be_refined_along_y))
14959 changed |= cell->flag_for_face_refinement(
14966 cell->neighbor_or_periodic_neighbor(i)
14967 ->clear_coarsen_flag();
14968 const unsigned int nb_nb =
14969 has_periodic_neighbor ?
14971 ->periodic_neighbor_of_periodic_neighbor(
14973 cell->neighbor_of_neighbor(i);
14975 cell->neighbor_or_periodic_neighbor(i);
14978 neighbor->refine_flag_set(),
14980 neighbor->face_orientation(nb_nb),
14981 neighbor->face_flip(nb_nb),
14982 neighbor->face_rotation(nb_nb));
14985 cell->refine_flag_set(),
14987 cell->face_orientation(i),
14988 cell->face_flip(i),
14989 cell->face_rotation(i));
14994 if ((face_ref_case ==
14996 needed_face_ref_case ==
15000 needed_face_ref_case ==
15003 changed = cell->flag_for_face_refinement(
15005 neighbor->flag_for_face_refinement(
15006 nb_nb, needed_face_ref_case);
15013 face_ref_case = cell->face(i)->refinement_case(),
15014 needed_face_ref_case =
15016 cell->refine_flag_set(),
15018 cell->face_orientation(i),
15019 cell->face_flip(i),
15020 cell->face_rotation(i));
15025 needed_face_ref_case ==
15028 needed_face_ref_case ==
15031 cell->flag_for_face_refinement(i,
15052 fix_coarsen_flags();
15055 std::vector<bool> flags_after_loop[2];
15056 save_coarsen_flags(flags_after_loop[0]);
15057 save_refine_flags(flags_after_loop[1]);
15061 mesh_changed_in_this_loop =
15062 ((flags_before_loop[0] != flags_after_loop[0]) ||
15063 (flags_before_loop[1] != flags_after_loop[1]));
15067 flags_before_loop[0].swap(flags_after_loop[0]);
15068 flags_before_loop[1].swap(flags_after_loop[1]);
15070 while (mesh_changed_in_this_loop);
15076 return ((flags_before[0] != flags_before_loop[0]) ||
15077 (flags_before[1] != flags_before_loop[1]));
15082 template <
int dim,
int spacedim>
15085 const unsigned int magic_number1,
15086 const std::vector<bool> &v,
15087 const unsigned int magic_number2,
15088 std::ostream & out)
15090 const unsigned int N = v.size();
15091 unsigned char * flags =
new unsigned char[
N / 8 + 1];
15092 for (
unsigned int i = 0; i <
N / 8 + 1; ++i)
15095 for (
unsigned int position = 0; position <
N; ++position)
15096 flags[position / 8] |= (v[position] ? (1 << (position % 8)) : 0);
15105 out << magic_number1 <<
' ' <<
N << std::endl;
15106 for (
unsigned int i = 0; i <
N / 8 + 1; ++i)
15107 out <<
static_cast<unsigned int>(flags[i]) <<
' ';
15109 out << std::endl << magic_number2 << std::endl;
15117 template <
int dim,
int spacedim>
15120 std::vector<bool> &v,
15121 const unsigned int magic_number2,
15126 unsigned int magic_number;
15127 in >> magic_number;
15128 AssertThrow(magic_number == magic_number1, ExcGridReadError());
15134 unsigned char * flags =
new unsigned char[
N / 8 + 1];
15135 unsigned short int tmp;
15136 for (
unsigned int i = 0; i <
N / 8 + 1; ++i)
15142 for (
unsigned int position = 0; position !=
N; ++position)
15143 v[position] = (flags[position / 8] & (1 << (position % 8)));
15145 in >> magic_number;
15146 AssertThrow(magic_number == magic_number2, ExcGridReadError());
15155 template <
int dim,
int spacedim>
15159 std::size_t mem = 0;
15161 for (
unsigned int i = 0; i < levels.size(); ++i)
15165 mem +=
sizeof(manifold);
15166 mem +=
sizeof(smooth_grid);
15168 mem +=
sizeof(faces);
15177 #include "tria.inst"