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())
398 for (
const unsigned int v : cell->vertex_indices())
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)
441 std::swap(cell.vertices[2], cell.vertices[3]);
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)
457 tmp[i] = cell.vertices[i];
459 cell.vertices[local_vertex_numbering[i]] = tmp[i];
465 std::swap(boundary_quad.vertices[2], boundary_quad.vertices[3]);
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))
600 if (determinants[i] <= 1
e-9 * std::pow(cell->diameter(), 1. * dim))
602 distorted_cells.distorted_cells.push_back(cell);
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] <=
635 1
e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
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 =
1772 if (level_objects > 0)
1805 for (; line != endc; ++line)
1808 if (line->has_children() ==
false)
1826 for (; line != endc; ++line)
1829 if (line->has_children() ==
false)
1849 template <
int dim,
int spacedim>
1853 const unsigned int level_objects,
1864 &compute_number_cache_dim<dim, spacedim>),
1870 using quad_iterator =
1887 unsigned int n_levels = 0;
1888 if (level_objects > 0)
1892 n_levels =
level + 1;
1905 (
level == n_levels - 1 ?
1908 for (; quad != endc; ++quad)
1911 if (quad->has_children() ==
false)
1929 for (; quad != endc; ++quad)
1932 if (quad->has_children() ==
false)
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 =
1995 unsigned int n_levels = 0;
1996 if (level_objects > 0)
2000 n_levels =
level + 1;
2012 endc = (
level == n_levels - 1 ?
2015 for (; hex != endc; ++hex)
2018 if (hex->has_children() ==
false)
2036 for (; hex != endc; ++hex)
2039 if (hex->has_children() ==
false)
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 =
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)
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)
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] <<
')';
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)
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);
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();
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())
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();
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))
4462 cells_with_distorted_children.distorted_cells.push_back(
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();
4602 second_child = next_unused_cell;
4603 second_child->set_used_flag();
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)});
5278 for (
const auto &child : children)
5280 child->set_used_flag();
5281 child->clear_children();
5283 child->clear_user_flag();
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());
5711 for (
const unsigned int i : hex->vertex_indices())
5714 const std::array<unsigned int, 12> line_indices =
5715 TriaAccessorImplementation::Implementation::
5716 get_line_indices_of_cell(*hex);
5720 for (
unsigned int l = 0;
l < n_lines; ++
l)
5730 for (
const unsigned int i : hex->face_indices())
5732 hex->face(i)->child(0)->vertex_index(3);
5736 get_next_unused_vertex(current_vertex,
5741 hex->center(
true,
true);
5747 static constexpr ::ndarray<unsigned int, 6, 2>
5748 new_line_vertices_hex = {{{{22, 26}},
5755 static constexpr ::ndarray<unsigned int, 6, 2>
5756 new_line_vertices_tet = {{{{6, 8}},
5763 const auto &new_line_vertices =
5765 new_line_vertices_hex :
5766 new_line_vertices_tet;
5768 for (
unsigned int i = 0; i < n_new_lines; ++i)
5769 new_lines[i]->set_bounding_object_indices(
5776 boost::container::small_vector<raw_line_iterator, 30>
5781 relevant_lines.resize(30);
5782 for (
unsigned int f = 0, k = 0; f < 6; ++f)
5783 for (
unsigned int c = 0; c < 4; ++c, ++k)
5786 ndarray<unsigned int, 4, 2>
5788 {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}};
5794 standard_to_real_face_vertex(
5796 hex->face_orientation(f),
5798 hex->face_rotation(f)))
5800 standard_to_real_face_line(
5802 hex->face_orientation(f),
5804 hex->face_rotation(f)));
5807 for (
unsigned int i = 0, k = 24; i < 6; ++i, ++k)
5808 relevant_lines[k] = new_lines[i];
5812 relevant_lines.resize(13);
5815 for (
unsigned int f = 0; f < 4; ++f)
5816 for (
unsigned int l = 0;
l < 3; ++
l, ++k)
5820 array<std::array<unsigned int, 3>, 6>
5821 table = {{{{1, 0, 2}},
5834 .get_combined_orientation(
5840 relevant_lines[k++] = new_lines[0];
5847 boost::container::small_vector<unsigned int, 30>
5848 relevant_line_indices(relevant_lines.size());
5849 for (
unsigned int i = 0; i < relevant_line_indices.size();
5851 relevant_line_indices[i] = relevant_lines[i]->
index();
5853 static constexpr ::ndarray<unsigned int, 12, 4>
5854 new_quad_lines_hex = {{{{10, 28, 16, 24}},
5865 {{25, 7, 27, 13}}}};
5867 static constexpr ::ndarray<unsigned int, 12, 4>
5868 new_quad_lines_tet = {{{{2, 3, 8, X}},
5881 static constexpr ::ndarray<unsigned int, 12, 4, 2>
5883 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
5884 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
5885 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
5886 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
5887 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
5888 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
5889 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
5890 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
5891 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
5892 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
5893 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
5894 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
5896 static constexpr ::ndarray<unsigned int, 12, 4, 2>
5898 {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
5899 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
5900 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
5901 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
5902 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
5903 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
5904 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
5905 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
5906 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5907 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5908 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5909 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
5911 const auto &new_quad_lines =
5913 new_quad_lines_hex :
5921 static constexpr ::ndarray<unsigned int, 4, 2>
5922 representative_lines{
5923 {{{0, 2}}, {{2, 0}}, {{3, 3}}, {{1, 1}}}};
5925 for (
unsigned int q = 0; q < n_new_quads; ++q)
5927 auto &new_quad = new_quads[q];
5929 if (new_quad->n_lines() == 3)
5930 new_quad->set_bounding_object_indices(
5931 {relevant_line_indices[new_quad_lines[q][0]],
5932 relevant_line_indices[new_quad_lines[q][1]],
5933 relevant_line_indices[new_quad_lines[q][2]]});
5934 else if (new_quad->n_lines() == 4)
5935 new_quad->set_bounding_object_indices(
5936 {relevant_line_indices[new_quad_lines[q][0]],
5937 relevant_line_indices[new_quad_lines[q][1]],
5938 relevant_line_indices[new_quad_lines[q][2]],
5939 relevant_line_indices[new_quad_lines[q][3]]});
5947 const unsigned int n_compute_lines =
5951 for (
unsigned int line = 0; line < n_compute_lines;
5954 const unsigned int l =
5955 (reference_cell_type ==
5957 representative_lines[q % 4][0] :
5960 const std::array<unsigned int, 2> vertices_0 = {
5961 {relevant_lines[new_quad_lines[q][
l]]
5963 relevant_lines[new_quad_lines[q][
l]]
5964 ->vertex_index(1)}};
5966 const std::array<unsigned int, 2> vertices_1 = {
5970 const auto orientation =
5975 new_quad->set_line_orientation(
l, orientation);
5980 if (reference_cell_type ==
5982 new_quads[representative_lines[q % 4][1] + q -
5984 ->set_line_orientation(
l, orientation);
5991 std::array<int, 36> quad_indices;
5995 for (
unsigned int i = 0; i < n_new_quads; ++i)
5996 quad_indices[i] = new_quads[i]->
index();
5998 for (
unsigned int f = 0, k = n_new_quads; f < 6; ++f)
5999 for (
unsigned int c = 0; c < 4; ++c, ++k)
6001 hex->face(f)->isotropic_child_index(
6004 hex->face_orientation(f),
6006 hex->face_rotation(f)));
6010 for (
unsigned int i = 0; i < n_new_quads; ++i)
6011 quad_indices[i] = new_quads[i]->
index();
6013 for (
unsigned int f = 0, k = n_new_quads; f < 4; ++f)
6014 for (
unsigned int c = 0; c < 4; ++c, ++k)
6016 quad_indices[k] = hex->face(f)->child_index(
6020 .standard_to_real_face_vertex(
6025 .get_combined_orientation(
6036 static constexpr ::ndarray<unsigned int, 8, 6>
6038 {{12, 0, 20, 4, 28, 8}},
6039 {{0, 16, 22, 6, 29, 9}},
6040 {{13, 1, 4, 24, 30, 10}},
6041 {{1, 17, 6, 26, 31, 11}},
6042 {{14, 2, 21, 5, 8, 32}},
6043 {{2, 18, 23, 7, 9, 33}},
6044 {{15, 3, 5, 25, 10, 34}},
6045 {{3, 19, 7, 27, 11, 35}}
6048 static constexpr ::ndarray<unsigned int, 8, 6>
6049 cell_quads_tet{{{{8, 13, 16, 0, X, X}},
6050 {{9, 12, 1, 21, X, X}},
6051 {{10, 2, 17, 20, X, X}},
6052 {{3, 14, 18, 22, X, X}},
6053 {{11, 1, 4, 5, X, X}},
6054 {{15, 0, 4, 6, X, X}},
6055 {{19, 7, 6, 3, X, X}},
6056 {{23, 5, 2, 7, X, X}}}};
6058 static constexpr ::ndarray<unsigned int, 8, 6, 4>
6059 cell_face_vertices_tet{{{{{{0, 4, 6, X}},
6108 const auto &cell_quads =
6113 for (
unsigned int c = 0;
6114 c < GeometryInfo<dim>::max_children_per_cell;
6117 auto &new_hex = new_hexes[c];
6119 if (new_hex->n_faces() == 4)
6121 new_hex->set_bounding_object_indices(
6122 {quad_indices[cell_quads[c][0]],
6123 quad_indices[cell_quads[c][1]],
6124 quad_indices[cell_quads[c][2]],
6125 quad_indices[cell_quads[c][3]]});
6129 for (
const auto f : new_hex->face_indices())
6131 const auto &face = new_hex->face(f);
6133 Assert(face->n_vertices() == 3,
6136 const std::array<unsigned int, 3> vertices_0 = {
6137 {face->vertex_index(0),
6138 face->vertex_index(1),
6139 face->vertex_index(2)}};
6141 const std::array<unsigned int, 3> vertices_1 = {
6151 new_hex->set_combined_face_orientation(
6153 face->reference_cell()
6154 .get_combined_orientation(
6159 else if (new_hex->n_faces() == 6)
6160 new_hex->set_bounding_object_indices(
6161 {quad_indices[cell_quads[c][0]],
6162 quad_indices[cell_quads[c][1]],
6163 quad_indices[cell_quads[c][2]],
6164 quad_indices[cell_quads[c][3]],
6165 quad_indices[cell_quads[c][4]],
6166 quad_indices[cell_quads[c][5]]});
6175 static constexpr ::ndarray<unsigned int, 6, 4>
6176 face_to_child_indices_hex{{{{0, 2, 4, 6}},
6183 for (
const auto f : hex->face_indices())
6185 const unsigned char combined_orientation =
6186 hex->combined_face_orientation(f);
6187 for (
unsigned int c = 0; c < 4; ++c)
6188 new_hexes[face_to_child_indices_hex[f][c]]
6189 ->set_combined_face_orientation(
6190 f, combined_orientation);
6195 if (check_for_distorted_cells &&
6196 has_distorted_children<dim, spacedim>(hex))
6205 return cells_with_distorted_children;
6212 template <
int spacedim>
6215 const bool check_for_distorted_cells)
6217 const unsigned int dim = 3;
6220 bool flag_isotropic_mesh =
true;
6224 for (; cell != endc; ++cell)
6234 flag_isotropic_mesh =
false;
6238 if (flag_isotropic_mesh)
6240 check_for_distorted_cells);
6251 for (
const auto &cell :
triangulation.active_cell_iterators_on_level(
6253 if (cell->refine_flag_set())
6270 line->clear_user_flag();
6275 quad->clear_user_flag();
6298 unsigned int needed_vertices = 0;
6299 unsigned int needed_lines_single = 0;
6300 unsigned int needed_quads_single = 0;
6301 unsigned int needed_lines_pair = 0;
6302 unsigned int needed_quads_pair = 0;
6307 unsigned int new_cells = 0;
6309 for (
const auto &acell :
6311 if (acell->refine_flag_set())
6321 ++needed_quads_single;
6329 ++needed_lines_single;
6330 needed_quads_single += 4;
6337 needed_lines_single += 6;
6338 needed_quads_single += 12;
6352 for (
const unsigned int face :
6356 aface = acell->face(face);
6363 acell->face_orientation(face),
6364 acell->face_flip(face),
6365 acell->face_rotation(face));
6370 if (face_ref_case ==
6373 if (aface->n_active_descendants() < 4)
6376 aface->set_user_flag();
6378 else if (aface->refinement_case() != face_ref_case)
6389 Assert(aface->refinement_case() ==
6391 dim - 1>::isotropic_refinement ||
6392 aface->refinement_case() ==
6395 aface->set_user_index(face_ref_case);
6401 for (
unsigned int line = 0;
6402 line < GeometryInfo<dim>::lines_per_cell;
6406 !acell->line(line)->has_children())
6407 acell->line(line)->set_user_flag();
6413 const unsigned int used_cells =
6423 used_cells + new_cells,