46 namespace TriangulationImplementation
68 if (active_cell_index_partitioner)
69 mem += active_cell_index_partitioner->memory_consumption();
71 for (
const auto &partitioner : level_cell_index_partitioners)
73 mem += partitioner->memory_consumption();
130 template <
int dim,
int spacedim>
132 cell_is_patch_level_1(
137 unsigned int n_active_children = 0;
138 for (
unsigned int i = 0; i < cell->n_children(); ++i)
139 if (cell->child(i)->is_active())
142 return (n_active_children == 0) ||
143 (n_active_children == cell->n_children());
153 template <
int dim,
int spacedim>
155 cell_will_be_coarsened(
161 if (cell->has_children())
163 unsigned int children_to_coarsen = 0;
164 const unsigned int n_children = cell->n_children();
166 for (
unsigned int c = 0; c < n_children; ++c)
167 if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
168 ++children_to_coarsen;
169 if (children_to_coarsen == n_children)
172 for (
unsigned int c = 0; c < n_children; ++c)
173 if (cell->child(c)->is_active())
174 cell->child(c)->clear_coarsen_flag();
210 template <
int dim,
int spacedim>
212 face_will_be_refined_by_neighbor_internal(
214 const unsigned int face_no,
223 cell->neighbor(face_no);
230 if (neighbor->has_children())
235 if (cell_will_be_coarsened(neighbor))
244 expected_face_ref_case = cell->face(face_no)->refinement_case();
256 const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
264 neighbor->face_orientation(neighbor_neighbor),
265 neighbor->face_flip(neighbor_neighbor),
266 neighbor->face_rotation(neighbor_neighbor));
270 neighbor_face = neighbor->face(neighbor_neighbor);
271 const int this_face_index = cell->face_index(face_no);
277 if (neighbor_face->index() == this_face_index)
283 expected_face_ref_case = face_ref_case;
301 for (
unsigned int c = 0; c < neighbor_face->n_children(); ++c)
302 if (neighbor_face->child_index(c) == this_face_index)
311 if ((neighbor_face->refinement_case() | face_ref_case) ==
312 neighbor_face->refinement_case())
331 expected_face_ref_case =
332 ~neighbor_face->refinement_case();
359 template <
int dim,
int spacedim>
361 face_will_be_refined_by_neighbor(
363 const unsigned int face_no)
366 return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
373 template <
int dim,
int spacedim>
375 face_will_be_refined_by_neighbor(
377 const unsigned int face_no,
380 return face_will_be_refined_by_neighbor_internal(cell,
382 expected_face_ref_case);
387 template <
int dim,
int spacedim>
389 satisfies_level1_at_vertex_rule(
392 std::vector<unsigned int> min_adjacent_cell_level(
394 std::vector<unsigned int> max_adjacent_cell_level(
397 for (
const auto &cell :
triangulation.active_cell_iterators())
400 min_adjacent_cell_level[cell->vertex_index(v)] =
401 std::min<unsigned int>(
402 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
403 max_adjacent_cell_level[cell->vertex_index(v)] =
404 std::max<unsigned int>(
405 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
408 for (
unsigned int k = 0; k <
triangulation.n_vertices(); ++k)
410 if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
439 for (
auto &cell : cells)
446 reorder_compatibility(std::vector<
CellData<3>> &cells,
450 static constexpr std::array<
unsigned int,
452 local_vertex_numbering{{0, 1, 5, 4, 2, 3, 7, 6}};
453 for (
auto &cell : cells)
459 cell.
vertices[local_vertex_numbering[i]] = tmp[i];
463 for (
auto &boundary_quad : subcelldata.boundary_quads)
487 template <
int dim,
int spacedim>
492 if (line->has_children())
493 return line->child(0)->vertex_index(1);
498 template <
int dim,
int spacedim>
503 switch (
static_cast<unsigned char>(quad->refinement_case()))
506 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
509 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
512 return quad->child(0)->vertex_index(3);
521 template <
int dim,
int spacedim>
526 switch (
static_cast<unsigned char>(hex->refinement_case()))
529 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
532 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
535 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
538 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
541 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
544 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
547 return hex->child(0)->vertex_index(7);
568 template <
class TRIANGULATION>
569 inline typename TRIANGULATION::DistortedCellList
570 collect_distorted_coarse_cells(
const TRIANGULATION &)
572 return typename TRIANGULATION::DistortedCellList();
590 for (
const auto &cell :
triangulation.cell_iterators_on_level(0))
607 return distorted_cells;
619 has_distorted_children(
624 for (
unsigned int c = 0; c < cell->n_children(); ++c)
628 vertices[i] = cell->child(c)->vertex(i);
634 if (determinants[i] <=
650 template <
int dim,
int spacedim>
652 has_distorted_children(
659 template <
int dim,
int spacedim>
661 update_periodic_face_map_recursively(
664 unsigned int n_face_1,
665 unsigned int n_face_2,
666 const std::bitset<3> & orientation,
672 std::bitset<3>>> &periodic_face_map)
675 const FaceIterator face_1 = cell_1->face(n_face_1);
676 const FaceIterator face_2 = cell_2->face(n_face_2);
678 const bool face_orientation = orientation[0];
679 const bool face_flip = orientation[1];
680 const bool face_rotation = orientation[2];
682 Assert((dim != 1) || (face_orientation ==
true && face_flip ==
false &&
683 face_rotation ==
false),
685 "(face_orientation, face_flip, face_rotation) "
686 "is invalid for 1d"));
688 Assert((dim != 2) || (face_orientation ==
true && face_rotation ==
false),
690 "(face_orientation, face_flip, face_rotation) "
691 "is invalid for 2d"));
695 Assert(face_1->at_boundary() && face_2->at_boundary(),
696 ExcMessage(
"Periodic faces must be on the boundary"));
702 Assert(dim == 1 ||
std::abs(cell_1->level() - cell_2->level()) < 2,
707 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
709 const CellFace cell_face_1(cell_1, n_face_1);
710 const CellFace cell_face_2(cell_2, n_face_2);
711 const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
712 cell_face_2, orientation);
714 const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
715 periodic_faces(cell_face_1, cell_face_orientation_2);
719 periodic_face_map.insert(periodic_faces);
723 if (cell_1->has_children())
725 if (cell_2->has_children())
727 update_periodic_face_map_recursively<dim, spacedim>(
728 cell_1->child(n_face_1),
729 cell_2->child(n_face_2),
737 update_periodic_face_map_recursively<dim, spacedim>(
738 cell_1->child(n_face_1),
753 static const int lookup_table_2d[2][2] =
760 static const int lookup_table_3d[2][2][2][4] =
779 if (cell_1->has_children())
781 if (cell_2->has_children())
787 Assert(face_1->n_children() ==
793 for (
unsigned int i = 0;
794 i < GeometryInfo<dim>::max_children_per_face;
802 j = lookup_table_2d[face_flip][i];
805 j = lookup_table_3d[face_orientation][face_flip]
813 unsigned int child_cell_1 =
815 cell_1->refinement_case(),
818 cell_1->face_orientation(n_face_1),
819 cell_1->face_flip(n_face_1),
820 cell_1->face_rotation(n_face_1),
821 face_1->refinement_case());
822 unsigned int child_cell_2 =
824 cell_2->refinement_case(),
827 cell_2->face_orientation(n_face_2),
828 cell_2->face_flip(n_face_2),
829 cell_2->face_rotation(n_face_2),
830 face_2->refinement_case());
832 Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
835 Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
841 update_periodic_face_map_recursively<dim, spacedim>(
842 cell_1->child(child_cell_1),
843 cell_2->child(child_cell_2),
852 for (
unsigned int i = 0;
853 i < GeometryInfo<dim>::max_children_per_face;
857 unsigned int child_cell_1 =
859 cell_1->refinement_case(),
862 cell_1->face_orientation(n_face_1),
863 cell_1->face_flip(n_face_1),
864 cell_1->face_rotation(n_face_1),
865 face_1->refinement_case());
868 update_periodic_face_map_recursively<dim, spacedim>(
869 cell_1->child(child_cell_1),
887 namespace TriangulationImplementation
895 using ::Triangulation;
903 <<
"Something went wrong upon construction of cell "
917 <<
" has negative measure. This typically "
918 <<
"indicates some distortion in the cell, or a mistakenly "
919 <<
"swapped pair of vertices in the input to "
920 <<
"Triangulation::create_triangulation().");
932 <<
"Error while creating cell " << arg1
933 <<
": the vertex index " << arg2 <<
" must be between 0 and "
942 <<
"While trying to assign a boundary indicator to a line: "
943 <<
"the line with end vertices " << arg1 <<
" and " << arg2
944 <<
" does not exist.");
954 <<
"While trying to assign a boundary indicator to a quad: "
955 <<
"the quad with bounding lines " << arg1 <<
", " << arg2
956 <<
", " << arg3 <<
", " << arg4 <<
" does not exist.");
966 <<
"The input data for creating a triangulation contained "
967 <<
"information about a line with indices " << arg1 <<
" and " << arg2
968 <<
" that is described to have boundary indicator "
969 <<
static_cast<int>(arg3)
970 <<
". However, this is an internal line not located on the "
971 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
973 <<
"If this happened at a place where you call "
974 <<
"Triangulation::create_triangulation() yourself, you need "
975 <<
"to check the SubCellData object you pass to this function."
978 <<
"If this happened in a place where you are reading a mesh "
979 <<
"from a file, then you need to investigate why such a line "
980 <<
"ended up in the input file. A typical case is a geometry "
981 <<
"that consisted of multiple parts and for which the mesh "
982 <<
"generator program assumes that the interface between "
983 <<
"two parts is a boundary when that isn't supposed to be "
984 <<
"the case, or where the mesh generator simply assigns "
985 <<
"'geometry indicators' to lines at the perimeter of "
986 <<
"a part that are not supposed to be interpreted as "
987 <<
"'boundary indicators'.");
999 <<
"The input data for creating a triangulation contained "
1000 <<
"information about a quad with indices " << arg1 <<
", " << arg2
1001 <<
", " << arg3 <<
", and " << arg4
1002 <<
" that is described to have boundary indicator "
1003 <<
static_cast<int>(arg5)
1004 <<
". However, this is an internal quad not located on the "
1005 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1007 <<
"If this happened at a place where you call "
1008 <<
"Triangulation::create_triangulation() yourself, you need "
1009 <<
"to check the SubCellData object you pass to this function."
1012 <<
"If this happened in a place where you are reading a mesh "
1013 <<
"from a file, then you need to investigate why such a quad "
1014 <<
"ended up in the input file. A typical case is a geometry "
1015 <<
"that consisted of multiple parts and for which the mesh "
1016 <<
"generator program assumes that the interface between "
1017 <<
"two parts is a boundary when that isn't supposed to be "
1018 <<
"the case, or where the mesh generator simply assigns "
1019 <<
"'geometry indicators' to quads at the surface of "
1020 <<
"a part that are not supposed to be interpreted as "
1021 <<
"'boundary indicators'.");
1030 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1031 <<
" and " << arg2 <<
" appears more than once. "
1032 <<
"This is not allowed.");
1042 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1043 <<
" and " << arg2 <<
" appears multiple times with different (valid) "
1044 << arg3 <<
". This is not allowed.");
1056 <<
"In SubCellData the quad info of the quad with line indices " << arg1
1057 <<
", " << arg2 <<
", " << arg3 <<
" and " << arg4
1058 <<
" appears multiple times with different (valid) " << arg5
1059 <<
". This is not allowed.");
1070 const unsigned int new_quads_in_pairs,
1071 const unsigned int new_quads_single)
1077 unsigned int next_free_single = 0;
1078 unsigned int next_free_pair = 0;
1082 unsigned int n_quads = 0;
1083 unsigned int n_unused_pairs = 0;
1084 unsigned int n_unused_singles = 0;
1085 for (
unsigned int i = 0; i < tria_faces.
quads.
used.size(); ++i)
1089 else if (i + 1 < tria_faces.
quads.
used.size())
1094 if (next_free_single == 0)
1095 next_free_single = i;
1100 if (next_free_pair == 0)
1108 Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
1114 const int additional_single_quads = new_quads_single - n_unused_singles;
1116 unsigned int new_size =
1117 tria_faces.
quads.
used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
1118 if (additional_single_quads > 0)
1119 new_size += additional_single_quads;
1129 q_is_q.reserve(new_size);
1130 q_is_q.insert(q_is_q.end(), new_size - q_is_q.size(),
true);
1151 const unsigned int total_cells,
1152 const unsigned int dimension,
1153 const unsigned int space_dimension)
1205 if (dimension < space_dimension)
1216 tria_level.
parents.reserve((total_cells + 1) / 2);
1218 (total_cells + 1) / 2 -
1222 tria_level.
neighbors.reserve(total_cells * (2 * dimension));
1224 total_cells * (2 * dimension) -
1226 std::make_pair(-1, -1));
1228 if (tria_level.
dim == 2 || tria_level.
dim == 3)
1230 const unsigned int max_faces_per_cell = 2 * dimension;
1232 max_faces_per_cell);
1252 <<
"The containers have sizes " << arg1 <<
" and " << arg2
1253 <<
", which is not as expected.");
1262 const unsigned int true_dimension)
1265 (void)true_dimension;
1292 const unsigned int new_objects_in_pairs,
1293 const unsigned int new_objects_single = 0)
1305 unsigned int n_objects = 0;
1306 unsigned int n_unused_pairs = 0;
1307 unsigned int n_unused_singles = 0;
1308 for (
unsigned int i = 0; i < tria_objects.
used.size(); ++i)
1310 if (tria_objects.
used[i])
1312 else if (i + 1 < tria_objects.
used.size())
1314 if (tria_objects.
used[i + 1])
1331 Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
1332 tria_objects.
used.size(),
1338 const int additional_single_objects =
1339 new_objects_single - n_unused_singles;
1341 unsigned int new_size = tria_objects.
used.size() +
1342 new_objects_in_pairs - 2 * n_unused_pairs;
1343 if (additional_single_objects > 0)
1344 new_size += additional_single_objects;
1347 if (new_size > tria_objects.
n_objects())
1349 const unsigned int max_faces_per_cell =
1351 const unsigned int max_children_per_cell =
1354 tria_objects.
cells.reserve(new_size * max_faces_per_cell);
1355 tria_objects.
cells.insert(tria_objects.
cells.end(),
1360 tria_objects.
used.reserve(new_size);
1361 tria_objects.
used.insert(tria_objects.
used.end(),
1362 new_size - tria_objects.
used.size(),
1371 const unsigned int factor = max_children_per_cell / 2;
1372 tria_objects.
children.reserve(factor * new_size);
1392 tria_objects.
user_data.reserve(new_size);
1393 tria_objects.
user_data.resize(new_size);
1402 if (n_unused_singles == 0)
1410 const unsigned int new_hexes = new_objects_in_pairs;
1412 const unsigned int new_size =
1413 new_hexes + std::count(tria_objects.
used.begin(),
1414 tria_objects.
used.end(),
1418 if (new_size > tria_objects.
n_objects())
1420 const unsigned int max_faces_per_cell =
1423 tria_objects.
cells.reserve(new_size * max_faces_per_cell);
1424 tria_objects.
cells.insert(tria_objects.
cells.end(),
1429 tria_objects.
used.reserve(new_size);
1430 tria_objects.
used.insert(tria_objects.
used.end(),
1431 new_size - tria_objects.
used.size(),
1440 tria_objects.
children.reserve(4 * new_size);
1459 tria_objects.
user_data.reserve(new_size);
1460 tria_objects.
user_data.resize(new_size);
1484 tria_object.
used.size()));
1525 template <
int dim,
int spacedim>
1547 std::vector<unsigned int> & line_cell_count,
1548 std::vector<unsigned int> &quad_cell_count) = 0;
1555 const bool check_for_distorted_cells) = 0;
1584 virtual std::unique_ptr<Policy<dim, spacedim>>
1595 template <
int dim,
int spacedim,
typename T>
1602 T::update_neighbors(
tria);
1609 std::vector<unsigned int> & line_cell_count,
1610 std::vector<unsigned int> &quad_cell_count)
override
1612 T::delete_children(
tria, cell, line_cell_count, quad_cell_count);
1617 const bool check_for_distorted_cells)
override
1619 return T::execute_refinement(
triangulation, check_for_distorted_cells);
1641 return T::template coarsening_allowed<dim, spacedim>(cell);
1644 std::unique_ptr<Policy<dim, spacedim>>
1647 return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
1761 template <
int dim,
int spacedim>
1765 const unsigned int level_objects,
1768 using line_iterator =
1771 number_cache.n_levels = 0;
1772 if (level_objects > 0)
1776 number_cache.n_levels =
level + 1;
1784 number_cache.n_lines = 0;
1785 number_cache.n_active_lines = 0;
1791 number_cache.n_lines_level.resize(number_cache.n_levels);
1792 number_cache.n_active_lines_level.resize(number_cache.n_levels);
1797 number_cache.n_lines_level[
level] = 0;
1798 number_cache.n_active_lines_level[
level] = 0;
1802 (
level == number_cache.n_levels - 1 ?
1805 for (; line != endc; ++line)
1807 ++number_cache.n_lines_level[
level];
1808 if (line->has_children() ==
false)
1809 ++number_cache.n_active_lines_level[
level];
1813 number_cache.n_lines += number_cache.n_lines_level[
level];
1814 number_cache.n_active_lines +=
1815 number_cache.n_active_lines_level[
level];
1821 number_cache.n_lines_level.clear();
1822 number_cache.n_active_lines_level.clear();
1826 for (; line != endc; ++line)
1828 ++number_cache.n_lines;
1829 if (line->has_children() ==
false)
1830 ++number_cache.n_active_lines;
1849 template <
int dim,
int spacedim>
1853 const unsigned int level_objects,
1864 &compute_number_cache_dim<dim, spacedim>),
1870 using quad_iterator =
1876 number_cache.n_quads = 0;
1877 number_cache.n_active_quads = 0;
1887 unsigned int n_levels = 0;
1888 if (level_objects > 0)
1892 n_levels =
level + 1;
1894 number_cache.n_quads_level.resize(n_levels);
1895 number_cache.n_active_quads_level.resize(n_levels);
1900 number_cache.n_quads_level[
level] = 0;
1901 number_cache.n_active_quads_level[
level] = 0;
1905 (
level == n_levels - 1 ?
1908 for (; quad != endc; ++quad)
1910 ++number_cache.n_quads_level[
level];
1911 if (quad->has_children() ==
false)
1912 ++number_cache.n_active_quads_level[
level];
1916 number_cache.n_quads += number_cache.n_quads_level[
level];
1917 number_cache.n_active_quads +=
1918 number_cache.n_active_quads_level[
level];
1924 number_cache.n_quads_level.clear();
1925 number_cache.n_active_quads_level.clear();
1929 for (; quad != endc; ++quad)
1931 ++number_cache.n_quads;
1932 if (quad->has_children() ==
false)
1933 ++number_cache.n_active_quads;
1938 update_lines.
join();
1956 template <
int dim,
int spacedim>
1960 const unsigned int level_objects,
1971 &compute_number_cache_dim<dim, spacedim>),
1977 using hex_iterator =
1983 number_cache.n_hexes = 0;
1984 number_cache.n_active_hexes = 0;
1995 unsigned int n_levels = 0;
1996 if (level_objects > 0)
2000 n_levels =
level + 1;
2002 number_cache.n_hexes_level.resize(n_levels);
2003 number_cache.n_active_hexes_level.resize(n_levels);
2008 number_cache.n_hexes_level[
level] = 0;
2009 number_cache.n_active_hexes_level[
level] = 0;
2012 endc = (
level == n_levels - 1 ?
2015 for (; hex != endc; ++hex)
2017 ++number_cache.n_hexes_level[
level];
2018 if (hex->has_children() ==
false)
2019 ++number_cache.n_active_hexes_level[
level];
2023 number_cache.n_hexes += number_cache.n_hexes_level[
level];
2024 number_cache.n_active_hexes +=
2025 number_cache.n_active_hexes_level[
level];
2031 number_cache.n_hexes_level.clear();
2032 number_cache.n_active_hexes_level.clear();
2036 for (; hex != endc; ++hex)
2038 ++number_cache.n_hexes;
2039 if (hex->has_children() ==
false)
2040 ++number_cache.n_active_hexes;
2045 update_quads_and_lines.
join();
2049 template <
int dim,
int spacedim>
2053 const unsigned int level_objects,
2058 number_cache.active_cell_index_partitioner =
2059 std::make_shared<const Utilities::MPI::Partitioner>(
2062 number_cache.level_cell_index_partitioners.resize(
2065 number_cache.level_cell_index_partitioners[
level] =
2066 std::make_shared<const Utilities::MPI::Partitioner>(
2071 template <
int spacedim>
2077 template <
int dim,
int spacedim>
2132 static const unsigned int left_right_offset[2][6][2] = {
2141 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
2152 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
2156 for (
auto f : cell->face_indices())
2161 const unsigned int offset =
2162 (cell->direction_flag() ?
2163 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2165 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2177 if (cell->is_active() && face->has_children())
2207 if (face->has_children() &&
2208 (cell->is_active() ||
2210 cell->refinement_case(), f) ==
2213 for (
unsigned int c = 0; c < face->n_children(); ++c)
2216 if (face->child(0)->has_children())
2223 if (face->child(1)->has_children())
2240 for (
auto f : cell->face_indices())
2242 const unsigned int offset =
2243 (cell->direction_flag() ?
2244 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2246 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2256 template <
int dim,
int spacedim>
2270 if (dim == spacedim)
2271 for (
unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2279 const double cell_measure = GridTools::cell_measure<spacedim>(
2302 const auto connectivity = build_connectivity<unsigned int>(cells);
2303 const unsigned int n_cell = cells.size();
2311 const auto & crs = connectivity.entity_to_entities(1, 0);
2312 const unsigned int n_lines = crs.ptr.size() - 1;
2318 for (
unsigned int line = 0; line < n_lines; ++line)
2319 for (
unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
2332 const auto & crs = connectivity.entity_to_entities(2, 1);
2333 const unsigned int n_quads = crs.ptr.size() - 1;
2340 for (
unsigned int q = 0, k = 0; q < n_quads; ++q)
2343 faces.set_quad_type(q, connectivity.entity_types(2)[q]);
2346 for (
unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
2354 const unsigned char combined_orientation =
2355 connectivity.entity_orientations(1)
2356 .get_combined_orientation(k);
2360 combined_orientation ==
2362 combined_orientation ==
2365 faces.quads_line_orientations
2367 combined_orientation ==
2379 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2380 const auto &nei = connectivity.entity_to_entities(dim, dim);
2384 bool orientation_needed =
false;
2386 orientation_needed =
true;
2389 const auto &orientations = connectivity.entity_orientations(1);
2390 for (
unsigned int i = 0; i < orientations.n_objects(); ++i)
2391 if (orientations.get_combined_orientation(i) !=
2394 orientation_needed =
true;
2404 for (
unsigned int cell = 0; cell < n_cell; ++cell)
2407 cells_0.boundary_or_material_id[cell].material_id =
2408 cells[cell].material_id;
2411 cells_0.manifold_id[cell] = cells[cell].manifold_id;
2414 level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
2417 for (
unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2421 if (nei.col[i] !=
static_cast<unsigned int>(-1))
2423 j] = {0, nei.col[i]};
2430 if (orientation_needed)
2432 level.face_orientations.set_combined_orientation(
2434 connectivity.entity_orientations(dim - 1)
2435 .get_combined_orientation(i));
2444 auto &bids_face = dim == 3 ?
2445 tria.
faces->quads.boundary_or_material_id :
2446 tria.
faces->lines.boundary_or_material_id;
2449 std::vector<unsigned int> count(bids_face.size(), 0);
2452 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2455 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
2456 for (
unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
2457 count[crs.col[i]]++;
2460 for (
unsigned int face = 0; face < count.size(); ++face)
2462 if (count[face] != 1)
2466 bids_face[face].boundary_id = 0;
2472 const auto &crs = connectivity.entity_to_entities(2, 1);
2473 for (
unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
2474 tria.
faces->lines.boundary_or_material_id[crs.col[i]]
2480 static const unsigned int t_tba =
static_cast<unsigned int>(-1);
2481 static const unsigned int t_inner =
static_cast<unsigned int>(-2);
2483 std::vector<unsigned int> type(
vertices.size(), t_tba);
2485 const auto &crs = connectivity.entity_to_entities(1, 0);
2487 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
2488 for (
unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2490 if (type[crs.col[i]] != t_inner)
2491 type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
2493 for (
unsigned int face = 0; face < type.size(); ++face)
2498 if (type[face] != t_inner && type[face] != t_tba)
2519 template <
int structdim,
int spacedim,
typename T>
2529 if (boundary_objects_in.size() == 0)
2533 auto boundary_objects = boundary_objects_in;
2536 for (
auto &boundary_object : boundary_objects)
2537 std::sort(boundary_object.vertices.begin(),
2538 boundary_object.vertices.end());
2541 std::sort(boundary_objects.begin(),
2542 boundary_objects.end(),
2543 [](
const auto &a,
const auto &b) {
2544 return a.vertices < b.vertices;
2547 unsigned int counter = 0;
2549 std::vector<unsigned int> key;
2552 for (
unsigned int o = 0; o < obj.
n_objects(); ++o)
2566 key.assign(crs.
col.data() + crs.
ptr[o],
2567 crs.
col.data() + crs.
ptr[o + 1]);
2568 std::sort(key.begin(), key.end());
2571 const auto subcell_object =
2572 std::lower_bound(boundary_objects.begin(),
2573 boundary_objects.end(),
2575 [&](
const auto &cell,
const auto &key) {
2576 return cell.vertices < key;
2580 if (subcell_object == boundary_objects.end() ||
2581 subcell_object->vertices != key)
2587 manifold_id = subcell_object->manifold_id;
2590 if (subcell_object->boundary_id !=
2593 (void)vertex_locations;
2597 "The input arguments for creating a triangulation "
2598 "specified a boundary id for an internal face. This "
2601 "The object in question has vertex indices " +
2602 [subcell_object]() {
2604 for (
const auto v : subcell_object->vertices)
2605 s += std::to_string(v) +
',';
2608 " which are located at positions " +
2609 [vertex_locations, subcell_object]() {
2610 std::ostringstream s;
2611 for (
const auto v : subcell_object->vertices)
2612 s <<
'(' << vertex_locations[v] <<
')';
2616 boundary_id = subcell_object->boundary_id;
2630 const unsigned structdim,
2631 const unsigned int size)
2633 const unsigned int dim = faces.
dim;
2635 const unsigned int max_lines_per_face = 2 * structdim;
2637 if (dim == 3 && structdim == 2)
2652 const unsigned int spacedim,
2653 const unsigned int size,
2654 const bool orientation_needed)
2656 const unsigned int dim =
level.dim;
2658 const unsigned int max_faces_per_cell = 2 * dim;
2660 level.active_cell_indices.assign(size, -1);
2661 level.subdomain_ids.assign(size, 0);
2662 level.level_subdomain_ids.assign(size, 0);
2664 level.refine_flags.assign(size, 0u);
2665 level.coarsen_flags.assign(size,
false);
2667 level.parents.assign((size + 1) / 2, -1);
2670 level.direction_flags.assign(size,
true);
2672 level.neighbors.assign(size * max_faces_per_cell, {-1, -1});
2676 if (orientation_needed)
2677 level.face_orientations.reinit(size * max_faces_per_cell);
2680 level.global_active_cell_indices.assign(size,
2682 level.global_level_cell_indices.assign(size,
2691 const unsigned int structdim = obj.
structdim;
2693 const unsigned int max_children_per_cell = 1 << structdim;
2694 const unsigned int max_faces_per_cell = 2 * structdim;
2696 obj.
used.assign(size,
true);
2700 BoundaryOrMaterialId());
2708 obj.
children.assign(max_children_per_cell / 2 * size, -1);
2710 obj.
cells.assign(max_faces_per_cell * size, -1);
2740 template <
int spacedim>
2744 std::vector<unsigned int> &,
2745 std::vector<unsigned int> &)
2747 const unsigned int dim = 1;
2759 Assert(!cell->child(0)->has_children() &&
2760 !cell->child(1)->has_children(),
2766 if (cell->neighbor(0)->has_children())
2773 neighbor = neighbor->child(1);
2776 Assert(neighbor->neighbor(1) == cell->child(0),
2778 neighbor->set_neighbor(1, cell);
2784 if (neighbor->has_children())
2785 neighbor = neighbor->child(1);
2794 if (cell->neighbor(1)->has_children())
2801 neighbor = neighbor->child(0);
2804 Assert(neighbor->neighbor(0) == cell->child(1),
2806 neighbor->set_neighbor(0, cell);
2812 if (neighbor->has_children())
2813 neighbor = neighbor->child(0);
2823 triangulation.vertices_used[cell->child(0)->vertex_index(1)] =
false;
2829 for (
unsigned int child = 0; child < cell->n_children(); ++child)
2831 cell->child(child)->clear_user_data();
2832 cell->child(child)->clear_user_flag();
2833 cell->child(child)->clear_used_flag();
2838 cell->clear_children();
2839 cell->clear_user_flag();
2844 template <
int spacedim>
2848 std::vector<unsigned int> &line_cell_count,
2849 std::vector<unsigned int> &)
2851 const unsigned int dim = 2;
2859 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2862 lines_to_delete.reserve(4 * 2 + 4);
2867 for (
unsigned int c = 0; c < cell->n_children(); ++c)
2871 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2872 --line_cell_count[child->line_index(l)];
2887 .vertices_used[cell->child(0)->line(1)->vertex_index(1)] =
false;
2889 lines_to_delete.push_back(cell->child(0)->line(1));
2890 lines_to_delete.push_back(cell->child(0)->line(3));
2891 lines_to_delete.push_back(cell->child(3)->line(0));
2892 lines_to_delete.push_back(cell->child(3)->line(2));
2896 unsigned int inner_face_no =
2901 lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2905 for (
unsigned int child = 0; child < cell->n_children(); ++child)
2907 cell->child(child)->clear_user_data();
2908 cell->child(child)->clear_user_flag();
2909 cell->child(child)->clear_used_flag();
2914 cell->clear_children();
2915 cell->clear_refinement_case();
2916 cell->clear_user_flag();
2922 for (
unsigned int line_no = 0;
2923 line_no < GeometryInfo<dim>::lines_per_cell;
2927 cell->line(line_no);
2929 if (line->has_children())
2934 Assert((line_cell_count[line->child_index(0)] == 0 &&
2935 line_cell_count[line->child_index(1)] == 0) ||
2936 (line_cell_count[line->child_index(0)] > 0 &&
2937 line_cell_count[line->child_index(1)] > 0),
2940 if (line_cell_count[line->child_index(0)] == 0)
2942 for (
unsigned int c = 0; c < 2; ++c)
2943 Assert(!line->child(c)->has_children(),
2951 .vertices_used[line->child(0)->vertex_index(1)] =
false;
2953 lines_to_delete.push_back(line->child(0));
2954 lines_to_delete.push_back(line->child(1));
2956 line->clear_children();
2968 typename std::vector<
2970 line = lines_to_delete.
begin(),
2971 endline = lines_to_delete.end();
2972 for (; line != endline; ++line)
2974 (*line)->clear_user_data();
2975 (*line)->clear_user_flag();
2976 (*line)->clear_used_flag();
2982 template <
int spacedim>
2986 std::vector<unsigned int> &line_cell_count,
2987 std::vector<unsigned int> &quad_cell_count)
2989 const unsigned int dim = 3;
3001 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3003 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3006 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3007 quads_to_delete.reserve(6 * 4 + 12);
3011 for (
unsigned int c = 0; c < cell->n_children(); ++c)
3015 const auto line_indices = TriaAccessorImplementation::
3016 Implementation::get_line_indices_of_cell(*child);
3017 for (
const unsigned int l : cell->line_indices())
3018 --line_cell_count[line_indices[l]];
3020 --quad_cell_count[child->quad_index(f)];
3034 quads_to_delete.push_back(cell->child(0)->face(1));
3037 quads_to_delete.push_back(cell->child(0)->face(3));
3040 quads_to_delete.push_back(cell->child(0)->face(5));
3043 quads_to_delete.push_back(cell->child(0)->face(1));
3044 quads_to_delete.push_back(cell->child(0)->face(3));
3045 quads_to_delete.push_back(cell->child(3)->face(0));
3046 quads_to_delete.push_back(cell->child(3)->face(2));
3048 lines_to_delete.push_back(cell->child(0)->line(11));
3051 quads_to_delete.push_back(cell->child(0)->face(1));
3052 quads_to_delete.push_back(cell->child(0)->face(5));
3053 quads_to_delete.push_back(cell->child(3)->face(0));
3054 quads_to_delete.push_back(cell->child(3)->face(4));
3056 lines_to_delete.push_back(cell->child(0)->line(5));
3059 quads_to_delete.push_back(cell->child(0)->face(3));
3060 quads_to_delete.push_back(cell->child(0)->face(5));
3061 quads_to_delete.push_back(cell->child(3)->face(2));
3062 quads_to_delete.push_back(cell->child(3)->face(4));
3064 lines_to_delete.push_back(cell->child(0)->line(7));
3067 quads_to_delete.push_back(cell->child(0)->face(1));
3068 quads_to_delete.push_back(cell->child(2)->face(1));
3069 quads_to_delete.push_back(cell->child(4)->face(1));
3070 quads_to_delete.push_back(cell->child(6)->face(1));
3072 quads_to_delete.push_back(cell->child(0)->face(3));
3073 quads_to_delete.push_back(cell->child(1)->face(3));
3074 quads_to_delete.push_back(cell->child(4)->face(3));
3075 quads_to_delete.push_back(cell->child(5)->face(3));
3077 quads_to_delete.push_back(cell->child(0)->face(5));
3078 quads_to_delete.push_back(cell->child(1)->face(5));
3079 quads_to_delete.push_back(cell->child(2)->face(5));
3080 quads_to_delete.push_back(cell->child(3)->face(5));
3082 lines_to_delete.push_back(cell->child(0)->line(5));
3083 lines_to_delete.push_back(cell->child(0)->line(7));
3084 lines_to_delete.push_back(cell->child(0)->line(11));
3085 lines_to_delete.push_back(cell->child(7)->line(0));
3086 lines_to_delete.push_back(cell->child(7)->line(2));
3087 lines_to_delete.push_back(cell->child(7)->line(8));
3093 triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3105 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3107 cell->child(child)->clear_user_data();
3108 cell->child(child)->clear_user_flag();
3113 cell->child(child)->set_combined_face_orientation(
3116 cell->child(child)->clear_used_flag();
3121 cell->clear_children();
3122 cell->clear_refinement_case();
3123 cell->clear_user_flag();
3136 cell->face(quad_no);
3140 quad->has_children()) ||
3145 switch (quad->refinement_case())
3157 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3158 quad_cell_count[quad->child_index(1)] == 0) ||
3159 (quad_cell_count[quad->child_index(0)] > 0 &&
3160 quad_cell_count[quad->child_index(1)] > 0),
3166 unsigned int deleted_grandchildren = 0;
3167 unsigned int number_of_child_refinements = 0;
3169 for (
unsigned int c = 0; c < 2; ++c)
3170 if (quad->child(c)->has_children())
3172 ++number_of_child_refinements;
3177 (quad_cell_count[quad->child(c)->child_index(0)] ==
3179 quad_cell_count[quad->child(c)->child_index(1)] ==
3181 (quad_cell_count[quad->child(c)->child_index(0)] >
3183 quad_cell_count[quad->child(c)->child_index(1)] >
3186 if (quad_cell_count[quad->child(c)->child_index(0)] ==
3193 Assert(quad->refinement_case() +
3194 quad->child(c)->refinement_case() ==
3202 quads_to_delete.push_back(
3203 quad->child(c)->child(0));
3204 quads_to_delete.push_back(
3205 quad->child(c)->child(1));
3206 if (quad->child(c)->refinement_case() ==
3208 lines_to_delete.push_back(
3209 quad->child(c)->child(0)->line(1));
3211 lines_to_delete.push_back(
3212 quad->child(c)->child(0)->line(3));
3213 quad->child(c)->clear_children();
3214 quad->child(c)->clear_refinement_case();
3215 ++deleted_grandchildren;
3223 if (number_of_child_refinements > 0 &&
3224 deleted_grandchildren == number_of_child_refinements)
3229 middle_line = quad->child(0)->line(1);
3231 middle_line = quad->child(0)->line(3);
3233 lines_to_delete.push_back(middle_line->child(0));
3234 lines_to_delete.push_back(middle_line->child(1));
3236 .vertices_used[middle_vertex_index<dim, spacedim>(
3237 middle_line)] =
false;
3238 middle_line->clear_children();
3243 if (quad_cell_count[quad->child_index(0)] == 0)
3249 quads_to_delete.push_back(quad->child(0));
3250 quads_to_delete.push_back(quad->child(1));
3252 lines_to_delete.push_back(quad->child(0)->line(1));
3254 lines_to_delete.push_back(quad->child(0)->line(3));
3274 if (quad->child(0)->has_children())
3275 if (quad->refinement_case() ==
3309 quad_iterator switch_1 =
3310 quad->child(0)->child(1),
3312 quad->child(1)->child(0);
3314 Assert(!switch_1->has_children(),
3316 Assert(!switch_2->has_children(),
3319 const int switch_1_index = switch_1->index();
3320 const int switch_2_index = switch_2->index();
3321 for (
unsigned int l = 0;
3324 for (
unsigned int h = 0;
3328 for (
const unsigned int q :
3333 ->cells.get_bounding_object_indices(
3335 if (
index == switch_1_index)
3337 ->cells.get_bounding_object_indices(
3338 h)[q] = switch_2_index;
3339 else if (
index == switch_2_index)
3341 ->cells.get_bounding_object_indices(
3342 h)[q] = switch_1_index;
3347 const int switch_1_lines[4] = {
3348 static_cast<signed int>(
3349 switch_1->line_index(0)),
3350 static_cast<signed int>(
3351 switch_1->line_index(1)),
3352 static_cast<signed int>(
3353 switch_1->line_index(2)),
3354 static_cast<signed int>(
3355 switch_1->line_index(3))};
3356 const bool switch_1_line_orientations[4] = {
3357 switch_1->line_orientation(0),
3358 switch_1->line_orientation(1),
3359 switch_1->line_orientation(2),
3360 switch_1->line_orientation(3)};
3362 switch_1->boundary_id();
3363 const unsigned int switch_1_user_index =
3364 switch_1->user_index();
3365 const bool switch_1_user_flag =
3366 switch_1->user_flag_set();
3368 switch_1->set_bounding_object_indices(
3369 {switch_2->line_index(0),
3370 switch_2->line_index(1),
3371 switch_2->line_index(2),
3372 switch_2->line_index(3)});
3373 switch_1->set_line_orientation(
3374 0, switch_2->line_orientation(0));
3375 switch_1->set_line_orientation(
3376 1, switch_2->line_orientation(1));
3377 switch_1->set_line_orientation(
3378 2, switch_2->line_orientation(2));
3379 switch_1->set_line_orientation(
3380 3, switch_2->line_orientation(3));
3381 switch_1->set_boundary_id_internal(
3382 switch_2->boundary_id());
3383 switch_1->set_manifold_id(
3384 switch_2->manifold_id());
3385 switch_1->set_user_index(switch_2->user_index());
3386 if (switch_2->user_flag_set())
3387 switch_1->set_user_flag();
3389 switch_1->clear_user_flag();
3391 switch_2->set_bounding_object_indices(
3395 switch_1_lines[3]});
3396 switch_2->set_line_orientation(
3397 0, switch_1_line_orientations[0]);
3398 switch_2->set_line_orientation(
3399 1, switch_1_line_orientations[1]);
3400 switch_2->set_line_orientation(
3401 2, switch_1_line_orientations[2]);
3402 switch_2->set_line_orientation(
3403 3, switch_1_line_orientations[3]);
3404 switch_2->set_boundary_id_internal(
3405 switch_1_boundary_id);
3406 switch_2->set_manifold_id(
3407 switch_1->manifold_id());
3408 switch_2->set_user_index(switch_1_user_index);
3409 if (switch_1_user_flag)
3410 switch_2->set_user_flag();
3412 switch_2->clear_user_flag();
3414 const unsigned int child_0 =
3415 quad->child(0)->child_index(0);
3416 const unsigned int child_2 =
3417 quad->child(1)->child_index(0);
3418 quad->clear_children();
3419 quad->clear_refinement_case();
3420 quad->set_refinement_case(
3422 quad->set_children(0, child_0);
3423 quad->set_children(2, child_2);
3424 std::swap(quad_cell_count[child_0 + 1],
3425 quad_cell_count[child_2]);
3440 const unsigned int child_0 =
3441 quad->child(0)->child_index(0);
3442 const unsigned int child_2 =
3443 quad->child(1)->child_index(0);
3444 quad->clear_children();
3445 quad->clear_refinement_case();
3446 quad->set_refinement_case(
3448 quad->set_children(0, child_0);
3449 quad->set_children(2, child_2);
3453 quad->clear_children();
3454 quad->clear_refinement_case();
3465 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3466 quad_cell_count[quad->child_index(1)] == 0 &&
3467 quad_cell_count[quad->child_index(2)] == 0 &&
3468 quad_cell_count[quad->child_index(3)] == 0) ||
3469 (quad_cell_count[quad->child_index(0)] > 0 &&
3470 quad_cell_count[quad->child_index(1)] > 0 &&
3471 quad_cell_count[quad->child_index(2)] > 0 &&
3472 quad_cell_count[quad->child_index(3)] > 0),
3475 if (quad_cell_count[quad->child_index(0)] == 0)
3481 lines_to_delete.push_back(quad->child(0)->line(1));
3482 lines_to_delete.push_back(quad->child(3)->line(0));
3483 lines_to_delete.push_back(quad->child(0)->line(3));
3484 lines_to_delete.push_back(quad->child(3)->line(2));
3486 for (
unsigned int child = 0; child < quad->n_children();
3488 quads_to_delete.push_back(quad->child(child));
3491 .vertices_used[quad->child(0)->vertex_index(3)] =
3494 quad->clear_children();
3495 quad->clear_refinement_case();
3514 for (
unsigned int line_no = 0;
3515 line_no < GeometryInfo<dim>::lines_per_cell;
3519 cell->line(line_no);
3523 line->has_children()) ||
3528 if (line->has_children())
3533 Assert((line_cell_count[line->child_index(0)] == 0 &&
3534 line_cell_count[line->child_index(1)] == 0) ||
3535 (line_cell_count[line->child_index(0)] > 0 &&
3536 line_cell_count[line->child_index(1)] > 0),
3539 if (line_cell_count[line->child_index(0)] == 0)
3541 for (
unsigned int c = 0; c < 2; ++c)
3542 Assert(!line->child(c)->has_children(),
3550 .vertices_used[line->child(0)->vertex_index(1)] =
false;
3552 lines_to_delete.push_back(line->child(0));
3553 lines_to_delete.push_back(line->child(1));
3555 line->clear_children();
3567 typename std::vector<
3569 line = lines_to_delete.
begin(),
3570 endline = lines_to_delete.end();
3571 for (; line != endline; ++line)
3573 (*line)->clear_user_data();
3574 (*line)->clear_user_flag();
3575 (*line)->clear_used_flag();
3578 typename std::vector<
3580 quad = quads_to_delete.
begin(),
3581 endquad = quads_to_delete.end();
3582 for (; quad != endquad; ++quad)
3584 (*quad)->clear_user_data();
3585 (*quad)->clear_children();
3586 (*quad)->clear_refinement_case();
3587 (*quad)->clear_user_flag();
3588 (*quad)->clear_used_flag();
3610 template <
int spacedim>
3614 unsigned int & next_unused_vertex,
3621 const unsigned int dim = 2;
3624 cell->clear_refine_flag();
3717 int new_vertices[9];
3718 for (
unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3719 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3720 for (
unsigned int line_no = 0; line_no < 4; ++line_no)
3721 if (cell->line(line_no)->has_children())
3722 new_vertices[4 + line_no] =
3723 cell->line(line_no)->child(0)->vertex_index(1);
3732 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
3733 ++next_unused_vertex;
3736 "Internal error: During refinement, the triangulation "
3737 "wants to access an element of the 'vertices' array "
3738 "but it turns out that the array is not large enough."));
3741 new_vertices[8] = next_unused_vertex;
3747 cell->center(
true,
true);
3753 unsigned int lmin = 8;
3754 unsigned int lmax = 12;
3761 for (
unsigned int l = lmin; l < lmax; ++l)
3763 while (next_unused_line->used() ==
true)
3765 new_lines[l] = next_unused_line;
3783 for (
unsigned int c = 0; c < 2; ++c, ++l)
3784 new_lines[l] = cell->line(face_no)->child(c);
3787 new_lines[8]->set_bounding_object_indices(
3788 {new_vertices[6], new_vertices[8]});
3789 new_lines[9]->set_bounding_object_indices(
3790 {new_vertices[8], new_vertices[7]});
3791 new_lines[10]->set_bounding_object_indices(
3792 {new_vertices[4], new_vertices[8]});
3793 new_lines[11]->set_bounding_object_indices(
3794 {new_vertices[8], new_vertices[5]});
3803 new_lines[0] = cell->line(0);
3804 new_lines[1] = cell->line(1);
3805 new_lines[2] = cell->line(2)->child(0);
3806 new_lines[3] = cell->line(2)->child(1);
3807 new_lines[4] = cell->line(3)->child(0);
3808 new_lines[5] = cell->line(3)->child(1);
3809 new_lines[6]->set_bounding_object_indices(
3810 {new_vertices[6], new_vertices[7]});
3820 new_lines[0] = cell->line(0)->child(0);
3821 new_lines[1] = cell->line(0)->child(1);
3822 new_lines[2] = cell->line(1)->child(0);
3823 new_lines[3] = cell->line(1)->child(1);
3824 new_lines[4] = cell->line(2);
3825 new_lines[5] = cell->line(3);
3826 new_lines[6]->set_bounding_object_indices(
3827 {new_vertices[4], new_vertices[5]});
3830 for (
unsigned int l = lmin; l < lmax; ++l)
3832 new_lines[l]->set_used_flag();
3833 new_lines[l]->clear_user_flag();
3835 new_lines[l]->clear_children();
3837 new_lines[l]->set_boundary_id_internal(
3839 new_lines[l]->set_manifold_id(cell->manifold_id());
3846 while (next_unused_cell->used() ==
true)
3850 for (
unsigned int i = 0; i < n_children; ++i)
3853 subcells[i] = next_unused_cell;
3855 if (i % 2 == 1 && i < n_children - 1)
3856 while (next_unused_cell->used() ==
true)
3874 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3875 new_lines[8]->
index(),
3876 new_lines[4]->
index(),
3877 new_lines[10]->
index()});
3878 subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3879 new_lines[2]->
index(),
3880 new_lines[5]->
index(),
3881 new_lines[11]->
index()});
3882 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3883 new_lines[9]->
index(),
3884 new_lines[10]->
index(),
3885 new_lines[6]->
index()});
3886 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3887 new_lines[3]->
index(),
3888 new_lines[11]->
index(),
3889 new_lines[7]->
index()});
3905 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3906 new_lines[6]->
index(),
3907 new_lines[2]->
index(),
3908 new_lines[4]->
index()});
3909 subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3910 new_lines[1]->
index(),
3911 new_lines[3]->
index(),
3912 new_lines[5]->
index()});
3929 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3930 new_lines[2]->
index(),
3931 new_lines[4]->
index(),
3932 new_lines[6]->
index()});
3933 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
3934 new_lines[3]->
index(),
3935 new_lines[6]->
index(),
3936 new_lines[5]->
index()});
3941 for (
unsigned int i = 0; i < n_children; ++i)
3943 subcells[i]->set_used_flag();
3944 subcells[i]->clear_refine_flag();
3945 subcells[i]->clear_user_flag();
3946 subcells[i]->clear_user_data();
3947 subcells[i]->clear_children();
3949 subcells[i]->set_material_id(cell->material_id());
3950 subcells[i]->set_manifold_id(cell->manifold_id());
3951 subcells[i]->set_subdomain_id(subdomainid);
3954 subcells[i]->set_parent(cell->index());
3960 for (
unsigned int i = 0; i < n_children / 2; ++i)
3961 cell->set_children(2 * i, subcells[2 * i]->index());
3963 cell->set_refinement_case(ref_case);
3971 for (
unsigned int c = 0; c < n_children; ++c)
3972 cell->child(c)->set_direction_flag(cell->direction_flag());
3977 template <
int dim,
int spacedim>
3980 const bool check_for_distorted_cells)
3986 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
3988 if (cell->refine_flag_set())
4001 line->clear_user_flag();
4005 unsigned int n_single_lines = 0;
4006 unsigned int n_lines_in_pairs = 0;
4007 unsigned int needed_vertices = 0;
4013 unsigned int needed_cells = 0;
4015 for (
const auto &cell :
4017 if (cell->refine_flag_set())
4022 needed_vertices += 0;
4023 n_single_lines += 3;
4025 else if (cell->reference_cell() ==
4029 needed_vertices += 1;
4030 n_single_lines += 4;
4037 for (
const auto line_no : cell->face_indices())
4039 auto line = cell->line(line_no);
4040 if (line->has_children() ==
false)
4041 line->set_user_flag();
4046 const unsigned int used_cells =
4053 used_cells + needed_cells,
4065 if (line->user_flag_set())
4068 n_lines_in_pairs += 2;
4069 needed_vertices += 1;
4074 needed_vertices += std::count(
triangulation.vertices_used.begin(),
4084 unsigned int next_unused_vertex = 0;
4093 for (; line != endl; ++line)
4094 if (line->user_flag_set())
4100 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
4101 ++next_unused_vertex;
4105 "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."));
4108 triangulation.vertices[next_unused_vertex] = line->center(
true);
4110 bool pair_found =
false;
4112 for (; next_unused_line != endl; ++next_unused_line)
4113 if (!next_unused_line->used() &&
4114 !(++next_unused_line)->used())
4122 line->set_children(0, next_unused_line->index());
4125 children[2] = {next_unused_line, ++next_unused_line};
4130 children[0]->set_bounding_object_indices(
4131 {line->vertex_index(0), next_unused_vertex});
4132 children[1]->set_bounding_object_indices(
4133 {next_unused_vertex, line->vertex_index(1)});
4135 children[0]->set_used_flag();
4136 children[1]->set_used_flag();
4137 children[0]->clear_children();
4138 children[1]->clear_children();
4141 children[0]->clear_user_flag();
4142 children[1]->clear_user_flag();
4145 children[0]->set_boundary_id_internal(line->boundary_id());
4146 children[1]->set_boundary_id_internal(line->boundary_id());
4148 children[0]->set_manifold_id(line->manifold_id());
4149 children[1]->set_manifold_id(line->manifold_id());
4151 line->clear_user_flag();
4158 cells_with_distorted_children;
4164 unsigned int &next_unused_vertex,
4165 auto & next_unused_line,
4166 auto & next_unused_cell,
4167 const auto & cell) {
4168 const auto ref_case = cell->refine_flag_set();
4169 cell->clear_refine_flag();
4171 unsigned int n_new_vertices = 0;
4180 std::vector<unsigned int> new_vertices(n_new_vertices,
4182 for (
unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
4184 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4185 for (
unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
4186 if (cell->line(line_no)->has_children())
4187 new_vertices[cell->n_vertices() + line_no] =
4188 cell->line(line_no)->child(0)->vertex_index(1);
4192 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
4193 ++next_unused_vertex;
4197 "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."));
4200 new_vertices[8] = next_unused_vertex;
4203 cell->center(
true,
true);
4206 std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
4209 unsigned int lmin = 0;
4210 unsigned int lmax = 0;
4227 for (
unsigned int l = lmin; l < lmax; ++l)
4229 while (next_unused_line->used() ==
true)
4231 new_lines[l] = next_unused_line;
4242 const auto ref = [&](
const unsigned int face_no,
4243 const unsigned int vertex_no) {
4244 auto l = cell->line(face_no);
4247 if (l->child(0)->vertex_index(0) == new_vertices[vertex_no] ||
4248 l->child(0)->vertex_index(1) == new_vertices[vertex_no])
4250 new_lines[2 * face_no + 0] = l->child(0);
4251 new_lines[2 * face_no + 1] = l->child(1);
4255 new_lines[2 * face_no + 0] = l->child(1);
4256 new_lines[2 * face_no + 1] = l->child(0);
4265 new_lines[6]->set_bounding_object_indices(
4266 {new_vertices[3], new_vertices[4]});
4267 new_lines[7]->set_bounding_object_indices(
4268 {new_vertices[4], new_vertices[5]});
4269 new_lines[8]->set_bounding_object_indices(
4270 {new_vertices[5], new_vertices[3]});
4275 for (
const unsigned int face_no : cell->face_indices())
4276 for (
unsigned int c = 0; c < 2; ++c, ++l)
4277 new_lines[l] = cell->line(face_no)->child(c);
4279 new_lines[8]->set_bounding_object_indices(
4280 {new_vertices[6], new_vertices[8]});
4281 new_lines[9]->set_bounding_object_indices(
4282 {new_vertices[8], new_vertices[7]});
4283 new_lines[10]->set_bounding_object_indices(
4284 {new_vertices[4], new_vertices[8]});
4285 new_lines[11]->set_bounding_object_indices(
4286 {new_vertices[8], new_vertices[5]});
4293 for (
unsigned int l = lmin; l < lmax; ++l)
4295 new_lines[l]->set_used_flag();
4296 new_lines[l]->clear_user_flag();
4297 new_lines[l]->clear_user_data();
4298 new_lines[l]->clear_children();
4300 new_lines[l]->set_boundary_id_internal(
4302 new_lines[l]->set_manifold_id(cell->manifold_id());
4307 while (next_unused_cell->used() ==
true)
4310 unsigned int n_children = 0;
4319 for (
unsigned int i = 0; i < n_children; ++i)
4322 subcells[i] = next_unused_cell;
4324 if (i % 2 == 1 && i < n_children - 1)
4325 while (next_unused_cell->used() ==
true)
4332 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4333 new_lines[8]->
index(),
4334 new_lines[5]->
index()});
4335 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4336 new_lines[2]->
index(),
4337 new_lines[6]->
index()});
4338 subcells[2]->set_bounding_object_indices({new_lines[7]->index(),
4339 new_lines[3]->
index(),
4340 new_lines[4]->
index()});
4341 subcells[3]->set_bounding_object_indices({new_lines[6]->index(),
4342 new_lines[7]->
index(),
4343 new_lines[8]->
index()});
4348 const auto fix_line_orientation =
4349 [&](
const unsigned int line_no,
4350 const unsigned int vertex_no,
4351 const unsigned int subcell_no,
4352 const unsigned int subcell_line_no) {
4353 if (new_lines[line_no]->vertex_index(1) !=
4354 new_vertices[vertex_no])
4356 ->face_orientations.set_combined_orientation(
4357 subcells[subcell_no]->
index() *
4363 fix_line_orientation(0, 3, 0, 0);
4364 fix_line_orientation(8, 5, 0, 1);
4365 fix_line_orientation(5, 0, 0, 2);
4367 fix_line_orientation(1, 1, 1, 0);
4368 fix_line_orientation(2, 4, 1, 1);
4369 fix_line_orientation(6, 3, 1, 2);
4371 fix_line_orientation(7, 4, 2, 0);
4372 fix_line_orientation(3, 2, 2, 1);
4373 fix_line_orientation(4, 5, 2, 2);
4377 fix_line_orientation(6, 4, 3, 0);
4378 fix_line_orientation(7, 5, 3, 1);
4379 fix_line_orientation(8, 3, 3, 2);
4381 else if ((dim == 2) &&
4384 subcells[0]->set_bounding_object_indices(
4385 {new_lines[0]->index(),
4386 new_lines[8]->
index(),
4387 new_lines[4]->
index(),
4388 new_lines[10]->
index()});
4389 subcells[1]->set_bounding_object_indices(
4390 {new_lines[8]->index(),
4391 new_lines[2]->
index(),
4392 new_lines[5]->
index(),
4393 new_lines[11]->
index()});
4394 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
4395 new_lines[9]->
index(),
4396 new_lines[10]->
index(),
4397 new_lines[6]->
index()});
4398 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
4399 new_lines[3]->
index(),
4400 new_lines[11]->
index(),
4401 new_lines[7]->
index()});
4410 for (
unsigned int i = 0; i < n_children; ++i)
4412 subcells[i]->set_used_flag();
4413 subcells[i]->clear_refine_flag();
4414 subcells[i]->clear_user_flag();
4415 subcells[i]->clear_user_data();
4416 subcells[i]->clear_children();
4419 subcells[i]->set_material_id(cell->material_id());
4420 subcells[i]->set_manifold_id(cell->manifold_id());
4421 subcells[i]->set_subdomain_id(subdomainid);
4426 ->reference_cell[subcells[i]->index()] = cell->reference_cell();
4429 subcells[i]->set_parent(cell->index());
4432 for (
unsigned int i = 0; i < n_children / 2; ++i)
4433 cell->set_children(2 * i, subcells[2 * i]->index());
4435 cell->set_refinement_case(ref_case);
4438 for (
unsigned int c = 0; c < n_children; ++c)
4439 cell->child(c)->set_direction_flag(cell->direction_flag());
4449 for (
const auto &cell :
4451 if (cell->refine_flag_set())
4460 check_for_distorted_cells &&
4461 has_distorted_children<dim, spacedim>(cell))
4469 return cells_with_distorted_children;
4478 template <
int spacedim>
4483 const unsigned int dim = 1;
4487 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
4489 if (cell->refine_flag_set())
4502 unsigned int needed_vertices = 0;
4507 unsigned int flagged_cells = 0;
4509 for (
const auto &acell :
4511 if (acell->refine_flag_set())
4516 const unsigned int used_cells =
4536 needed_vertices += flagged_cells;
4541 needed_vertices += std::count(
triangulation.vertices_used.begin(),
4558 unsigned int next_unused_vertex = 0;
4565 for (
const auto &cell :
4567 if (cell->refine_flag_set())
4570 cell->clear_refine_flag();
4576 ++next_unused_vertex;
4580 "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."));
4595 while (next_unused_cell->used() ==
true)
4597 first_child = next_unused_cell;
4598 first_child->set_used_flag();
4599 first_child->clear_user_data();
4602 second_child = next_unused_cell;
4603 second_child->set_used_flag();
4604 second_child->clear_user_data();
4609 cell->set_children(0, first_child->index());
4610 first_child->clear_children();
4611 first_child->set_bounding_object_indices(
4612 {cell->vertex_index(0), next_unused_vertex});
4613 first_child->set_material_id(cell->material_id());
4614 first_child->set_manifold_id(cell->manifold_id());
4615 first_child->set_subdomain_id(subdomainid);
4616 first_child->set_direction_flag(cell->direction_flag());
4618 first_child->set_parent(cell->index());
4622 first_child->face(1)->set_manifold_id(cell->manifold_id());
4627 first_child->set_neighbor(1, second_child);
4629 first_child->set_neighbor(0, cell->neighbor(0));
4630 else if (cell->neighbor(0)->is_active())
4636 Assert(cell->neighbor(0)->level() <= cell->level(),
4638 first_child->set_neighbor(0, cell->neighbor(0));
4644 const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4645 first_child->set_neighbor(0,
4646 cell->neighbor(0)->child(nbnb));
4651 left_neighbor = cell->neighbor(0);
4652 while (left_neighbor->has_children())
4654 left_neighbor = left_neighbor->child(nbnb);
4655 left_neighbor->set_neighbor(nbnb, first_child);
4660 second_child->clear_children();
4661 second_child->set_bounding_object_indices(
4662 {next_unused_vertex, cell->vertex_index(1)});
4663 second_child->set_neighbor(0, first_child);
4664 second_child->set_material_id(cell->material_id());
4665 second_child->set_manifold_id(cell->manifold_id());
4666 second_child->set_subdomain_id(subdomainid);
4667 second_child->set_direction_flag(cell->direction_flag());
4670 second_child->set_neighbor(1, cell->neighbor(1));
4671 else if (cell->neighbor(1)->is_active())
4673 Assert(cell->neighbor(1)->level() <= cell->level(),
4675 second_child->set_neighbor(1, cell->neighbor(1));
4680 const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4681 second_child->set_neighbor(
4682 1, cell->neighbor(1)->child(nbnb));
4685 right_neighbor = cell->neighbor(1);
4686 while (right_neighbor->has_children())
4688 right_neighbor = right_neighbor->child(nbnb);
4689 right_neighbor->set_neighbor(nbnb, second_child);
4709 template <
int spacedim>
4712 const bool check_for_distorted_cells)
4714 const unsigned int dim = 2;
4720 bool do_isotropic_refinement =
true;
4721 for (
const auto &cell :
triangulation.active_cell_iterators())
4725 do_isotropic_refinement =
false;
4729 if (do_isotropic_refinement)
4731 check_for_distorted_cells);
4736 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
4738 if (cell->refine_flag_set())
4753 line->clear_user_flag();
4759 unsigned int n_single_lines = 0;
4763 unsigned int n_lines_in_pairs = 0;
4769 unsigned int needed_vertices = 0;
4774 unsigned int needed_cells = 0;
4776 for (
const auto &cell :
4778 if (cell->refine_flag_set())
4789 n_single_lines += 4;
4801 n_single_lines += 1;
4809 for (
const unsigned int line_no :
4813 cell->refine_flag_set(), line_no) ==
4817 line = cell->line(line_no);
4818 if (line->has_children() ==
false)
4819 line->set_user_flag();
4826 const unsigned int used_cells =
4836 used_cells + needed_cells,
4852 if (line->user_flag_set())
4855 n_lines_in_pairs += 2;
4856 needed_vertices += 1;
4868 needed_vertices += std::count(
triangulation.vertices_used.begin(),
4885 unsigned int next_unused_vertex = 0;
4897 for (; line != endl; ++line)
4898 if (line->user_flag_set())
4904 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
4905 ++next_unused_vertex;
4909 "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."));
4912 triangulation.vertices[next_unused_vertex] = line->center(
true);
4917 bool pair_found =
false;
4919 for (; next_unused_line != endl; ++next_unused_line)
4920 if (!next_unused_line->used() &&
4921 !(++next_unused_line)->used())
4934 line->set_children(0, next_unused_line->index());
4938 children[2] = {next_unused_line, ++next_unused_line};
4944 children[0]->set_bounding_object_indices(
4945 {line->vertex_index(0), next_unused_vertex});
4946 children[1]->set_bounding_object_indices(
4947 {next_unused_vertex, line->vertex_index(1)});
4949 children[0]->set_used_flag();
4950 children[1]->set_used_flag();
4951 children[0]->clear_children();
4952 children[1]->clear_children();
4955 children[0]->clear_user_flag();
4956 children[1]->clear_user_flag();
4959 children[0]->set_boundary_id_internal(line->boundary_id());
4960 children[1]->set_boundary_id_internal(line->boundary_id());
4962 children[0]->set_manifold_id(line->manifold_id());
4963 children[1]->set_manifold_id(line->manifold_id());
4967 line->clear_user_flag();
4979 cells_with_distorted_children;
4993 for (
const auto &cell :
4995 if (cell->refine_flag_set())
5005 if (check_for_distorted_cells &&
5006 has_distorted_children<dim, spacedim>(cell))
5014 return cells_with_distorted_children;
5018 template <
int spacedim>
5021 const bool check_for_distorted_cells)
5023 static const int dim = 3;
5025 using raw_line_iterator =
5027 using raw_quad_iterator =
5038 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
5040 if (cell->refine_flag_set())
5066 unsigned int needed_vertices = 0;
5067 unsigned int needed_lines_single = 0;
5068 unsigned int needed_quads_single = 0;
5069 unsigned int needed_lines_pair = 0;
5070 unsigned int needed_quads_pair = 0;
5073 unsigned int new_cells = 0;
5075 for (
const auto &cell :
5077 if (cell->refine_flag_set())
5080 Assert(cell->refine_flag_set() ==
5086 new_cells += cell->reference_cell().n_isotropic_children();
5091 needed_lines_single += 6;
5092 needed_quads_single += 12;
5094 else if (cell->reference_cell() ==
5097 needed_lines_single += 1;
5098 needed_quads_single += 8;
5109 for (
const auto face : cell->face_indices())
5110 if (cell->face(face)->n_children() == 0)
5111 cell->face(face)->set_user_flag();
5113 Assert(cell->face(face)->n_children() ==
5114 cell->reference_cell()
5115 .face_reference_cell(face)
5116 .n_isotropic_children(),
5119 for (
const auto line : cell->line_indices())
5120 if (cell->line(line)->has_children() ==
false)
5121 cell->line(line)->set_user_flag();
5123 Assert(cell->line(line)->n_children() == 2,
5127 const unsigned int used_cells =
5133 used_cells + new_cells,
5147 if (quad->user_flag_set() ==
false)
5152 needed_quads_pair += 4;
5153 needed_lines_pair += 4;
5154 needed_vertices += 1;
5158 needed_quads_pair += 4;
5159 needed_lines_single += 3;
5172 if (line->user_flag_set() ==
false)
5175 needed_lines_pair += 2;
5176 needed_vertices += 1;
5181 needed_lines_single);
5184 needed_quads_single);
5187 needed_quads_single);
5191 needed_vertices += std::count(
triangulation.vertices_used.begin(),
5212 for (
const auto &cell :
triangulation.active_cell_iterators())
5213 if (!cell->refine_flag_set())
5214 for (
unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
5215 if (cell->line(line_n)->has_children())
5216 for (
unsigned int c = 0; c < 2; ++c)
5217 Assert(cell->line(line_n)->child(c)->user_flag_set() ==
false,
5221 unsigned int current_vertex = 0;
5225 auto get_next_unused_vertex = [](
const unsigned int current_vertex,
5226 std::vector<bool> &vertices_used) {
5227 unsigned int next_vertex = current_vertex;
5228 while (next_vertex < vertices_used.size() &&
5229 vertices_used[next_vertex] ==
true)
5232 vertices_used[next_vertex] =
true;
5242 raw_line_iterator next_unused_line =
triangulation.begin_raw_line();
5244 for (; line != endl; ++line)
5246 if (line->user_flag_set() ==
false)
5250 triangulation.faces->lines.template next_free_pair_object<1>(
5258 line->set_children(0, next_unused_line->index());
5260 const std::array<raw_line_iterator, 2> children{
5261 {next_unused_line, ++next_unused_line}};
5267 get_next_unused_vertex(current_vertex,
5269 triangulation.vertices[current_vertex] = line->center(
true);
5271 children[0]->set_bounding_object_indices(
5272 {line->vertex_index(0), current_vertex});
5273 children[1]->set_bounding_object_indices(
5274 {current_vertex, line->vertex_index(1)});
5276 const auto manifold_id = line->manifold_id();
5277 const auto boundary_id = line->boundary_id();
5278 for (
const auto &child : children)
5280 child->set_used_flag();
5281 child->clear_children();
5283 child->clear_user_flag();
5284 child->set_boundary_id_internal(boundary_id);
5285 child->set_manifold_id(manifold_id);
5288 line->clear_user_flag();
5298 for (; quad != endq; ++quad)
5300 if (quad->user_flag_set() ==
false)
5303 const auto reference_face_type = quad->reference_cell();
5307 std::array<raw_line_iterator, 4> new_lines;
5310 for (
unsigned int l = 0; l < 2; ++l)
5312 auto next_unused_line =
5315 new_lines[2 * l] = next_unused_line;
5316 new_lines[2 * l + 1] = ++next_unused_line;
5321 for (
unsigned int l = 0; l < 3; ++l)
5331 for (
const unsigned int line : quad->line_indices())
5339 std::array<raw_quad_iterator, 4> new_quads;
5340 for (
unsigned int q = 0; q < 2; ++q)
5342 auto next_unused_quad =
5346 new_quads[2 * q] = next_unused_quad;
5347 new_quads[2 * q + 1] = ++next_unused_quad;
5349 quad->set_children(2 * q, new_quads[2 * q]->
index());
5353 for (
const auto &quad : new_quads)
5365 for (
const auto i : quad->vertex_indices())
5368 for (
const auto i : quad->line_indices())
5374 get_next_unused_vertex(current_vertex,
5379 quad->center(
true,
true);
5383 std::array<raw_line_iterator, 12> lines;
5384 unsigned int n_lines = 0;
5385 for (
unsigned int l = 0; l < quad->
n_lines(); ++l)
5386 for (
unsigned int c = 0; c < 2; ++c)
5388 static constexpr ::ndarray<unsigned int, 2, 2>
index =
5395 quad->line(l)->child(
index[c][quad->line_orientation(l)]);
5398 for (
unsigned int l = 0; l < quad->
n_lines(); ++l)
5399 lines[n_lines++] = new_lines[l];
5401 std::array<int, 12> line_indices;
5402 for (
unsigned int i = 0; i < n_lines; ++i)
5403 line_indices[i] = lines[i]->
index();
5405 static constexpr ::ndarray<unsigned int, 12, 2>
5406 line_vertices_quad{{{{0, 4}},
5419 static constexpr ::ndarray<unsigned int, 4, 4>
5420 quad_lines_quad{{{{0, 8, 4, 10}},
5425 static constexpr ::ndarray<unsigned int, 12, 2>
5426 line_vertices_tri{{{{0, 3}},
5439 static constexpr ::ndarray<unsigned int, 4, 4>
5440 quad_lines_tri{{{{0, 8, 5, X}},
5445 static constexpr ::ndarray<unsigned int, 4, 4, 2>
5446 quad_line_vertices_tri{
5447 {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
5448 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
5449 {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
5450 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
5452 const auto &line_vertices =
5454 line_vertices_quad :
5456 const auto &quad_lines =
5461 for (
unsigned int i = 0, j = 2 * quad->
n_lines();
5465 auto &new_line = new_lines[i];
5466 new_line->set_bounding_object_indices(
5469 new_line->set_used_flag();
5470 new_line->clear_user_flag();
5471 new_line->clear_user_data();
5472 new_line->clear_children();
5473 new_line->set_boundary_id_internal(quad->boundary_id());
5474 new_line->set_manifold_id(quad->manifold_id());
5478 for (
unsigned int i = 0; i < new_quads.size(); ++i)
5480 auto &new_quad = new_quads[i];
5485 reference_face_type);
5488 new_quad->set_bounding_object_indices(
5489 {line_indices[quad_lines[i][0]],
5490 line_indices[quad_lines[i][1]],
5491 line_indices[quad_lines[i][2]]});
5493 new_quad->set_bounding_object_indices(
5494 {line_indices[quad_lines[i][0]],
5495 line_indices[quad_lines[i][1]],
5496 line_indices[quad_lines[i][2]],
5497 line_indices[quad_lines[i][3]]});
5501 new_quad->set_used_flag();
5502 new_quad->clear_user_flag();
5503 new_quad->clear_user_data();
5504 new_quad->clear_children();
5505 new_quad->set_boundary_id_internal(quad->boundary_id());
5506 new_quad->set_manifold_id(quad->manifold_id());
5509 std::set<unsigned int> s;
5517 for (
const auto f : new_quad->line_indices())
5519 const std::array<unsigned int, 2> vertices_0 = {
5520 {lines[quad_lines[i][f]]->vertex_index(0),
5521 lines[quad_lines[i][f]]->vertex_index(1)}};
5523 const std::array<unsigned int, 2> vertices_1 = {
5527 const auto orientation =
5533 for (
const auto i : vertices_0)
5535 for (
const auto i : vertices_1)
5539 new_quad->set_line_orientation(f, orientation);
5551 static constexpr ::ndarray<unsigned int, 4, 2>
5552 quad_child_boundary_lines{
5553 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
5555 for (
unsigned int i = 0; i < 4; ++i)
5556 for (
unsigned int j = 0; j < 2; ++j)
5557 new_quads[quad_child_boundary_lines[i][j]]
5558 ->set_line_orientation(i, quad->line_orientation(i));
5561 quad->clear_user_flag();
5566 cells_with_distorted_children;
5576 hex->level() >=
static_cast<int>(
level),
5580 hex->level() ==
static_cast<int>(
level);
5583 if (hex->refine_flag_set() ==
5587 const auto &reference_cell_type = hex->reference_cell();
5590 hex->clear_refine_flag();
5591 hex->set_refinement_case(ref_case);
5593 unsigned int n_new_lines = 0;
5594 unsigned int n_new_quads = 0;
5595 unsigned int n_new_hexes = 0;
5612 std::array<raw_line_iterator, 6> new_lines;
5613 for (
unsigned int i = 0; i < n_new_lines; ++i)
5620 new_lines[i]->set_used_flag();
5621 new_lines[i]->clear_user_flag();
5622 new_lines[i]->clear_user_data();
5623 new_lines[i]->clear_children();
5624 new_lines[i]->set_boundary_id_internal(
5626 new_lines[i]->set_manifold_id(hex->manifold_id());
5629 std::array<raw_quad_iterator, 12> new_quads;
5630 for (
unsigned int i = 0; i < n_new_quads; ++i)
5636 auto &new_quad = new_quads[i];
5642 reference_cell_type.face_reference_cell(0));
5645 new_quad->set_used_flag();
5646 new_quad->clear_user_flag();
5647 new_quad->clear_user_data();
5648 new_quad->clear_children();
5649 new_quad->set_boundary_id_internal(
5651 new_quad->set_manifold_id(hex->manifold_id());
5652 for (
const auto j : new_quads[i]->line_indices())
5653 new_quad->set_line_orientation(j,
true);
5662 for (
unsigned int i = 0; i < n_new_hexes; ++i)
5671 new_hexes[i] = next_unused_hex;
5673 auto &new_hex = new_hexes[i];
5677 ->reference_cell[new_hex->index()] =
5678 reference_cell_type;
5681 new_hex->set_used_flag();
5682 new_hex->clear_user_flag();
5683 new_hex->clear_user_data();
5684 new_hex->clear_children();
5685 new_hex->set_material_id(hex->material_id());
5686 new_hex->set_manifold_id(hex->manifold_id());
5687 new_hex->set_subdomain_id(hex->subdomain_id());
5690 new_hex->set_parent(hex->index());
5695 for (
const auto f : new_hex->face_indices())
5696 new_hex->set_combined_face_orientation(
5700 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
5701 hex->set_children(2 * i, new_hexes[2 * i]->index());
5713 const unsigned int n_vertices =
5715 for (
unsigned int i = 0; i < n_vertices; ++i)
5718 const std::array<unsigned int, 12> line_indices =
5719 TriaAccessorImplementation::Implementation::
5720 get_line_indices_of_cell(*hex);
5724 for (
unsigned int l = 0; l < n_lines; ++l)
5734 for (
const unsigned int i : hex->face_indices())
5736 hex->face(i)->child(0)->vertex_index(3);
5740 get_next_unused_vertex(current_vertex,
5745 hex->center(
true,
true);
5751 static constexpr ::ndarray<unsigned int, 6, 2>
5752 new_line_vertices_hex = {{{{22, 26}},
5759 static constexpr ::ndarray<unsigned int, 6, 2>
5760 new_line_vertices_tet = {{{{6, 8}},
5767 const auto &new_line_vertices =
5769 new_line_vertices_hex :
5770 new_line_vertices_tet;
5772 for (
unsigned int i = 0; i < n_new_lines; ++i)
5773 new_lines[i]->set_bounding_object_indices(
5780 boost::container::small_vector<raw_line_iterator, 30>
5785 relevant_lines.resize(30);
5786 for (
unsigned int f = 0, k = 0; f < 6; ++f)
5787 for (
unsigned int c = 0; c < 4; ++c, ++k)
5790 ndarray<unsigned int, 4, 2>
5792 {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}};
5798 standard_to_real_face_vertex(
5800 hex->face_orientation(f),
5802 hex->face_rotation(f)))
5804 standard_to_real_face_line(
5806 hex->face_orientation(f),
5808 hex->face_rotation(f)));
5811 for (
unsigned int i = 0, k = 24; i < 6; ++i, ++k)
5812 relevant_lines[k] = new_lines[i];
5816 relevant_lines.resize(13);
5819 for (
unsigned int f = 0; f < 4; ++f)
5820 for (
unsigned int l = 0; l < 3; ++l, ++k)
5824 array<std::array<unsigned int, 3>, 6>
5825 table = {{{{1, 0, 2}},
5838 .get_combined_orientation(
5844 relevant_lines[k++] = new_lines[0];
5851 boost::container::small_vector<unsigned int, 30>
5852 relevant_line_indices(relevant_lines.size());
5853 for (
unsigned int i = 0; i < relevant_line_indices.size();
5855 relevant_line_indices[i] = relevant_lines[i]->
index();
5857 static constexpr ::ndarray<unsigned int, 12, 4>
5858 new_quad_lines_hex = {{{{10, 28, 16, 24}},
5869 {{25, 7, 27, 13}}}};
5871 static constexpr ::ndarray<unsigned int, 12, 4>
5872 new_quad_lines_tet = {{{{2, 3, 8, X}},
5885 static constexpr ::ndarray<unsigned int, 12, 4, 2>
5887 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
5888 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
5889 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
5890 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
5891 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
5892 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
5893 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
5894 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
5895 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
5896 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
5897 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
5898 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
5900 static constexpr ::ndarray<unsigned int, 12, 4, 2>
5902 {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
5903 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
5904 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
5905 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
5906 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
5907 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
5908 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
5909 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
5910 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5911 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5912 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5913 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
5915 const auto &new_quad_lines =
5917 new_quad_lines_hex :
5925 static constexpr ::ndarray<unsigned int, 4, 2>
5926 representative_lines{
5927 {{{0, 2}}, {{2, 0}}, {{3, 3}}, {{1, 1}}}};
5929 for (
unsigned int q = 0; q < n_new_quads; ++q)
5931 auto &new_quad = new_quads[q];
5933 if (new_quad->n_lines() == 3)
5934 new_quad->set_bounding_object_indices(
5935 {relevant_line_indices[new_quad_lines[q][0]],
5936 relevant_line_indices[new_quad_lines[q][1]],
5937 relevant_line_indices[new_quad_lines[q][2]]});
5938 else if (new_quad->n_lines() == 4)
5939 new_quad->set_bounding_object_indices(
5940 {relevant_line_indices[new_quad_lines[q][0]],
5941 relevant_line_indices[new_quad_lines[q][1]],
5942 relevant_line_indices[new_quad_lines[q][2]],
5943 relevant_line_indices[new_quad_lines[q][3]]});
5951 const unsigned int n_compute_lines =
5955 for (
unsigned int line = 0; line < n_compute_lines;
5958 const unsigned int l =
5959 (reference_cell_type ==
5961 representative_lines[q % 4][0] :
5964 const std::array<unsigned int, 2> vertices_0 = {
5965 {relevant_lines[new_quad_lines[q][l]]
5967 relevant_lines[new_quad_lines[q][l]]
5968 ->vertex_index(1)}};
5970 const std::array<unsigned int, 2> vertices_1 = {
5974 const auto orientation =
5979 new_quad->set_line_orientation(l, orientation);
5984 if (reference_cell_type ==
5986 new_quads[representative_lines[q % 4][1] + q -
5988 ->set_line_orientation(l, orientation);
5995 std::array<int, 36> quad_indices;
5999 for (
unsigned int i = 0; i < n_new_quads; ++i)
6000 quad_indices[i] = new_quads[i]->
index();
6002 for (
unsigned int f = 0, k = n_new_quads; f < 6; ++f)
6003 for (
unsigned int c = 0; c < 4; ++c, ++k)
6005 hex->face(f)->isotropic_child_index(
6008 hex->face_orientation(f),
6010 hex->face_rotation(f)));
6014 for (
unsigned int i = 0; i < n_new_quads; ++i)
6015 quad_indices[i] = new_quads[i]->
index();
6017 for (
unsigned int f = 0, k = n_new_quads; f < 4; ++f)
6018 for (
unsigned int c = 0; c < 4; ++c, ++k)
6020 quad_indices[k] = hex->face(f)->child_index(
6024 .standard_to_real_face_vertex(
6029 .get_combined_orientation(
6040 static constexpr ::ndarray<unsigned int, 8, 6>
6042 {{12, 0, 20, 4, 28, 8}},
6043 {{0, 16, 22, 6, 29, 9}},
6044 {{13, 1, 4, 24, 30, 10}},
6045 {{1, 17, 6, 26, 31, 11}},
6046 {{14, 2, 21, 5, 8, 32}},
6047 {{2, 18, 23, 7, 9, 33}},
6048 {{15, 3, 5, 25, 10, 34}},
6049 {{3, 19, 7, 27, 11, 35}}
6052 static constexpr ::ndarray<unsigned int, 8, 6>
6053 cell_quads_tet{{{{8, 13, 16, 0, X, X}},
6054 {{9, 12, 1, 21, X, X}},
6055 {{10, 2, 17, 20, X, X}},
6056 {{3, 14, 18, 22, X, X}},
6057 {{11, 1, 4, 5, X, X}},
6058 {{15, 0, 4, 6, X, X}},
6059 {{19, 7, 6, 3, X, X}},
6060 {{23, 5, 2, 7, X, X}}}};
6062 static constexpr ::ndarray<unsigned int, 8, 6, 4>
6063 cell_face_vertices_tet{{{{{{0, 4, 6, X}},
6112 const auto &cell_quads =
6117 for (
unsigned int c = 0;
6118 c < GeometryInfo<dim>::max_children_per_cell;
6121 auto &new_hex = new_hexes[c];
6123 if (new_hex->n_faces() == 4)
6125 new_hex->set_bounding_object_indices(
6126 {quad_indices[cell_quads[c][0]],
6127 quad_indices[cell_quads[c][1]],
6128 quad_indices[cell_quads[c][2]],
6129 quad_indices[cell_quads[c][3]]});
6133 for (
const auto f : new_hex->face_indices())
6135 const auto &face = new_hex->face(f);
6137 Assert(face->n_vertices() == 3,
6140 const std::array<unsigned int, 3> vertices_0 = {
6141 {face->vertex_index(0),
6142 face->vertex_index(1),
6143 face->vertex_index(2)}};
6145 const std::array<unsigned int, 3> vertices_1 = {
6155 new_hex->set_combined_face_orientation(
6157 face->reference_cell()
6158 .get_combined_orientation(
6163 else if (new_hex->n_faces() == 6)
6164 new_hex->set_bounding_object_indices(
6165 {quad_indices[cell_quads[c][0]],
6166 quad_indices[cell_quads[c][1]],
6167 quad_indices[cell_quads[c][2]],
6168 quad_indices[cell_quads[c][3]],
6169 quad_indices[cell_quads[c][4]],
6170 quad_indices[cell_quads[c][5]]});
6179 static constexpr ::ndarray<unsigned int, 6, 4>
6180 face_to_child_indices_hex{{{{0, 2, 4, 6}},
6187 for (
const auto f : hex->face_indices())
6189 const unsigned char combined_orientation =
6190 hex->combined_face_orientation(f);
6191 for (
unsigned int c = 0; c < 4; ++c)
6192 new_hexes[face_to_child_indices_hex[f][c]]
6193 ->set_combined_face_orientation(
6194 f, combined_orientation);
6199 if (check_for_distorted_cells &&
6200 has_distorted_children<dim, spacedim>(hex))
6209 return cells_with_distorted_children;
6216 template <
int spacedim>
6219 const bool check_for_distorted_cells)
6221 const unsigned int dim = 3;
6224 bool flag_isotropic_mesh =
true;
6228 for (; cell != endc; ++cell)
6238 flag_isotropic_mesh =
false;
6242 if (flag_isotropic_mesh)
6244 check_for_distorted_cells);
6255 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
6257 if (cell->refine_flag_set())
6274 line->clear_user_flag();
6279 quad->clear_user_flag();
6302 unsigned int needed_vertices = 0;
6303 unsigned int needed_lines_single = 0;
6304 unsigned int needed_quads_single = 0;
6305 unsigned int needed_lines_pair = 0;
6306 unsigned int needed_quads_pair = 0;
6311 unsigned int new_cells = 0;
6313 for (
const auto &acell :
6315 if (acell->refine_flag_set())
6325 ++needed_quads_single;
6333 ++needed_lines_single;
6334 needed_quads_single += 4;
6341 needed_lines_single += 6;
6342 needed_quads_single += 12;
6356 for (
const unsigned int face :
6360 aface = acell->face(face);
6367 acell->face_orientation(face),
6368 acell->face_flip(face),
6369 acell->face_rotation(face));
6374 if (face_ref_case ==
6377 if (aface->n_active_descendants() < 4)
6380 aface->set_user_flag();
6382 else if (aface->refinement_case() != face_ref_case)
6393 Assert(aface->refinement_case() ==
6395 dim - 1>::isotropic_refinement ||
6396 aface->refinement_case() ==
6399 aface->set_user_index(face_ref_case);
6405 for (
unsigned int line = 0;
6406 line < GeometryInfo<dim>::lines_per_cell;
6410 !acell->line(line)->has_children())
6411 acell->line(line)->set_user_flag();
6417 const unsigned int used_cells =
6427 used_cells + new_cells,
6441 if (quad->user_flag_set())
6447 needed_quads_pair += 4;
6448 needed_lines_pair += 4;
6449 needed_vertices += 1;
6451 if (quad->user_index())
6455 needed_quads_pair += 2;
6456 needed_lines_single += 1;
6469 if (quad->has_children())
6471 Assert(quad->refinement_case() ==
6474 if ((face_refinement_cases[quad->user_index()] ==
6476 (quad->child(0)->line_index(1) + 1 !=
6477 quad->child(2)->line_index(1))) ||
6478 (face_refinement_cases[quad->user_index()] ==
6480 (quad->child(0)->line_index(3) + 1 !=
6481 quad->child(1)->line_index(3))))
6482 needed_lines_pair += 2;
6491 if (line->user_flag_set())
6493 needed_lines_pair += 2;
6494 needed_vertices += 1;
6500 needed_lines_single);
6504 needed_quads_single);
6507 needed_quads_single);
6511 needed_vertices += std::count(
triangulation.vertices_used.begin(),
6535 for (
const auto &cell :
triangulation.active_cell_iterators())
6536 if (!cell->refine_flag_set())
6537 for (
unsigned int line = 0;
6538 line < GeometryInfo<dim>::lines_per_cell;
6540 if (cell->line(line)->has_children())
6541 for (
unsigned int c = 0; c < 2; ++c)
6542 Assert(cell->line(line)->child(c)->user_flag_set() ==
false,
6554 unsigned int next_unused_vertex = 0;
6565 for (; line != endl; ++line)
6566 if (line->user_flag_set())
6572 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
6573 ++next_unused_vertex;
6577 "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."));
6580 triangulation.vertices[next_unused_vertex] = line->center(
true);
6586 triangulation.faces->lines.template next_free_pair_object<1>(
6594 line->set_children(0, next_unused_line->index());
6598 children[2] = {next_unused_line, ++next_unused_line};
6605 children[0]->set_bounding_object_indices(
6606 {line->vertex_index(0), next_unused_vertex});
6607 children[1]->set_bounding_object_indices(
6608 {next_unused_vertex, line->vertex_index(1)});
6610 children[0]->set_used_flag();
6611 children[1]->set_used_flag();
6612 children[0]->clear_children();
6613 children[1]->clear_children();
6616 children[0]->clear_user_flag();
6617 children[1]->clear_user_flag();
6619 children[0]->set_boundary_id_internal(line->boundary_id());
6620 children[1]->set_boundary_id_internal(line->boundary_id());
6622 children[0]->set_manifold_id(line->manifold_id());
6623 children[1]->set_manifold_id(line->manifold_id());
6628 line->clear_user_flag();
6658 for (
unsigned int loop = 0; loop < 2; ++loop)
6671 for (; quad != endq; ++quad)
6673 if (quad->user_index())
6676 face_refinement_cases[quad->user_index()];
6685 if (aniso_quad_ref_case == quad->refinement_case())
6688 Assert(quad->refinement_case() ==
6690 quad->refinement_case() ==
6695 Assert(quad->user_index() ==
6697 quad->user_index() ==
6727 quad->line(2)->child(0)->vertex_index(1);
6729 quad->line(3)->child(0)->vertex_index(1);
6734 quad->line(0)->child(0)->vertex_index(1);
6736 quad->line(1)->child(0)->vertex_index(1);
6739 new_line->set_bounding_object_indices(
6741 new_line->set_used_flag();
6742 new_line->clear_user_flag();
6744 new_line->clear_children();
6745 new_line->set_boundary_id_internal(quad->boundary_id());
6746 new_line->set_manifold_id(quad->manifold_id());
6754 const unsigned int index[2][2] = {
6766 new_quads[0] = next_unused_quad;
6770 new_quads[1] = next_unused_quad;
6775 new_quads[0]->set_bounding_object_indices(
6776 {
static_cast<int>(quad->line_index(0)),
6779 ->child(
index[0][quad->line_orientation(2)])
6782 ->child(
index[0][quad->line_orientation(3)])
6784 new_quads[1]->set_bounding_object_indices(
6786 static_cast<int>(quad->line_index(1)),
6788 ->child(
index[1][quad->line_orientation(2)])
6791 ->child(
index[1][quad->line_orientation(3)])
6796 new_quads[0]->set_bounding_object_indices(
6798 ->child(
index[0][quad->line_orientation(0)])
6801 ->child(
index[0][quad->line_orientation(1)])
6803 static_cast<int>(quad->line_index(2)),
6804 new_line->index()});
6805 new_quads[1]->set_bounding_object_indices(
6807 ->child(
index[1][quad->line_orientation(0)])
6810 ->child(
index[1][quad->line_orientation(1)])
6813 static_cast<int>(quad->line_index(3))});
6816 for (
const auto &new_quad : new_quads)
6818 new_quad->set_used_flag();
6819 new_quad->clear_user_flag();
6821 new_quad->clear_children();
6822 new_quad->set_boundary_id_internal(quad->boundary_id());
6823 new_quad->set_manifold_id(quad->manifold_id());
6827 for (
unsigned int j = 0;
6828 j < GeometryInfo<dim>::lines_per_face;
6830 new_quad->set_line_orientation(j,
true);
6836 new_quads[0]->set_line_orientation(
6837 0, quad->line_orientation(0));
6838 new_quads[0]->set_line_orientation(
6839 2, quad->line_orientation(2));
6840 new_quads[1]->set_line_orientation(
6841 1, quad->line_orientation(1));
6842 new_quads[1]->set_line_orientation(
6843 3, quad->line_orientation(3));
6846 new_quads[0]->set_line_orientation(
6847 3, quad->line_orientation(3));
6848 new_quads[1]->set_line_orientation(
6849 2, quad->line_orientation(2));
6853 new_quads[0]->set_line_orientation(
6854 1, quad->line_orientation(1));
6855 new_quads[1]->set_line_orientation(
6856 0, quad->line_orientation(0));
6862 if (quad->refinement_case() ==
6885 if (aniso_quad_ref_case ==
6888 old_child[0] = quad->child(0)->line(1);
6889 old_child[1] = quad->child(2)->line(1);
6893 Assert(aniso_quad_ref_case ==
6897 old_child[0] = quad->child(0)->line(3);
6898 old_child[1] = quad->child(1)->line(3);
6901 if (old_child[0]->
index() + 1 != old_child[1]->
index())
6907 spacedim>::raw_line_iterator
6910 new_child[0] = new_child[1] =
6912 .template next_free_pair_object<1>(
6916 new_child[0]->set_used_flag();
6917 new_child[1]->set_used_flag();
6919 const int old_index_0 = old_child[0]->index(),
6920 old_index_1 = old_child[1]->index(),
6921 new_index_0 = new_child[0]->index(),
6922 new_index_1 = new_child[1]->index();
6926 for (
unsigned int q = 0;
6929 for (
unsigned int l = 0;
6930 l < GeometryInfo<dim>::lines_per_face;
6933 const int this_index =
6935 .get_bounding_object_indices(q)[l];
6936 if (this_index == old_index_0)
6938 .get_bounding_object_indices(q)[l] =
6940 else if (this_index == old_index_1)
6942 .get_bounding_object_indices(q)[l] =
6947 for (
unsigned int i = 0; i < 2; ++i)
6949 Assert(!old_child[i]->has_children(),
6952 new_child[i]->set_bounding_object_indices(
6953 {old_child[i]->vertex_index(0),
6954 old_child[i]->vertex_index(1)});
6955 new_child[i]->set_boundary_id_internal(
6956 old_child[i]->boundary_id());
6957 new_child[i]->set_manifold_id(
6958 old_child[i]->manifold_id());
6959 new_child[i]->set_user_index(
6960 old_child[i]->user_index());
6961 if (old_child[i]->user_flag_set())
6962 new_child[i]->set_user_flag();
6964 new_child[i]->clear_user_flag();
6966 new_child[i]->clear_children();
6968 old_child[i]->clear_user_flag();
6969 old_child[i]->clear_user_index();
6970 old_child[i]->clear_used_flag();
6976 if (aniso_quad_ref_case ==
6979 new_line->set_children(
6980 0, quad->child(0)->line_index(1));
6981 Assert(new_line->child(1) ==
6982 quad->child(2)->line(1),
7019 quad_iterator switch_1 = quad->child(1),
7020 switch_2 = quad->child(2);
7021 const int switch_1_index = switch_1->index();
7022 const int switch_2_index = switch_2->index();
7023 for (
unsigned int l = 0;
7026 for (
unsigned int h = 0;
7030 for (
const unsigned int q :
7033 const int face_index =
7035 ->cells.get_bounding_object_indices(
7037 if (face_index == switch_1_index)
7039 ->cells.get_bounding_object_indices(
7040 h)[q] = switch_2_index;
7041 else if (face_index == switch_2_index)
7043 ->cells.get_bounding_object_indices(
7044 h)[q] = switch_1_index;
7048 const unsigned int switch_1_lines[4] = {
7049 switch_1->line_index(0),
7050 switch_1->line_index(1),
7051 switch_1->line_index(2),
7052 switch_1->line_index(3)};
7053 const bool switch_1_line_orientations[4] = {
7054 switch_1->line_orientation(0),
7055 switch_1->line_orientation(1),
7056 switch_1->line_orientation(2),
7057 switch_1->line_orientation(3)};
7059 switch_1->boundary_id();
7060 const unsigned int switch_1_user_index =
7061 switch_1->user_index();
7062 const bool switch_1_user_flag =
7063 switch_1->user_flag_set();
7065 switch_1_refinement_case =
7066 switch_1->refinement_case();
7067 const int switch_1_first_child_pair =
7068 (switch_1_refinement_case ?
7069 switch_1->child_index(0) :
7071 const int switch_1_second_child_pair =
7072 (switch_1_refinement_case ==
7074 switch_1->child_index(2) :
7077 switch_1->set_bounding_object_indices(
7078 {switch_2->line_index(0),
7079 switch_2->line_index(1),
7080 switch_2->line_index(2),
7081 switch_2->line_index(3)});
7082 switch_1->set_line_orientation(
7083 0, switch_2->line_orientation(0));
7084 switch_1->set_line_orientation(
7085 1, switch_2->line_orientation(1));
7086 switch_1->set_line_orientation(
7087 2, switch_2->line_orientation(2));
7088 switch_1->set_line_orientation(
7089 3, switch_2->line_orientation(3));
7090 switch_1->set_boundary_id_internal(
7091 switch_2->boundary_id());
7092 switch_1->set_manifold_id(switch_2->manifold_id());
7093 switch_1->set_user_index(switch_2->user_index());
7094 if (switch_2->user_flag_set())
7095 switch_1->set_user_flag();
7097 switch_1->clear_user_flag();
7098 switch_1->clear_refinement_case();
7099 switch_1->set_refinement_case(
7100 switch_2->refinement_case());
7101 switch_1->clear_children();
7102 if (switch_2->refinement_case())
7103 switch_1->set_children(0,
7104 switch_2->child_index(0));
7105 if (switch_2->refinement_case() ==
7107 switch_1->set_children(2,
7108 switch_2->child_index(2));
7110 switch_2->set_bounding_object_indices(
7114 switch_1_lines[3]});
7115 switch_2->set_line_orientation(
7116 0, switch_1_line_orientations[0]);
7117 switch_2->set_line_orientation(
7118 1, switch_1_line_orientations[1]);
7119 switch_2->set_line_orientation(
7120 2, switch_1_line_orientations[2]);
7121 switch_2->set_line_orientation(
7122 3, switch_1_line_orientations[3]);
7123 switch_2->set_boundary_id_internal(
7124 switch_1_boundary_id);
7125 switch_2->set_manifold_id(switch_1->manifold_id());
7126 switch_2->set_user_index(switch_1_user_index);
7127 if (switch_1_user_flag)
7128 switch_2->set_user_flag();
7130 switch_2->clear_user_flag();
7131 switch_2->clear_refinement_case();
7132 switch_2->set_refinement_case(
7133 switch_1_refinement_case);
7134 switch_2->clear_children();
7135 switch_2->set_children(0,
7136 switch_1_first_child_pair);
7137 switch_2->set_children(2,
7138 switch_1_second_child_pair);
7140 new_quads[0]->set_refinement_case(
7142 new_quads[0]->set_children(0, quad->child_index(0));
7143 new_quads[1]->set_refinement_case(
7145 new_quads[1]->set_children(0, quad->child_index(2));
7149 new_quads[0]->set_refinement_case(
7151 new_quads[0]->set_children(0, quad->child_index(0));
7152 new_quads[1]->set_refinement_case(
7154 new_quads[1]->set_children(0, quad->child_index(2));
7155 new_line->set_children(
7156 0, quad->child(0)->line_index(3));
7157 Assert(new_line->child(1) ==
7158 quad->child(1)->line(3),
7161 quad->clear_children();
7165 quad->set_children(0, new_quads[0]->
index());
7167 quad->set_refinement_case(aniso_quad_ref_case);
7174 if (quad->user_flag_set())
7186 ++next_unused_vertex;
7190 "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."));
7198 quad->refinement_case();
7204 quad->child(0)->set_user_index(
7206 quad->child(1)->set_user_index(
7212 middle_line = quad->child(0)->line(1);
7214 middle_line = quad->child(0)->line(3);
7221 if (!middle_line->has_children())
7229 middle_line->center(
true);
7237 .template next_free_pair_object<1>(
7242 middle_line->set_children(
7243 0, next_unused_line->index());
7247 raw_line_iterator children[2] = {
7248 next_unused_line, ++next_unused_line};
7256 children[0]->set_bounding_object_indices(
7257 {middle_line->vertex_index(0),
7258 next_unused_vertex});
7259 children[1]->set_bounding_object_indices(
7260 {next_unused_vertex,
7261 middle_line->vertex_index(1)});
7263 children[0]->set_used_flag();
7264 children[1]->set_used_flag();
7265 children[0]->clear_children();
7266 children[1]->clear_children();
7269 children[0]->clear_user_flag();
7270 children[1]->clear_user_flag();
7272 children[0]->set_boundary_id_internal(
7273 middle_line->boundary_id());
7274 children[1]->set_boundary_id_internal(
7275 middle_line->boundary_id());
7277 children[0]->set_manifold_id(
7278 middle_line->manifold_id());
7279 children[1]->set_manifold_id(
7280 middle_line->manifold_id());
7286 quad->clear_user_flag();
7316 quad->center(
true,
true);
7325 for (
unsigned int i = 0; i < 4; ++i)
7338 new_lines[i] = next_unused_line;
7362 quad->line(0)->child(0)->vertex_index(1),
7363 quad->line(1)->child(0)->vertex_index(1),
7364 quad->line(2)->child(0)->vertex_index(1),
7365 quad->line(3)->child(0)->vertex_index(1),
7366 next_unused_vertex};
7368 new_lines[0]->set_bounding_object_indices(
7370 new_lines[1]->set_bounding_object_indices(
7372 new_lines[2]->set_bounding_object_indices(
7374 new_lines[3]->set_bounding_object_indices(
7377 for (
const auto &new_line : new_lines)
7379 new_line->set_used_flag();
7380 new_line->clear_user_flag();
7382 new_line->clear_children();
7383 new_line->set_boundary_id_internal(quad->boundary_id());
7384 new_line->set_manifold_id(quad->manifold_id());
7403 const unsigned int index[2][2] = {
7407 const int line_indices[12] = {
7409 ->child(
index[0][quad->line_orientation(0)])
7412 ->child(
index[1][quad->line_orientation(0)])
7415 ->child(
index[0][quad->line_orientation(1)])
7418 ->child(
index[1][quad->line_orientation(1)])
7421 ->child(
index[0][quad->line_orientation(2)])
7424 ->child(
index[1][quad->line_orientation(2)])
7427 ->child(
index[0][quad->line_orientation(3)])
7430 ->child(
index[1][quad->line_orientation(3)])
7432 new_lines[0]->index(),
7433 new_lines[1]->index(),
7434 new_lines[2]->index(),
7435 new_lines[3]->index()};
7447 new_quads[0] = next_unused_quad;
7451 new_quads[1] = next_unused_quad;
7457 new_quads[2] = next_unused_quad;
7461 new_quads[3] = next_unused_quad;
7465 quad->set_children(0, new_quads[0]->
index());
7466 quad->set_children(2, new_quads[2]->
index());
7469 new_quads[0]->set_bounding_object_indices(
7474 new_quads[1]->set_bounding_object_indices(
7479 new_quads[2]->set_bounding_object_indices(
7484 new_quads[3]->set_bounding_object_indices(
7489 for (
const auto &new_quad : new_quads)
7491 new_quad->set_used_flag();
7492 new_quad->clear_user_flag();
7494 new_quad->clear_children();
7495 new_quad->set_boundary_id_internal(quad->boundary_id());
7496 new_quad->set_manifold_id(quad->manifold_id());
7500 for (
unsigned int j = 0;
7501 j < GeometryInfo<dim>::lines_per_face;
7503 new_quad->set_line_orientation(j,
true);
7509 new_quads[0]->set_line_orientation(
7510 0, quad->line_orientation(0));
7511 new_quads[0]->set_line_orientation(
7512 2, quad->line_orientation(2));
7513 new_quads[1]->set_line_orientation(
7514 1, quad->line_orientation(1));
7515 new_quads[1]->set_line_orientation(
7516 2, quad->line_orientation(2));
7517 new_quads[2]->set_line_orientation(
7518 0, quad->line_orientation(0));
7519 new_quads[2]->set_line_orientation(
7520 3, quad->line_orientation(3));
7521 new_quads[3]->set_line_orientation(
7522 1, quad->line_orientation(1));
7523 new_quads[3]->set_line_orientation(
7524 3, quad->line_orientation(3));
7528 quad->clear_user_flag();
7539 cells_with_distorted_children;
7553 for (; hex != endh; ++hex)
7554 if (hex->refine_flag_set())
7562 hex->clear_refine_flag();
7563 hex->set_refinement_case(ref_case);
7574 unsigned int n_new_lines = 0;
7575 unsigned int n_new_quads = 0;
7576 unsigned int n_new_hexes = 0;
7607 new_lines(n_new_lines);
7608 for (
unsigned int i = 0; i < n_new_lines; ++i)
7615 new_lines[i]->set_used_flag();
7616 new_lines[i]->clear_user_flag();
7617 new_lines[i]->clear_user_data();
7618 new_lines[i]->clear_children();
7620 new_lines[i]->set_boundary_id_internal(
7624 new_lines[i]->set_manifold_id(hex->manifold_id());
7631 new_quads(n_new_quads);
7632 for (
unsigned int i = 0; i < n_new_quads; ++i)
7639 new_quads[i]->set_used_flag();
7640 new_quads[i]->clear_user_flag();
7641 new_quads[i]->clear_user_data();
7642 new_quads[i]->clear_children();
7644 new_quads[i]->set_boundary_id_internal(
7648 new_quads[i]->set_manifold_id(hex->manifold_id());
7651 for (
unsigned int j = 0;
7652 j < GeometryInfo<dim>::lines_per_face;
7654 new_quads[i]->set_line_orientation(j,
true);
7663 new_hexes(n_new_hexes);
7664 for (
unsigned int i = 0; i < n_new_hexes; ++i)
7673 new_hexes[i] = next_unused_hex;
7676 new_hexes[i]->set_used_flag();
7677 new_hexes[i]->clear_user_flag();
7678 new_hexes[i]->clear_user_data();
7679 new_hexes[i]->clear_children();
7682 new_hexes[i]->set_material_id(hex->material_id());
7683 new_hexes[i]->set_manifold_id(hex->manifold_id());
7684 new_hexes[i]->set_subdomain_id(subdomainid);
7687 new_hexes[i]->set_parent(hex->index());
7699 for (
const unsigned int f :
7701 new_hexes[i]->set_combined_face_orientation(
7706 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
7707 hex->set_children(2 * i, new_hexes[2 * i]->index());
7714 const bool f_or[6] = {hex->face_orientation(0),
7715 hex->face_orientation(1),
7716 hex->face_orientation(2),
7717 hex->face_orientation(3),
7718 hex->face_orientation(4),
7719 hex->face_orientation(5)};
7722 const bool f_fl[6] = {hex->face_flip(0),
7730 const bool f_ro[6] = {hex->face_rotation(0),
7731 hex->face_rotation(1),
7732 hex->face_rotation(2),
7733 hex->face_rotation(3),
7734 hex->face_rotation(4),
7735 hex->face_rotation(5)};
7738 const unsigned char f_co[6] = {
7739 hex->combined_face_orientation(0),
7740 hex->combined_face_orientation(1),
7741 hex->combined_face_orientation(2),
7742 hex->combined_face_orientation(3),
7743 hex->combined_face_orientation(4),
7744 hex->combined_face_orientation(5)};
7754 const unsigned int child_at_origin[2][2][2] = {
7841 raw_line_iterator lines[4] = {
7842 hex->face(2)->child(0)->line(
7843 (hex->face(2)->refinement_case() ==
7847 hex->face(3)->child(0)->line(
7848 (hex->face(3)->refinement_case() ==
7852 hex->face(4)->child(0)->line(
7853 (hex->face(4)->refinement_case() ==
7857 hex->face(5)->child(0)->line(
7858 (hex->face(5)->refinement_case() ==
7864 unsigned int line_indices[4];
7865 for (
unsigned int i = 0; i < 4; ++i)
7866 line_indices[i] = lines[i]->
index();
7873 bool line_orientation[4];
7879 const unsigned int middle_vertices[2] = {
7880 hex->line(2)->child(0)->vertex_index(1),
7881 hex->line(7)->child(0)->vertex_index(1)};
7883 for (
unsigned int i = 0; i < 4; ++i)
7884 if (lines[i]->vertex_index(i % 2) ==
7885 middle_vertices[i % 2])
7886 line_orientation[i] =
true;
7891 Assert(lines[i]->vertex_index((i + 1) % 2) ==
7892 middle_vertices[i % 2],
7894 line_orientation[i] =
false;
7899 new_quads[0]->set_bounding_object_indices(
7905 new_quads[0]->set_line_orientation(
7906 0, line_orientation[0]);
7907 new_quads[0]->set_line_orientation(
7908 1, line_orientation[1]);
7909 new_quads[0]->set_line_orientation(
7910 2, line_orientation[2]);
7911 new_quads[0]->set_line_orientation(
7912 3, line_orientation[3]);
7945 const int quad_indices[11] = {
7946 new_quads[0]->index(),
7948 hex->face(0)->index(),
7950 hex->face(1)->index(),
7952 hex->face(2)->child_index(
7953 child_at_origin[hex->face(2)->refinement_case() -
7954 1][f_fl[2]][f_ro[2]]),
7955 hex->face(2)->child_index(
7957 child_at_origin[hex->face(2)->refinement_case() -
7958 1][f_fl[2]][f_ro[2]]),
7960 hex->face(3)->child_index(
7961 child_at_origin[hex->face(3)->refinement_case() -
7962 1][f_fl[3]][f_ro[3]]),
7963 hex->face(3)->child_index(
7965 child_at_origin[hex->face(3)->refinement_case() -
7966 1][f_fl[3]][f_ro[3]]),
7968 hex->face(4)->child_index(
7969 child_at_origin[hex->face(4)->refinement_case() -
7970 1][f_fl[4]][f_ro[4]]),
7971 hex->face(4)->child_index(
7973 child_at_origin[hex->face(4)->refinement_case() -
7974 1][f_fl[4]][f_ro[4]]),
7976 hex->face(5)->child_index(
7977 child_at_origin[hex->face(5)->refinement_case() -
7978 1][f_fl[5]][f_ro[5]]),
7979 hex->face(5)->child_index(
7981 child_at_origin[hex->face(5)->refinement_case() -
7982 1][f_fl[5]][f_ro[5]])
7986 new_hexes[0]->set_bounding_object_indices(
7993 new_hexes[1]->set_bounding_object_indices(
8068 raw_line_iterator lines[4] = {
8069 hex->face(0)->child(0)->line(
8070 (hex->face(0)->refinement_case() ==
8074 hex->face(1)->child(0)->line(
8075 (hex->face(1)->refinement_case() ==
8079 hex->face(4)->child(0)->line(
8080 (hex->face(4)->refinement_case() ==
8084 hex->face(5)->child(0)->line(
8085 (hex->face(5)->refinement_case() ==
8091 unsigned int line_indices[4];
8092 for (
unsigned int i = 0; i < 4; ++i)
8093 line_indices[i] = lines[i]->
index();
8100 bool line_orientation[4];
8106 const unsigned int middle_vertices[2] = {
8107 hex->line(0)->child(0)->vertex_index(1),
8108 hex->line(5)->child(0)->vertex_index(1)};
8110 for (
unsigned int i = 0; i < 4; ++i)
8111 if (lines[i]->vertex_index(i % 2) ==
8112 middle_vertices[i % 2])
8113 line_orientation[i] =
true;
8117 Assert(lines[i]->vertex_index((i + 1) % 2) ==
8118 middle_vertices[i % 2],
8120 line_orientation[i] =
false;
8125 new_quads[0]->set_bounding_object_indices(
8131 new_quads[0]->set_line_orientation(
8132 0, line_orientation[2]);
8133 new_quads[0]->set_line_orientation(
8134 1, line_orientation[3]);
8135 new_quads[0]->set_line_orientation(
8136 2, line_orientation[0]);
8137 new_quads[0]->set_line_orientation(
8138 3, line_orientation[1]);
8171 const int quad_indices[11] = {
8172 new_quads[0]->index(),
8174 hex->face(0)->child_index(
8175 child_at_origin[hex->face(0)->refinement_case() -
8176 1][f_fl[0]][f_ro[0]]),
8177 hex->face(0)->child_index(
8179 child_at_origin[hex->face(0)->refinement_case() -
8180 1][f_fl[0]][f_ro[0]]),
8182 hex->face(1)->child_index(
8183 child_at_origin[hex->face(1)->refinement_case() -
8184 1][f_fl[1]][f_ro[1]]),
8185 hex->face(1)->child_index(
8187 child_at_origin[hex->face(1)->refinement_case() -
8188 1][f_fl[1]][f_ro[1]]),
8190 hex->face(2)->index(),
8192 hex->face(3)->index(),
8194 hex->face(4)->child_index(
8195 child_at_origin[hex->face(4)->refinement_case() -
8196 1][f_fl[4]][f_ro[4]]),
8197 hex->face(4)->child_index(
8199 child_at_origin[hex->face(4)->refinement_case() -
8200 1][f_fl[4]][f_ro[4]]),
8202 hex->face(5)->child_index(
8203 child_at_origin[hex->face(5)->refinement_case() -
8204 1][f_fl[5]][f_ro[5]]),
8205 hex->face(5)->child_index(
8207 child_at_origin[hex->face(5)->refinement_case() -
8208 1][f_fl[5]][f_ro[5]])
8212 new_hexes[0]->set_bounding_object_indices(
8219 new_hexes[1]->set_bounding_object_indices(
8296 raw_line_iterator lines[4] = {
8297 hex->face(0)->child(0)->line(
8298 (hex->face(0)->refinement_case() ==
8302 hex->face(1)->child(0)->line(
8303 (hex->face(1)->refinement_case() ==
8307 hex->face(2)->child(0)->line(
8308 (hex->face(2)->refinement_case() ==
8312 hex->face(3)->child(0)->line(
8313 (hex->face(3)->refinement_case() ==
8319 unsigned int line_indices[4];
8320 for (
unsigned int i = 0; i < 4; ++i)
8321 line_indices[i] = lines[i]->
index();
8328 bool line_orientation[4];
8334 const unsigned int middle_vertices[2] = {
8335 middle_vertex_index<dim, spacedim>(hex->line(8)),
8336 middle_vertex_index<dim, spacedim>(hex->line(11))};
8338 for (
unsigned int i = 0; i < 4; ++i)
8339 if (lines[i]->vertex_index(i % 2) ==
8340 middle_vertices[i % 2])
8341 line_orientation[i] =
true;
8345 Assert(lines[i]->vertex_index((i + 1) % 2) ==
8346 middle_vertices[i % 2],
8348 line_orientation[i] =
false;
8353 new_quads[0]->set_bounding_object_indices(
8359 new_quads[0]->set_line_orientation(
8360 0, line_orientation[0]);
8361 new_quads[0]->set_line_orientation(
8362 1, line_orientation[1]);
8363 new_quads[0]->set_line_orientation(
8364 2, line_orientation[2]);
8365 new_quads[0]->set_line_orientation(
8366 3, line_orientation[3]);
8400 const int quad_indices[11] = {
8401 new_quads[0]->index(),
8403 hex->face(0)->child_index(
8404 child_at_origin[hex->face(0)->refinement_case() -
8405 1][f_fl[0]][f_ro[0]]),
8406 hex->face(0)->child_index(
8408 child_at_origin[hex->face(0)->refinement_case() -
8409 1][f_fl[0]][f_ro[0]]),
8411 hex->face(1)->child_index(
8412 child_at_origin[hex->face(1)->refinement_case() -
8413 1][f_fl[1]][f_ro[1]]),
8414 hex->face(1)->child_index(
8416 child_at_origin[hex->face(1)->refinement_case() -
8417 1][f_fl[1]][f_ro[1]]),
8419 hex->face(2)->child_index(
8420 child_at_origin[hex->face(2)->refinement_case() -
8421 1][f_fl[2]][f_ro[2]]),
8422 hex->face(2)->child_index(
8424 child_at_origin[hex->face(2)->refinement_case() -
8425 1][f_fl[2]][f_ro[2]]),
8427 hex->face(3)->child_index(
8428 child_at_origin[hex->face(3)->refinement_case() -
8429 1][f_fl[3]][f_ro[3]]),
8430 hex->face(3)->child_index(
8432 child_at_origin[hex->face(3)->refinement_case() -
8433 1][f_fl[3]][f_ro[3]]),
8435 hex->face(4)->index(),
8437 hex->face(5)->index()
8440 new_hexes[0]->set_bounding_object_indices(
8447 new_hexes[1]->set_bounding_object_indices(
8479 new_lines[0]->set_bounding_object_indices(
8480 {middle_vertex_index<dim, spacedim>(hex->face(4)),
8481 middle_vertex_index<dim, spacedim>(hex->face(5))});
8549 spacedim>::raw_line_iterator lines[13] = {
8550 hex->face(0)->child(0)->line(
8551 (hex->face(0)->refinement_case() ==
8555 hex->face(1)->child(0)->line(
8556 (hex->face(1)->refinement_case() ==
8560 hex->face(2)->child(0)->line(
8561 (hex->face(2)->refinement_case() ==
8565 hex->face(3)->child(0)->line(
8566 (hex->face(3)->refinement_case() ==
8574 0, f_or[4], f_fl[4], f_ro[4]))
8577 1, f_or[4], f_fl[4], f_ro[4])),
8581 3, f_or[4], f_fl[4], f_ro[4]))
8584 0, f_or[4], f_fl[4], f_ro[4])),
8588 0, f_or[4], f_fl[4], f_ro[4]))
8591 3, f_or[4], f_fl[4], f_ro[4])),
8595 3, f_or[4], f_fl[4], f_ro[4]))
8598 2, f_or[4], f_fl[4], f_ro[4])),
8603 0, f_or[5], f_fl[5], f_ro[5]))
8606 1, f_or[5], f_fl[5], f_ro[5])),
8610 3, f_or[5], f_fl[5], f_ro[5]))
8613 0, f_or[5], f_fl[5], f_ro[5])),
8617 0, f_or[5], f_fl[5], f_ro[5]))
8620 3, f_or[5], f_fl[5], f_ro[5])),
8624 3, f_or[5], f_fl[5], f_ro[5]))
8627 2, f_or[5], f_fl[5], f_ro[5])),
8632 unsigned int line_indices[13];
8633 for (
unsigned int i = 0; i < 13; ++i)
8634 line_indices[i] = lines[i]->
index();
8641 bool line_orientation[13];
8645 const unsigned int middle_vertices[4] = {
8646 hex->line(0)->child(0)->vertex_index(1),
8647 hex->line(1)->child(0)->vertex_index(1),
8648 hex->line(2)->child(0)->vertex_index(1),
8649 hex->line(3)->child(0)->vertex_index(1),
8655 for (
unsigned int i = 0; i < 4; ++i)
8656 if (lines[i]->vertex_index(0) == middle_vertices[i])
8657 line_orientation[i] =
true;
8661 Assert(lines[i]->vertex_index(1) ==
8664 line_orientation[i] =
false;
8673 for (
unsigned int i = 4; i < 12; ++i)
8674 if (lines[i]->vertex_index((i + 1) % 2) ==
8675 middle_vertex_index<dim, spacedim>(
8676 hex->face(3 + i / 4)))
8677 line_orientation[i] =
true;
8682 Assert(lines[i]->vertex_index(i % 2) ==
8683 (middle_vertex_index<dim, spacedim>(
8684 hex->face(3 + i / 4))),
8686 line_orientation[i] =
false;
8691 line_orientation[12] =
true;
8716 new_quads[0]->set_bounding_object_indices(
8721 new_quads[1]->set_bounding_object_indices(
8726 new_quads[2]->set_bounding_object_indices(
8731 new_quads[3]->set_bounding_object_indices(
8737 new_quads[0]->set_line_orientation(
8738 0, line_orientation[2]);
8739 new_quads[0]->set_line_orientation(
8740 2, line_orientation[4]);
8741 new_quads[0]->set_line_orientation(
8742 3, line_orientation[8]);
8744 new_quads[1]->set_line_orientation(
8745 1, line_orientation[3]);
8746 new_quads[1]->set_line_orientation(
8747 2, line_orientation[5]);
8748 new_quads[1]->set_line_orientation(
8749 3, line_orientation[9]);
8751 new_quads[2]->set_line_orientation(
8752 0, line_orientation[6]);
8753 new_quads[2]->set_line_orientation(
8754 1, line_orientation[10]);
8755 new_quads[2]->set_line_orientation(
8756 2, line_orientation[0]);
8758 new_quads[3]->set_line_orientation(
8759 0, line_orientation[7]);
8760 new_quads[3]->set_line_orientation(
8761 1, line_orientation[11]);
8762 new_quads[3]->set_line_orientation(
8763 3, line_orientation[1]);
8796 const int quad_indices[20] = {
8797 new_quads[0]->index(),
8798 new_quads[1]->index(),
8799 new_quads[2]->index(),
8800 new_quads[3]->index(),
8802 hex->face(0)->child_index(
8803 child_at_origin[hex->face(0)->refinement_case() -
8804 1][f_fl[0]][f_ro[0]]),
8805 hex->face(0)->child_index(
8807 child_at_origin[hex->face(0)->refinement_case() -
8808 1][f_fl[0]][f_ro[0]]),
8810 hex->face(1)->child_index(
8811 child_at_origin[hex->face(1)->refinement_case() -
8812 1][f_fl[1]][f_ro[1]]),
8813 hex->face(1)->child_index(
8815 child_at_origin[hex->face(1)->refinement_case() -
8816 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)->isotropic_child_index(
8836 0, f_or[4], f_fl[4], f_ro[4])),
8837 hex->face(4)->isotropic_child_index(
8839 1, f_or[4], f_fl[4], f_ro[4])),
8840 hex->face(4)->isotropic_child_index(
8842 2, f_or[4], f_fl[4], f_ro[4])),
8843 hex->face(4)->isotropic_child_index(
8845 3, f_or[4], f_fl[4], f_ro[4])),
8847 hex->face(5)->isotropic_child_index(
8849 0, f_or[5], f_fl[5], f_ro[5])),
8850 hex->face(5)->isotropic_child_index(
8852 1, f_or[5], f_fl[5], f_ro[5])),
8853 hex->face(5)->isotropic_child_index(
8855 2, f_or[5], f_fl[5], f_ro[5])),
8856 hex->face(5)->isotropic_child_index(
8858 3, f_or[5], f_fl[5], f_ro[5]))};
8860 new_hexes[0]->set_bounding_object_indices(
8867 new_hexes[1]->set_bounding_object_indices(
8874 new_hexes[2]->set_bounding_object_indices(
8881 new_hexes[3]->set_bounding_object_indices(
8913 new_lines[0]->set_bounding_object_indices(
8914 {middle_vertex_index<dim, spacedim>(hex->face(2)),
8915 middle_vertex_index<dim, spacedim>(hex->face(3))});
8983 spacedim>::raw_line_iterator lines[13] = {
8984 hex->face(0)->child(0)->line(
8985 (hex->face(0)->refinement_case() ==
8989 hex->face(1)->child(0)->line(
8990 (hex->face(1)->refinement_case() ==
8994 hex->face(4)->child(0)->line(
8995 (hex->face(4)->refinement_case() ==
8999 hex->face(5)->child(0)->line(
9000 (hex->face(5)->refinement_case() ==
9008 0, f_or[2], f_fl[2], f_ro[2]))
9011 3, f_or[2], f_fl[2], f_ro[2])),
9015 3, f_or[2], f_fl[2], f_ro[2]))
9018 2, f_or[2], f_fl[2], f_ro[2])),
9022 0, f_or[2], f_fl[2], f_ro[2]))
9025 1, f_or[2], f_fl[2], f_ro[2])),
9029 3, f_or[2], f_fl[2], f_ro[2]))
9032 0, f_or[2], f_fl[2], f_ro[2])),
9037 0, f_or[3], f_fl[3], f_ro[3]))
9040 3, f_or[3], f_fl[3], f_ro[3])),
9044 3, f_or[3], f_fl[3], f_ro[3]))
9047 2, f_or[3], f_fl[3], f_ro[3])),
9051 0, f_or[3], f_fl[3], f_ro[3]))
9054 1, f_or[3], f_fl[3], f_ro[3])),
9058 3, f_or[3], f_fl[3], f_ro[3]))
9061 0, f_or[3], f_fl[3], f_ro[3])),
9066 unsigned int line_indices[13];
9067 for (
unsigned int i = 0; i < 13; ++i)
9068 line_indices[i] = lines[i]->
index();
9075 bool line_orientation[13];
9079 const unsigned int middle_vertices[4] = {
9080 hex->line(8)->child(0)->vertex_index(1),
9081 hex->line(9)->child(0)->vertex_index(1),
9082 hex->line(2)->child(0)->vertex_index(1),
9083 hex->line(6)->child(0)->vertex_index(1),
9088 for (
unsigned int i = 0; i < 4; ++i)
9089 if (lines[i]->vertex_index(0) == middle_vertices[i])
9090 line_orientation[i] =
true;
9094 Assert(lines[i]->vertex_index(1) ==
9097 line_orientation[i] =
false;
9106 for (
unsigned int i = 4; i < 12; ++i)
9107 if (lines[i]->vertex_index((i + 1) % 2) ==
9108 middle_vertex_index<dim, spacedim>(
9109 hex->face(1 + i / 4)))
9110 line_orientation[i] =
true;
9115 Assert(lines[i]->vertex_index(i % 2) ==
9116 (middle_vertex_index<dim, spacedim>(
9117 hex->face(1 + i / 4))),
9119 line_orientation[i] =
false;
9124 line_orientation[12] =
true;
9150 new_quads[0]->set_bounding_object_indices(
9155 new_quads[1]->set_bounding_object_indices(
9160 new_quads[2]->set_bounding_object_indices(
9165 new_quads[3]->set_bounding_object_indices(
9171 new_quads[0]->set_line_orientation(
9172 0, line_orientation[0]);
9173 new_quads[0]->set_line_orientation(
9174 2, line_orientation[6]);
9175 new_quads[0]->set_line_orientation(
9176 3, line_orientation[10]);
9178 new_quads[1]->set_line_orientation(
9179 1, line_orientation[1]);
9180 new_quads[1]->set_line_orientation(
9181 2, line_orientation[7]);
9182 new_quads[1]->set_line_orientation(
9183 3, line_orientation[11]);
9185 new_quads[2]->set_line_orientation(
9186 0, line_orientation[4]);
9187 new_quads[2]->set_line_orientation(
9188 1, line_orientation[8]);
9189 new_quads[2]->set_line_orientation(
9190 2, line_orientation[2]);
9192 new_quads[3]->set_line_orientation(
9193 0, line_orientation[5]);
9194 new_quads[3]->set_line_orientation(
9195 1, line_orientation[9]);
9196 new_quads[3]->set_line_orientation(
9197 3, line_orientation[3]);
9231 const int quad_indices[20] = {
9232 new_quads[0]->index(),
9233 new_quads[1]->index(),
9234 new_quads[2]->index(),
9235 new_quads[3]->index(),
9237 hex->face(0)->child_index(
9238 child_at_origin[hex->face(0)->refinement_case() -
9239 1][f_fl[0]][f_ro[0]]),
9240 hex->face(0)->child_index(
9242 child_at_origin[hex->face(0)->refinement_case() -
9243 1][f_fl[0]][f_ro[0]]),
9245 hex->face(1)->child_index(
9246 child_at_origin[hex->face(1)->refinement_case() -
9247 1][f_fl[1]][f_ro[1]]),
9248 hex->face(1)->child_index(
9250 child_at_origin[hex->face(1)->refinement_case() -
9251 1][f_fl[1]][f_ro[1]]),
9253 hex->face(2)->isotropic_child_index(
9255 0, f_or[2], f_fl[2], f_ro[2])),
9256 hex->face(2)->isotropic_child_index(
9258 1, f_or[2], f_fl[2], f_ro[2])),
9259 hex->face(2)->isotropic_child_index(
9261 2, f_or[2], f_fl[2], f_ro[2])),
9262 hex->face(2)->isotropic_child_index(
9264 3, f_or[2], f_fl[2], f_ro[2])),
9266 hex->face(3)->isotropic_child_index(
9268 0, f_or[3], f_fl[3], f_ro[3])),
9269 hex->face(3)->isotropic_child_index(
9271 1, f_or[3], f_fl[3], f_ro[3])),
9272 hex->face(3)->isotropic_child_index(
9274 2, f_or[3], f_fl[3], f_ro[3])),
9275 hex->face(3)->isotropic_child_index(
9277 3, f_or[3], f_fl[3], f_ro[3])),
9279 hex->face(4)->child_index(
9280 child_at_origin[hex->face(4)->refinement_case() -
9281 1][f_fl[4]][f_ro[4]]),
9282 hex->face(4)->child_index(
9284 child_at_origin[hex->face(4)->refinement_case() -
9285 1][f_fl[4]][f_ro[4]]),
9287 hex->face(5)->child_index(
9288 child_at_origin[hex->face(5)->refinement_case() -
9289 1][f_fl[5]][f_ro[5]]),
9290 hex->face(5)->child_index(
9292 child_at_origin[hex->face(5)->refinement_case() -
9293 1][f_fl[5]][f_ro[5]])};
9304 new_hexes[0]->set_bounding_object_indices(
9311 new_hexes[1]->set_bounding_object_indices(
9318 new_hexes[2]->set_bounding_object_indices(
9325 new_hexes[3]->set_bounding_object_indices(
9358 new_lines[0]->set_bounding_object_indices(
9360 {middle_vertex_index<dim, spacedim>(hex->face(0)),
9361 middle_vertex_index<dim, spacedim>(hex->face(1))});
9430 spacedim>::raw_line_iterator lines[13] = {
9431 hex->face(2)->child(0)->line(
9432 (hex->face(2)->refinement_case() ==
9436 hex->face(3)->child(0)->line(
9437 (hex->face(3)->refinement_case() ==
9441 hex->face(4)->child(0)->line(
9442 (hex->face(4)->refinement_case() ==
9446 hex->face(5)->child(0)->line(
9447 (hex->face(5)->refinement_case() ==
9455 0, f_or[0], f_fl[0], f_ro[0]))
9458 1, f_or[0], f_fl[0], f_ro[0])),
9462 3, f_or[0], f_fl[0], f_ro[0]))
9465 0, f_or[0], f_fl[0], f_ro[0])),
9469 0, f_or[0], f_fl[0], f_ro[0]))
9472 3, f_or[0], f_fl[0], f_ro[0])),
9476 3, f_or[0], f_fl[0], f_ro[0]))
9479 2, f_or[0], f_fl[0], f_ro[0])),
9484 0, f_or[1], f_fl[1], f_ro[1]))
9487 1, f_or[1], f_fl[1], f_ro[1])),
9491 3, f_or[1], f_fl[1], f_ro[1]))
9494 0, f_or[1], f_fl[1], f_ro[1])),
9498 0, f_or[1], f_fl[1], f_ro[1]))
9501 3, f_or[1], f_fl[1], f_ro[1])),
9505 3, f_or[1], f_fl[1], f_ro[1]))
9508 2, f_or[1], f_fl[1], f_ro[1])),
9513 unsigned int line_indices[13];
9515 for (
unsigned int i = 0; i < 13; ++i)
9516 line_indices[i] = lines[i]->
index();
9523 bool line_orientation[13];
9527 const unsigned int middle_vertices[4] = {
9528 hex->line(8)->child(0)->vertex_index(1),
9529 hex->line(10)->child(0)->vertex_index(1),
9530 hex->line(0)->child(0)->vertex_index(1),
9531 hex->line(4)->child(0)->vertex_index(1),
9536 for (
unsigned int i = 0; i < 4; ++i)
9537 if (lines[i]->vertex_index(0) == middle_vertices[i])
9538 line_orientation[i] =
true;
9542 Assert(lines[i]->vertex_index(1) ==
9545 line_orientation[i] =
false;
9554 for (
unsigned int i = 4; i < 12; ++i)
9555 if (lines[i]->vertex_index((i + 1) % 2) ==
9556 middle_vertex_index<dim, spacedim>(
9557 hex->face(i / 4 - 1)))
9558 line_orientation[i] =
true;
9563 Assert(lines[i]->vertex_index(i % 2) ==
9564 (middle_vertex_index<dim, spacedim>(
9565 hex->face(i / 4 - 1))),
9567 line_orientation[i] =
false;
9572 line_orientation[12] =
true;
9592 new_quads[0]->set_bounding_object_indices(
9597 new_quads[1]->set_bounding_object_indices(
9602 new_quads[2]->set_bounding_object_indices(
9607 new_quads[3]->set_bounding_object_indices(
9613 new_quads[0]->set_line_orientation(
9614 0, line_orientation[6]);
9615 new_quads[0]->set_line_orientation(
9616 1, line_orientation[10]);
9617 new_quads[0]->set_line_orientation(
9618 2, line_orientation[0]);
9620 new_quads[1]->set_line_orientation(
9621 0, line_orientation[7]);
9622 new_quads[1]->set_line_orientation(
9623 1, line_orientation[11]);
9624 new_quads[1]->set_line_orientation(
9625 3, line_orientation[1]);
9627 new_quads[2]->set_line_orientation(
9628 0, line_orientation[2]);
9629 new_quads[2]->set_line_orientation(
9630 2, line_orientation[4]);
9631 new_quads[2]->set_line_orientation(
9632 3, line_orientation[8]);
9634 new_quads[3]->set_line_orientation(
9635 1, line_orientation[3]);
9636 new_quads[3]->set_line_orientation(
9637 2, line_orientation[5]);
9638 new_quads[3]->set_line_orientation(
9639 3, line_orientation[9]);
9673 const int quad_indices[20] = {
9674 new_quads[0]->index(),
9675 new_quads[1]->index(),
9676 new_quads[2]->index(),
9677 new_quads[3]->index(),
9679 hex->face(0)->isotropic_child_index(
9681 0, f_or[0], f_fl[0], f_ro[0])),
9682 hex->face(0)->isotropic_child_index(
9684 1, f_or[0], f_fl[0], f_ro[0])),
9685 hex->face(0)->isotropic_child_index(
9687 2, f_or[0], f_fl[0], f_ro[0])),
9688 hex->face(0)->isotropic_child_index(
9690 3, f_or[0], f_fl[0], f_ro[0])),
9692 hex->face(1)->isotropic_child_index(
9694 0, f_or[1], f_fl[1], f_ro[1])),
9695 hex->face(1)->isotropic_child_index(
9697 1, f_or[1], f_fl[1], f_ro[1])),
9698 hex->face(1)->isotropic_child_index(
9700 2, f_or[1], f_fl[1], f_ro[1])),
9701 hex->face(1)->isotropic_child_index(
9703 3, f_or[1], f_fl[1], f_ro[1])),
9705 hex->face(2)->child_index(
9706 child_at_origin[hex->face(2)->refinement_case() -
9707 1][f_fl[2]][f_ro[2]]),
9708 hex->face(2)->child_index(
9710 child_at_origin[hex->face(2)->refinement_case() -
9711 1][f_fl[2]][f_ro[2]]),
9713 hex->face(3)->child_index(
9714 child_at_origin[hex->face(3)->refinement_case() -
9715 1][f_fl[3]][f_ro[3]]),
9716 hex->face(3)->child_index(
9718 child_at_origin[hex->face(3)->refinement_case() -
9719 1][f_fl[3]][f_ro[3]]),
9721 hex->face(4)->child_index(
9722 child_at_origin[hex->face(4)->refinement_case() -
9723 1][f_fl[4]][f_ro[4]]),
9724 hex->face(4)->child_index(
9726 child_at_origin[hex->face(4)->refinement_case() -
9727 1][f_fl[4]][f_ro[4]]),
9729 hex->face(5)->child_index(
9730 child_at_origin[hex->face(5)->refinement_case() -
9731 1][f_fl[5]][f_ro[5]]),
9732 hex->face(5)->child_index(
9734 child_at_origin[hex->face(5)->refinement_case() -
9735 1][f_fl[5]][f_ro[5]])};
9737 new_hexes[0]->set_bounding_object_indices(
9744 new_hexes[1]->set_bounding_object_indices(
9751 new_hexes[2]->set_bounding_object_indices(
9758 new_hexes[3]->set_bounding_object_indices(
9796 ++next_unused_vertex;
9800 "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."));
9811 hex->center(
true,
true);
9833 middle_vertex_index<dim, spacedim>(hex->face(0)),
9834 middle_vertex_index<dim, spacedim>(hex->face(1)),
9835 middle_vertex_index<dim, spacedim>(hex->face(2)),
9836 middle_vertex_index<dim, spacedim>(hex->face(3)),
9837 middle_vertex_index<dim, spacedim>(hex->face(4)),
9838 middle_vertex_index<dim, spacedim>(hex->face(5)),
9839 next_unused_vertex};
9841 new_lines[0]->set_bounding_object_indices(
9843 new_lines[1]->set_bounding_object_indices(
9845 new_lines[2]->set_bounding_object_indices(
9847 new_lines[3]->set_bounding_object_indices(
9849 new_lines[4]->set_bounding_object_indices(
9851 new_lines[5]->set_bounding_object_indices(
9921 spacedim>::raw_line_iterator lines[30] = {
9925 0, f_or[0], f_fl[0], f_ro[0]))
9928 1, f_or[0], f_fl[0], f_ro[0])),
9932 3, f_or[0], f_fl[0], f_ro[0]))
9935 0, f_or[0], f_fl[0], f_ro[0])),
9939 0, f_or[0], f_fl[0], f_ro[0]))
9942 3, f_or[0], f_fl[0], f_ro[0])),
9946 3, f_or[0], f_fl[0], f_ro[0]))
9949 2, f_or[0], f_fl[0], f_ro[0])),
9954 0, f_or[1], f_fl[1], f_ro[1]))
9957 1, f_or[1], f_fl[1], f_ro[1])),
9961 3, f_or[1], f_fl[1], f_ro[1]))
9964 0, f_or[1], f_fl[1], f_ro[1])),
9968 0, f_or[1], f_fl[1], f_ro[1]))
9971 3, f_or[1], f_fl[1], f_ro[1])),
9975 3, f_or[1], f_fl[1], f_ro[1]))
9978 2, f_or[1], f_fl[1], f_ro[1])),
9983 0, f_or[2], f_fl[2], f_ro[2]))
9986 1, f_or[2], f_fl[2], f_ro[2])),
9990 3, f_or[2], f_fl[2], f_ro[2]))
9993 0, f_or[2], f_fl[2], f_ro[2])),
9997 0, f_or[2], f_fl[2], f_ro[2]))
10000 3, f_or[2], f_fl[2], f_ro[2])),
10004 3, f_or[2], f_fl[2], f_ro[2]))
10007 2, f_or[2], f_fl[2], f_ro[2])),
10012 0, f_or[3], f_fl[3], f_ro[3]))
10015 1, f_or[3], f_fl[3], f_ro[3])),
10019 3, f_or[3], f_fl[3], f_ro[3]))
10022 0, f_or[3], f_fl[3], f_ro[3])),
10026 0, f_or[3], f_fl[3], f_ro[3]))
10029 3, f_or[3], f_fl[3], f_ro[3])),
10033 3, f_or[3], f_fl[3], f_ro[3]))
10036 2, f_or[3], f_fl[3], f_ro[3])),
10041 0, f_or[4], f_fl[4], f_ro[4]))
10044 1, f_or[4], f_fl[4], f_ro[4])),
10048 3, f_or[4], f_fl[4], f_ro[4]))
10051 0, f_or[4], f_fl[4], f_ro[4])),
10055 0, f_or[4], f_fl[4], f_ro[4]))
10058 3, f_or[4], f_fl[4], f_ro[4])),
10062 3, f_or[4], f_fl[4], f_ro[4]))
10065 2, f_or[4], f_fl[4], f_ro[4])),
10070 0, f_or[5], f_fl[5], f_ro[5]))
10073 1, f_or[5], f_fl[5], f_ro[5])),
10077 3, f_or[5], f_fl[5], f_ro[5]))
10080 0, f_or[5], f_fl[5], f_ro[5])),
10084 0, f_or[5], f_fl[5], f_ro[5]))
10087 3, f_or[5], f_fl[5], f_ro[5])),
10091 3, f_or[5], f_fl[5], f_ro[5]))
10094 2, f_or[5], f_fl[5], f_ro[5])),
10104 unsigned int line_indices[30];
10105 for (
unsigned int i = 0; i < 30; ++i)
10106 line_indices[i] = lines[i]->
index();
10113 bool line_orientation[30];
10121 for (
unsigned int i = 0; i < 24; ++i)
10122 if (lines[i]->vertex_index((i + 1) % 2) ==
10124 line_orientation[i] =
true;
10129 Assert(lines[i]->vertex_index(i % 2) ==
10132 line_orientation[i] =
false;
10137 for (
unsigned int i = 24; i < 30; ++i)
10138 line_orientation[i] =
true;
10170 new_quads[0]->set_bounding_object_indices(
10174 line_indices[24]});
10175 new_quads[1]->set_bounding_object_indices(
10179 line_indices[25]});
10180 new_quads[2]->set_bounding_object_indices(
10184 line_indices[20]});
10185 new_quads[3]->set_bounding_object_indices(
10189 line_indices[21]});
10190 new_quads[4]->set_bounding_object_indices(
10194 line_indices[28]});
10195 new_quads[5]->set_bounding_object_indices(
10199 line_indices[29]});
10200 new_quads[6]->set_bounding_object_indices(
10205 new_quads[7]->set_bounding_object_indices(
10210 new_quads[8]->set_bounding_object_indices(
10214 line_indices[26]});
10215 new_quads[9]->set_bounding_object_indices(
10219 line_indices[27]});
10220 new_quads[10]->set_bounding_object_indices(
10224 line_indices[12]});
10225 new_quads[11]->set_bounding_object_indices(
10229 line_indices[13]});
10234 new_quads[0]->set_line_orientation(
10235 0, line_orientation[10]);
10236 new_quads[0]->set_line_orientation(
10237 2, line_orientation[16]);
10239 new_quads[1]->set_line_orientation(
10240 1, line_orientation[14]);
10241 new_quads[1]->set_line_orientation(
10242 2, line_orientation[17]);
10244 new_quads[2]->set_line_orientation(
10245 0, line_orientation[11]);
10246 new_quads[2]->set_line_orientation(
10247 3, line_orientation[20]);
10249 new_quads[3]->set_line_orientation(
10250 1, line_orientation[15]);
10251 new_quads[3]->set_line_orientation(
10252 3, line_orientation[21]);
10254 new_quads[4]->set_line_orientation(
10255 0, line_orientation[18]);
10256 new_quads[4]->set_line_orientation(
10257 2, line_orientation[0]);
10259 new_quads[5]->set_line_orientation(
10260 1, line_orientation[22]);
10261 new_quads[5]->set_line_orientation(
10262 2, line_orientation[1]);
10264 new_quads[6]->set_line_orientation(
10265 0, line_orientation[19]);
10266 new_quads[6]->set_line_orientation(
10267 3, line_orientation[4]);
10269 new_quads[7]->set_line_orientation(
10270 1, line_orientation[23]);
10271 new_quads[7]->set_line_orientation(
10272 3, line_orientation[5]);
10274 new_quads[8]->set_line_orientation(
10275 0, line_orientation[2]);
10276 new_quads[8]->set_line_orientation(
10277 2, line_orientation[8]);
10279 new_quads[9]->set_line_orientation(
10280 1, line_orientation[6]);
10281 new_quads[9]->set_line_orientation(
10282 2, line_orientation[9]);
10284 new_quads[10]->set_line_orientation(
10285 0, line_orientation[3]);
10286 new_quads[10]->set_line_orientation(
10287 3, line_orientation[12]);
10289 new_quads[11]->set_line_orientation(
10290 1, line_orientation[7]);
10291 new_quads[11]->set_line_orientation(
10292 3, line_orientation[13]);
10333 const int quad_indices[36] = {
10334 new_quads[0]->index(),
10335 new_quads[1]->index(),
10336 new_quads[2]->index(),
10337 new_quads[3]->index(),
10338 new_quads[4]->index(),
10339 new_quads[5]->index(),
10340 new_quads[6]->index(),
10341 new_quads[7]->index(),
10342 new_quads[8]->index(),
10343 new_quads[9]->index(),
10344 new_quads[10]->index(),
10345 new_quads[11]->index(),
10347 hex->face(0)->isotropic_child_index(
10349 0, f_or[0], f_fl[0], f_ro[0])),
10350 hex->face(0)->isotropic_child_index(
10352 1, f_or[0], f_fl[0], f_ro[0])),
10353 hex->face(0)->isotropic_child_index(
10355 2, f_or[0], f_fl[0], f_ro[0])),
10356 hex->face(0)->isotropic_child_index(
10358 3, f_or[0], f_fl[0], f_ro[0])),
10360 hex->face(1)->isotropic_child_index(
10362 0, f_or[1], f_fl[1], f_ro[1])),
10363 hex->face(1)->isotropic_child_index(
10365 1, f_or[1], f_fl[1], f_ro[1])),
10366 hex->face(1)->isotropic_child_index(
10368 2, f_or[1], f_fl[1], f_ro[1])),
10369 hex->face(1)->isotropic_child_index(
10371 3, f_or[1], f_fl[1], f_ro[1])),
10373 hex->face(2)->isotropic_child_index(
10375 0, f_or[2], f_fl[2], f_ro[2])),
10376 hex->face(2)->isotropic_child_index(
10378 1, f_or[2], f_fl[2], f_ro[2])),
10379 hex->face(2)->isotropic_child_index(
10381 2, f_or[2], f_fl[2], f_ro[2])),
10382 hex->face(2)->isotropic_child_index(
10384 3, f_or[2], f_fl[2], f_ro[2])),
10386 hex->face(3)->isotropic_child_index(
10388 0, f_or[3], f_fl[3], f_ro[3])),
10389 hex->face(3)->isotropic_child_index(
10391 1, f_or[3], f_fl[3], f_ro[3])),
10392 hex->face(3)->isotropic_child_index(
10394 2, f_or[3], f_fl[3], f_ro[3])),
10395 hex->face(3)->isotropic_child_index(
10397 3, f_or[3], f_fl[3], f_ro[3])),
10399 hex->face(4)->isotropic_child_index(
10401 0, f_or[4], f_fl[4], f_ro[4])),
10402 hex->face(4)->isotropic_child_index(
10404 1, f_or[4], f_fl[4], f_ro[4])),
10405 hex->face(4)->isotropic_child_index(
10407 2, f_or[4], f_fl[4], f_ro[4])),
10408 hex->face(4)->isotropic_child_index(
10410 3, f_or[4], f_fl[4], f_ro[4])),
10412 hex->face(5)->isotropic_child_index(
10414 0, f_or[5], f_fl[5], f_ro[5])),
10415 hex->face(5)->isotropic_child_index(
10417 1, f_or[5], f_fl[5], f_ro[5])),
10418 hex->face(5)->isotropic_child_index(
10420 2, f_or[5], f_fl[5], f_ro[5])),
10421 hex->face(5)->isotropic_child_index(
10423 3, f_or[5], f_fl[5], f_ro[5]))};
10426 new_hexes[0]->set_bounding_object_indices(
10433 new_hexes[1]->set_bounding_object_indices(
10440 new_hexes[2]->set_bounding_object_indices(
10446 quad_indices[10]});
10447 new_hexes[3]->set_bounding_object_indices(
10453 quad_indices[11]});
10456 new_hexes[4]->set_bounding_object_indices(
10462 quad_indices[32]});
10463 new_hexes[5]->set_bounding_object_indices(
10469 quad_indices[33]});
10470 new_hexes[6]->set_bounding_object_indices(
10476 quad_indices[34]});
10477 new_hexes[7]->set_bounding_object_indices(
10483 quad_indices[35]});
10517 for (
unsigned int s = 0;
10524 const unsigned int current_child =
10533 ref_case, f, f_or[f], f_fl[f], f_ro[f]));
10534 new_hexes[current_child]->set_combined_face_orientation(
10540 if (check_for_distorted_cells &&
10541 has_distorted_children<dim, spacedim>(hex))
10561 return cells_with_distorted_children;
10577 template <
int spacedim>
10584 template <
int dim,
int spacedim>
10591 if (spacedim > dim)
10595 if (cell->at_boundary() && cell->refine_flag_set() &&
10596 cell->refine_flag_set() !=
10603 for (
const unsigned int face_no :
10605 if (cell->face(face_no)->at_boundary())
10616 spacedim>::face_iterator
10617 face = cell->face(face_no);
10624 cell->reference_cell()
10627 .transform_real_to_unit_cell(cell, new_bound);
10637 const double dist =
10638 std::fabs(new_unit[face_no / 2] - face_no % 2);
10643 const double allowed = 0.25;
10645 if (dist > allowed)
10646 cell->flag_for_face_refinement(face_no);
10665 template <
int dim,
int spacedim>
10670 ExcMessage(
"Wrong function called -- there should "
10671 "be a specialization."));
10675 template <
int spacedim>
10680 const unsigned int dim = 3;
10681 using raw_line_iterator =
10686 bool mesh_changed =
false;
10690 mesh_changed =
false;
10704 if (cell->refine_flag_set())
10706 const std::array<unsigned int, 12> line_indices =
10707 TriaAccessorImplementation::Implementation::
10708 get_line_indices_of_cell(*cell);
10709 for (
unsigned int l = 0; l < cell->n_lines(); ++l)
10711 cell->refine_flag_set(), l) ==
10718 line->set_user_flag();
10721 else if (cell->has_children() &&
10722 !cell->child(0)->coarsen_flag_set())
10724 const std::array<unsigned int, 12> line_indices =
10725 TriaAccessorImplementation::Implementation::
10726 get_line_indices_of_cell(*cell);
10727 for (
unsigned int l = 0; l < cell->n_lines(); ++l)
10729 cell->refinement_case(), l) ==
10736 line->set_user_flag();
10739 else if (cell->has_children() &&
10740 cell->child(0)->coarsen_flag_set())
10741 cell->set_user_flag();
10753 const std::array<unsigned int, 12> line_indices =
10754 TriaAccessorImplementation::Implementation::
10755 get_line_indices_of_cell(*cell);
10756 for (
unsigned int l = 0; l < cell->n_lines(); ++l)
10758 raw_line_iterator line(&
triangulation, 0, line_indices[l]);
10759 if (line->has_children())
10768 bool offending_line_found =
false;
10770 for (
unsigned int c = 0; c < 2; ++c)
10772 Assert(line->child(c)->has_children() ==
false,
10775 if (line->child(c)->user_flag_set() &&
10777 cell->refine_flag_set(), l) ==
10781 cell->clear_coarsen_flag();
10788 allow_anisotropic_smoothing)
10789 cell->flag_for_line_refinement(l);
10791 cell->set_refine_flag();
10793 for (
unsigned int k = 0; k < cell->n_lines();
10796 cell->refine_flag_set(), l) ==
10805 offending_line_found =
true;
10811 for (
unsigned int k = 0; k < cell->n_lines();
10818 if (!line->has_children() &&
10820 line_refinement_case(
10821 cell->refine_flag_set(), k) !=
10823 line->set_user_flag();
10830 if (offending_line_found)
10832 mesh_changed =
true;
10854 if (cell->user_flag_set())
10856 const std::array<unsigned int, 12> line_indices =
10857 TriaAccessorImplementation::Implementation::
10858 get_line_indices_of_cell(*cell);
10859 for (
unsigned int l = 0; l < cell->n_lines(); ++l)
10864 if (line->has_children() &&
10865 (line->child(0)->user_flag_set() ||
10866 line->child(1)->user_flag_set()))
10868 for (
unsigned int c = 0; c < cell->n_children(); ++c)
10869 cell->child(c)->clear_coarsen_flag();
10870 cell->clear_user_flag();
10871 for (
unsigned int k = 0; k < cell->n_lines(); ++k)
10873 cell->refinement_case(), k) ==
10881 mesh_changed =
true;
10887 while (mesh_changed ==
true);
10898 template <
int dim,
int spacedim>
10918 const unsigned int n_subfaces =
10921 if (n_subfaces == 0 || cell->at_boundary(n))
10923 for (
unsigned int c = 0; c < n_subfaces; ++c)
10926 child = cell->child(
10930 child_neighbor = child->neighbor(n);
10931 if (!child->neighbor_is_coarser(n))
10948 if (child_neighbor->has_children() &&
10949 !(child_neighbor->child(0)->is_active() &&
10950 child_neighbor->child(0)->coarsen_flag_set()))
10955 if (child_neighbor->refine_flag_set())
10970 template <
int spacedim>
10975 template <
int dim,
int spacedim>
10978 std::vector<std::pair<unsigned int, unsigned int>>
adjacent_cells(
10980 {numbers::invalid_unsigned_int, numbers::invalid_unsigned_int});
10982 const auto set_entry = [&](
const auto &face_index,
const auto &cell) {
10983 const std::pair<unsigned int, unsigned int> cell_pair = {
10984 cell->level(), cell->index()};
10985 unsigned int index;
10992 index = 2 * face_index + 0;
11001 index = 2 * face_index + 1;
11007 const auto get_entry =
11008 [&](
const auto &face_index,
11012 if (test == std::pair<unsigned int, unsigned int>(cell->level(),
11026 for (
const auto &face : cell->face_iterators())
11028 set_entry(face->index(), cell);
11030 if (cell->is_active() && face->has_children())
11031 for (
unsigned int c = 0; c < face->n_children(); ++c)
11032 set_entry(face->child(c)->index(), cell);
11036 for (
auto f : cell->face_indices())
11037 cell->set_neighbor(f, get_entry(cell->face(f)->index(), cell));
11040 template <
int dim,
int spacedim>
11045 std::vector<unsigned int> & line_cell_count,
11046 std::vector<unsigned int> & quad_cell_count)
11051 (void)line_cell_count;
11052 (void)quad_cell_count;
11055 template <
int dim,
int spacedim>
11058 const bool check_for_distorted_cells)
11060 return Implementation::execute_refinement_isotropic(
11064 template <
int dim,
int spacedim>
11073 template <
int dim,
int spacedim>
11078 Implementation::prepare_refinement_dim_dependent(
triangulation);
11081 template <
int dim,
int spacedim>
11093 template <
int dim,
int spacedim>
11098 return flat_manifold;
11105template <
int dim,
int spacedim>
11111template <
int dim,
int spacedim>
11114 const MeshSmoothing smooth_grid,
11115 const bool check_for_distorted_cells)
11116 : smooth_grid(smooth_grid)
11117 , anisotropic_refinement(false)
11118 , check_for_distorted_cells(check_for_distorted_cells)
11122 vertex_to_boundary_id_map_1d =
11123 std::make_unique<std::map<unsigned int, types::boundary_id>>();
11124 vertex_to_manifold_id_map_1d =
11125 std::make_unique<std::map<unsigned int, types::manifold_id>>();
11129 signals.create.connect(signals.any_change);
11130 signals.post_refinement.connect(signals.any_change);
11131 signals.clear.connect(signals.any_change);
11132 signals.mesh_movement.connect(signals.any_change);
11137template <
int dim,
int spacedim>
11164template <
int dim,
int spacedim>
11195template <
int dim,
int spacedim>
11213 AssertNothrow((dim == 1) || (vertex_to_boundary_id_map_1d ==
nullptr),
11218 AssertNothrow((dim == 1) || (vertex_to_manifold_id_map_1d ==
nullptr),
11224template <
int dim,
int spacedim>
11232 clear_despite_subscriptions();
11233 periodic_face_pairs_level_0.clear();
11234 periodic_face_map.clear();
11235 reference_cells.clear();
11239template <
int dim,
int spacedim>
11243 return MPI_COMM_SELF;
11248template <
int dim,
int spacedim>
11250const std::weak_ptr<const Utilities::MPI::Partitioner>
Triangulation<
11252 spacedim>::global_active_cell_index_partitioner()
const
11254 return number_cache.active_cell_index_partitioner;
11259template <
int dim,
int spacedim>
11261const std::weak_ptr<const Utilities::MPI::Partitioner>
Triangulation<
11263 spacedim>::global_level_cell_index_partitioner(
const unsigned int level)
const
11267 return number_cache.level_cell_index_partitioners[
level];
11272template <
int dim,
int spacedim>
11275 const MeshSmoothing mesh_smoothing)
11278 ExcTriangulationNotEmpty(
vertices.size(), levels.size()));
11279 smooth_grid = mesh_smoothing;
11284template <
int dim,
int spacedim>
11289 return smooth_grid;
11294template <
int dim,
int spacedim>
11302 manifolds[m_number] = manifold_object.
clone();
11307template <
int dim,
int spacedim>
11315 manifolds.erase(m_number);
11319template <
int dim,
int spacedim>
11327template <
int dim,
int spacedim>
11335 "Error: set_all_manifold_ids() can not be called on an empty Triangulation."));
11337 for (
const auto &cell : this->active_cell_iterators())
11338 cell->set_all_manifold_ids(m_number);
11342template <
int dim,
int spacedim>
11350 "Error: set_all_manifold_ids_on_boundary() can not be called on an empty Triangulation."));
11352 for (
const auto &cell : this->active_cell_iterators())
11354 if (cell->face(f)->at_boundary())
11355 cell->face(f)->set_all_manifold_ids(m_number);
11359template <
int dim,
int spacedim>
11368 "Error: set_all_manifold_ids_on_boundary() can not be called on an empty Triangulation."));
11370 bool boundary_found =
false;
11372 for (
const auto &cell : this->active_cell_iterators())
11376 if (cell->face(f)->at_boundary() &&
11379 boundary_found =
true;
11380 cell->face(f)->set_manifold_id(m_number);
11385 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++
e)
11386 if (cell->line(e)->at_boundary() &&
11387 cell->line(e)->boundary_id() == b_id)
11389 boundary_found =
true;
11390 cell->line(e)->set_manifold_id(m_number);
11394 (void)boundary_found;
11395 Assert(boundary_found, ExcBoundaryIdNotFound(b_id));
11400template <
int dim,
int spacedim>
11407 const auto it = manifolds.find(m_number);
11409 if (it != manifolds.end())
11412 return *(it->second);
11417 return internal::TriangulationImplementation::
11418 get_default_flat_manifold<dim, spacedim>();
11423template <
int dim,
int spacedim>
11432 std::vector<types::boundary_id> boundary_ids;
11433 for (std::map<unsigned int, types::boundary_id>::const_iterator p =
11434 vertex_to_boundary_id_map_1d->begin();
11435 p != vertex_to_boundary_id_map_1d->end();
11437 boundary_ids.push_back(p->second);
11439 return boundary_ids;
11443 std::set<types::boundary_id> b_ids;
11444 for (
auto cell : active_cell_iterators())
11445 if (cell->is_locally_owned())
11446 for (const unsigned
int face : cell->face_indices())
11447 if (cell->at_boundary(face))
11449 std::vector<types::boundary_id> boundary_ids(b_ids.begin(), b_ids.end());
11450 return boundary_ids;
11456template <
int dim,
int spacedim>
11461 std::set<types::manifold_id> m_ids;
11462 for (
const auto &cell : active_cell_iterators())
11463 if (cell->is_locally_owned())
11465 m_ids.insert(cell->manifold_id());
11466 for (
const auto &face : cell->face_iterators())
11470 const auto line_indices = internal::TriaAccessorImplementation::
11471 Implementation::get_line_indices_of_cell(*cell);
11472 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
11474 raw_line_iterator line(
this, 0, line_indices[l]);
11475 m_ids.insert(line->manifold_id());
11479 return {m_ids.begin(), m_ids.end()};
11485template <
int dim,
int spacedim>
11490 Assert((
vertices.size() == 0) && (levels.size() == 0) && (faces ==
nullptr),
11491 ExcTriangulationNotEmpty(
vertices.size(), levels.size()));
11493 (dim == 1 || other_tria.
faces !=
nullptr),
11495 "When calling Triangulation::copy_triangulation(), "
11496 "the target triangulation must be empty but the source "
11497 "triangulation (the argument to this function) must contain "
11498 "something. Here, it seems like the source does not "
11499 "contain anything at all."));
11510 faces = std::make_unique<internal::TriangulationImplementation::TriaFaces>(
11511 *other_tria.
faces);
11513 for (
const auto &p : other_tria.manifolds)
11517 levels.reserve(other_tria.
levels.size());
11520 std::make_unique<internal::TriangulationImplementation::TriaLevel>(
11527 vertex_to_boundary_id_map_1d =
11528 std::make_unique<std::map<unsigned int, types::boundary_id>>(
11531 vertex_to_manifold_id_map_1d =
11532 std::make_unique<std::map<unsigned int, types::manifold_id>>(
11537 this->policy = other_tria.
policy->clone();
11551template <
int dim,
int spacedim>
11558 std::vector<CellData<dim>> reordered_cells(cells);
11562 reorder_compatibility(reordered_cells, reordered_subcelldata);
11570template <
int dim,
int spacedim>
11574 this->update_reference_cells();
11576 if (this->all_reference_cells_are_hyper_cube())
11596template <
int dim,
int spacedim>
11603 Assert((
vertices.size() == 0) && (levels.size() == 0) && (faces ==
nullptr),
11604 ExcTriangulationNotEmpty(
vertices.size(), levels.size()));
11619 clear_despite_subscriptions();
11627 reset_cell_vertex_indices_cache();
11629 *
this, levels.size(), number_cache);
11630 reset_active_cell_indices();
11631 reset_global_cell_indices();
11636 if (check_for_distorted_cells)
11638 DistortedCellList distorted_cells = collect_distorted_coarse_cells(*
this);
11642 AssertThrow(distorted_cells.distorted_cells.size() == 0, distorted_cells);
11675 if (dim < spacedim && all_reference_cells_are_hyper_cube())
11683 bool values[][2] = {{
false,
true}, {
true,
false}};
11684 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
11686 correct(i, j) = (
values[i][j]);
11691 bool values[][4] = {{
false,
true,
true,
false},
11692 {
true,
false,
false,
true},
11693 {
true,
false,
false,
true},
11694 {
false,
true,
true,
false}};
11695 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
11697 correct(i, j) = (
values[i][j]);
11705 std::list<active_cell_iterator> this_round, next_round;
11706 active_cell_iterator neighbor;
11708 this_round.push_back(begin_active());
11709 begin_active()->set_direction_flag(
true);
11710 begin_active()->set_user_flag();
11712 while (this_round.size() > 0)
11714 for (
typename std::list<active_cell_iterator>::iterator cell =
11715 this_round.begin();
11716 cell != this_round.end();
11719 for (
const unsigned int i : (*cell)->face_indices())
11721 if (!((*cell)->face(i)->at_boundary()))
11723 neighbor = (*cell)->neighbor(i);
11725 unsigned int cf = (*cell)->face_index(i);
11726 unsigned int j = 0;
11727 while (neighbor->face_index(j) != cf)
11735 if (neighbor->user_flag_set())
11739 Assert(!(correct(i, j) ^
11740 (neighbor->direction_flag() ==
11741 (*cell)->direction_flag())),
11742 ExcNonOrientableTriangulation());
11746 next_round.push_back(neighbor);
11747 neighbor->set_user_flag();
11748 if ((correct(i, j) ^ (neighbor->direction_flag() ==
11749 (*cell)->direction_flag())))
11750 neighbor->set_direction_flag(
11751 !neighbor->direction_flag());
11759 if (next_round.size() == 0)
11760 for (
const auto &cell : this->active_cell_iterators())
11761 if (cell->user_flag_set() == false)
11763 next_round.push_back(cell);
11764 cell->set_direction_flag(
true);
11765 cell->set_user_flag();
11769 this_round = next_round;
11770 next_round.clear();
11772 clear_user_flags();
11781template <
int dim,
int spacedim>
11801 const CellId a_id(a.id);
11804 const auto a_coarse_cell_index =
11805 this->coarse_cell_id_to_coarse_cell_index(a_id.get_coarse_cell_id());
11806 const auto b_coarse_cell_index =
11807 this->coarse_cell_id_to_coarse_cell_index(b_id.get_coarse_cell_id());
11814 if (a_coarse_cell_index != b_coarse_cell_index)
11815 return a_coarse_cell_index < b_coarse_cell_index;
11817 return a_id < b_id;
11823 for (
unsigned int level = 0;
11834 while (cell_info->id != cell->id().template to_binary<dim>())
11837 for (
const auto face : cell->face_indices())
11838 cell->face(face)->set_manifold_id(
11839 cell_info->manifold_line_ids[face]);
11842 for (
const auto face : cell->face_indices())
11843 cell->face(face)->set_manifold_id(
11844 cell_info->manifold_quad_ids[face]);
11846 const auto line_indices = internal::TriaAccessorImplementation::
11847 Implementation::get_line_indices_of_cell(*cell);
11848 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
11850 raw_line_iterator line(
this, 0, line_indices[l]);
11851 line->set_manifold_id(cell_info->manifold_line_ids[l]);
11855 cell->set_manifold_id(cell_info->manifold_id);
11873 !coarse_cell->id().is_parent_of(
CellId(fine_cell_info->id)))
11877 coarse_cell->set_refine_flag();
11882 spacedim>::execute_coarsening_and_refinement();
11887 for (
unsigned int level = 0;
11896 while (cell_info->id != cell->id().template to_binary<dim>())
11900 for (
auto pair : cell_info->boundary_ids)
11901 if (cell->face(pair.
first)->at_boundary())
11902 cell->face(pair.
first)->set_boundary_id(pair.
second);
11911template <
int dim,
int spacedim>
11916 ExcMessage(
"Only works for dim == spacedim-1"));
11917 for (
const auto &cell : this->active_cell_iterators())
11918 cell->set_direction_flag(!cell->direction_flag());
11923template <
int dim,
int spacedim>
11928 ExcMessage(
"Error: An empty Triangulation can not be refined."));
11930 for (
const auto &cell : this->active_cell_iterators())
11932 cell->clear_coarsen_flag();
11933 cell->set_refine_flag();
11939template <
int dim,
int spacedim>
11943 for (
unsigned int i = 0; i < times; ++i)
11945 set_all_refine_flags();
11946 execute_coarsening_and_refinement();
11952template <
int dim,
int spacedim>
11956 for (
unsigned int i = 0; i < times; ++i)
11958 for (
const auto &cell : this->active_cell_iterators())
11960 cell->clear_refine_flag();
11961 cell->set_coarsen_flag();
11963 execute_coarsening_and_refinement();
11972template <
int dim,
int spacedim>
11977 std::vector<bool>::iterator i = v.begin();
11979 for (
const auto &cell : this->active_cell_iterators())
11980 for (unsigned
int j = 0; j < dim; ++j, ++i)
11981 if (cell->refine_flag_set() & (1 << j))
11989template <
int dim,
int spacedim>
11993 std::vector<bool> v;
11994 save_refine_flags(v);
12003template <
int dim,
int spacedim>
12007 std::vector<bool> v;
12009 load_refine_flags(v);
12014template <
int dim,
int spacedim>
12020 std::vector<bool>::const_iterator i = v.begin();
12021 for (
const auto &cell : this->active_cell_iterators())
12023 unsigned int ref_case = 0;
12025 for (
unsigned int j = 0; j < dim; ++j, ++i)
12027 ref_case += 1 << j;
12029 ExcGridReadError());
12033 cell->clear_refine_flag();
12041template <
int dim,
int spacedim>
12044 std::vector<bool> &v)
const
12047 std::vector<bool>::iterator i = v.begin();
12048 for (
const auto &cell : this->active_cell_iterators())
12050 *i = cell->coarsen_flag_set();
12059template <
int dim,
int spacedim>
12063 std::vector<bool> v;
12064 save_coarsen_flags(v);
12073template <
int dim,
int spacedim>
12077 std::vector<bool> v;
12082 load_coarsen_flags(v);
12087template <
int dim,
int spacedim>
12090 const std::vector<bool> &v)
12094 std::vector<bool>::const_iterator i = v.begin();
12095 for (
const auto &cell : this->active_cell_iterators())
12098 cell->set_coarsen_flag();
12100 cell->clear_coarsen_flag();
12108template <
int dim,
int spacedim>
12112 return anisotropic_refinement;
12121 std::vector<std::vector<bool>>
12122 extract_raw_coarsen_flags(
12123 const std::vector<std::unique_ptr<
12126 std::vector<std::vector<bool>> coarsen_flags(levels.size());
12128 coarsen_flags[
level] = levels[
level]->coarsen_flags;
12129 return coarsen_flags;
12132 std::vector<std::vector<std::uint8_t>>
12133 extract_raw_refine_flags(
12134 const std::vector<std::unique_ptr<
12137 std::vector<std::vector<std::uint8_t>> refine_flags(levels.size());
12139 refine_flags[
level] = levels[
level]->refine_flags;
12140 return refine_flags;
12153 clear_user_data(std::vector<std::unique_ptr<
12156 for (
auto &
level : levels)
12157 level->cells.clear_user_data();
12165 if (faces->
dim == 2)
12171 if (faces->
dim == 3)
12180template <
int dim,
int spacedim>
12185 ::clear_user_data(levels);
12187 ::clear_user_data(faces.get());
12195 clear_user_flags_line(
12198 std::unique_ptr<internal::TriangulationImplementation::TriaLevel>>
12204 for (
const auto &
level : levels)
12205 level->cells.clear_user_flags();
12207 else if (dim == 2 || dim == 3)
12219template <
int dim,
int spacedim>
12223 ::clear_user_flags_line(dim, levels, faces.get());
12231 clear_user_flags_quad(
12234 std::unique_ptr<internal::TriangulationImplementation::TriaLevel>>
12244 for (
const auto &
level : levels)
12245 level->cells.clear_user_flags();
12259template <
int dim,
int spacedim>
12263 ::clear_user_flags_quad(dim, levels, faces.get());
12271 clear_user_flags_hex(
12274 std::unique_ptr<internal::TriangulationImplementation::TriaLevel>>
12288 for (
const auto &
level : levels)
12289 level->cells.clear_user_flags();
12299template <
int dim,
int spacedim>
12303 ::clear_user_flags_hex(dim, levels, faces.get());
12308template <
int dim,
int spacedim>
12312 clear_user_flags_line();
12313 clear_user_flags_quad();
12314 clear_user_flags_hex();
12319template <
int dim,
int spacedim>
12323 save_user_flags_line(out);
12326 save_user_flags_quad(out);
12329 save_user_flags_hex(out);
12337template <
int dim,
int spacedim>
12345 std::vector<bool> tmp;
12347 save_user_flags_line(tmp);
12348 v.insert(v.end(), tmp.begin(), tmp.end());
12352 save_user_flags_quad(tmp);
12353 v.insert(v.end(), tmp.begin(), tmp.end());
12358 save_user_flags_hex(tmp);
12359 v.insert(v.end(), tmp.begin(), tmp.end());
12368template <
int dim,
int spacedim>
12372 load_user_flags_line(in);
12375 load_user_flags_quad(in);
12378 load_user_flags_hex(in);
12386template <
int dim,
int spacedim>
12391 std::vector<bool> tmp;
12395 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
12397 load_user_flags_line(tmp);
12402 tmp.insert(tmp.end(),
12403 v.begin() + n_lines(),
12404 v.begin() + n_lines() + n_quads());
12405 load_user_flags_quad(tmp);
12411 tmp.insert(tmp.end(),
12412 v.begin() + n_lines() + n_quads(),
12413 v.begin() + n_lines() + n_quads() + n_hexs());
12414 load_user_flags_hex(tmp);
12423template <
int dim,
int spacedim>
12426 std::vector<bool> &v)
const
12428 v.resize(n_lines(),
false);
12429 std::vector<bool>::iterator i = v.begin();
12430 line_iterator line = begin_line(), endl = end_line();
12431 for (; line != endl; ++line, ++i)
12432 *i = line->user_flag_set();
12439template <
int dim,
int spacedim>
12443 std::vector<bool> v;
12444 save_user_flags_line(v);
12453template <
int dim,
int spacedim>
12457 std::vector<bool> v;
12462 load_user_flags_line(v);
12467template <
int dim,
int spacedim>
12470 const std::vector<bool> &v)
12472 Assert(v.size() == n_lines(), ExcGridReadError());
12474 line_iterator line = begin_line(), endl = end_line();
12475 std::vector<bool>::const_iterator i = v.begin();
12476 for (; line != endl; ++line, ++i)
12478 line->set_user_flag();
12480 line->clear_user_flag();
12488 template <
typename Iterator>
12490 get_user_flag(
const Iterator &i)
12492 return i->user_flag_set();
12497 template <
int structdim,
int dim,
int spacedim>
12507 template <
typename Iterator>
12509 set_user_flag(
const Iterator &i)
12511 i->set_user_flag();
12516 template <
int structdim,
int dim,
int spacedim>
12525 template <
typename Iterator>
12527 clear_user_flag(
const Iterator &i)
12529 i->clear_user_flag();
12534 template <
int structdim,
int dim,
int spacedim>
12544template <
int dim,
int spacedim>
12547 std::vector<bool> &v)
const
12549 v.resize(n_quads(),
false);
12553 std::vector<bool>::iterator i = v.begin();
12554 quad_iterator quad = begin_quad(), endq = end_quad();
12555 for (; quad != endq; ++quad, ++i)
12556 *i = get_user_flag(quad);
12564template <
int dim,
int spacedim>
12568 std::vector<bool> v;
12569 save_user_flags_quad(v);
12578template <
int dim,
int spacedim>
12582 std::vector<bool> v;
12587 load_user_flags_quad(v);
12592template <
int dim,
int spacedim>
12595 const std::vector<bool> &v)
12597 Assert(v.size() == n_quads(), ExcGridReadError());
12601 quad_iterator quad = begin_quad(), endq = end_quad();
12602 std::vector<bool>::const_iterator i = v.begin();
12603 for (; quad != endq; ++quad, ++i)
12605 set_user_flag(quad);
12607 clear_user_flag(quad);
12615template <
int dim,
int spacedim>
12618 std::vector<bool> &v)
const
12620 v.resize(n_hexs(),
false);
12624 std::vector<bool>::iterator i = v.begin();
12625 hex_iterator hex = begin_hex(), endh = end_hex();
12626 for (; hex != endh; ++hex, ++i)
12627 *i = get_user_flag(hex);
12635template <
int dim,
int spacedim>
12639 std::vector<bool> v;
12640 save_user_flags_hex(v);
12649template <
int dim,
int spacedim>
12653 std::vector<bool> v;
12658 load_user_flags_hex(v);
12663template <
int dim,
int spacedim>
12666 const std::vector<bool> &v)
12668 Assert(v.size() == n_hexs(), ExcGridReadError());
12672 hex_iterator hex = begin_hex(), endh = end_hex();
12673 std::vector<bool>::const_iterator i = v.begin();
12674 for (; hex != endh; ++hex, ++i)
12676 set_user_flag(hex);
12678 clear_user_flag(hex);
12686template <
int dim,
int spacedim>
12689 std::vector<unsigned int> &v)
const
12695 std::vector<unsigned int> tmp;
12697 save_user_indices_line(tmp);
12698 v.insert(v.end(), tmp.begin(), tmp.end());
12702 save_user_indices_quad(tmp);
12703 v.insert(v.end(), tmp.begin(), tmp.end());
12708 save_user_indices_hex(tmp);
12709 v.insert(v.end(), tmp.begin(), tmp.end());
12718template <
int dim,
int spacedim>
12721 const std::vector<unsigned int> &v)
12724 std::vector<unsigned int> tmp;
12728 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
12730 load_user_indices_line(tmp);
12735 tmp.insert(tmp.end(),
12736 v.begin() + n_lines(),
12737 v.begin() + n_lines() + n_quads());
12738 load_user_indices_quad(tmp);
12744 tmp.insert(tmp.end(),
12745 v.begin() + n_lines() + n_quads(),
12746 v.begin() + n_lines() + n_quads() + n_hexs());
12747 load_user_indices_hex(tmp);
12758 template <
typename Iterator>
12760 get_user_index(
const Iterator &i)
12762 return i->user_index();
12767 template <
int structdim,
int dim,
int spacedim>
12778 template <
typename Iterator>
12780 set_user_index(
const Iterator &i,
const unsigned int x)
12782 i->set_user_index(x);
12787 template <
int structdim,
int dim,
int spacedim>
12791 const unsigned int)
12798template <
int dim,
int spacedim>
12801 std::vector<unsigned int> &v)
const
12803 v.resize(n_lines(), 0);
12804 std::vector<unsigned int>::iterator i = v.begin();
12805 line_iterator line = begin_line(), endl = end_line();
12806 for (; line != endl; ++line, ++i)
12807 *i = line->user_index();
12812template <
int dim,
int spacedim>
12815 const std::vector<unsigned int> &v)
12817 Assert(v.size() == n_lines(), ExcGridReadError());
12819 line_iterator line = begin_line(), endl = end_line();
12820 std::vector<unsigned int>::const_iterator i = v.begin();
12821 for (; line != endl; ++line, ++i)
12822 line->set_user_index(*i);
12826template <
int dim,
int spacedim>
12829 std::vector<unsigned int> &v)
const
12831 v.resize(n_quads(), 0);
12835 std::vector<unsigned int>::iterator i = v.begin();
12836 quad_iterator quad = begin_quad(), endq = end_quad();
12837 for (; quad != endq; ++quad, ++i)
12838 *i = get_user_index(quad);
12844template <
int dim,
int spacedim>
12847 const std::vector<unsigned int> &v)
12849 Assert(v.size() == n_quads(), ExcGridReadError());
12853 quad_iterator quad = begin_quad(), endq = end_quad();
12854 std::vector<unsigned int>::const_iterator i = v.begin();
12855 for (; quad != endq; ++quad, ++i)
12856 set_user_index(quad, *i);
12861template <
int dim,
int spacedim>
12864 std::vector<unsigned int> &v)
const
12866 v.resize(n_hexs(), 0);
12870 std::vector<unsigned int>::iterator i = v.begin();
12871 hex_iterator hex = begin_hex(), endh = end_hex();
12872 for (; hex != endh; ++hex, ++i)
12873 *i = get_user_index(hex);
12879template <
int dim,
int spacedim>
12882 const std::vector<unsigned int> &v)
12884 Assert(v.size() == n_hexs(), ExcGridReadError());
12888 hex_iterator hex = begin_hex(), endh = end_hex();
12889 std::vector<unsigned int>::const_iterator i = v.begin();
12890 for (; hex != endh; ++hex, ++i)
12891 set_user_index(hex, *i);
12902 template <
typename Iterator>
12904 get_user_pointer(
const Iterator &i)
12906 return i->user_pointer();
12911 template <
int structdim,
int dim,
int spacedim>
12922 template <
typename Iterator>
12924 set_user_pointer(
const Iterator &i,
void *x)
12926 i->set_user_pointer(x);
12931 template <
int structdim,
int dim,
int spacedim>
12942template <
int dim,
int spacedim>
12945 std::vector<void *> &v)
const
12951 std::vector<void *> tmp;
12953 save_user_pointers_line(tmp);
12954 v.insert(v.end(), tmp.begin(), tmp.end());
12958 save_user_pointers_quad(tmp);
12959 v.insert(v.end(), tmp.begin(), tmp.end());
12964 save_user_pointers_hex(tmp);
12965 v.insert(v.end(), tmp.begin(), tmp.end());
12974template <
int dim,
int spacedim>
12977 const std::vector<void *> &v)
12980 std::vector<void *> tmp;
12984 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
12986 load_user_pointers_line(tmp);
12991 tmp.insert(tmp.end(),
12992 v.begin() + n_lines(),
12993 v.begin() + n_lines() + n_quads());
12994 load_user_pointers_quad(tmp);
13000 tmp.insert(tmp.end(),
13001 v.begin() + n_lines() + n_quads(),
13002 v.begin() + n_lines() + n_quads() + n_hexs());
13003 load_user_pointers_hex(tmp);
13012template <
int dim,
int spacedim>
13015 std::vector<void *> &v)
const
13017 v.resize(n_lines(),
nullptr);
13018 std::vector<void *>::iterator i = v.begin();
13019 line_iterator line = begin_line(), endl = end_line();
13020 for (; line != endl; ++line, ++i)
13021 *i = line->user_pointer();
13026template <
int dim,
int spacedim>
13029 const std::vector<void *> &v)
13031 Assert(v.size() == n_lines(), ExcGridReadError());
13033 line_iterator line = begin_line(), endl = end_line();
13034 std::vector<void *>::const_iterator i = v.begin();
13035 for (; line != endl; ++line, ++i)
13036 line->set_user_pointer(*i);
13041template <
int dim,
int spacedim>
13044 std::vector<void *> &v)
const
13046 v.resize(n_quads(),
nullptr);
13050 std::vector<void *>::iterator i = v.begin();
13051 quad_iterator quad = begin_quad(), endq = end_quad();
13052 for (; quad != endq; ++quad, ++i)
13053 *i = get_user_pointer(quad);
13059template <
int dim,
int spacedim>
13062 const std::vector<void *> &v)
13064 Assert(v.size() == n_quads(), ExcGridReadError());
13068 quad_iterator quad = begin_quad(), endq = end_quad();
13069 std::vector<void *>::const_iterator i = v.begin();
13070 for (; quad != endq; ++quad, ++i)
13071 set_user_pointer(quad, *i);
13076template <
int dim,
int spacedim>
13079 std::vector<void *> &v)
const
13081 v.resize(n_hexs(),
nullptr);
13085 std::vector<void *>::iterator i = v.begin();
13086 hex_iterator hex = begin_hex(), endh = end_hex();
13087 for (; hex != endh; ++hex, ++i)
13088 *i = get_user_pointer(hex);
13094template <
int dim,
int spacedim>
13097 const std::vector<void *> &v)
13099 Assert(v.size() == n_hexs(), ExcGridReadError());
13103 hex_iterator hex = begin_hex(), endh = end_hex();
13104 std::vector<void *>::const_iterator i = v.begin();
13105 for (; hex != endh; ++hex, ++i)
13106 set_user_pointer(hex, *i);
13115template <
int dim,
int spacedim>
13123 return begin_raw_line(
level);
13125 return begin_raw_quad(
level);
13127 return begin_raw_hex(
level);
13130 return raw_cell_iterator();
13136template <
int dim,
int spacedim>
13144 return begin_line(
level);
13146 return begin_quad(
level);
13148 return begin_hex(
level);
13151 return cell_iterator();
13157template <
int dim,
int spacedim>
13165 return begin_active_line(
level);
13167 return begin_active_quad(
level);
13169 return begin_active_hex(
level);
13172 return active_cell_iterator();
13178template <
int dim,
int spacedim>
13183 const unsigned int level = levels.size() - 1;
13184 if (levels[
level]->cells.n_objects() == 0)
13191 levels[
level]->cells.n_objects() - 1);
13194 if (ri->used() ==
true)
13197 if (ri->used() ==
true)
13204template <
int dim,
int spacedim>
13210 cell_iterator cell = last();
13215 if (cell->is_active() ==
true)
13218 if (cell->is_active() ==
true)
13226template <
int dim,
int spacedim>
13230 const CellId &cell_id)
const
13232 cell_iterator cell(
13235 for (
const auto &child_index : cell_id.get_child_indices())
13238 cell->has_children(),
13240 "CellId is invalid for this triangulation.\n"
13241 "Either the provided CellId does not correspond to a cell in this "
13242 "triangulation object, or, in case you are using a parallel "
13243 "triangulation, may correspond to an artificial cell that is less "
13244 "refined on this processor."));
13245 cell = cell->child(
static_cast<unsigned int>(child_index));
13253template <
int dim,
int spacedim>
13265template <
int dim,
int spacedim>
13278 if (
level >= levels.size())
13281 ExcInvalidLevel(
level, n_global_levels()));
13288 if (
level < levels.size() - 1)
13289 return begin_raw(
level + 1);
13295template <
int dim,
int spacedim>
13308 if (
level >= levels.size())
13311 ExcInvalidLevel(
level, n_global_levels()));
13318 if (
level < levels.size() - 1)
13325template <
int dim,
int spacedim>
13338 if (
level >= levels.size())
13341 ExcInvalidLevel(
level, n_global_levels()));
13348 return (
level >= levels.size() - 1 ? active_cell_iterator(
end()) :
13349 begin_active(
level + 1));
13354template <
int dim,
int spacedim>
13365template <
int dim,
int spacedim>
13378template <
int dim,
int spacedim>
13390template <
int dim,
int spacedim>
13405template <
int dim,
int spacedim>
13414 return raw_face_iterator();
13416 return begin_line();
13418 return begin_quad();
13421 return face_iterator();
13427template <
int dim,
int spacedim>
13436 return raw_face_iterator();
13438 return begin_active_line();
13440 return begin_active_quad();
13443 return active_face_iterator();
13449template <
int dim,
int spacedim>
13458 return raw_face_iterator();
13465 return raw_face_iterator();
13471template <
int dim,
int spacedim>
13479 begin_active_face(), end_face());
13485template <
int dim,
int spacedim>
13490 vertex_iterator i =
13495 while (i->used() ==
false)
13503template <
int dim,
int spacedim>
13509 return begin_vertex();
13514template <
int dim,
int spacedim>
13530template <
int dim,
int spacedim>
13543 if (
level >= levels.size())
13546 ExcInvalidLevel(
level, n_global_levels()));
13557 if (
level >= levels.size() || levels[
level]->cells.n_objects() == 0)
13560 return raw_line_iterator(
13565 return raw_line_iterator(
13571template <
int dim,
int spacedim>
13577 raw_line_iterator ri = begin_raw_line(
level);
13580 while (ri->used() ==
false)
13588template <
int dim,
int spacedim>
13592 const unsigned int level)
const
13595 line_iterator i = begin_line(
level);
13598 while (i->has_children())
13606template <
int dim,
int spacedim>
13621template <
int dim,
int spacedim>
13634 if (
level >= levels.size())
13637 ExcInvalidLevel(
level, n_global_levels()));
13645 return raw_hex_iterator();
13652 if (
level >= levels.size() || levels[
level]->cells.n_objects() == 0)
13655 return raw_quad_iterator(
13663 return raw_quad_iterator(
13670 return raw_hex_iterator();
13676template <
int dim,
int spacedim>
13682 raw_quad_iterator ri = begin_raw_quad(
level);
13685 while (ri->used() ==
false)
13693template <
int dim,
int spacedim>
13697 const unsigned int level)
const
13700 quad_iterator i = begin_quad(
level);
13703 while (i->has_children())
13711template <
int dim,
int spacedim>
13725template <
int dim,
int spacedim>
13738 if (
level >= levels.size())
13741 ExcInvalidLevel(
level, n_global_levels()));
13750 return raw_hex_iterator();
13757 if (
level >= levels.size() || levels[
level]->cells.n_objects() == 0)
13760 return raw_hex_iterator(
13766 return raw_hex_iterator();
13772template <
int dim,
int spacedim>
13778 raw_hex_iterator ri = begin_raw_hex(
level);
13781 while (ri->used() ==
false)
13789template <
int dim,
int spacedim>
13795 hex_iterator i = begin_hex(
level);
13798 while (i->has_children())
13806template <
int dim,
int spacedim>
13823 namespace TriangulationImplementation
13825 inline unsigned int
13832 inline unsigned int
13836 return c.n_active_lines;
13840 inline unsigned int
13847 inline unsigned int
13851 return c.n_active_quads;
13855 inline unsigned int
13862 inline unsigned int
13866 return c.n_active_hexes;
13873template <
int dim,
int spacedim>
13881template <
int dim,
int spacedim>
13888template <
int dim,
int spacedim>
13896template <
int dim,
int spacedim>
13904template <
int dim,
int spacedim>
13911 return n_used_vertices();
13923template <
int dim,
int spacedim>
13930 return n_vertices();
13932 return n_raw_lines();
13934 return n_raw_quads();
13942template <
int dim,
int spacedim>
13949 return n_used_vertices();
13951 return n_active_lines();
13953 return n_active_quads();
13961template <
int dim,
int spacedim>
13964 const unsigned int level)
const
13969 return n_raw_lines(
level);
13971 return n_raw_quads(
level);
13973 return n_raw_hexs(
level);
13982template <
int dim,
int spacedim>
13985 const unsigned int level)
const
13990 return n_lines(
level);
13992 return n_quads(
level);
13994 return n_hexs(
level);
14003template <
int dim,
int spacedim>
14006 const unsigned int level)
const
14011 return n_active_lines(
level);
14013 return n_active_quads(
level);
14015 return n_active_hexs(
level);
14023template <
int dim,
int spacedim>
14027 if (anisotropic_refinement ==
false)
14029 for (
unsigned int lvl = 0; lvl < n_global_levels() - 1; ++lvl)
14035 for (
const auto &cell : active_cell_iterators())
14036 for (const auto &i : cell->face_indices())
14037 if (cell->face(i)->has_children())
14044template <
int dim,
int spacedim>
14048 return number_cache.n_lines;
14053template <
int dim,
int spacedim>
14056 const unsigned int level)
const
14061 return levels[
level]->cells.n_objects();
14064 Assert(
false, ExcFacesHaveNoLevel());
14069template <
int dim,
int spacedim>
14083template <
int dim,
int spacedim>
14086 const unsigned int level)
const
14089 Assert(dim == 1, ExcFacesHaveNoLevel());
14090 return number_cache.n_lines_level[
level];
14094template <
int dim,
int spacedim>
14098 return number_cache.n_active_lines;
14102template <
int dim,
int spacedim>
14105 const unsigned int level)
const
14108 Assert(dim == 1, ExcFacesHaveNoLevel());
14110 return number_cache.n_active_lines_level[
level];
14260template <
int dim,
int spacedim>
14264 return number_cache.n_quads;
14268template <
int dim,
int spacedim>
14271 const unsigned int level)
const
14273 Assert(dim == 2, ExcFacesHaveNoLevel());
14275 return number_cache.n_quads_level[
level];
14285 return levels[
level]->cells.n_objects();
14295 return levels[
level]->cells.n_objects();
14303 Assert(
false, ExcFacesHaveNoLevel());
14309template <
int dim,
int spacedim>
14328template <
int dim,
int spacedim>
14332 return number_cache.n_active_quads;
14336template <
int dim,
int spacedim>
14339 const unsigned int level)
const
14342 Assert(dim == 2, ExcFacesHaveNoLevel());
14344 return number_cache.n_active_quads_level[
level];
14348template <
int dim,
int spacedim>
14357template <
int dim,
int spacedim>
14366template <
int dim,
int spacedim>
14374template <
int dim,
int spacedim>
14383template <
int dim,
int spacedim>
14386 const unsigned int)
const
14396 return number_cache.n_hexes;
14407 return number_cache.n_hexes_level[
level];
14417 return levels[
level]->cells.n_objects();
14425 return number_cache.n_active_hexes;
14436 return number_cache.n_active_hexes_level[
level];
14441template <
int dim,
int spacedim>
14445 return std::count(vertices_used.begin(), vertices_used.end(),
true);
14450template <
int dim,
int spacedim>
14454 return vertices_used;
14484template <
int dim,
int spacedim>
14488 cell_iterator cell = begin(0),
14489 endc = (n_levels() > 1 ? begin(1) : cell_iterator(end()));
14492 unsigned int max_vertex_index = 0;
14493 for (; cell != endc; ++cell)
14495 if (cell->vertex_index(vertex) > max_vertex_index)
14496 max_vertex_index = cell->vertex_index(vertex);
14502 std::vector<unsigned short int> usage_count(max_vertex_index + 1, 0);
14506 for (cell = begin(); cell != endc; ++cell)
14508 ++usage_count[cell->vertex_index(vertex)];
14511 static_cast<unsigned int>(
14512 *std::max_element(usage_count.begin(), usage_count.end())));
14517template <
int dim,
int spacedim>
14527template <
int dim,
int spacedim>
14536template <
int dim,
int spacedim>
14546template <
int dim,
int spacedim>
14550 &periodicity_vector)
14552 periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(),
14553 periodicity_vector.begin(),
14554 periodicity_vector.end());
14557 update_periodic_face_map();
14562template <
int dim,
int spacedim>
14564const typename std::map<
14565 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
unsigned int>,
14566 std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
14571 return periodic_face_map;
14576template <
int dim,
int spacedim>
14590 if (smooth_grid & limit_level_difference_at_vertices)
14594 signals.pre_refinement();
14596 execute_coarsening();
14598 const DistortedCellList cells_with_distorted_children = execute_refinement();
14600 reset_cell_vertex_indices_cache();
14605 if (smooth_grid & limit_level_difference_at_vertices)
14610 this->policy->update_neighbors(*
this);
14611 reset_active_cell_indices();
14613 reset_global_cell_indices();
14616 signals.post_refinement();
14618 AssertThrow(cells_with_distorted_children.distorted_cells.size() == 0,
14619 cells_with_distorted_children);
14621 update_periodic_face_map();
14626template <
int dim,
int spacedim>
14630 unsigned int active_cell_index = 0;
14631 for (raw_cell_iterator cell = begin_raw(); cell !=
end(); ++cell)
14632 if ((cell->used() ==
false) || cell->has_children())
14636 cell->set_active_cell_index(active_cell_index);
14637 ++active_cell_index;
14645template <
int dim,
int spacedim>
14651 for (
const auto &cell : active_cell_iterators())
14652 cell->set_global_active_cell_index(
cell_index++);
14655 for (
unsigned int l = 0;
l < levels.size(); ++
l)
14658 for (
const auto &cell : cell_iterators_on_level(
l))
14659 cell->set_global_level_cell_index(
cell_index++);
14665template <
int dim,
int spacedim>
14669 for (
unsigned int l = 0;
l < levels.size(); ++
l)
14671 constexpr unsigned int max_vertices_per_cell = 1 << dim;
14672 std::vector<unsigned int> &cache = levels[
l]->cell_vertex_indices_cache;
14674 cache.resize(levels[l]->refine_flags.size() * max_vertices_per_cell,
14676 for (
const auto &cell : cell_iterators_on_level(
l))
14678 const unsigned int my_index = cell->index() * max_vertices_per_cell;
14685 for (
unsigned int face = 4; face < 6; ++face)
14687 const auto face_iter = cell->face(face);
14688 const std::array<bool, 2> line_orientations{
14689 {face_iter->line_orientation(0),
14690 face_iter->line_orientation(1)}};
14691 std::array<unsigned int, 4> raw_vertex_indices{
14692 {face_iter->line(0)->vertex_index(1 - line_orientations[0]),
14693 face_iter->line(1)->vertex_index(1 - line_orientations[1]),
14694 face_iter->line(0)->vertex_index(line_orientations[0]),
14695 face_iter->line(1)->vertex_index(line_orientations[1])}};
14697 const unsigned char combined_orientation =
14700 std::array<unsigned int, 4> vertex_order{
14703 combined_orientation),
14706 combined_orientation),
14709 combined_orientation),
14711 3, face, combined_orientation)}};
14713 const unsigned int index = my_index + 4 * (face - 4);
14714 for (
unsigned int i = 0; i < 4; ++i)
14715 cache[
index + i] = raw_vertex_indices[vertex_order[i]];
14719 const std::array<bool, 2> line_orientations{
14720 {cell->line_orientation(0), cell->line_orientation(1)}};
14721 std::array<unsigned int, 4> raw_vertex_indices{
14722 {cell->line(0)->vertex_index(1 - line_orientations[0]),
14723 cell->line(1)->vertex_index(1 - line_orientations[1]),
14724 cell->line(0)->vertex_index(line_orientations[0]),
14725 cell->line(1)->vertex_index(line_orientations[1])}};
14726 for (
unsigned int i = 0; i < 4; ++i)
14727 cache[my_index + i] = raw_vertex_indices[i];
14731 cache[my_index + i] =
internal::TriaAccessorImplementation::
14732 Implementation::vertex_index(*cell, i);
14739template <
int dim,
int spacedim>
14744 periodic_face_map.clear();
14746 typename std::vector<
14748 for (it = periodic_face_pairs_level_0.begin();
14749 it != periodic_face_pairs_level_0.end();
14752 update_periodic_face_map_recursively<dim, spacedim>(it->cell[0],
14757 periodic_face_map);
14760 std::bitset<3> inverted_orientation;
14762 bool orientation, flip, rotation;
14763 orientation = it->orientation[0];
14764 rotation = it->orientation[2];
14765 flip = orientation ? rotation ^ it->orientation[1] : it->orientation[1];
14766 inverted_orientation[0] = orientation;
14767 inverted_orientation[1] = flip;
14768 inverted_orientation[2] = rotation;
14770 update_periodic_face_map_recursively<dim, spacedim>(it->cell[1],
14774 inverted_orientation,
14775 periodic_face_map);
14779 typename std::map<std::pair<cell_iterator, unsigned int>,
14780 std::pair<std::pair<cell_iterator, unsigned int>,
14781 std::bitset<3>>>::const_iterator it_test;
14782 for (it_test = periodic_face_map.begin(); it_test != periodic_face_map.end();
14786 it_test->first.first;
14788 it_test->second.first.first;
14789 if (cell_1->level() == cell_2->level())
14793 Assert(periodic_face_map[it_test->second.first].first ==
14802template <
int dim,
int spacedim>
14806 std::set<ReferenceCell> reference_cells_set;
14807 for (
auto cell : active_cell_iterators())
14808 if (cell->is_locally_owned())
14811 this->reference_cells =
14812 std::vector<ReferenceCell>(reference_cells_set.begin(),
14813 reference_cells_set.end());
14818template <
int dim,
int spacedim>
14820const std::vector<ReferenceCell>
14823 return this->reference_cells;
14828template <
int dim,
int spacedim>
14832 Assert(this->reference_cells.size() > 0,
14833 ExcMessage(
"You can't ask about the kinds of reference "
14834 "cells used by this triangulation if the "
14835 "triangulation doesn't yet have any cells in it."));
14836 return (this->reference_cells.size() == 1 &&
14837 this->reference_cells[0].is_hyper_cube());
14842template <
int dim,
int spacedim>
14846 Assert(this->reference_cells.size() > 0,
14847 ExcMessage(
"You can't ask about the kinds of reference "
14848 "cells used by this triangulation if the "
14849 "triangulation doesn't yet have any cells in it."));
14850 return (this->reference_cells.size() == 1 &&
14851 this->reference_cells[0].is_simplex());
14856template <
int dim,
int spacedim>
14860 Assert(this->reference_cells.size() > 0,
14861 ExcMessage(
"You can't ask about the kinds of reference "
14862 "cells used by this triangulation if the "
14863 "triangulation doesn't yet have any cells in it."));
14864 return reference_cells.size() > 1 ||
14865 ((reference_cells[0].is_hyper_cube() ==
false) &&
14866 (reference_cells[0].is_simplex() ==
false));
14871template <
int dim,
int spacedim>
14879 vertices_used.clear();
14888template <
int dim,
int spacedim>
14893 const DistortedCellList cells_with_distorted_children =
14894 this->policy->execute_refinement(*
this, check_for_distorted_cells);
14900 *
this, levels.size(), number_cache);
14903 for (
const auto &
level : levels)
14910 for (
const auto &cell : this->cell_iterators())
14914 return cells_with_distorted_children;
14919template <
int dim,
int spacedim>
14925 const cell_iterator endc =
end();
14926 bool do_coarsen =
false;
14927 if (levels.size() >= 2)
14928 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
14929 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
14941 std::vector<unsigned int> line_cell_count(dim > 1 ? this->n_raw_lines() : 0);
14942 std::vector<unsigned int> quad_cell_count(dim > 2 ? this->n_raw_quads() : 0);
14944 for (
const auto &cell : this->cell_iterators())
14948 const auto line_indices = internal::TriaAccessorImplementation::
14949 Implementation::get_line_indices_of_cell(*cell);
14952 const unsigned int n_lines =
std::min(cell->n_lines(), 12u);
14953 for (
unsigned int l = 0;
l < n_lines; ++
l)
14954 ++line_cell_count[line_indices[l]];
14955 for (
const unsigned int q : cell->face_indices())
14956 ++quad_cell_count[cell->face_index(q)];
14959 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
14960 ++line_cell_count[cell->line(l)->index()];
14977 if (levels.size() >= 2)
14978 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
14979 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
14981 for (
unsigned int child = 0; child < cell->n_children(); ++child)
14983 Assert(cell->child(child)->coarsen_flag_set(),
14985 cell->child(child)->clear_coarsen_flag();
14988 signals.pre_coarsening_on_cell(cell);
14990 this->policy->delete_children(*
this,
14998 *
this, levels.size(), number_cache);
15003template <
int dim,
int spacedim>
15020 auto previous_coarsen_flags = internal::extract_raw_coarsen_flags(levels);
15022 bool continue_iterating =
true;
15026 if (smooth_grid & limit_level_difference_at_vertices)
15028 Assert(!anisotropic_refinement,
15029 ExcMessage(
"In case of anisotropic refinement the "
15030 "limit_level_difference_at_vertices flag for "
15031 "mesh smoothing must not be set!"));
15035 std::vector<int> vertex_level(
vertices.size(), 0);
15036 for (
const auto &cell : this->active_cell_iterators())
15038 if (cell->refine_flag_set())
15039 for (
const unsigned int vertex :
15041 vertex_level[cell->vertex_index(vertex)] =
15042 std::
max(vertex_level[cell->vertex_index(vertex)],
15043 cell->
level() + 1);
15044 else if (!cell->coarsen_flag_set())
15045 for (
const unsigned int vertex :
15047 vertex_level[cell->vertex_index(vertex)] =
15048 std::
max(vertex_level[cell->vertex_index(vertex)],
15059 for (
const unsigned int vertex :
15061 vertex_level[cell->vertex_index(vertex)] =
15062 std::
max(vertex_level[cell->vertex_index(vertex)],
15063 cell->
level() - 1);
15077 active_cell_iterator endc =
end();
15078 for (active_cell_iterator cell = last_active(); cell != endc; --cell)
15079 if (cell->refine_flag_set() ==
false)
15081 for (
const unsigned int vertex :
15083 if (vertex_level[cell->vertex_index(vertex)] >=
15087 cell->clear_coarsen_flag();
15092 if (vertex_level[cell->vertex_index(vertex)] >
15095 cell->set_refine_flag();
15097 for (
const unsigned int v :
15099 vertex_level[cell->vertex_index(v)] =
15100 std::
max(vertex_level[cell->vertex_index(v)],
15101 cell->
level() + 1);
15117 for (
const auto &acell : this->active_cell_iterators_on_level(0))
15118 acell->clear_coarsen_flag();
15120 const cell_iterator endc =
end();
15121 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
15124 if (cell->is_active())
15127 const unsigned int n_children = cell->n_children();
15128 unsigned int flagged_children = 0;
15129 for (
unsigned int child = 0; child < n_children; ++child)
15131 const auto child_cell = cell->child(child);
15132 if (child_cell->is_active() && child_cell->coarsen_flag_set())
15134 ++flagged_children;
15136 child_cell->clear_coarsen_flag();
15142 if (flagged_children == n_children &&
15143 this->policy->coarsening_allowed(cell))
15144 for (
unsigned int c = 0; c < n_children; ++c)
15146 Assert(cell->child(c)->refine_flag_set() ==
false,
15149 cell->child(c)->set_coarsen_flag();
15155 auto current_coarsen_flags = internal::extract_raw_coarsen_flags(levels);
15157 continue_iterating = (current_coarsen_flags != previous_coarsen_flags);
15158 previous_coarsen_flags.swap(current_coarsen_flags);
15160 while (continue_iterating ==
true);
15172 const auto flags_before = internal::extract_raw_coarsen_flags(levels);
15175 fix_coarsen_flags();
15177 const auto flags_after = internal::extract_raw_coarsen_flags(levels);
15179 return (flags_before != flags_after);
15190 const auto flags_before = internal::extract_raw_coarsen_flags(levels);
15193 fix_coarsen_flags();
15195 const auto flags_after = internal::extract_raw_coarsen_flags(levels);
15197 return (flags_before != flags_after);
15208 const auto flags_before = internal::extract_raw_coarsen_flags(levels);
15211 fix_coarsen_flags();
15213 const auto flags_after = internal::extract_raw_coarsen_flags(levels);
15215 return (flags_before != flags_after);
15226 template <
int dim,
int spacedim>
15228 possibly_do_not_produce_unrefined_islands(
15233 unsigned int n_neighbors = 0;
15236 unsigned int count = 0;
15237 for (
const unsigned int n :
GeometryInfo<dim>::face_indices())
15244 if (face_will_be_refined_by_neighbor(cell, n))
15251 if (count == n_neighbors ||
15252 (count >= n_neighbors - 1 &&
15255 for (
unsigned int c = 0; c < cell->n_children(); ++c)
15256 cell->child(c)->clear_coarsen_flag();
15258 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
15259 if (!cell->at_boundary(face) &&
15260 (!cell->neighbor(face)->is_active()) &&
15261 (cell_will_be_coarsened(cell->neighbor(face))))
15262 possibly_do_not_produce_unrefined_islands<dim, spacedim>(
15263 cell->neighbor(face));
15278 template <
int dim,
int spacedim>
15280 possibly_refine_unrefined_island(
15282 const bool allow_anisotropic_smoothing)
15297 if (allow_anisotropic_smoothing ==
false)
15300 unsigned int refined_neighbors = 0, unrefined_neighbors = 0;
15301 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
15302 if (!cell->at_boundary(face))
15304 if (face_will_be_refined_by_neighbor(cell, face))
15305 ++refined_neighbors;
15307 ++unrefined_neighbors;
15310 if (unrefined_neighbors < refined_neighbors)
15312 cell->clear_coarsen_flag();
15313 cell->set_refine_flag();
15318 if (unrefined_neighbors > 0)
15319 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
15320 if (!cell->at_boundary(face) &&
15321 (face_will_be_refined_by_neighbor(cell, face) == false) &&
15322 (cell->neighbor(face)->has_children() == false) &&
15323 (cell->neighbor(face)->refine_flag_set() == false))
15324 possibly_refine_unrefined_island<dim, spacedim>(
15325 cell->neighbor(face), allow_anisotropic_smoothing);
15337 for (
unsigned int face_pair = 0;
15338 face_pair < GeometryInfo<dim>::faces_per_cell / 2;
15347 for (
unsigned int face_index = 0; face_index < 2; ++face_index)
15349 unsigned int face = 2 * face_pair + face_index;
15356 face_will_be_refined_by_neighbor<dim, spacedim>(
15357 cell, face, expected_face_ref_case);
15371 directional_cell_refinement_case =
15372 (directional_cell_refinement_case &
15375 expected_face_ref_case,
15377 cell->face_orientation(face),
15378 cell->face_flip(face),
15379 cell->face_rotation(face)));
15385 Assert(directional_cell_refinement_case <
15388 smoothing_cell_refinement_case =
15389 smoothing_cell_refinement_case | directional_cell_refinement_case;
15394 if (smoothing_cell_refinement_case)
15396 cell->clear_coarsen_flag();
15397 cell->set_refine_flag(cell->refine_flag_set() |
15398 smoothing_cell_refinement_case);
15405template <
int dim,
int spacedim>
15411 const auto coarsen_flags_before = internal::extract_raw_coarsen_flags(levels);
15412 const auto refine_flags_before = internal::extract_raw_refine_flags(levels);
15430 auto coarsen_flags_before_loop = coarsen_flags_before;
15431 auto refine_flags_before_loop = refine_flags_before;
15508 if (((smooth_grid & coarsest_level_1) || (smooth_grid & patch_level_1)) &&
15511 for (
const auto &cell : active_cell_iterators_on_level(1))
15512 cell->clear_coarsen_flag();
15515 bool mesh_changed_in_this_loop =
false;
15524 if (smooth_grid & do_not_produce_unrefined_islands &&
15525 !(smooth_grid & patch_level_1))
15527 for (
const auto &cell : cell_iterators())
15531 if (!cell->is_active() && cell_will_be_coarsened(cell))
15532 possibly_do_not_produce_unrefined_islands<dim, spacedim>(cell);
15556 if (smooth_grid & (eliminate_refined_inner_islands |
15557 eliminate_refined_boundary_islands) &&
15558 !(smooth_grid & patch_level_1))
15560 for (
const auto &cell : cell_iterators())
15561 if (!cell->is_active() ||
15562 (cell->is_active() && cell->refine_flag_set() &&
15563 cell->is_locally_owned()))
15570 bool all_children_active =
true;
15571 if (!cell->is_active())
15572 for (
unsigned int c = 0; c < cell->n_children(); ++c)
15573 if (!cell->child(c)->is_active() ||
15574 cell->child(c)->is_ghost() ||
15575 cell->child(c)->is_artificial())
15577 all_children_active =
false;
15581 if (all_children_active)
15588 unsigned int unrefined_neighbors = 0, total_neighbors = 0;
15599 bool at_periodic_boundary =
false;
15601 for (
const unsigned int n :
15604 const cell_iterator neighbor = cell->neighbor(n);
15609 if (!face_will_be_refined_by_neighbor(cell, n))
15610 ++unrefined_neighbors;
15612 else if (cell->has_periodic_neighbor(n))
15615 at_periodic_boundary =
true;
15630 if ((unrefined_neighbors == total_neighbors) &&
15631 ((!cell->at_boundary() &&
15632 (smooth_grid & eliminate_refined_inner_islands)) ||
15633 (cell->at_boundary() && !at_periodic_boundary &&
15635 eliminate_refined_boundary_islands))) &&
15636 (total_neighbors != 0))
15638 if (!cell->is_active())
15639 for (
unsigned int c = 0; c < cell->n_children(); ++c)
15641 cell->child(c)->clear_refine_flag();
15642 cell->child(c)->set_coarsen_flag();
15645 cell->clear_refine_flag();
15660 if (smooth_grid & limit_level_difference_at_vertices)
15662 Assert(!anisotropic_refinement,
15663 ExcMessage(
"In case of anisotropic refinement the "
15664 "limit_level_difference_at_vertices flag for "
15665 "mesh smoothing must not be set!"));
15669 std::vector<int> vertex_level(
vertices.size(), 0);
15670 for (
const auto &cell : active_cell_iterators())
15672 if (cell->refine_flag_set())
15673 for (
const unsigned int vertex :
15675 vertex_level[cell->vertex_index(vertex)] =
15676 std::
max(vertex_level[cell->vertex_index(vertex)],
15677 cell->
level() + 1);
15678 else if (!cell->coarsen_flag_set())
15679 for (
const unsigned int vertex :
15681 vertex_level[cell->vertex_index(vertex)] =
15682 std::
max(vertex_level[cell->vertex_index(vertex)],
15691 for (
const unsigned int vertex :
15693 vertex_level[cell->vertex_index(vertex)] =
15694 std::
max(vertex_level[cell->vertex_index(vertex)],
15695 cell->
level() - 1);
15709 for (active_cell_iterator cell = last_active(); cell !=
end(); --cell)
15710 if (cell->refine_flag_set() ==
false)
15712 for (
const unsigned int vertex :
15714 if (vertex_level[cell->vertex_index(vertex)] >=
15718 cell->clear_coarsen_flag();
15723 if (vertex_level[cell->vertex_index(vertex)] >
15726 cell->set_refine_flag();
15728 for (
const unsigned int v :
15730 vertex_level[cell->vertex_index(v)] =
15731 std::
max(vertex_level[cell->vertex_index(v)],
15732 cell->
level() + 1);
15753 if (smooth_grid & eliminate_unrefined_islands)
15755 for (active_cell_iterator cell = last_active(); cell !=
end(); --cell)
15758 if (cell->refine_flag_set() !=
15760 possibly_refine_unrefined_island<dim, spacedim>(
15761 cell, (smooth_grid & allow_anisotropic_smoothing) != 0);
15789 if (smooth_grid & patch_level_1)
15799 for (
const auto &cell : cell_iterators())
15800 if (!cell->is_active())
15806 if (cell->child(0)->has_children() ==
true)
15813 for (
unsigned int i = 0; i < cell->n_children(); ++i)
15814 combined_ref_case =
15815 combined_ref_case | cell->child(i)->refine_flag_set();
15817 for (
unsigned int i = 0; i < cell->n_children(); ++i)
15819 cell_iterator child = cell->child(i);
15821 child->clear_coarsen_flag();
15822 child->set_refine_flag(combined_ref_case);
15835 for (
const auto &cell : cell_iterators())
15841 if (cell->is_active() || cell->child(0)->is_active())
15848 const unsigned int n_children = cell->n_children();
15849 bool has_active_grandchildren =
false;
15851 for (
unsigned int i = 0; i < n_children; ++i)
15852 if (cell->child(i)->child(0)->is_active())
15854 has_active_grandchildren =
true;
15858 if (has_active_grandchildren ==
false)
15864 unsigned int n_grandchildren = 0;
15867 unsigned int n_coarsen_flags = 0;
15872 for (
unsigned int c = 0; c < n_children; ++c)
15877 cell_iterator child = cell->child(c);
15879 const unsigned int nn_children = child->n_children();
15880 n_grandchildren += nn_children;
15885 if (child->child(0)->is_active())
15886 for (
unsigned int cc = 0; cc < nn_children; ++cc)
15887 if (child->child(cc)->coarsen_flag_set())
15902 if ((n_coarsen_flags != n_grandchildren) && (n_coarsen_flags > 0))
15903 for (
unsigned int c = 0; c < n_children; ++c)
15905 const cell_iterator child = cell->child(c);
15906 if (child->child(0)->is_active())
15907 for (
unsigned int cc = 0; cc < child->n_children(); ++cc)
15908 child->child(cc)->clear_coarsen_flag();
15919 this->policy->prevent_distorted_boundary_cells(*
this);
15935 bool changed =
true;
15939 active_cell_iterator cell = last_active(), endc =
end();
15941 for (; cell != endc; --cell)
15942 if (cell->refine_flag_set())
15945 for (
const auto i : cell->face_indices())
15950 const bool has_periodic_neighbor =
15951 cell->has_periodic_neighbor(i);
15952 const bool has_neighbor_or_periodic_neighbor =
15953 !cell->at_boundary(i) || has_periodic_neighbor;
15954 if (has_neighbor_or_periodic_neighbor &&
15956 cell->refine_flag_set(), i) !=
15969 if (cell->neighbor_or_periodic_neighbor(i)->is_active())
15971 if ((!has_periodic_neighbor &&
15972 cell->neighbor_is_coarser(i)) ||
15973 (has_periodic_neighbor &&
15974 cell->periodic_neighbor_is_coarser(i)))
15976 if (cell->neighbor_or_periodic_neighbor(i)
15977 ->coarsen_flag_set())
15978 cell->neighbor_or_periodic_neighbor(i)
15979 ->clear_coarsen_flag();
15996 allow_anisotropic_smoothing)
15998 has_periodic_neighbor ?
15999 cell->periodic_neighbor(i)
16000 ->flag_for_face_refinement(
16002 ->periodic_neighbor_of_coarser_periodic_neighbor(
16007 ->flag_for_face_refinement(
16009 ->neighbor_of_coarser_neighbor(
16016 ->neighbor_or_periodic_neighbor(
16018 ->refine_flag_set())
16020 cell->neighbor_or_periodic_neighbor(i)
16021 ->set_refine_flag();
16153 std::pair<unsigned int, unsigned int>
16155 has_periodic_neighbor ?
16157 ->periodic_neighbor_of_coarser_periodic_neighbor(
16159 cell->neighbor_of_coarser_neighbor(i);
16160 unsigned int refined_along_x = 0,
16161 refined_along_y = 0,
16162 to_be_refined_along_x = 0,
16163 to_be_refined_along_y = 0;
16165 const int this_face_index =
16166 cell->face_index(i);
16175 const auto parent_face = [&]() {
16176 if (has_periodic_neighbor)
16178 const auto neighbor =
16179 cell->periodic_neighbor(i);
16180 const auto parent_face_no =
16182 ->periodic_neighbor_of_periodic_neighbor(
16185 neighbor->periodic_neighbor(
16187 return parent->face(parent_face_no);
16190 return cell->neighbor(i)->face(
16194 if ((this_face_index ==
16195 parent_face->child_index(0)) ||
16196 (this_face_index ==
16197 parent_face->child_index(1)))
16206 parent_face->refinement_case();
16224 cell->refine_flag_set(),
16226 cell->face_orientation(i),
16227 cell->face_flip(i),
16228 cell->face_rotation(i));
16231 ++to_be_refined_along_x;
16234 ++to_be_refined_along_y;
16239 allow_anisotropic_smoothing) ||
16240 cell->neighbor_or_periodic_neighbor(i)
16241 ->refine_flag_set())
16243 if (refined_along_x +
16244 to_be_refined_along_x >
16248 ->neighbor_or_periodic_neighbor(i)
16249 ->flag_for_face_refinement(
16253 if (refined_along_y +
16254 to_be_refined_along_y >
16258 ->neighbor_or_periodic_neighbor(i)
16259 ->flag_for_face_refinement(
16267 ->neighbor_or_periodic_neighbor(i)
16268 ->refine_flag_set() !=
16270 dim>::isotropic_refinement)
16272 cell->neighbor_or_periodic_neighbor(i)
16273 ->set_refine_flag();
16279 cell->neighbor_or_periodic_neighbor(i);
16282 nb->refine_flag_set(),
16284 nb->face_orientation(nb_indices.first),
16285 nb->face_flip(nb_indices.first),
16286 nb->face_rotation(nb_indices.first));
16288 !((refined_along_x != 0u) ||
16289 (to_be_refined_along_x != 0u)))
16290 changed |= cell->flag_for_face_refinement(
16294 !((refined_along_y != 0u) ||
16295 (to_be_refined_along_y != 0u)))
16296 changed |= cell->flag_for_face_refinement(
16303 cell->neighbor_or_periodic_neighbor(i)
16304 ->clear_coarsen_flag();
16305 const unsigned int nb_nb =
16306 has_periodic_neighbor ?
16308 ->periodic_neighbor_of_periodic_neighbor(
16310 cell->neighbor_of_neighbor(i);
16311 const cell_iterator neighbor =
16312 cell->neighbor_or_periodic_neighbor(i);
16315 neighbor->refine_flag_set(),
16317 neighbor->face_orientation(nb_nb),
16318 neighbor->face_flip(nb_nb),
16319 neighbor->face_rotation(nb_nb));
16322 cell->refine_flag_set(),
16324 cell->face_orientation(i),
16325 cell->face_flip(i),
16326 cell->face_rotation(i));
16331 if ((face_ref_case ==
16333 needed_face_ref_case ==
16337 needed_face_ref_case ==
16340 changed = cell->flag_for_face_refinement(
16342 neighbor->flag_for_face_refinement(
16343 nb_nb, needed_face_ref_case);
16350 face_ref_case = cell->face(i)->refinement_case(),
16351 needed_face_ref_case =
16353 cell->refine_flag_set(),
16355 cell->face_orientation(i),
16356 cell->face_flip(i),
16357 cell->face_rotation(i));
16362 needed_face_ref_case ==
16365 needed_face_ref_case ==
16368 cell->flag_for_face_refinement(i,
16380 this->policy->prepare_refinement_dim_dependent(*
this);
16386 fix_coarsen_flags();
16389 auto coarsen_flags_after_loop =
16390 internal::extract_raw_coarsen_flags(levels);
16391 auto refine_flags_after_loop = internal::extract_raw_refine_flags(levels);
16394 mesh_changed_in_this_loop =
16395 ((coarsen_flags_before_loop != coarsen_flags_after_loop) ||
16396 (refine_flags_before_loop != refine_flags_after_loop));
16399 coarsen_flags_before_loop.swap(coarsen_flags_after_loop);
16400 refine_flags_before_loop.swap(refine_flags_after_loop);
16402 while (mesh_changed_in_this_loop);
16408 return ((coarsen_flags_before != coarsen_flags_before_loop) ||
16409 (refine_flags_before != refine_flags_before_loop));
16414template <
int dim,
int spacedim>
16417 const unsigned int magic_number1,
16418 const std::vector<bool> &v,
16419 const unsigned int magic_number2,
16420 std::ostream & out)
16422 const unsigned int N = v.size();
16423 unsigned char * flags =
new unsigned char[N / 8 + 1];
16424 for (
unsigned int i = 0; i < N / 8 + 1; ++i)
16427 for (
unsigned int position = 0; position < N; ++position)
16428 flags[position / 8] |= (v[position] ? (1 << (position % 8)) : 0);
16437 out << magic_number1 <<
' ' << N << std::endl;
16438 for (
unsigned int i = 0; i < N / 8 + 1; ++i)
16439 out <<
static_cast<unsigned int>(flags[i]) <<
' ';
16441 out << std::endl << magic_number2 << std::endl;
16449template <
int dim,
int spacedim>
16452 const unsigned int magic_number1,
16453 std::vector<bool> &v,
16454 const unsigned int magic_number2,
16459 unsigned int magic_number;
16460 in >> magic_number;
16461 AssertThrow(magic_number == magic_number1, ExcGridReadError());
16467 unsigned char * flags =
new unsigned char[N / 8 + 1];
16468 unsigned short int tmp;
16469 for (
unsigned int i = 0; i < N / 8 + 1; ++i)
16475 for (
unsigned int position = 0; position != N; ++position)
16476 v[position] = ((flags[position / 8] & (1 << (position % 8))) != 0);
16478 in >> magic_number;
16479 AssertThrow(magic_number == magic_number2, ExcGridReadError());
16488template <
int dim,
int spacedim>
16492 std::size_t mem = 0;
16493 mem +=
sizeof(MeshSmoothing);
16497 for (
const auto &
level : levels)
16501 mem +=
sizeof(manifolds);
16502 mem +=
sizeof(smooth_grid);
16504 mem +=
sizeof(faces);
16513template <
int dim,
int spacedim>
16520#include "tria.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
types::coarse_cell_id get_coarse_cell_id() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const =0
static constexpr unsigned char default_combined_face_orientation()
unsigned int standard_to_real_face_vertex(const unsigned int vertex, const unsigned int face, const unsigned char face_orientation) const
static constexpr unsigned char reversed_combined_line_orientation()
unsigned int n_lines() const
unsigned char get_combined_orientation(const ArrayView< const T > &vertices_0, const ArrayView< const T > &vertices_1) const
Subscriptor & operator=(const Subscriptor &)
IteratorState::IteratorStates state() const
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator > > &)
virtual types::global_cell_index n_global_active_cells() const
quad_iterator begin_quad(const unsigned int level=0) const
typename IteratorSelector::raw_line_iterator raw_line_iterator
active_vertex_iterator begin_active_vertex() const
virtual MPI_Comm get_communicator() const
void load_user_indices_quad(const std::vector< unsigned int > &v)
unsigned int n_quads() const
Triangulation & operator=(Triangulation< dim, spacedim > &&tria) noexcept
void load_user_indices(const std::vector< unsigned int > &v)
std::vector< bool > vertices_used
bool anisotropic_refinement
active_quad_iterator begin_active_quad(const unsigned int level=0) const
bool get_anisotropic_refinement_flag() const
virtual const MeshSmoothing & get_mesh_smoothing() const
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
virtual types::coarse_cell_id n_global_coarse_cells() const
std::unique_ptr< std::map< unsigned int, types::manifold_id > > vertex_to_manifold_id_map_1d
void save_user_pointers_quad(std::vector< void * > &v) const
void save_user_flags_hex(std::ostream &out) const
void clear_user_flags_quad()
unsigned int n_faces() const
active_hex_iterator begin_active_hex(const unsigned int level=0) const
static void read_bool_vector(const unsigned int magic_number1, std::vector< bool > &v, const unsigned int magic_number2, std::istream &in)
bool all_reference_cells_are_hyper_cube() const
void load_user_flags_line(std::istream &in)
raw_hex_iterator begin_raw_hex(const unsigned int level=0) const
void save_user_flags_line(std::ostream &out) const
active_cell_iterator last_active() const
void reset_global_cell_indices()
face_iterator end_face() const
void reset_active_cell_indices()
cell_iterator create_cell_iterator(const CellId &cell_id) const
cell_iterator begin(const unsigned int level=0) const
std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, std::bitset< 3 > > > periodic_face_map
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, std::bitset< 3 > > > & get_periodic_face_map() const
void save_user_pointers_line(std::vector< void * > &v) const
void load_refine_flags(std::istream &in)
void save_user_indices_line(std::vector< unsigned int > &v) const
raw_cell_iterator begin_raw(const unsigned int level=0) const
unsigned int n_lines() const
virtual void set_mesh_smoothing(const MeshSmoothing mesh_smoothing)
unsigned int n_raw_lines() const
virtual std::size_t memory_consumption() const
std::vector< Point< spacedim > > vertices
raw_quad_iterator begin_raw_quad(const unsigned int level=0) const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_raw_faces() const
unsigned int n_active_faces() const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
const bool check_for_distorted_cells
raw_cell_iterator end_raw(const unsigned int level) const
line_iterator end_line() const
std::unique_ptr< std::map< unsigned int, types::boundary_id > > vertex_to_boundary_id_map_1d
void load_user_flags_quad(std::istream &in)
unsigned int n_active_cells() const
virtual void update_reference_cells()
std::vector< ReferenceCell > reference_cells
void update_periodic_face_map()
void clear_despite_subscriptions()
void coarsen_global(const unsigned int times=1)
Triangulation(const MeshSmoothing smooth_grid=none, const bool check_for_distorted_cells=false)
void save_user_flags(std::ostream &out) const
void refine_global(const unsigned int times=1)
void load_user_flags_hex(std::istream &in)
void load_user_pointers_quad(const std::vector< void * > &v)
virtual void create_triangulation_compatibility(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
std::unique_ptr<::internal::TriangulationImplementation::TriaFaces > faces
unsigned int n_used_vertices() const
void reset_cell_vertex_indices_cache()
unsigned int n_active_lines() const
void load_user_indices_line(const std::vector< unsigned int > &v)
void clear_user_flags_hex()
void save_user_pointers_hex(std::vector< void * > &v) const
const std::vector< ReferenceCell > & get_reference_cells() const
typename IteratorSelector::raw_quad_iterator raw_quad_iterator
void load_user_pointers(const std::vector< void * > &v)
::internal::TriangulationImplementation::NumberCache< dim > number_cache
void save_user_indices_hex(std::vector< unsigned int > &v) const
DistortedCellList execute_refinement()
active_line_iterator begin_active_line(const unsigned int level=0) const
void save_user_indices_quad(std::vector< unsigned int > &v) const
void load_user_pointers_hex(const std::vector< void * > &v)
cell_iterator end() const
virtual bool has_hanging_nodes() const
std::vector< GridTools::PeriodicFacePair< cell_iterator > > periodic_face_pairs_level_0
unsigned int n_raw_cells(const unsigned int level) const
void load_coarsen_flags(std::istream &out)
quad_iterator end_quad() const
line_iterator begin_line(const unsigned int level=0) const
unsigned int max_adjacent_cells() const
vertex_iterator begin_vertex() const
unsigned int n_hexs() const
vertex_iterator end_vertex() const
void load_user_pointers_line(const std::vector< void * > &v)
hex_iterator end_hex() const
hex_iterator begin_hex(const unsigned int level=0) const
virtual void execute_coarsening_and_refinement()
active_cell_iterator end_active(const unsigned int level) const
bool is_mixed_mesh() const
cell_iterator last() const
unsigned int n_active_quads() const
void load_user_indices_hex(const std::vector< unsigned int > &v)
unsigned int n_raw_quads() const
void save_user_pointers(std::vector< void * > &v) const
face_iterator begin_face() const
unsigned int n_cells() const
virtual bool prepare_coarsening_and_refinement()
const std::vector< bool > & get_used_vertices() const
typename IteratorSelector::raw_hex_iterator raw_hex_iterator
std::map< types::manifold_id, std::unique_ptr< const Manifold< dim, spacedim > > > manifolds
MeshSmoothing smooth_grid
void save_refine_flags(std::ostream &out) const
std::unique_ptr< ::internal::TriangulationImplementation::Policy< dim, spacedim > > policy
Triangulation< dim, spacedim > & get_triangulation()
void save_user_flags_quad(std::ostream &out) const
virtual ~Triangulation() override
unsigned int n_vertices() const
void save_user_indices(std::vector< unsigned int > &v) const
bool all_reference_cells_are_simplex() const
std::vector< std::unique_ptr<::internal::TriangulationImplementation::TriaLevel > > levels
unsigned int n_raw_hexs(const unsigned int level) const
void set_all_refine_flags()
unsigned int n_active_hexs() const
virtual std::vector< types::boundary_id > get_boundary_ids() const
void load_user_flags(std::istream &in)
void save_coarsen_flags(std::ostream &out) const
active_face_iterator begin_active_face() const
void clear_user_flags_line()
raw_line_iterator begin_raw_line(const unsigned int level=0) const
static void write_bool_vector(const unsigned int magic_number1, const std::vector< bool > &v, const unsigned int magic_number2, std::ostream &out)
void flip_all_direction_flags()
active_cell_iterator begin_active(const unsigned int level=0) const
void execute_coarsening()
void prevent_distorted_boundary_cells(Triangulation< dim, spacedim > &triangulation) override
void prepare_refinement_dim_dependent(Triangulation< dim, spacedim > &triangulation) override
void delete_children(Triangulation< dim, spacedim > &tria, typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< unsigned int > &line_cell_count, std::vector< unsigned int > &quad_cell_count) override
void update_neighbors(Triangulation< dim, spacedim > &tria) override
bool coarsening_allowed(const typename Triangulation< dim, spacedim >::cell_iterator &cell) override
Triangulation< dim, spacedim >::DistortedCellList execute_refinement(Triangulation< dim, spacedim > &triangulation, const bool check_for_distorted_cells) override
std::unique_ptr< Policy< dim, spacedim > > clone() override
virtual std::unique_ptr< Policy< dim, spacedim > > clone()=0
virtual void update_neighbors(Triangulation< dim, spacedim > &tria)=0
virtual void prepare_refinement_dim_dependent(Triangulation< dim, spacedim > &triangulation)=0
virtual void prevent_distorted_boundary_cells(Triangulation< dim, spacedim > &triangulation)=0
virtual Triangulation< dim, spacedim >::DistortedCellList execute_refinement(Triangulation< dim, spacedim > &triangulation, const bool check_for_distorted_cells)=0
virtual bool coarsening_allowed(const typename Triangulation< dim, spacedim >::cell_iterator &cell)=0
virtual void delete_children(Triangulation< dim, spacedim > &triangulation, typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< unsigned int > &line_cell_count, std::vector< unsigned int > &quad_cell_count)=0
virtual ~Policy()=default
std::vector< bool > quads_line_orientations
std::vector< bool > quad_is_quadrilateral
std::vector< std::pair< int, int > > neighbors
std::vector< types::global_cell_index > global_active_cell_indices
std::vector< types::global_cell_index > global_level_cell_indices
std::vector< ReferenceCell > reference_cell
std::vector< int > parents
std::vector< types::subdomain_id > level_subdomain_ids
std::vector< bool > coarsen_flags
std::vector< std::uint8_t > refine_flags
std::vector< types::subdomain_id > subdomain_ids
std::vector< unsigned int > active_cell_indices
std::vector< bool > direction_flags
TriaObjectsOrientations face_orientations
void resize(const unsigned int n_objects)
std::vector< types::manifold_id > manifold_id
void clear_user_data(const unsigned int i)
std::vector< int > children
unsigned int next_free_pair
std::vector< bool > user_flags
bool reverse_order_next_free_single
unsigned int next_free_single
unsigned int n_objects() const
std::vector< BoundaryOrMaterialId > boundary_or_material_id
std::vector< std::uint8_t > refinement_cases
std::vector< UserData > user_data
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
IteratorRange< active_face_iterator > active_face_iterators() const
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcInternalErrorOnCell(int arg1)
static ::ExceptionBase & ExcIO()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
static ::ExceptionBase & ExcInteriorQuadCantBeBoundary(int arg1, int arg2, int arg3, int arg4, types::boundary_id arg5)
static ::ExceptionBase & ExcQuadInexistant(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInconsistentLineInfoOfLine(int arg1, int arg2, std::string arg3)
static ::ExceptionBase & ExcCellHasNegativeMeasure(int arg1)
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcMemoryInexact(int arg1, int arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcGridHasInvalidCell(int arg1)
static ::ExceptionBase & ExcLineInexistant(int arg1, int arg2)
static ::ExceptionBase & ExcMultiplySetLineInfoOfLine(int arg1, int arg2)
#define AssertNothrow(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInteriorLineCantBeBoundary(int arg1, int arg2, types::boundary_id arg3)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcInvalidVertexIndex(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DeclException5( Exception5, type1, type2, type3, type4, type5, outsequence)
#define AssertThrow(cond, exc)
static ::ExceptionBase & ExcInconsistentQuadInfoOfQuad(int arg1, int arg2, int arg3, int arg4, std::string arg5)
typename IteratorSelector::hex_iterator hex_iterator
typename IteratorSelector::active_quad_iterator active_quad_iterator
typename IteratorSelector::active_hex_iterator active_hex_iterator
typename IteratorSelector::quad_iterator quad_iterator
typename IteratorSelector::line_iterator line_iterator
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
typename IteratorSelector::active_line_iterator active_line_iterator
void set_all_manifold_ids_on_boundary(const types::manifold_id number)
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
virtual std::vector< types::manifold_id > get_manifold_ids() const
void reset_manifold(const types::manifold_id manifold_number)
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
void reset_all_manifolds()
void set_all_manifold_ids(const types::manifold_id number)
Task< RT > new_task(const std::function< RT()> &function)
#define AssertIsNotUsed(obj)
const unsigned int mn_tria_refine_flags_end
const unsigned int mn_tria_coarsen_flags_end
const unsigned int mn_tria_refine_flags_begin
const unsigned int mn_tria_hex_user_flags_end
const unsigned int mn_tria_line_user_flags_begin
const unsigned int mn_tria_line_user_flags_end
const unsigned int mn_tria_quad_user_flags_end
const unsigned int mn_tria_coarsen_flags_begin
const unsigned int mn_tria_hex_user_flags_begin
const unsigned int mn_tria_quad_user_flags_begin
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
@ valid
Iterator points to a valid object.
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Invalid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Line
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
const Manifold< dim, spacedim > & get_default_flat_manifold()
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void reserve_space(TriaFaces &tria_faces, const unsigned int new_quads_in_pairs, const unsigned int new_quads_single)
void monitor_memory(const TriaLevel &tria_level, const unsigned int true_dimension)
const types::boundary_id internal_face_boundary_id
const types::subdomain_id invalid_subdomain_id
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
const types::global_dof_index invalid_dof_index
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
void swap(SmartPointer< T, P > &t1, SmartPointer< T, Q > &t2)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static RefinementCase< dim - 1 > face_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static RefinementCase< dim > min_cell_refinement_case_for_face_refinement(const RefinementCase< dim - 1 > &face_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim - dim, spacedim >(&forms)[vertices_per_cell])
std::vector< CellData< 2 > > boundary_quads
bool check_consistency(const unsigned int dim) const
std::vector< CellData< 1 > > boundary_lines
std::vector< std::vector< CellData< dim > > > cell_infos
std::vector<::CellData< dim > > coarse_cells
std::vector< Point< spacedim > > coarse_cell_vertices
virtual ~DistortedCellList() noexcept override
std::list< typename Triangulation< dim, spacedim >::cell_iterator > distorted_cells
boost::signals2::signal< void(const Triangulation< dim, spacedim > &destination_tria)> copy
std::vector< std::size_t > ptr
static void update_neighbors(Triangulation< 1, spacedim > &)
static Triangulation< dim, spacedim >::DistortedCellList execute_refinement(Triangulation< dim, spacedim > &triangulation, const bool check_for_distorted_cells)
static bool coarsening_allowed(const typename Triangulation< dim, spacedim >::cell_iterator &)
static void update_neighbors(Triangulation< dim, spacedim > &triangulation)
static void delete_children(Triangulation< dim, spacedim > &triangulation, typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< unsigned int > &line_cell_count, std::vector< unsigned int > &quad_cell_count)
static void prepare_refinement_dim_dependent(Triangulation< dim, spacedim > &triangulation)
static void prevent_distorted_boundary_cells(Triangulation< dim, spacedim > &triangulation)
static void reserve_space_(TriaObjects &obj, const unsigned int size)
static void reserve_space_(TriaFaces &faces, const unsigned structdim, const unsigned int size)
static void compute_number_cache_dim(const Triangulation< dim, spacedim > &triangulation, const unsigned int level_objects, internal::TriangulationImplementation::NumberCache< 2 > &number_cache)
static void prevent_distorted_boundary_cells(Triangulation< 1, spacedim > &)
static void update_neighbors(Triangulation< dim, spacedim > &triangulation)
static void prepare_refinement_dim_dependent(const Triangulation< dim, spacedim > &)
static void delete_children(Triangulation< 3, spacedim > &triangulation, typename Triangulation< 3, spacedim >::cell_iterator &cell, std::vector< unsigned int > &line_cell_count, std::vector< unsigned int > &quad_cell_count)
static void reserve_space_(TriaLevel &level, const unsigned int spacedim, const unsigned int size, const bool orientation_needed)
static void update_neighbors(Triangulation< 1, spacedim > &)
static Triangulation< 3, spacedim >::DistortedCellList execute_refinement(Triangulation< 3, spacedim > &triangulation, const bool check_for_distorted_cells)
static Triangulation< dim, spacedim >::DistortedCellList execute_refinement_isotropic(Triangulation< dim, spacedim > &triangulation, const bool check_for_distorted_cells)
static void compute_number_cache(const Triangulation< dim, spacedim > &triangulation, const unsigned int level_objects, internal::TriangulationImplementation::NumberCache< dim > &number_cache)
static void create_children(Triangulation< 2, spacedim > &triangulation, unsigned int &next_unused_vertex, typename Triangulation< 2, spacedim >::raw_line_iterator &next_unused_line, typename Triangulation< 2, spacedim >::raw_cell_iterator &next_unused_cell, const typename Triangulation< 2, spacedim >::cell_iterator &cell)
static void prevent_distorted_boundary_cells(Triangulation< dim, spacedim > &triangulation)
static bool coarsening_allowed(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
static void compute_number_cache_dim(const Triangulation< dim, spacedim > &triangulation, const unsigned int level_objects, internal::TriangulationImplementation::NumberCache< 3 > &number_cache)
static void delete_children(Triangulation< 1, spacedim > &triangulation, typename Triangulation< 1, spacedim >::cell_iterator &cell, std::vector< unsigned int > &, std::vector< unsigned int > &)
static void prepare_refinement_dim_dependent(Triangulation< 3, spacedim > &triangulation)
static void delete_children(Triangulation< 2, spacedim > &triangulation, typename Triangulation< 2, spacedim >::cell_iterator &cell, std::vector< unsigned int > &line_cell_count, std::vector< unsigned int > &)
static void compute_number_cache_dim(const Triangulation< dim, spacedim > &triangulation, const unsigned int level_objects, internal::TriangulationImplementation::NumberCache< 1 > &number_cache)
static Triangulation< 1, spacedim >::DistortedCellList execute_refinement(Triangulation< 1, spacedim > &triangulation, const bool)
static Triangulation< 3, spacedim >::DistortedCellList execute_refinement_isotropic(Triangulation< 3, spacedim > &triangulation, const bool check_for_distorted_cells)
static void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata, Triangulation< dim, spacedim > &tria)
static Triangulation< 2, spacedim >::DistortedCellList execute_refinement(Triangulation< 2, spacedim > &triangulation, const bool check_for_distorted_cells)
static void process_subcelldata(const CRS< T > &crs, TriaObjects &obj, const std::vector< CellData< structdim > > &boundary_objects_in, const std::vector< Point< spacedim > > &vertex_locations)
std::vector< std::vector< CellData< dim > > > cell_infos
const ::Triangulation< dim, spacedim > & tria