50 namespace TriangulationImplementation
124 template <
int dim,
int spacedim>
126 cell_is_patch_level_1(
131 unsigned int n_active_children = 0;
132 for (
unsigned int i = 0; i < cell->n_children(); ++i)
133 if (cell->child(i)->is_active())
136 return (n_active_children == 0) ||
137 (n_active_children == cell->n_children());
147 template <
int dim,
int spacedim>
149 cell_will_be_coarsened(
155 if (cell->has_children())
157 unsigned int children_to_coarsen = 0;
158 const unsigned int n_children = cell->n_children();
160 for (
unsigned int c = 0; c < n_children; ++c)
161 if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
162 ++children_to_coarsen;
163 if (children_to_coarsen == n_children)
166 for (
unsigned int c = 0; c < n_children; ++c)
167 if (cell->child(c)->is_active())
168 cell->child(c)->clear_coarsen_flag();
204 template <
int dim,
int spacedim>
206 face_will_be_refined_by_neighbor_internal(
208 const unsigned int face_no,
217 cell->neighbor(face_no);
224 if (neighbor->has_children())
229 if (cell_will_be_coarsened(neighbor))
238 expected_face_ref_case = cell->face(face_no)->refinement_case();
250 const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
258 neighbor->face_orientation(neighbor_neighbor),
259 neighbor->face_flip(neighbor_neighbor),
260 neighbor->face_rotation(neighbor_neighbor));
264 neighbor_face = neighbor->face(neighbor_neighbor);
265 const int this_face_index = cell->face_index(face_no);
271 if (neighbor_face->index() == this_face_index)
277 expected_face_ref_case = face_ref_case;
295 for (
unsigned int c = 0; c < neighbor_face->n_children(); ++c)
296 if (neighbor_face->child_index(c) == this_face_index)
305 if ((neighbor_face->refinement_case() | face_ref_case) ==
306 neighbor_face->refinement_case())
325 expected_face_ref_case =
326 ~neighbor_face->refinement_case();
353 template <
int dim,
int spacedim>
355 face_will_be_refined_by_neighbor(
357 const unsigned int face_no)
360 return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
367 template <
int dim,
int spacedim>
369 face_will_be_refined_by_neighbor(
371 const unsigned int face_no,
374 return face_will_be_refined_by_neighbor_internal(cell,
376 expected_face_ref_case);
381 template <
int dim,
int spacedim>
383 satisfies_level1_at_vertex_rule(
386 std::vector<unsigned int> min_adjacent_cell_level(
388 std::vector<unsigned int> max_adjacent_cell_level(
391 for (
const auto &cell :
triangulation.active_cell_iterators())
392 for (
const unsigned int v : cell->vertex_indices())
394 min_adjacent_cell_level[cell->vertex_index(v)] =
395 std::min<unsigned int>(
396 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
397 max_adjacent_cell_level[cell->vertex_index(v)] =
398 std::max<unsigned int>(
399 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
402 for (
unsigned int k = 0; k <
triangulation.n_vertices(); ++k)
404 if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
417 template <
int dim,
int spacedim>
418 std::vector<unsigned int>
423 std::vector<unsigned int> line_cell_count(
triangulation.n_raw_lines(),
426 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
427 ++line_cell_count[cell->line_index(
l)];
428 return line_cell_count;
431 return std::vector<unsigned int>();
442 template <
int dim,
int spacedim>
443 std::vector<unsigned int>
448 std::vector<unsigned int> quad_cell_count(
triangulation.n_raw_quads(),
451 for (
unsigned int q : cell->face_indices())
452 ++quad_cell_count[cell->quad_index(q)];
453 return quad_cell_count;
483 for (
auto &cell : cells)
485 std::swap(cell.vertices[2], cell.vertices[3]);
490 reorder_compatibility(std::vector<
CellData<3>> &cells,
494 for (
auto &cell : cells)
498 tmp[i] = cell.vertices[i];
506 std::swap(boundary_quad.vertices[2], boundary_quad.vertices[3]);
528 template <
int dim,
int spacedim>
533 if (line->has_children())
534 return line->child(0)->vertex_index(1);
539 template <
int dim,
int spacedim>
544 switch (
static_cast<unsigned char>(quad->refinement_case()))
547 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
550 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
553 return quad->child(0)->vertex_index(3);
562 template <
int dim,
int spacedim>
567 switch (
static_cast<unsigned char>(hex->refinement_case()))
570 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
573 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
576 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
579 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
582 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
585 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
588 return hex->child(0)->vertex_index(7);
609 template <
class TRIANGULATION>
610 inline typename TRIANGULATION::DistortedCellList
611 collect_distorted_coarse_cells(
const TRIANGULATION &)
613 return typename TRIANGULATION::DistortedCellList();
631 for (
const auto &cell :
triangulation.cell_iterators_on_level(0))
641 if (determinants[i] <= 1
e-9 *
std::pow(cell->diameter(), 1. * dim))
643 distorted_cells.distorted_cells.push_back(cell);
648 return distorted_cells;
660 has_distorted_children(
665 for (
unsigned int c = 0; c < cell->n_children(); ++c)
669 vertices[i] = cell->child(c)->vertex(i);
675 if (determinants[i] <=
676 1
e-9 *
std::pow(cell->child(c)->diameter(), 1. * dim))
691 template <
int dim,
int spacedim>
693 has_distorted_children(
700 template <
int dim,
int spacedim>
702 update_periodic_face_map_recursively(
705 unsigned int n_face_1,
706 unsigned int n_face_2,
707 const std::bitset<3> & orientation,
713 std::bitset<3>>> &periodic_face_map)
716 const FaceIterator face_1 = cell_1->face(n_face_1);
717 const FaceIterator face_2 = cell_2->face(n_face_2);
719 const bool face_orientation = orientation[0];
720 const bool face_flip = orientation[1];
721 const bool face_rotation = orientation[2];
723 Assert((dim != 1) || (face_orientation ==
true && face_flip ==
false &&
724 face_rotation ==
false),
726 "(face_orientation, face_flip, face_rotation) "
727 "is invalid for 1D"));
729 Assert((dim != 2) || (face_orientation ==
true && face_rotation ==
false),
731 "(face_orientation, face_flip, face_rotation) "
732 "is invalid for 2D"));
736 Assert(face_1->at_boundary() && face_2->at_boundary(),
737 ExcMessage(
"Periodic faces must be on the boundary"));
743 Assert(dim == 1 ||
std::abs(cell_1->level() - cell_2->level()) < 2,
748 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
750 const CellFace cell_face_1(cell_1, n_face_1);
751 const CellFace cell_face_2(cell_2, n_face_2);
752 const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
753 cell_face_2, orientation);
755 const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
756 periodic_faces(cell_face_1, cell_face_orientation_2);
760 periodic_face_map.insert(periodic_faces);
764 if (cell_1->has_children())
766 if (cell_2->has_children())
768 update_periodic_face_map_recursively<dim, spacedim>(
769 cell_1->child(n_face_1),
770 cell_2->child(n_face_2),
778 update_periodic_face_map_recursively<dim, spacedim>(
779 cell_1->child(n_face_1),
794 static const int lookup_table_2d[2][2] =
801 static const int lookup_table_3d[2][2][2][4] =
820 if (cell_1->has_children())
822 if (cell_2->has_children())
828 Assert(face_1->n_children() ==
834 for (
unsigned int i = 0;
835 i < GeometryInfo<dim>::max_children_per_face;
843 j = lookup_table_2d[face_flip][i];
846 j = lookup_table_3d[face_orientation][face_flip]
854 unsigned int child_cell_1 =
856 cell_1->refinement_case(),
859 cell_1->face_orientation(n_face_1),
860 cell_1->face_flip(n_face_1),
861 cell_1->face_rotation(n_face_1),
862 face_1->refinement_case());
863 unsigned int child_cell_2 =
865 cell_2->refinement_case(),
868 cell_2->face_orientation(n_face_2),
869 cell_2->face_flip(n_face_2),
870 cell_2->face_rotation(n_face_2),
871 face_2->refinement_case());
873 Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
876 Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
882 update_periodic_face_map_recursively<dim, spacedim>(
883 cell_1->child(child_cell_1),
884 cell_2->child(child_cell_2),
893 for (
unsigned int i = 0;
894 i < GeometryInfo<dim>::max_children_per_face;
898 unsigned int child_cell_1 =
900 cell_1->refinement_case(),
903 cell_1->face_orientation(n_face_1),
904 cell_1->face_flip(n_face_1),
905 cell_1->face_rotation(n_face_1),
906 face_1->refinement_case());
909 update_periodic_face_map_recursively<dim, spacedim>(
910 cell_1->child(child_cell_1),
928 namespace TriangulationImplementation
936 using ::Triangulation;
944 <<
"Something went wrong when making cell " << arg1
945 <<
". Read the docs and the source code "
946 <<
"for more information.");
953 <<
"Something went wrong upon construction of cell "
967 <<
" has negative measure. This typically "
968 <<
"indicates some distortion in the cell, or a mistakenly "
969 <<
"swapped pair of vertices in the input to "
970 <<
"Triangulation::create_triangulation().");
982 <<
"Error while creating cell " << arg1
983 <<
": the vertex index " << arg2 <<
" must be between 0 and "
992 <<
"While trying to assign a boundary indicator to a line: "
993 <<
"the line with end vertices " << arg1 <<
" and " << arg2
994 <<
" does not exist.");
1004 <<
"While trying to assign a boundary indicator to a quad: "
1005 <<
"the quad with bounding lines " << arg1 <<
", " << arg2
1006 <<
", " << arg3 <<
", " << arg4 <<
" does not exist.");
1016 <<
"The input data for creating a triangulation contained "
1017 <<
"information about a line with indices " << arg1 <<
" and " << arg2
1018 <<
" that is described to have boundary indicator "
1019 <<
static_cast<int>(arg3)
1020 <<
". However, this is an internal line not located on the "
1021 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1023 <<
"If this happened at a place where you call "
1024 <<
"Triangulation::create_triangulation() yourself, you need "
1025 <<
"to check the SubCellData object you pass to this function."
1028 <<
"If this happened in a place where you are reading a mesh "
1029 <<
"from a file, then you need to investigate why such a line "
1030 <<
"ended up in the input file. A typical case is a geometry "
1031 <<
"that consisted of multiple parts and for which the mesh "
1032 <<
"generator program assumes that the interface between "
1033 <<
"two parts is a boundary when that isn't supposed to be "
1034 <<
"the case, or where the mesh generator simply assigns "
1035 <<
"'geometry indicators' to lines at the perimeter of "
1036 <<
"a part that are not supposed to be interpreted as "
1037 <<
"'boundary indicators'.");
1049 <<
"The input data for creating a triangulation contained "
1050 <<
"information about a quad with indices " << arg1 <<
", " << arg2
1051 <<
", " << arg3 <<
", and " << arg4
1052 <<
" that is described to have boundary indicator "
1053 <<
static_cast<int>(arg5)
1054 <<
". However, this is an internal quad not located on the "
1055 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1057 <<
"If this happened at a place where you call "
1058 <<
"Triangulation::create_triangulation() yourself, you need "
1059 <<
"to check the SubCellData object you pass to this function."
1062 <<
"If this happened in a place where you are reading a mesh "
1063 <<
"from a file, then you need to investigate why such a quad "
1064 <<
"ended up in the input file. A typical case is a geometry "
1065 <<
"that consisted of multiple parts and for which the mesh "
1066 <<
"generator program assumes that the interface between "
1067 <<
"two parts is a boundary when that isn't supposed to be "
1068 <<
"the case, or where the mesh generator simply assigns "
1069 <<
"'geometry indicators' to quads at the surface of "
1070 <<
"a part that are not supposed to be interpreted as "
1071 <<
"'boundary indicators'.");
1080 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1081 <<
" and " << arg2 <<
" appears more than once. "
1082 <<
"This is not allowed.");
1092 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1093 <<
" and " << arg2 <<
" appears multiple times with different (valid) "
1094 << arg3 <<
". This is not allowed.");
1106 <<
"In SubCellData the quad info of the quad with line indices " << arg1
1107 <<
", " << arg2 <<
", " << arg3 <<
" and " << arg4
1108 <<
" appears multiple times with different (valid) " << arg5
1109 <<
". This is not allowed.");
1120 const unsigned int new_quads_in_pairs,
1121 const unsigned int new_quads_single)
1127 unsigned int next_free_single = 0;
1128 unsigned int next_free_pair = 0;
1132 unsigned int n_quads = 0;
1133 unsigned int n_unused_pairs = 0;
1134 unsigned int n_unused_singles = 0;
1135 for (
unsigned int i = 0; i < tria_faces.
quads.
used.size(); ++i)
1139 else if (i + 1 < tria_faces.
quads.
used.size())
1144 if (next_free_single == 0)
1145 next_free_single = i;
1150 if (next_free_pair == 0)
1158 Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
1163 const int additional_single_quads = new_quads_single - n_unused_singles;
1165 unsigned int new_size =
1166 tria_faces.
quads.
used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
1167 if (additional_single_quads > 0)
1168 new_size += additional_single_quads;
1207 const unsigned int total_cells,
1208 const unsigned int dimension,
1209 const unsigned int space_dimension)
1261 if (dimension < space_dimension)
1272 tria_level.
parents.reserve((total_cells + 1) / 2);
1274 (total_cells + 1) / 2 -
1278 tria_level.
neighbors.reserve(total_cells * (2 * dimension));
1280 total_cells * (2 * dimension) -
1282 std::make_pair(-1, -1));
1284 if (tria_level.
dim == 2 || tria_level.
dim == 3)
1286 const unsigned int max_faces_per_cell = 2 * dimension;
1288 max_faces_per_cell);
1291 total_cells * max_faces_per_cell -
1313 <<
"The containers have sizes " << arg1 <<
" and " << arg2
1314 <<
", which is not as expected.");
1323 const unsigned int true_dimension)
1326 (void)true_dimension;
1353 const unsigned int new_objects_in_pairs,
1354 const unsigned int new_objects_single = 0)
1366 unsigned int n_objects = 0;
1367 unsigned int n_unused_pairs = 0;
1368 unsigned int n_unused_singles = 0;
1369 for (
unsigned int i = 0; i < tria_objects.
used.size(); ++i)
1371 if (tria_objects.
used[i])
1373 else if (i + 1 < tria_objects.
used.size())
1375 if (tria_objects.
used[i + 1])
1392 Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
1393 tria_objects.
used.size(),
1398 const int additional_single_objects =
1399 new_objects_single - n_unused_singles;
1401 unsigned int new_size = tria_objects.
used.size() +
1402 new_objects_in_pairs - 2 * n_unused_pairs;
1403 if (additional_single_objects > 0)
1404 new_size += additional_single_objects;
1407 if (new_size > tria_objects.
n_objects())
1409 const unsigned int max_faces_per_cell =
1411 const unsigned int max_children_per_cell =
1414 tria_objects.
cells.reserve(new_size * max_faces_per_cell);
1415 tria_objects.
cells.insert(tria_objects.
cells.end(),
1420 tria_objects.
used.reserve(new_size);
1421 tria_objects.
used.insert(tria_objects.
used.end(),
1422 new_size - tria_objects.
used.size(),
1431 const unsigned int factor = max_children_per_cell / 2;
1432 tria_objects.
children.reserve(factor * new_size);
1452 tria_objects.
user_data.reserve(new_size);
1453 tria_objects.
user_data.resize(new_size);
1462 if (n_unused_singles == 0)
1470 const unsigned int new_hexes = new_objects_in_pairs;
1472 const unsigned int new_size =
1473 new_hexes + std::count(tria_objects.
used.begin(),
1474 tria_objects.
used.end(),
1478 if (new_size > tria_objects.
n_objects())
1480 const unsigned int max_faces_per_cell =
1483 tria_objects.
cells.reserve(new_size * max_faces_per_cell);
1484 tria_objects.
cells.insert(tria_objects.
cells.end(),
1489 tria_objects.
used.reserve(new_size);
1490 tria_objects.
used.insert(tria_objects.
used.end(),
1491 new_size - tria_objects.
used.size(),
1500 tria_objects.
children.reserve(4 * new_size);
1519 tria_objects.
user_data.reserve(new_size);
1520 tria_objects.
user_data.resize(new_size);
1544 tria_object.
used.size()));
1585 template <
int dim,
int spacedim>
1607 std::vector<unsigned int> & line_cell_count,
1608 std::vector<unsigned int> &quad_cell_count) = 0;
1615 const bool check_for_distorted_cells) = 0;
1644 virtual std::unique_ptr<Policy<dim, spacedim>>
1655 template <
int dim,
int spacedim,
typename T>
1662 T::update_neighbors(
tria);
1669 std::vector<unsigned int> & line_cell_count,
1670 std::vector<unsigned int> &quad_cell_count)
override
1672 T::delete_children(
tria, cell, line_cell_count, quad_cell_count);
1677 const bool check_for_distorted_cells)
override
1679 return T::execute_refinement(
triangulation, check_for_distorted_cells);
1701 return T::template coarsening_allowed<dim, spacedim>(cell);
1704 std::unique_ptr<Policy<dim, spacedim>>
1707 return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
1821 template <
int dim,
int spacedim>
1825 const unsigned int level_objects,
1828 using line_iterator =
1832 if (level_objects > 0)
1865 for (; line != endc; ++line)
1868 if (line->has_children() ==
false)
1886 for (; line != endc; ++line)
1889 if (line->has_children() ==
false)
1909 template <
int dim,
int spacedim>
1913 const unsigned int level_objects,
1924 &compute_number_cache<dim, spacedim>),
1930 using quad_iterator =
1947 unsigned int n_levels = 0;
1948 if (level_objects > 0)
1952 n_levels =
level + 1;
1965 (
level == n_levels - 1 ?
1968 for (; quad != endc; ++quad)
1971 if (quad->has_children() ==
false)
1989 for (; quad != endc; ++quad)
1992 if (quad->has_children() ==
false)
1998 update_lines.
join();
2016 template <
int dim,
int spacedim>
2020 const unsigned int level_objects,
2031 &compute_number_cache<dim, spacedim>),
2037 using hex_iterator =
2055 unsigned int n_levels = 0;
2056 if (level_objects > 0)
2060 n_levels =
level + 1;
2072 endc = (
level == n_levels - 1 ?
2075 for (; hex != endc; ++hex)
2078 if (hex->has_children() ==
false)
2096 for (; hex != endc; ++hex)
2099 if (hex->has_children() ==
false)
2105 update_quads_and_lines.
join();
2110 template <
int spacedim>
2116 template <
int dim,
int spacedim>
2171 static const unsigned int left_right_offset[2][6][2] = {
2180 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
2191 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
2195 for (
auto f : cell->face_indices())
2200 const unsigned int offset =
2201 (cell->direction_flag() ?
2202 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2204 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2216 if (cell->is_active() && face->has_children())
2246 if (face->has_children() &&
2247 (cell->is_active() ||
2249 cell->refinement_case(), f) ==
2252 for (
unsigned int c = 0; c < face->n_children(); ++c)
2255 if (face->child(0)->has_children())
2262 if (face->child(1)->has_children())
2279 for (
auto f : cell->face_indices())
2281 const unsigned int offset =
2282 (cell->direction_flag() ?
2283 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2285 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2295 template <
int dim,
int spacedim>
2309 if (dim == spacedim)
2310 for (
unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2318 const double cell_measure = GridTools::cell_measure<spacedim>(
2341 const auto connectivity = build_connectivity<unsigned int>(cells);
2342 const unsigned int n_cell = cells.size();
2350 const auto & crs = connectivity.entity_to_entities(1, 0);
2351 const unsigned int n_lines = crs.ptr.size() - 1;
2357 for (
unsigned int line = 0; line < n_lines; ++line)
2358 for (
unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
2371 const auto & crs = connectivity.entity_to_entities(2, 1);
2372 const unsigned int n_quads = crs.ptr.size() - 1;
2379 for (
unsigned int q = 0, k = 0; q < n_quads; ++q)
2382 faces.quad_reference_cell[q] = connectivity.entity_types(2)[q];
2385 for (
unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
2393 faces.quads_line_orientations
2395 connectivity.entity_orientations(1)[k];
2406 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2407 const auto &nei = connectivity.entity_to_entities(dim, dim);
2411 const bool orientation_needed =
2414 std::any_of(connectivity.entity_orientations(1).begin(),
2415 connectivity.entity_orientations(1).end(),
2416 [](
const auto &i) { return i == 0; }));
2423 for (
unsigned int cell = 0; cell < n_cell; ++cell)
2426 cells_0.boundary_or_material_id[cell].material_id =
2427 cells[cell].material_id;
2430 cells_0.manifold_id[cell] = cells[cell].manifold_id;
2433 level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
2436 for (
unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2440 if (nei.col[i] !=
static_cast<unsigned int>(-1))
2442 j] = {0, nei.col[i]};
2449 if (orientation_needed)
2451 level.face_orientations
2453 connectivity.entity_orientations(dim - 1)[i];
2462 auto &bids_face = dim == 3 ?
2463 tria.
faces->quads.boundary_or_material_id :
2464 tria.
faces->lines.boundary_or_material_id;
2467 std::vector<unsigned int> count(bids_face.size(), 0);
2470 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2473 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
2474 for (
unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
2475 count[crs.col[i]]++;
2478 for (
unsigned int face = 0; face < count.size(); ++face)
2480 if (count[face] != 1)
2484 bids_face[face].boundary_id = 0;
2490 const auto &crs = connectivity.entity_to_entities(2, 1);
2491 for (
unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
2492 tria.
faces->lines.boundary_or_material_id[crs.col[i]]
2498 static const unsigned int t_tba =
static_cast<unsigned int>(-1);
2499 static const unsigned int t_inner =
static_cast<unsigned int>(-2);
2501 std::vector<unsigned int> type(
vertices.size(), t_tba);
2503 const auto &crs = connectivity.entity_to_entities(1, 0);
2505 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
2506 for (
unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2508 if (type[crs.col[i]] != t_inner)
2509 type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
2511 for (
unsigned int face = 0; face < type.size(); ++face)
2516 if (type[face] != t_inner && type[face] != t_tba)
2537 template <
int structdim,
int spacedim,
typename T>
2547 if (boundary_objects_in.size() == 0)
2551 auto boundary_objects = boundary_objects_in;
2554 for (
auto &boundary_object : boundary_objects)
2555 std::sort(boundary_object.vertices.begin(),
2556 boundary_object.vertices.end());
2559 std::sort(boundary_objects.begin(),
2560 boundary_objects.end(),
2561 [](
const auto &a,
const auto &
b) {
2562 return a.vertices < b.vertices;
2565 unsigned int counter = 0;
2567 std::vector<unsigned int> key;
2570 for (
unsigned int o = 0; o < obj.
n_objects(); ++o)
2584 key.assign(crs.
col.data() + crs.
ptr[o],
2585 crs.
col.data() + crs.
ptr[o + 1]);
2586 std::sort(key.begin(), key.end());
2589 const auto subcell_object =
2591 boundary_objects.end(),
2593 [&](
const auto &cell,
const auto &key) {
2594 return cell.vertices < key;
2598 if (subcell_object == boundary_objects.end() ||
2599 subcell_object->vertices != key)
2608 if (subcell_object->boundary_id !=
2611 (void)vertex_locations;
2615 "The input arguments for creating a triangulation "
2616 "specified a boundary id for an internal face. This "
2619 "The object in question has vertex indices " +
2620 [subcell_object]() {
2622 for (
const auto v : subcell_object->vertices)
2626 " which are located at positions " +
2627 [vertex_locations, subcell_object]() {
2628 std::ostringstream s;
2629 for (
const auto v : subcell_object->vertices)
2630 s <<
'(' << vertex_locations[v] <<
')';
2647 const unsigned structdim,
2648 const unsigned int size)
2650 const unsigned int dim = faces.
dim;
2652 const unsigned int max_faces_per_cell = 2 * structdim;
2654 if (dim == 3 && structdim == 2)
2669 const unsigned int spacedim,
2670 const unsigned int size,
2671 const bool orientation_needed)
2673 const unsigned int dim =
level.dim;
2675 const unsigned int max_faces_per_cell = 2 * dim;
2677 level.active_cell_indices.assign(size, -1);
2678 level.subdomain_ids.assign(size, 0);
2679 level.level_subdomain_ids.assign(size, 0);
2681 level.refine_flags.assign(size, 0u);
2682 level.coarsen_flags.assign(size,
false);
2684 level.parents.assign((size + 1) / 2, -1);
2687 level.direction_flags.assign(size,
true);
2689 level.neighbors.assign(size * max_faces_per_cell, {-1, -1});
2693 if (orientation_needed)
2694 level.face_orientations.assign(size * max_faces_per_cell, -1);
2696 level.global_active_cell_indices.assign(size,
2698 level.global_level_cell_indices.assign(size,
2707 const unsigned int structdim = obj.
structdim;
2709 const unsigned int max_children_per_cell = 1 << structdim;
2710 const unsigned int max_faces_per_cell = 2 * structdim;
2712 obj.
used.assign(size,
true);
2716 BoundaryOrMaterialId());
2724 obj.
children.assign(max_children_per_cell / 2 * size, -1);
2726 obj.
cells.assign(max_faces_per_cell * size, -1);
2756 template <
int spacedim>
2760 std::vector<unsigned int> &,
2761 std::vector<unsigned int> &)
2763 const unsigned int dim = 1;
2775 Assert(!cell->child(0)->has_children() &&
2776 !cell->child(1)->has_children(),
2782 if (cell->neighbor(0)->has_children())
2789 neighbor = neighbor->child(1);
2792 Assert(neighbor->neighbor(1) == cell->child(0),
2794 neighbor->set_neighbor(1, cell);
2800 if (neighbor->has_children())
2801 neighbor = neighbor->child(1);
2810 if (cell->neighbor(1)->has_children())
2817 neighbor = neighbor->child(0);
2820 Assert(neighbor->neighbor(0) == cell->child(1),
2822 neighbor->set_neighbor(0, cell);
2828 if (neighbor->has_children())
2829 neighbor = neighbor->child(0);
2839 triangulation.vertices_used[cell->child(0)->vertex_index(1)] =
false;
2845 for (
unsigned int child = 0; child < cell->n_children(); ++child)
2848 cell->child(child)->clear_user_flag();
2849 cell->child(child)->clear_used_flag();
2854 cell->clear_children();
2855 cell->clear_user_flag();
2860 template <
int spacedim>
2864 std::vector<unsigned int> &line_cell_count,
2865 std::vector<unsigned int> &)
2867 const unsigned int dim = 2;
2875 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2878 lines_to_delete.reserve(4 * 2 + 4);
2883 for (
unsigned int c = 0; c < cell->n_children(); ++c)
2887 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
2888 --line_cell_count[child->line_index(
l)];
2903 .vertices_used[cell->child(0)->line(1)->vertex_index(1)] =
false;
2905 lines_to_delete.push_back(cell->child(0)->line(1));
2906 lines_to_delete.push_back(cell->child(0)->line(3));
2907 lines_to_delete.push_back(cell->child(3)->line(0));
2908 lines_to_delete.push_back(cell->child(3)->line(2));
2912 unsigned int inner_face_no =
2917 lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2921 for (
unsigned int child = 0; child < cell->n_children(); ++child)
2923 cell->child(child)->clear_user_data();
2924 cell->child(child)->clear_user_flag();
2925 cell->child(child)->clear_used_flag();
2930 cell->clear_children();
2931 cell->clear_refinement_case();
2932 cell->clear_user_flag();
2938 for (
unsigned int line_no = 0;
2939 line_no < GeometryInfo<dim>::lines_per_cell;
2943 cell->line(line_no);
2945 if (line->has_children())
2950 Assert((line_cell_count[line->child_index(0)] == 0 &&
2951 line_cell_count[line->child_index(1)] == 0) ||
2952 (line_cell_count[line->child_index(0)] > 0 &&
2953 line_cell_count[line->child_index(1)] > 0),
2956 if (line_cell_count[line->child_index(0)] == 0)
2958 for (
unsigned int c = 0; c < 2; ++c)
2959 Assert(!line->child(c)->has_children(),
2967 .vertices_used[line->child(0)->vertex_index(1)] =
false;
2969 lines_to_delete.push_back(line->child(0));
2970 lines_to_delete.push_back(line->child(1));
2972 line->clear_children();
2984 typename std::vector<
2986 line = lines_to_delete.
begin(),
2987 endline = lines_to_delete.end();
2988 for (; line != endline; ++line)
2990 (*line)->clear_user_data();
2991 (*line)->clear_user_flag();
2992 (*line)->clear_used_flag();
2998 template <
int spacedim>
3002 std::vector<unsigned int> &line_cell_count,
3003 std::vector<unsigned int> &quad_cell_count)
3005 const unsigned int dim = 3;
3017 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3019 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3022 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3023 quads_to_delete.reserve(6 * 4 + 12);
3027 for (
unsigned int c = 0; c < cell->n_children(); ++c)
3031 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
3032 --line_cell_count[child->line_index(
l)];
3034 --quad_cell_count[child->quad_index(f)];
3048 quads_to_delete.push_back(cell->child(0)->face(1));
3051 quads_to_delete.push_back(cell->child(0)->face(3));
3054 quads_to_delete.push_back(cell->child(0)->face(5));
3057 quads_to_delete.push_back(cell->child(0)->face(1));
3058 quads_to_delete.push_back(cell->child(0)->face(3));
3059 quads_to_delete.push_back(cell->child(3)->face(0));
3060 quads_to_delete.push_back(cell->child(3)->face(2));
3062 lines_to_delete.push_back(cell->child(0)->line(11));
3065 quads_to_delete.push_back(cell->child(0)->face(1));
3066 quads_to_delete.push_back(cell->child(0)->face(5));
3067 quads_to_delete.push_back(cell->child(3)->face(0));
3068 quads_to_delete.push_back(cell->child(3)->face(4));
3070 lines_to_delete.push_back(cell->child(0)->line(5));
3073 quads_to_delete.push_back(cell->child(0)->face(3));
3074 quads_to_delete.push_back(cell->child(0)->face(5));
3075 quads_to_delete.push_back(cell->child(3)->face(2));
3076 quads_to_delete.push_back(cell->child(3)->face(4));
3078 lines_to_delete.push_back(cell->child(0)->line(7));
3081 quads_to_delete.push_back(cell->child(0)->face(1));
3082 quads_to_delete.push_back(cell->child(2)->face(1));
3083 quads_to_delete.push_back(cell->child(4)->face(1));
3084 quads_to_delete.push_back(cell->child(6)->face(1));
3086 quads_to_delete.push_back(cell->child(0)->face(3));
3087 quads_to_delete.push_back(cell->child(1)->face(3));
3088 quads_to_delete.push_back(cell->child(4)->face(3));
3089 quads_to_delete.push_back(cell->child(5)->face(3));
3091 quads_to_delete.push_back(cell->child(0)->face(5));
3092 quads_to_delete.push_back(cell->child(1)->face(5));
3093 quads_to_delete.push_back(cell->child(2)->face(5));
3094 quads_to_delete.push_back(cell->child(3)->face(5));
3096 lines_to_delete.push_back(cell->child(0)->line(5));
3097 lines_to_delete.push_back(cell->child(0)->line(7));
3098 lines_to_delete.push_back(cell->child(0)->line(11));
3099 lines_to_delete.push_back(cell->child(7)->line(0));
3100 lines_to_delete.push_back(cell->child(7)->line(2));
3101 lines_to_delete.push_back(cell->child(7)->line(8));
3107 triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3119 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3121 cell->child(child)->clear_user_data();
3122 cell->child(child)->clear_user_flag();
3129 cell->child(child)->set_face_orientation(f,
true);
3130 cell->child(child)->set_face_flip(f,
false);
3131 cell->child(child)->set_face_rotation(f,
false);
3134 cell->child(child)->clear_used_flag();
3139 cell->clear_children();
3140 cell->clear_refinement_case();
3141 cell->clear_user_flag();
3154 cell->face(quad_no);
3158 quad->has_children()) ||
3163 switch (quad->refinement_case())
3175 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3176 quad_cell_count[quad->child_index(1)] == 0) ||
3177 (quad_cell_count[quad->child_index(0)] > 0 &&
3178 quad_cell_count[quad->child_index(1)] > 0),
3184 unsigned int deleted_grandchildren = 0;
3185 unsigned int number_of_child_refinements = 0;
3187 for (
unsigned int c = 0; c < 2; ++c)
3188 if (quad->child(c)->has_children())
3190 ++number_of_child_refinements;
3195 (quad_cell_count[quad->child(c)->child_index(0)] ==
3197 quad_cell_count[quad->child(c)->child_index(1)] ==
3199 (quad_cell_count[quad->child(c)->child_index(0)] >
3201 quad_cell_count[quad->child(c)->child_index(1)] >
3204 if (quad_cell_count[quad->child(c)->child_index(0)] ==
3211 Assert(quad->refinement_case() +
3212 quad->child(c)->refinement_case() ==
3220 quads_to_delete.push_back(
3221 quad->child(c)->child(0));
3222 quads_to_delete.push_back(
3223 quad->child(c)->child(1));
3224 if (quad->child(c)->refinement_case() ==
3226 lines_to_delete.push_back(
3227 quad->child(c)->child(0)->line(1));
3229 lines_to_delete.push_back(
3230 quad->child(c)->child(0)->line(3));
3231 quad->child(c)->clear_children();
3232 quad->child(c)->clear_refinement_case();
3233 ++deleted_grandchildren;
3241 if (number_of_child_refinements > 0 &&
3242 deleted_grandchildren == number_of_child_refinements)
3247 middle_line = quad->child(0)->line(1);
3249 middle_line = quad->child(0)->line(3);
3251 lines_to_delete.push_back(middle_line->child(0));
3252 lines_to_delete.push_back(middle_line->child(1));
3254 .vertices_used[middle_vertex_index<dim, spacedim>(
3255 middle_line)] =
false;
3256 middle_line->clear_children();
3261 if (quad_cell_count[quad->child_index(0)] == 0)
3267 quads_to_delete.push_back(quad->child(0));
3268 quads_to_delete.push_back(quad->child(1));
3270 lines_to_delete.push_back(quad->child(0)->line(1));
3272 lines_to_delete.push_back(quad->child(0)->line(3));
3292 if (quad->child(0)->has_children())
3293 if (quad->refinement_case() ==
3327 quad_iterator switch_1 =
3328 quad->child(0)->child(1),
3330 quad->child(1)->child(0);
3332 Assert(!switch_1->has_children(),
3334 Assert(!switch_2->has_children(),
3337 const int switch_1_index = switch_1->index();
3338 const int switch_2_index = switch_2->index();
3339 for (
unsigned int l = 0;
3342 for (
unsigned int h = 0;
3346 for (
const unsigned int q :
3351 ->cells.get_bounding_object_indices(
3353 if (
index == switch_1_index)
3355 ->cells.get_bounding_object_indices(
3356 h)[q] = switch_2_index;
3357 else if (
index == switch_2_index)
3359 ->cells.get_bounding_object_indices(
3360 h)[q] = switch_1_index;
3365 const int switch_1_lines[4] = {
3366 static_cast<signed int>(
3367 switch_1->line_index(0)),
3368 static_cast<signed int>(
3369 switch_1->line_index(1)),
3370 static_cast<signed int>(
3371 switch_1->line_index(2)),
3372 static_cast<signed int>(
3373 switch_1->line_index(3))};
3374 const bool switch_1_line_orientations[4] = {
3375 switch_1->line_orientation(0),
3376 switch_1->line_orientation(1),
3377 switch_1->line_orientation(2),
3378 switch_1->line_orientation(3)};
3380 switch_1->boundary_id();
3381 const unsigned int switch_1_user_index =
3382 switch_1->user_index();
3383 const bool switch_1_user_flag =
3384 switch_1->user_flag_set();
3386 switch_1->set_bounding_object_indices(
3387 {switch_2->line_index(0),
3388 switch_2->line_index(1),
3389 switch_2->line_index(2),
3390 switch_2->line_index(3)});
3391 switch_1->set_line_orientation(
3392 0, switch_2->line_orientation(0));
3393 switch_1->set_line_orientation(
3394 1, switch_2->line_orientation(1));
3395 switch_1->set_line_orientation(
3396 2, switch_2->line_orientation(2));
3397 switch_1->set_line_orientation(
3398 3, switch_2->line_orientation(3));
3399 switch_1->set_boundary_id_internal(
3400 switch_2->boundary_id());
3401 switch_1->set_manifold_id(
3402 switch_2->manifold_id());
3403 switch_1->set_user_index(switch_2->user_index());
3404 if (switch_2->user_flag_set())
3405 switch_1->set_user_flag();
3407 switch_1->clear_user_flag();
3409 switch_2->set_bounding_object_indices(
3413 switch_1_lines[3]});
3414 switch_2->set_line_orientation(
3415 0, switch_1_line_orientations[0]);
3416 switch_2->set_line_orientation(
3417 1, switch_1_line_orientations[1]);
3418 switch_2->set_line_orientation(
3419 2, switch_1_line_orientations[2]);
3420 switch_2->set_line_orientation(
3421 3, switch_1_line_orientations[3]);
3422 switch_2->set_boundary_id_internal(
3423 switch_1_boundary_id);
3424 switch_2->set_manifold_id(
3425 switch_1->manifold_id());
3426 switch_2->set_user_index(switch_1_user_index);
3427 if (switch_1_user_flag)
3428 switch_2->set_user_flag();
3430 switch_2->clear_user_flag();
3432 const unsigned int child_0 =
3433 quad->child(0)->child_index(0);
3434 const unsigned int child_2 =
3435 quad->child(1)->child_index(0);
3436 quad->clear_children();
3437 quad->clear_refinement_case();
3438 quad->set_refinement_case(
3440 quad->set_children(0, child_0);
3441 quad->set_children(2, child_2);
3443 quad_cell_count[child_2]);
3458 const unsigned int child_0 =
3459 quad->child(0)->child_index(0);
3460 const unsigned int child_2 =
3461 quad->child(1)->child_index(0);
3462 quad->clear_children();
3463 quad->clear_refinement_case();
3464 quad->set_refinement_case(
3466 quad->set_children(0, child_0);
3467 quad->set_children(2, child_2);
3471 quad->clear_children();
3472 quad->clear_refinement_case();
3483 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3484 quad_cell_count[quad->child_index(1)] == 0 &&
3485 quad_cell_count[quad->child_index(2)] == 0 &&
3486 quad_cell_count[quad->child_index(3)] == 0) ||
3487 (quad_cell_count[quad->child_index(0)] > 0 &&
3488 quad_cell_count[quad->child_index(1)] > 0 &&
3489 quad_cell_count[quad->child_index(2)] > 0 &&
3490 quad_cell_count[quad->child_index(3)] > 0),
3493 if (quad_cell_count[quad->child_index(0)] == 0)
3499 lines_to_delete.push_back(quad->child(0)->line(1));
3500 lines_to_delete.push_back(quad->child(3)->line(0));
3501 lines_to_delete.push_back(quad->child(0)->line(3));
3502 lines_to_delete.push_back(quad->child(3)->line(2));
3504 for (
unsigned int child = 0; child < quad->n_children();
3506 quads_to_delete.push_back(quad->child(child));
3509 .vertices_used[quad->child(0)->vertex_index(3)] =
3512 quad->clear_children();
3513 quad->clear_refinement_case();
3532 for (
unsigned int line_no = 0;
3533 line_no < GeometryInfo<dim>::lines_per_cell;
3537 cell->line(line_no);
3541 line->has_children()) ||
3546 if (line->has_children())
3551 Assert((line_cell_count[line->child_index(0)] == 0 &&
3552 line_cell_count[line->child_index(1)] == 0) ||
3553 (line_cell_count[line->child_index(0)] > 0 &&
3554 line_cell_count[line->child_index(1)] > 0),
3557 if (line_cell_count[line->child_index(0)] == 0)
3559 for (
unsigned int c = 0; c < 2; ++c)
3560 Assert(!line->child(c)->has_children(),
3568 .vertices_used[line->child(0)->vertex_index(1)] =
false;
3570 lines_to_delete.push_back(line->child(0));
3571 lines_to_delete.push_back(line->child(1));
3573 line->clear_children();
3585 typename std::vector<
3587 line = lines_to_delete.
begin(),
3588 endline = lines_to_delete.end();
3589 for (; line != endline; ++line)
3591 (*line)->clear_user_data();
3592 (*line)->clear_user_flag();
3593 (*line)->clear_used_flag();
3596 typename std::vector<
3598 quad = quads_to_delete.
begin(),
3599 endquad = quads_to_delete.end();
3600 for (; quad != endquad; ++quad)
3602 (*quad)->clear_user_data();
3603 (*quad)->clear_children();
3604 (*quad)->clear_refinement_case();
3605 (*quad)->clear_user_flag();
3606 (*quad)->clear_used_flag();
3628 template <
int spacedim>
3632 unsigned int & next_unused_vertex,
3639 const unsigned int dim = 2;
3642 cell->clear_refine_flag();
3735 int new_vertices[9];
3736 for (
unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3737 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3738 for (
unsigned int line_no = 0; line_no < 4; ++line_no)
3739 if (cell->line(line_no)->has_children())
3740 new_vertices[4 + line_no] =
3741 cell->line(line_no)->child(0)->vertex_index(1);
3750 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
3751 ++next_unused_vertex;
3754 "Internal error: During refinement, the triangulation "
3755 "wants to access an element of the 'vertices' array "
3756 "but it turns out that the array is not large enough."));
3759 new_vertices[8] = next_unused_vertex;
3765 cell->center(
true,
true);
3771 unsigned int lmin = 8;
3772 unsigned int lmax = 12;
3779 for (
unsigned int l = lmin;
l < lmax; ++
l)
3781 while (next_unused_line->used() ==
true)
3783 new_lines[
l] = next_unused_line;
3801 for (
unsigned int c = 0; c < 2; ++c, ++
l)
3802 new_lines[
l] = cell->line(face_no)->child(c);
3805 new_lines[8]->set_bounding_object_indices(
3806 {new_vertices[6], new_vertices[8]});
3807 new_lines[9]->set_bounding_object_indices(
3808 {new_vertices[8], new_vertices[7]});
3809 new_lines[10]->set_bounding_object_indices(
3810 {new_vertices[4], new_vertices[8]});
3811 new_lines[11]->set_bounding_object_indices(
3812 {new_vertices[8], new_vertices[5]});
3821 new_lines[0] = cell->line(0);
3822 new_lines[1] = cell->line(1);
3823 new_lines[2] = cell->line(2)->child(0);
3824 new_lines[3] = cell->line(2)->child(1);
3825 new_lines[4] = cell->line(3)->child(0);
3826 new_lines[5] = cell->line(3)->child(1);
3827 new_lines[6]->set_bounding_object_indices(
3828 {new_vertices[6], new_vertices[7]});
3838 new_lines[0] = cell->line(0)->child(0);
3839 new_lines[1] = cell->line(0)->child(1);
3840 new_lines[2] = cell->line(1)->child(0);
3841 new_lines[3] = cell->line(1)->child(1);
3842 new_lines[4] = cell->line(2);
3843 new_lines[5] = cell->line(3);
3844 new_lines[6]->set_bounding_object_indices(
3845 {new_vertices[4], new_vertices[5]});
3848 for (
unsigned int l = lmin;
l < lmax; ++
l)
3850 new_lines[
l]->set_used_flag();
3851 new_lines[
l]->clear_user_flag();
3853 new_lines[
l]->clear_children();
3855 new_lines[
l]->set_boundary_id_internal(
3857 new_lines[
l]->set_manifold_id(cell->manifold_id());
3864 while (next_unused_cell->used() ==
true)
3868 for (
unsigned int i = 0; i < n_children; ++i)
3871 subcells[i] = next_unused_cell;
3873 if (i % 2 == 1 && i < n_children - 1)
3874 while (next_unused_cell->used() ==
true)
3892 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3893 new_lines[8]->
index(),
3894 new_lines[4]->
index(),
3895 new_lines[10]->
index()});
3896 subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3897 new_lines[2]->
index(),
3898 new_lines[5]->
index(),
3899 new_lines[11]->
index()});
3900 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3901 new_lines[9]->
index(),
3902 new_lines[10]->
index(),
3903 new_lines[6]->
index()});
3904 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3905 new_lines[3]->
index(),
3906 new_lines[11]->
index(),
3907 new_lines[7]->
index()});
3923 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3924 new_lines[6]->
index(),
3925 new_lines[2]->
index(),
3926 new_lines[4]->
index()});
3927 subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3928 new_lines[1]->
index(),
3929 new_lines[3]->
index(),
3930 new_lines[5]->
index()});
3947 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3948 new_lines[2]->
index(),
3949 new_lines[4]->
index(),
3950 new_lines[6]->
index()});
3951 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
3952 new_lines[3]->
index(),
3953 new_lines[6]->
index(),
3954 new_lines[5]->
index()});
3959 for (
unsigned int i = 0; i < n_children; ++i)
3961 subcells[i]->set_used_flag();
3962 subcells[i]->clear_refine_flag();
3963 subcells[i]->clear_user_flag();
3965 subcells[i]->clear_children();
3967 subcells[i]->set_material_id(cell->material_id());
3968 subcells[i]->set_manifold_id(cell->manifold_id());
3969 subcells[i]->set_subdomain_id(subdomainid);
3972 subcells[i]->set_parent(cell->index());
3978 for (
unsigned int i = 0; i < n_children / 2; ++i)
3979 cell->set_children(2 * i, subcells[2 * i]->index());
3981 cell->set_refinement_case(ref_case);
3989 for (
unsigned int c = 0; c < n_children; ++c)
3990 cell->child(c)->set_direction_flag(cell->direction_flag());
3995 template <
int dim,
int spacedim>
3998 const bool check_for_distorted_cells)
4004 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
4006 if (cell->refine_flag_set())
4019 line->clear_user_flag();
4023 unsigned int n_single_lines = 0;
4024 unsigned int n_lines_in_pairs = 0;
4025 unsigned int needed_vertices = 0;
4031 unsigned int needed_cells = 0;
4033 for (
const auto &cell :
4035 if (cell->refine_flag_set())
4037 if (cell->reference_cell() ==
4041 needed_vertices += 0;
4042 n_single_lines += 3;
4044 else if (cell->reference_cell() ==
4048 needed_vertices += 1;
4049 n_single_lines += 4;
4056 for (
const auto line_no : cell->face_indices())
4058 auto line = cell->line(line_no);
4059 if (line->has_children() ==
false)
4060 line->set_user_flag();
4065 const unsigned int used_cells =
4072 used_cells + needed_cells,
4084 if (line->user_flag_set())
4087 n_lines_in_pairs += 2;
4088 needed_vertices += 1;
4093 needed_vertices += std::count(
triangulation.vertices_used.begin(),
4103 unsigned int next_unused_vertex = 0;
4112 for (; line != endl; ++line)
4113 if (line->user_flag_set())
4119 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
4120 ++next_unused_vertex;
4124 "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."));
4127 triangulation.vertices[next_unused_vertex] = line->center(
true);
4129 bool pair_found =
false;
4131 for (; next_unused_line != endl; ++next_unused_line)
4132 if (!next_unused_line->used() &&
4133 !(++next_unused_line)->used())
4141 line->set_children(0, next_unused_line->index());
4144 children[2] = {next_unused_line, ++next_unused_line};
4149 children[0]->set_bounding_object_indices(
4150 {line->vertex_index(0), next_unused_vertex});
4151 children[1]->set_bounding_object_indices(
4152 {next_unused_vertex, line->vertex_index(1)});
4154 children[0]->set_used_flag();
4155 children[1]->set_used_flag();
4156 children[0]->clear_children();
4157 children[1]->clear_children();
4160 children[0]->clear_user_flag();
4161 children[1]->clear_user_flag();
4164 children[0]->set_boundary_id_internal(line->boundary_id());
4165 children[1]->set_boundary_id_internal(line->boundary_id());
4167 children[0]->set_manifold_id(line->manifold_id());
4168 children[1]->set_manifold_id(line->manifold_id());
4170 line->clear_user_flag();
4177 cells_with_distorted_children;
4183 unsigned int &next_unused_vertex,
4184 auto & next_unused_line,
4185 auto & next_unused_cell,
4186 const auto & cell) {
4187 const auto ref_case = cell->refine_flag_set();
4188 cell->clear_refine_flag();
4190 unsigned int n_new_vertices = 0;
4194 else if (cell->reference_cell() ==
4200 std::vector<int> new_vertices(n_new_vertices);
4201 for (
unsigned int vertex_no = 0; vertex_no < cell->
n_vertices();
4203 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4204 for (
unsigned int line_no = 0; line_no < cell->
n_lines(); ++line_no)
4205 if (cell->line(line_no)->has_children())
4207 cell->line(line_no)->child(0)->vertex_index(1);
4211 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
4212 ++next_unused_vertex;
4216 "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."));
4219 new_vertices[8] = next_unused_vertex;
4222 cell->center(
true,
true);
4225 std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
4228 unsigned int lmin = 0;
4229 unsigned int lmax = 0;
4236 else if (cell->reference_cell() ==
4247 for (
unsigned int l = lmin;
l < lmax; ++
l)
4249 while (next_unused_line->used() ==
true)
4251 new_lines[
l] = next_unused_line;
4262 const auto ref = [&](
const unsigned int face_no,
4263 const unsigned int vertex_no) {
4264 if (cell->line(face_no)->child(0)->vertex_index(0) ==
4265 static_cast<unsigned int>(new_vertices[vertex_no]) ||
4266 cell->line(face_no)->child(0)->vertex_index(1) ==
4267 static_cast<unsigned int>(new_vertices[vertex_no]))
4269 new_lines[2 * face_no + 0] =
4270 cell->line(face_no)->child(0);
4271 new_lines[2 * face_no + 1] =
4272 cell->line(face_no)->child(1);
4276 new_lines[2 * face_no + 0] =
4277 cell->line(face_no)->child(1);
4278 new_lines[2 * face_no + 1] =
4279 cell->line(face_no)->child(0);
4287 new_lines[6]->set_bounding_object_indices(
4288 {new_vertices[3], new_vertices[4]});
4289 new_lines[7]->set_bounding_object_indices(
4290 {new_vertices[4], new_vertices[5]});
4291 new_lines[8]->set_bounding_object_indices(
4292 {new_vertices[5], new_vertices[3]});
4294 else if (cell->reference_cell() ==
4298 for (
const unsigned int face_no : cell->face_indices())
4299 for (
unsigned int c = 0; c < 2; ++c, ++
l)
4300 new_lines[
l] = cell->line(face_no)->child(c);
4302 new_lines[8]->set_bounding_object_indices(
4303 {new_vertices[6], new_vertices[8]});
4304 new_lines[9]->set_bounding_object_indices(
4305 {new_vertices[8], new_vertices[7]});
4306 new_lines[10]->set_bounding_object_indices(
4307 {new_vertices[4], new_vertices[8]});
4308 new_lines[11]->set_bounding_object_indices(
4309 {new_vertices[8], new_vertices[5]});
4318 for (
unsigned int l = lmin;
l < lmax; ++
l)
4320 new_lines[
l]->set_used_flag();
4321 new_lines[
l]->clear_user_flag();
4322 new_lines[
l]->clear_user_data();
4323 new_lines[
l]->clear_children();
4325 new_lines[
l]->set_boundary_id_internal(
4327 new_lines[
l]->set_manifold_id(cell->manifold_id());
4332 while (next_unused_cell->used() ==
true)
4335 unsigned int n_children = 0;
4339 else if (cell->reference_cell() ==
4345 for (
unsigned int i = 0; i < n_children; ++i)
4348 subcells[i] = next_unused_cell;
4350 if (i % 2 == 1 && i < n_children - 1)
4351 while (next_unused_cell->used() ==
true)
4358 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4359 new_lines[8]->
index(),
4360 new_lines[5]->
index()});
4361 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4362 new_lines[2]->
index(),
4363 new_lines[6]->
index()});
4364 subcells[2]->set_bounding_object_indices({new_lines[7]->index(),
4365 new_lines[3]->
index(),
4366 new_lines[4]->
index()});
4367 subcells[3]->set_bounding_object_indices({new_lines[6]->index(),
4368 new_lines[7]->
index(),
4369 new_lines[8]->
index()});
4373 const auto ref = [&](
const unsigned int line_no,
4374 const unsigned int vertex_no,
4375 const unsigned int subcell_no,
4376 const unsigned int subcell_line_no) {
4377 if (new_lines[line_no]->vertex_index(1) !=
4378 static_cast<unsigned int>(new_vertices[vertex_no]))
4380 ->face_orientations[subcells[subcell_no]->index() *
4382 subcell_line_no] = 0;
4406 else if ((dim == 2) && (cell->reference_cell() ==
4409 subcells[0]->set_bounding_object_indices(
4410 {new_lines[0]->index(),
4411 new_lines[8]->
index(),
4412 new_lines[4]->
index(),
4413 new_lines[10]->
index()});
4414 subcells[1]->set_bounding_object_indices(
4415 {new_lines[8]->index(),
4416 new_lines[2]->
index(),
4417 new_lines[5]->
index(),
4418 new_lines[11]->
index()});
4419 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
4420 new_lines[9]->
index(),
4421 new_lines[10]->
index(),
4422 new_lines[6]->
index()});
4423 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
4424 new_lines[3]->
index(),
4425 new_lines[11]->
index(),
4426 new_lines[7]->
index()});
4435 for (
unsigned int i = 0; i < n_children; ++i)
4437 subcells[i]->set_used_flag();
4438 subcells[i]->clear_refine_flag();
4439 subcells[i]->clear_user_flag();
4441 subcells[i]->clear_children();
4444 subcells[i]->set_material_id(cell->material_id());
4445 subcells[i]->set_manifold_id(cell->manifold_id());
4446 subcells[i]->set_subdomain_id(subdomainid);
4451 ->reference_cell[subcells[i]->index()] = cell->reference_cell();
4454 subcells[i]->set_parent(cell->index());
4457 for (
unsigned int i = 0; i < n_children / 2; ++i)
4458 cell->set_children(2 * i, subcells[2 * i]->index());
4460 cell->set_refinement_case(ref_case);
4463 for (
unsigned int c = 0; c < n_children; ++c)
4464 cell->child(c)->set_direction_flag(cell->direction_flag());
4474 for (
const auto &cell :
4476 if (cell->refine_flag_set())
4484 if (cell->reference_cell() ==
4486 check_for_distorted_cells &&
4487 has_distorted_children<dim, spacedim>(cell))
4488 cells_with_distorted_children.distorted_cells.push_back(
4495 return cells_with_distorted_children;
4504 template <
int spacedim>
4509 const unsigned int dim = 1;
4513 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
4515 if (cell->refine_flag_set())
4528 unsigned int needed_vertices = 0;
4533 unsigned int flagged_cells = 0;
4535 for (
const auto &acell :
4537 if (acell->refine_flag_set())
4542 const unsigned int used_cells =
4562 needed_vertices += flagged_cells;
4567 needed_vertices += std::count(
triangulation.vertices_used.begin(),
4584 unsigned int next_unused_vertex = 0;
4591 for (
const auto &cell :
4593 if (cell->refine_flag_set())
4596 cell->clear_refine_flag();
4602 ++next_unused_vertex;
4606 "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."));
4621 while (next_unused_cell->used() ==
true)
4623 first_child = next_unused_cell;
4624 first_child->set_used_flag();
4628 second_child = next_unused_cell;
4629 second_child->set_used_flag();
4635 cell->set_children(0, first_child->index());
4636 first_child->clear_children();
4637 first_child->set_bounding_object_indices(
4638 {cell->vertex_index(0), next_unused_vertex});
4639 first_child->set_material_id(cell->material_id());
4640 first_child->set_manifold_id(cell->manifold_id());
4641 first_child->set_subdomain_id(subdomainid);
4642 first_child->set_direction_flag(cell->direction_flag());
4644 first_child->set_parent(cell->index());
4648 first_child->face(1)->set_manifold_id(cell->manifold_id());
4653 first_child->set_neighbor(1, second_child);
4655 first_child->set_neighbor(0, cell->neighbor(0));
4656 else if (cell->neighbor(0)->is_active())
4662 Assert(cell->neighbor(0)->level() <= cell->level(),
4664 first_child->set_neighbor(0, cell->neighbor(0));
4670 const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4671 first_child->set_neighbor(0,
4672 cell->neighbor(0)->child(nbnb));
4677 left_neighbor = cell->neighbor(0);
4678 while (left_neighbor->has_children())
4680 left_neighbor = left_neighbor->child(nbnb);
4681 left_neighbor->set_neighbor(nbnb, first_child);
4686 second_child->clear_children();
4687 second_child->set_bounding_object_indices(
4688 {next_unused_vertex, cell->vertex_index(1)});
4689 second_child->set_neighbor(0, first_child);
4690 second_child->set_material_id(cell->material_id());
4691 second_child->set_manifold_id(cell->manifold_id());
4692 second_child->set_subdomain_id(subdomainid);
4693 second_child->set_direction_flag(cell->direction_flag());
4696 second_child->set_neighbor(1, cell->neighbor(1));
4697 else if (cell->neighbor(1)->is_active())
4699 Assert(cell->neighbor(1)->level() <= cell->level(),
4701 second_child->set_neighbor(1, cell->neighbor(1));
4706 const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4707 second_child->set_neighbor(
4708 1, cell->neighbor(1)->child(nbnb));
4711 right_neighbor = cell->neighbor(1);
4712 while (right_neighbor->has_children())
4714 right_neighbor = right_neighbor->child(nbnb);
4715 right_neighbor->set_neighbor(nbnb, second_child);
4735 template <
int spacedim>
4738 const bool check_for_distorted_cells)
4740 const unsigned int dim = 2;
4746 bool do_isotropic_refinement =
true;
4747 for (
const auto &cell :
triangulation.active_cell_iterators())
4751 do_isotropic_refinement =
false;
4755 if (do_isotropic_refinement)
4757 check_for_distorted_cells);
4762 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
4764 if (cell->refine_flag_set())
4779 line->clear_user_flag();
4785 unsigned int n_single_lines = 0;
4789 unsigned int n_lines_in_pairs = 0;
4795 unsigned int needed_vertices = 0;
4800 unsigned int needed_cells = 0;
4802 for (
const auto &cell :
4804 if (cell->refine_flag_set())
4815 n_single_lines += 4;
4827 n_single_lines += 1;
4835 for (
const unsigned int line_no :
4839 cell->refine_flag_set(), line_no) ==
4843 line = cell->line(line_no);
4844 if (line->has_children() ==
false)
4845 line->set_user_flag();
4852 const unsigned int used_cells =
4862 used_cells + needed_cells,
4878 if (line->user_flag_set())
4881 n_lines_in_pairs += 2;
4882 needed_vertices += 1;
4894 needed_vertices += std::count(
triangulation.vertices_used.begin(),
4911 unsigned int next_unused_vertex = 0;
4923 for (; line != endl; ++line)
4924 if (line->user_flag_set())
4930 while (
triangulation.vertices_used[next_unused_vertex] ==
true)
4931 ++next_unused_vertex;
4935 "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."));
4938 triangulation.vertices[next_unused_vertex] = line->center(
true);
4943 bool pair_found =
false;
4945 for (; next_unused_line != endl; ++next_unused_line)
4946 if (!next_unused_line->used() &&
4947 !(++next_unused_line)->used())
4960 line->set_children(0, next_unused_line->index());
4964 children[2] = {next_unused_line, ++next_unused_line};
4970 children[0]->set_bounding_object_indices(
4971 {line->vertex_index(0), next_unused_vertex});
4972 children[1]->set_bounding_object_indices(
4973 {next_unused_vertex, line->vertex_index(1)});
4975 children[0]->set_used_flag();
4976 children[1]->set_used_flag();
4977 children[0]->clear_children();
4978 children[1]->clear_children();
4981 children[0]->clear_user_flag();
4982 children[1]->clear_user_flag();
4985 children[0]->set_boundary_id_internal(line->boundary_id());
4986 children[1]->set_boundary_id_internal(line->boundary_id());
4988 children[0]->set_manifold_id(line->manifold_id());
4989 children[1]->set_manifold_id(line->manifold_id());
4993 line->clear_user_flag();
5005 cells_with_distorted_children;
5019 for (
const auto &cell :
5021 if (cell->refine_flag_set())
5031 if (check_for_distorted_cells &&
5032 has_distorted_children<dim, spacedim>(cell))
5040 return cells_with_distorted_children;
5044 template <
int spacedim>
5047 const bool check_for_distorted_cells)
5049 static const int dim = 3;
5060 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
5062 if (cell->refine_flag_set())
5078 line->clear_user_flag();
5084 quad->clear_user_flag();
5098 unsigned int needed_vertices = 0;
5099 unsigned int needed_lines_single = 0;
5100 unsigned int needed_quads_single = 0;
5101 unsigned int needed_lines_pair = 0;
5102 unsigned int needed_quads_pair = 0;
5105 unsigned int new_cells = 0;
5107 for (
const auto &cell :
5109 if (cell->refine_flag_set())
5112 Assert(cell->refine_flag_set() ==
5118 new_cells += cell->reference_cell().n_isotropic_children();
5123 needed_lines_single += 6;
5124 needed_quads_single += 12;
5126 else if (cell->reference_cell() ==
5129 needed_lines_single += 1;
5130 needed_quads_single += 8;
5141 for (
const auto face : cell->face_indices())
5142 if (cell->face(face)->n_children() == 0)
5143 cell->face(face)->set_user_flag();
5145 Assert(cell->face(face)->n_children() ==
5146 cell->reference_cell()
5147 .face_reference_cell(face)
5148 .n_isotropic_children(),
5151 for (
const auto line : cell->line_indices())
5152 if (cell->line(line)->has_children() ==
false)
5153 cell->line(line)->set_user_flag();
5155 Assert(cell->line(line)->n_children() == 2,
5159 const unsigned int used_cells =
5165 used_cells + new_cells,
5179 if (quad->user_flag_set() ==
false)
5184 needed_quads_pair += 4;
5185 needed_lines_pair += 4;
5186 needed_vertices += 1;
5190 needed_quads_pair += 4;
5191 needed_lines_single += 3;
5204 if (line->user_flag_set() ==
false)
5207 needed_lines_pair += 2;
5208 needed_vertices += 1;
5213 needed_lines_single);
5216 needed_quads_single);
5219 needed_quads_single);
5223 needed_vertices += std::count(
triangulation.vertices_used.begin(),
5244 for (
const auto &cell :
triangulation.active_cell_iterators())
5245 if (!cell->refine_flag_set())
5246 for (
unsigned int line_n = 0; line_n < cell->
n_lines(); ++line_n)
5247 if (cell->line(line_n)->has_children())
5248 for (
unsigned int c = 0; c < 2; ++c)
5249 Assert(cell->line(line_n)->child(c)->user_flag_set() ==
false,
5253 unsigned int current_vertex = 0;
5257 auto get_next_unused_vertex = [](
const unsigned int current_vertex,
5258 std::vector<bool> &vertices_used) {
5259 unsigned int next_vertex = current_vertex;
5260 while (next_vertex < vertices_used.size() &&
5261 vertices_used[next_vertex] ==
true)
5264 vertices_used[next_vertex] =
true;
5277 for (; line != endl; ++line)
5279 if (line->user_flag_set() ==
false)
5283 get_next_unused_vertex(current_vertex,
5285 triangulation.vertices[current_vertex] = line->center(
true);
5288 triangulation.faces->lines.template next_free_pair_object<1>(
5296 line->set_children(0, next_unused_line->index());
5299 children[2] = {next_unused_line, ++next_unused_line};
5304 children[0]->set_bounding_object_indices(
5305 {line->vertex_index(0), current_vertex});
5306 children[1]->set_bounding_object_indices(
5307 {current_vertex, line->vertex_index(1)});
5309 children[0]->set_used_flag();
5310 children[1]->set_used_flag();
5311 children[0]->clear_children();
5312 children[1]->clear_children();
5315 children[0]->clear_user_flag();
5316 children[1]->clear_user_flag();
5318 children[0]->set_boundary_id_internal(line->boundary_id());
5319 children[1]->set_boundary_id_internal(line->boundary_id());
5321 children[0]->set_manifold_id(line->manifold_id());
5322 children[1]->set_manifold_id(line->manifold_id());
5324 line->clear_user_flag();
5338 for (; quad != endq; ++quad)
5340 if (quad->user_flag_set() ==
false)
5343 const auto reference_face_type = quad->reference_cell();
5349 get_next_unused_vertex(current_vertex,
5352 quad->center(
true,
true);
5356 boost::container::small_vector<
5361 for (
unsigned int i = 0; i < new_lines.size(); ++i)
5381 new_lines[i] = next_unused_line;
5395 triangulation.faces->quads.template next_free_pair_object<2>(
5398 new_quads[0] = next_unused_quad;
5402 new_quads[1] = next_unused_quad;
5406 triangulation.faces->quads.template next_free_pair_object<2>(
5408 new_quads[2] = next_unused_quad;
5412 new_quads[3] = next_unused_quad;
5415 quad->set_children(0, new_quads[0]->
index());
5416 quad->set_children(2, new_quads[2]->
index());
5425 for (
const auto i : quad->vertex_indices())
5428 for (
const auto i : quad->line_indices())
5430 quad->line(i)->child(0)->vertex_index(1);
5435 boost::container::small_vector<
5444 for (
unsigned int l = 0;
l < quad->
n_lines(); ++
l)
5445 for (
unsigned int c = 0; c < 2; ++c)
5447 static constexpr std::array<std::array<unsigned int, 2>,
5454 lines[k++] = quad->line(
l)->child(
5455 index[c][quad->line_orientation(
l)]);
5458 for (
unsigned int l = 0;
l < new_lines.size(); ++
l)
5459 lines[k++] = new_lines[
l];
5462 boost::container::small_vector<int, 12> line_indices(
5464 for (
unsigned int i = 0; i < line_indices.size(); ++i)
5465 line_indices[i] = lines[i]->
index();
5467 static constexpr std::array<std::array<unsigned int, 2>, 12>
5468 line_vertices_quad{{{{0, 4}},
5481 static constexpr std::array<std::array<unsigned int, 4>, 4>
5482 quad_lines_quad{{{{0, 8, 4, 10}},
5487 static constexpr std::
5488 array<std::array<std::array<unsigned int, 2>, 4>, 4>
5489 quad_line_vertices_quad{
5490 {{{{{0, 4}}, {{6, 8}}, {{0, 6}}, {{4, 8}}}},
5491 {{{{6, 8}}, {{1, 5}}, {{6, 1}}, {{8, 5}}}},
5492 {{{{4, 2}}, {{8, 7}}, {{4, 8}}, {{2, 7}}}},
5493 {{{{8, 7}}, {{5, 3}}, {{8, 5}}, {{7, 3}}}}}};
5495 static constexpr std::array<std::array<unsigned int, 2>, 12>
5496 line_vertices_tri{{{{0, 3}},
5509 static constexpr std::array<std::array<unsigned int, 4>, 4>
5510 quad_lines_tri{{{{0, 8, 5, X}},
5515 static constexpr std::
5516 array<std::array<std::array<unsigned int, 2>, 4>, 4>
5517 quad_line_vertices_tri{
5518 {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
5519 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
5520 {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
5521 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
5523 const auto &line_vertices =
5525 line_vertices_quad :
5527 const auto &quad_lines =
5531 const auto &quad_line_vertices =
5533 quad_line_vertices_quad :
5534 quad_line_vertices_tri;
5537 for (
unsigned int i = 0, j = lines.size() - new_lines.size();
5538 i < new_lines.size();
5541 auto &new_line = new_lines[i];
5542 new_line->set_bounding_object_indices(
5545 new_line->set_used_flag();
5546 new_line->clear_user_flag();
5547 new_line->clear_user_data();
5548 new_line->clear_children();
5549 new_line->set_boundary_id_internal(quad->boundary_id());
5550 new_line->set_manifold_id(quad->manifold_id());
5554 for (
unsigned int i = 0; i < new_quads.size(); ++i)
5556 auto &new_quad = new_quads[i];
5560 triangulation.faces->quad_reference_cell[new_quad->index()] =
5561 reference_face_type;
5563 if (new_quad->n_lines() == 3)
5564 new_quad->set_bounding_object_indices(
5565 {line_indices[quad_lines[i][0]],
5566 line_indices[quad_lines[i][1]],
5567 line_indices[quad_lines[i][2]]});
5568 else if (new_quad->n_lines() == 4)
5569 new_quad->set_bounding_object_indices(
5570 {line_indices[quad_lines[i][0]],
5571 line_indices[quad_lines[i][1]],
5572 line_indices[quad_lines[i][2]],
5573 line_indices[quad_lines[i][3]]});
5577 new_quad->set_used_flag();
5578 new_quad->clear_user_flag();
5579 new_quad->clear_user_data();
5580 new_quad->clear_children();
5581 new_quad->set_boundary_id_internal(quad->boundary_id());
5582 new_quad->set_manifold_id(quad->manifold_id());
5585 std::set<unsigned int> s;
5589 for (
const auto f : new_quad->line_indices())
5591 std::array<unsigned int, 2> vertices_0, vertices_1;
5593 for (
unsigned int v = 0; v < 2; ++v)
5595 lines[quad_lines[i][f]]->vertex_index(v);
5597 for (
unsigned int v = 0; v < 2; ++v)
5601 const auto orientation =
5606 for (
const auto i : vertices_0)
5608 for (
const auto i : vertices_1)
5612 new_quad->set_line_orientation(f, orientation);
5622 quad->clear_user_flag();
5627 cells_with_distorted_children;
5638 for (; hex != endh; ++hex)
5640 if (hex->refine_flag_set() ==
5644 const auto &reference_cell_type = hex->reference_cell();
5647 hex->clear_refine_flag();
5648 hex->set_refinement_case(ref_case);
5650 unsigned int n_new_lines = 0;
5651 unsigned int n_new_quads = 0;
5652 unsigned int n_new_hexes = 0;
5673 get_next_unused_vertex(current_vertex,
5676 hex->center(
true,
true);
5679 boost::container::small_vector<
5682 new_lines(n_new_lines);
5683 for (
unsigned int i = 0; i < n_new_lines; ++i)
5690 new_lines[i]->set_used_flag();
5691 new_lines[i]->clear_user_flag();
5692 new_lines[i]->clear_user_data();
5693 new_lines[i]->clear_children();
5694 new_lines[i]->set_boundary_id_internal(
5696 new_lines[i]->set_manifold_id(hex->manifold_id());
5699 boost::container::small_vector<
5702 new_quads(n_new_quads);
5703 for (
unsigned int i = 0; i < n_new_quads; ++i)
5709 auto &new_quad = new_quads[i];
5714 ->quad_reference_cell[new_quad->index()] =
5720 new_quad->set_used_flag();
5721 new_quad->clear_user_flag();
5722 new_quad->clear_user_data();
5723 new_quad->clear_children();
5724 new_quad->set_boundary_id_internal(
5726 new_quad->set_manifold_id(hex->manifold_id());
5727 for (
const auto j : new_quads[i]->line_indices())
5728 new_quad->set_line_orientation(j,
true);
5737 for (
unsigned int i = 0; i < n_new_hexes; ++i)
5746 new_hexes[i] = next_unused_hex;
5748 auto &new_hex = new_hexes[i];
5752 ->reference_cell[new_hex->index()] =
5753 reference_cell_type;
5756 new_hex->set_used_flag();
5757 new_hex->clear_user_flag();
5758 new_hex->clear_user_data();
5759 new_hex->clear_children();
5760 new_hex->set_material_id(hex->material_id());
5761 new_hex->set_manifold_id(hex->manifold_id());
5762 new_hex->set_subdomain_id(hex->subdomain_id());
5765 new_hex->set_parent(hex->index());
5777 for (
const auto f : new_hex->face_indices())
5779 new_hex->set_face_orientation(f,
true);
5780 new_hex->set_face_flip(f,
false);
5781 new_hex->set_face_rotation(f,
false);
5784 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
5785 hex->set_children(2 * i, new_hexes[2 * i]->index());
5795 for (
const unsigned int i : hex->vertex_indices())
5798 for (
const unsigned int i : hex->line_indices())
5800 hex->line(i)->child(0)->vertex_index(1);
5804 for (
const unsigned int i : hex->face_indices())
5806 middle_vertex_index<dim, spacedim>(hex->face(i));
5814 static constexpr std::array<std::array<unsigned int, 2>, 6>
5815 new_line_vertices_hex = {{{{22, 26}},
5822 static constexpr std::array<std::array<unsigned int, 2>, 6>
5823 new_line_vertices_tet = {{{{6, 8}},
5830 const auto &new_line_vertices =
5832 new_line_vertices_hex :
5833 new_line_vertices_tet;
5835 for (
unsigned int i = 0; i < new_lines.size(); ++i)
5836 new_lines[i]->set_bounding_object_indices(
5843 boost::container::small_vector<
5850 relevant_lines.resize(30);
5851 for (
unsigned int f = 0, k = 0; f < 6; ++f)
5852 for (
unsigned int c = 0; c < 4; ++c, ++k)
5854 static constexpr std::
5855 array<std::array<unsigned int, 2>, 4>
5857 {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}};
5863 standard_to_real_face_vertex(
5865 hex->face_orientation(f),
5867 hex->face_rotation(f)))
5869 standard_to_real_face_line(
5871 hex->face_orientation(f),
5873 hex->face_rotation(f)));
5876 for (
unsigned int i = 0, k = 24; i < 6; ++i, ++k)
5877 relevant_lines[k] = new_lines[i];
5881 relevant_lines.resize(13);
5884 for (
unsigned int f = 0; f < 4; ++f)
5885 for (
unsigned int l = 0;
l < 3; ++
l, ++k)
5889 array<std::array<unsigned int, 3>, 6>
5890 table = {{{{1, 0, 2}},
5905 dim>::faces_per_cell +
5909 relevant_lines[k++] = new_lines[0];
5916 boost::container::small_vector<unsigned int, 30>
5917 relevant_line_indices(relevant_lines.size());
5918 for (
unsigned int i = 0; i < relevant_line_indices.size();
5920 relevant_line_indices[i] = relevant_lines[i]->
index();
5922 static constexpr std::array<std::array<unsigned int, 4>, 12>
5923 new_quad_lines_hex = {{{{10, 28, 16, 24}},
5934 {{25, 7, 27, 13}}}};
5936 static constexpr std::array<std::array<unsigned int, 4>, 12>
5937 new_quad_lines_tet = {{{{2, 3, 8, X}},
5950 static constexpr std::
5951 array<std::array<std::array<unsigned int, 2>, 4>, 12>
5953 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
5954 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
5955 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
5956 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
5957 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
5958 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
5959 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
5960 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
5961 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
5962 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
5963 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
5964 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
5966 static constexpr std::
5967 array<std::array<std::array<unsigned int, 2>, 4>, 12>
5969 {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
5970 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
5971 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
5972 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
5973 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
5974 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
5975 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
5976 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
5977 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5978 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5979 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5980 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
5982 const auto &new_quad_lines =
5984 new_quad_lines_hex :
5992 for (
unsigned int q = 0; q < new_quads.size(); ++q)
5994 for (
unsigned int l = 0;
l < 3; ++
l)
5996 std::array<unsigned int, 2> vertices_0, vertices_1;
5998 for (
unsigned int v = 0; v < 2; ++v)
6000 relevant_lines[new_quad_lines[q][
l]]
6003 for (
unsigned int v = 0; v < 2; ++v)
6008 for (
unsigned int q = 0; q < new_quads.size(); ++q)
6010 auto &new_quad = new_quads[q];
6012 if (new_quad->n_lines() == 3)
6013 new_quad->set_bounding_object_indices(
6014 {relevant_line_indices[new_quad_lines[q][0]],
6015 relevant_line_indices[new_quad_lines[q][1]],
6016 relevant_line_indices[new_quad_lines[q][2]]});
6017 else if (new_quad->n_lines() == 4)
6018 new_quad->set_bounding_object_indices(
6019 {relevant_line_indices[new_quad_lines[q][0]],
6020 relevant_line_indices[new_quad_lines[q][1]],
6021 relevant_line_indices[new_quad_lines[q][2]],
6022 relevant_line_indices[new_quad_lines[q][3]]});
6026 for (
const auto l : new_quad->line_indices())
6028 std::array<unsigned int, 2> vertices_0, vertices_1;
6030 for (
unsigned int v = 0; v < 2; ++v)
6032 relevant_lines[new_quad_lines[q][
l]]
6035 for (
unsigned int v = 0; v < 2; ++v)
6038 const auto orientation =
6040 vertices_0, vertices_1);
6042 new_quad->set_line_orientation(
l, orientation);
6049 std::array<int, 36> quad_indices;
6053 for (
unsigned int i = 0; i < new_quads.size(); ++i)
6054 quad_indices[i] = new_quads[i]->
index();
6056 for (
unsigned int f = 0, k = new_quads.size(); f < 6;
6058 for (
unsigned int c = 0; c < 4; ++c, ++k)
6060 hex->face(f)->isotropic_child_index(
6063 hex->face_orientation(f),
6065 hex->face_rotation(f)));
6069 for (
unsigned int i = 0; i < new_quads.size(); ++i)
6070 quad_indices[i] = new_quads[i]->
index();
6072 for (
unsigned int f = 0, k = new_quads.size(); f < 4;
6074 for (
unsigned int c = 0; c < 4; ++c, ++k)
6076 quad_indices[k] = hex->face(f)->child_index(
6080 .standard_to_real_face_vertex(
6095 static constexpr std::array<std::array<unsigned int, 6>, 8>
6097 {{12, 0, 20, 4, 28, 8}},
6098 {{0, 16, 22, 6, 29, 9}},
6099 {{13, 1, 4, 24, 30, 10}},
6100 {{1, 17, 6, 26, 31, 11}},
6101 {{14, 2, 21, 5, 8, 32}},
6102 {{2, 18, 23, 7, 9, 33}},
6103 {{15, 3, 5, 25, 10, 34}},
6104 {{3, 19, 7, 27, 11, 35}}
6107 static constexpr std::array<std::array<unsigned int, 6>, 8>
6108 cell_quads_tet{{{{8, 13, 16, 0, X, X}},
6109 {{9, 12, 1, 21, X, X}},
6110 {{10, 2, 17, 20, X, X}},
6111 {{3, 14, 18, 22, X, X}},
6112 {{11, 1, 4, 5, X, X}},
6113 {{15, 0, 4, 6, X, X}},
6114 {{19, 7, 6, 3, X, X}},
6115 {{23, 5, 2, 7, X, X}}}};
6117 static constexpr std::
6118 array<std::array<std::array<unsigned int, 4>, 6>, 8>
6119 cell_face_vertices_hex{{{{{{0, 8, 16, 20}},
6124 {{16, 22, 20, 26}}}},
6125 {{{{10, 24, 22, 26}},
6130 {{22, 17, 26, 21}}}},
6136 {{20, 26, 18, 23}}}},
6137 {{{{24, 11, 26, 23}},
6142 {{26, 21, 23, 19}}}},
6143 {{{{16, 20, 4, 12}},
6148 {{4, 14, 12, 25}}}},
6149 {{{{22, 26, 14, 25}},
6154 {{14, 5, 25, 13}}}},
6155 {{{{20, 18, 12, 6}},