489 auto cell_rel_it = cell_relations.begin();
491 for (; cell_rel_it != cell_relations.end();
494 cell_rel_it->second =
496 dest_fixed_it + size,
503 template <
int dim,
int spacedim>
509 const
unsigned int handle,
513 const
boost::iterator_range<
std::vector<
char>::const_iterator> &)>
514 &unpack_callback)
const
520 const bool callback_variable_transfer = (handle % 2 == 0);
521 const unsigned int callback_index = handle / 2;
529 if (cell_relations.size() > 0)
535 std::vector<char>::const_iterator dest_data_it;
536 std::vector<char>::const_iterator dest_sizes_cell_it;
546 if (callback_variable_transfer)
559 const unsigned int offset_variable_data_sizes =
568 offset_variable_data_sizes +
569 callback_index *
sizeof(
unsigned int);
587 if (cell_relations.begin() != cell_relations.end())
592 auto cell_rel_it = cell_relations.begin();
594 for (; cell_rel_it != cell_relations.end(); ++cell_rel_it)
596 const auto &dealii_cell = cell_rel_it->first;
597 const auto &cell_status = cell_rel_it->second;
599 if (callback_variable_transfer)
603 data_increment = *dest_sizes_it;
610 if (callback_index == 0)
614 &(*(dest_sizes_cell_it -
sizeof(
unsigned int))),
615 sizeof(
unsigned int));
618 &(*dest_sizes_cell_it),
619 sizeof(
unsigned int));
626 dest_data_it += offset;
627 data_increment -= offset;
632 if (cell_rel_it != cell_relations.end() - 1)
641 unpack_callback(dealii_cell,
643 boost::make_iterator_range(dest_data_it,
644 dest_data_it + size));
648 unpack_callback(dealii_cell->parent(),
650 boost::make_iterator_range(dest_data_it,
651 dest_data_it + size));
663 if (cell_rel_it != cell_relations.end() - 1)
664 dest_data_it += data_increment;
670 template <
int dim,
int spacedim>
673 const
unsigned int global_first_cell,
674 const
unsigned int global_num_cells,
675 const
std::
string &file_basename,
676 const
MPI_Comm &mpi_communicator)
const
681#ifdef DEAL_II_WITH_MPI
686 const unsigned int myrank =
688 const unsigned int mpisize =
699 const std::string fname_fixed =
700 std::string(file_basename) +
"_fixed.data";
703 int ierr = MPI_Info_create(&info);
707 ierr = MPI_File_open(mpi_communicator,
709 MPI_MODE_CREATE | MPI_MODE_WRONLY,
714 ierr = MPI_File_set_size(fh, 0);
718 ierr = MPI_Barrier(mpi_communicator);
720 ierr = MPI_Info_free(&info);
741 const MPI_Offset size_header =
746 const MPI_Offset my_global_file_position =
748 static_cast<MPI_Offset
>(global_first_cell) * bytes_per_cell;
752 my_global_file_position,
759 ierr = MPI_File_close(&fh);
770 const std::string fname_variable =
771 std::string(file_basename) +
"_variable.data";
774 int ierr = MPI_Info_create(&info);
778 ierr = MPI_File_open(mpi_communicator,
779 fname_variable.c_str(),
780 MPI_MODE_CREATE | MPI_MODE_WRONLY,
785 ierr = MPI_File_set_size(fh, 0);
789 ierr = MPI_Barrier(mpi_communicator);
791 ierr = MPI_Info_free(&info);
796 const MPI_Offset my_global_file_position =
797 static_cast<MPI_Offset
>(global_first_cell) *
798 sizeof(
unsigned int);
803 static_cast<std::size_t
>(
804 std::numeric_limits<int>::max()),
809 my_global_file_position,
821 std::uint64_t prefix_sum = 0;
822 ierr = MPI_Exscan(&size_on_proc,
830 const MPI_Offset my_global_file_position =
831 static_cast<MPI_Offset
>(global_num_cells) *
sizeof(
unsigned int) +
837 my_global_file_position,
845 ierr = MPI_File_close(&fh);
852 (void)global_first_cell;
853 (void)global_num_cells;
854 (void)mpi_communicator;
860 const std::string fname_fixed =
861 std::string(file_basename) +
"_fixed.data";
863 std::ofstream file(fname_fixed, std::ios::binary | std::ios::out);
867 file.write(
reinterpret_cast<const char *
>(
881 const std::string fname_variable =
882 std::string(file_basename) +
"_variable.data";
884 std::ofstream file(fname_variable,
885 std::ios::binary | std::ios::out);
889 file.write(
reinterpret_cast<const char *
>(
901 template <
int dim,
int spacedim>
904 const
unsigned int global_first_cell,
905 const
unsigned int global_num_cells,
906 const
unsigned int local_num_cells,
907 const
std::
string &file_basename,
908 const
unsigned int n_attached_deserialize_fixed,
909 const
unsigned int n_attached_deserialize_variable,
913 ExcMessage(
"Previously loaded data has not been released yet!"));
917#ifdef DEAL_II_WITH_MPI
922 const unsigned int mpisize =
931 const std::string fname_fixed =
932 std::string(file_basename) +
"_fixed.data";
935 int ierr = MPI_Info_create(&info);
939 ierr = MPI_File_open(
940 mpi_communicator, fname_fixed.c_str(), MPI_MODE_RDONLY, info, &fh);
943 ierr = MPI_Info_free(&info);
967 const MPI_Offset size_header =
972 const MPI_Offset my_global_file_position =
974 static_cast<MPI_Offset
>(global_first_cell) * bytes_per_cell;
978 my_global_file_position,
986 ierr = MPI_File_close(&fh);
995 const std::string fname_variable =
996 std::string(file_basename) +
"_variable.data";
999 int ierr = MPI_Info_create(&info);
1003 ierr = MPI_File_open(mpi_communicator,
1004 fname_variable.c_str(),
1010 ierr = MPI_Info_free(&info);
1016 const MPI_Offset my_global_file_position_sizes =
1017 static_cast<MPI_Offset
>(global_first_cell) *
sizeof(
unsigned int);
1021 my_global_file_position_sizes,
1031 const std::uint64_t size_on_proc =
1036 std::uint64_t prefix_sum = 0;
1037 ierr = MPI_Exscan(&size_on_proc,
1045 const MPI_Offset my_global_file_position =
1046 static_cast<MPI_Offset
>(global_num_cells) *
sizeof(
unsigned int) +
1053 my_global_file_position,
1060 ierr = MPI_File_close(&fh);
1067 (void)global_first_cell;
1068 (void)global_num_cells;
1069 (void)mpi_communicator;
1075 const std::string fname_fixed =
1076 std::string(file_basename) +
"_fixed.data";
1078 std::ifstream file(fname_fixed, std::ios::binary | std::ios::in);
1102 const std::string fname_variable =
1103 std::string(file_basename) +
"_variable.data";
1105 std::ifstream file(fname_variable, std::ios::binary | std::ios::in);
1114 const std::uint64_t size =
1126 template <
int dim,
int spacedim>
1168 template <
int dim,
int spacedim>
1170 cell_is_patch_level_1(
1175 unsigned int n_active_children = 0;
1176 for (
unsigned int i = 0; i < cell->n_children(); ++i)
1177 if (cell->child(i)->is_active())
1178 ++n_active_children;
1180 return (n_active_children == 0) ||
1181 (n_active_children == cell->n_children());
1191 template <
int dim,
int spacedim>
1193 cell_will_be_coarsened(
1199 if (cell->has_children())
1201 unsigned int children_to_coarsen = 0;
1202 const unsigned int n_children = cell->n_children();
1204 for (
unsigned int c = 0; c < n_children; ++c)
1205 if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
1206 ++children_to_coarsen;
1207 if (children_to_coarsen == n_children)
1210 for (
unsigned int c = 0; c < n_children; ++c)
1211 if (cell->child(c)->is_active())
1212 cell->child(c)->clear_coarsen_flag();
1248 template <
int dim,
int spacedim>
1250 face_will_be_refined_by_neighbor_internal(
1252 const unsigned int face_no,
1261 cell->neighbor(face_no);
1273 if (cell_will_be_coarsened(neighbor))
1282 expected_face_ref_case = cell->face(face_no)->refinement_case();
1294 const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
1308 neighbor_face = neighbor->
face(neighbor_neighbor);
1309 const int this_face_index = cell->face_index(face_no);
1315 if (neighbor_face->index() == this_face_index)
1321 expected_face_ref_case = face_ref_case;
1339 for (
unsigned int c = 0; c < neighbor_face->n_children(); ++c)
1340 if (neighbor_face->child_index(c) == this_face_index)
1349 if ((neighbor_face->refinement_case() | face_ref_case) ==
1350 neighbor_face->refinement_case())
1369 expected_face_ref_case =
1370 ~neighbor_face->refinement_case();
1397 template <
int dim,
int spacedim>
1399 face_will_be_refined_by_neighbor(
1401 const unsigned int face_no)
1404 return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
1411 template <
int dim,
int spacedim>
1413 face_will_be_refined_by_neighbor(
1415 const unsigned int face_no,
1418 return face_will_be_refined_by_neighbor_internal(cell,
1420 expected_face_ref_case);
1425 template <
int dim,
int spacedim>
1427 satisfies_level1_at_vertex_rule(
1430 std::vector<unsigned int> min_adjacent_cell_level(
1432 std::vector<unsigned int> max_adjacent_cell_level(
1436 for (
const unsigned int v : cell->vertex_indices())
1438 min_adjacent_cell_level[cell->vertex_index(v)] =
1440 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
1441 max_adjacent_cell_level[cell->vertex_index(v)] =
1443 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
1446 for (
unsigned int k = 0; k < triangulation.
n_vertices(); ++k)
1448 if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
1472 template <
int dim,
int spacedim>
1474 middle_vertex_index(
1477 if (line->has_children())
1478 return line->child(0)->vertex_index(1);
1483 template <
int dim,
int spacedim>
1485 middle_vertex_index(
1488 switch (
static_cast<unsigned char>(quad->refinement_case()))
1491 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
1494 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
1497 return quad->child(0)->vertex_index(3);
1506 template <
int dim,
int spacedim>
1508 middle_vertex_index(
1511 switch (
static_cast<unsigned char>(hex->refinement_case()))
1514 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
1517 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
1520 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
1523 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
1526 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
1529 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
1532 return hex->child(0)->vertex_index(7);
1553 template <
class TRIANGULATION>
1554 inline typename TRIANGULATION::DistortedCellList
1555 collect_distorted_coarse_cells(
const TRIANGULATION &)
1557 return typename TRIANGULATION::DistortedCellList();
1579 vertices[i] = cell->vertex(i);
1585 if (determinants[i] <=
1593 return distorted_cells;
1605 has_distorted_children(
1610 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
1620 if (determinants[i] <=
1636 template <
int dim,
int spacedim>
1638 has_distorted_children(
1645 template <
int dim,
int spacedim>
1647 update_periodic_face_map_recursively(
1650 unsigned int n_face_1,
1651 unsigned int n_face_2,
1652 const unsigned char orientation,
1658 unsigned char>> &periodic_face_map)
1661 const FaceIterator face_1 = cell_1->
face(n_face_1);
1662 const FaceIterator face_2 = cell_2->
face(n_face_2);
1664 const unsigned char inverse_orientation =
1665 face_1->reference_cell().get_inverse_combined_orientation(orientation);
1668 const auto [face_orientation, face_rotation, face_flip] =
1671 Assert((dim != 1) || (face_orientation ==
true && face_flip ==
false &&
1672 face_rotation ==
false),
1674 "(face_orientation, face_flip, face_rotation) "
1675 "is invalid for 1d"));
1677 Assert((dim != 2) || (face_flip ==
false && face_rotation ==
false),
1679 "(face_orientation, face_flip, face_rotation) "
1680 "is invalid for 2d"));
1685 Assert(face_1->at_boundary() && face_2->at_boundary(),
1686 ExcMessage(
"Periodic faces must be on the boundary"));
1697 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
1699 const CellFace cell_face_1(cell_1, n_face_1);
1700 const CellFace cell_face_2(cell_2, n_face_2);
1701 const std::pair<CellFace, unsigned char> cell_face_orientation_2(
1702 cell_face_2, orientation);
1704 const std::pair<CellFace, std::pair<CellFace, unsigned char>>
1705 periodic_faces(cell_face_1, cell_face_orientation_2);
1709 periodic_face_map.insert(periodic_faces);
1717 update_periodic_face_map_recursively<dim, spacedim>(
1718 cell_1->
child(n_face_1),
1719 cell_2->
child(n_face_2),
1727 update_periodic_face_map_recursively<dim, spacedim>(
1728 cell_1->
child(n_face_1),
1747 Assert(face_1->n_children() ==
1749 face_2->n_children() ==
1755 for (
unsigned int i = 0;
1756 i < GeometryInfo<dim>::max_children_per_face;
1760 const unsigned int j =
1762 i, n_face_1, inverse_orientation);
1765 unsigned int child_cell_1 =
1773 face_1->refinement_case());
1774 unsigned int child_cell_2 =
1782 face_2->refinement_case());
1793 update_periodic_face_map_recursively<dim, spacedim>(
1794 cell_1->
child(child_cell_1),
1795 cell_2->
child(child_cell_2),
1804 for (
unsigned int i = 0;
1805 i < GeometryInfo<dim>::max_children_per_face;
1809 unsigned int child_cell_1 =
1817 face_1->refinement_case());
1820 update_periodic_face_map_recursively<dim, spacedim>(
1821 cell_1->
child(child_cell_1),
1855 <<
"Something went wrong upon construction of cell "
1869 <<
" has negative measure. This typically "
1870 <<
"indicates some distortion in the cell, or a mistakenly "
1871 <<
"swapped pair of vertices in the input to "
1872 <<
"Triangulation::create_triangulation().");
1884 <<
"Error while creating cell " << arg1
1885 <<
": the vertex index " << arg2 <<
" must be between 0 and "
1896 <<
"The input data for creating a triangulation contained "
1897 <<
"information about a line with indices " << arg1 <<
" and " << arg2
1898 <<
" that is described to have boundary indicator "
1899 <<
static_cast<int>(arg3)
1900 <<
". However, this is an internal line not located on the "
1901 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1903 <<
"If this happened at a place where you call "
1904 <<
"Triangulation::create_triangulation() yourself, you need "
1905 <<
"to check the SubCellData object you pass to this function."
1908 <<
"If this happened in a place where you are reading a mesh "
1909 <<
"from a file, then you need to investigate why such a line "
1910 <<
"ended up in the input file. A typical case is a geometry "
1911 <<
"that consisted of multiple parts and for which the mesh "
1912 <<
"generator program assumes that the interface between "
1913 <<
"two parts is a boundary when that isn't supposed to be "
1914 <<
"the case, or where the mesh generator simply assigns "
1915 <<
"'geometry indicators' to lines at the perimeter of "
1916 <<
"a part that are not supposed to be interpreted as "
1917 <<
"'boundary indicators'.");
1929 <<
"The input data for creating a triangulation contained "
1930 <<
"information about a quad with indices " << arg1 <<
", " << arg2
1931 <<
", " << arg3 <<
", and " << arg4
1932 <<
" that is described to have boundary indicator "
1933 <<
static_cast<int>(arg5)
1934 <<
". However, this is an internal quad not located on the "
1935 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1937 <<
"If this happened at a place where you call "
1938 <<
"Triangulation::create_triangulation() yourself, you need "
1939 <<
"to check the SubCellData object you pass to this function."
1942 <<
"If this happened in a place where you are reading a mesh "
1943 <<
"from a file, then you need to investigate why such a quad "
1944 <<
"ended up in the input file. A typical case is a geometry "
1945 <<
"that consisted of multiple parts and for which the mesh "
1946 <<
"generator program assumes that the interface between "
1947 <<
"two parts is a boundary when that isn't supposed to be "
1948 <<
"the case, or where the mesh generator simply assigns "
1949 <<
"'geometry indicators' to quads at the surface of "
1950 <<
"a part that are not supposed to be interpreted as "
1951 <<
"'boundary indicators'.");
1960 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1961 <<
" and " << arg2 <<
" appears more than once. "
1962 <<
"This is not allowed.");
1972 <<
"In SubCellData the line info of the line with vertex indices " << arg1
1973 <<
" and " << arg2 <<
" appears multiple times with different (valid) "
1974 << arg3 <<
". This is not allowed.");
1986 <<
"In SubCellData the quad info of the quad with line indices " << arg1
1987 <<
", " << arg2 <<
", " << arg3 <<
" and " << arg4
1988 <<
" appears multiple times with different (valid) " << arg5
1989 <<
". This is not allowed.");
2000 const unsigned int new_quads_in_pairs,
2001 const unsigned int new_quads_single)
2007 unsigned int next_free_single = 0;
2008 unsigned int next_free_pair = 0;
2012 unsigned int n_quads = 0;
2013 unsigned int n_unused_pairs = 0;
2014 unsigned int n_unused_singles = 0;
2015 for (
unsigned int i = 0; i < tria_faces.
quads.
used.size(); ++i)
2019 else if (i + 1 < tria_faces.
quads.
used.size())
2024 if (next_free_single == 0)
2025 next_free_single = i;
2030 if (next_free_pair == 0)
2038 Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
2044 const int additional_single_quads = new_quads_single - n_unused_singles;
2046 unsigned int new_size =
2047 tria_faces.
quads.
used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
2048 if (additional_single_quads > 0)
2049 new_size += additional_single_quads;
2059 q_is_q.reserve(new_size);
2060 q_is_q.insert(q_is_q.end(), new_size - q_is_q.size(),
true);
2081 const unsigned int total_cells,
2084 const bool tetraheder_in_mesh =
false)
2100 if (tetraheder_in_mesh)
2157 tria_level.
parents.reserve((total_cells + 1) / 2);
2159 (total_cells + 1) / 2 -
2167 std::make_pair(-1, -1));
2169 if (tria_level.
dim == 2 || tria_level.
dim == 3)
2171 const unsigned int max_faces_per_cell = 2 *
dimension;
2173 max_faces_per_cell);
2193 <<
"The containers have sizes " << arg1 <<
" and " << arg2
2194 <<
", which is not as expected.");
2203 const unsigned int true_dimension)
2206 (void)true_dimension;
2233 const unsigned int new_objects_in_pairs,
2234 const unsigned int new_objects_single = 0)
2246 unsigned int n_objects = 0;
2247 unsigned int n_unused_pairs = 0;
2248 unsigned int n_unused_singles = 0;
2249 for (
unsigned int i = 0; i < tria_objects.
used.size(); ++i)
2251 if (tria_objects.
used[i])
2253 else if (i + 1 < tria_objects.
used.size())
2255 if (tria_objects.
used[i + 1])
2272 Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
2273 tria_objects.
used.size(),
2279 const int additional_single_objects =
2280 new_objects_single - n_unused_singles;
2282 unsigned int new_size = tria_objects.
used.size() +
2283 new_objects_in_pairs - 2 * n_unused_pairs;
2284 if (additional_single_objects > 0)
2285 new_size += additional_single_objects;
2288 if (new_size > tria_objects.
n_objects())
2290 const unsigned int max_faces_per_cell =
2295 tria_objects.
cells.reserve(new_size * max_faces_per_cell);
2296 tria_objects.
cells.insert(tria_objects.
cells.end(),
2301 tria_objects.
used.reserve(new_size);
2302 tria_objects.
used.insert(tria_objects.
used.end(),
2303 new_size - tria_objects.
used.size(),
2313 tria_objects.
children.reserve(factor * new_size);
2333 tria_objects.
user_data.reserve(new_size);
2334 tria_objects.
user_data.resize(new_size);
2343 if (n_unused_singles == 0)
2351 const unsigned int new_hexes = new_objects_in_pairs;
2353 const unsigned int new_size =
2354 new_hexes + std::count(tria_objects.
used.begin(),
2355 tria_objects.
used.end(),
2359 if (new_size > tria_objects.
n_objects())
2361 const unsigned int max_faces_per_cell =
2364 tria_objects.
cells.reserve(new_size * max_faces_per_cell);
2365 tria_objects.
cells.insert(tria_objects.
cells.end(),
2370 tria_objects.
used.reserve(new_size);
2371 tria_objects.
used.insert(tria_objects.
used.end(),
2372 new_size - tria_objects.
used.size(),
2381 tria_objects.
children.reserve(4 * new_size);
2400 tria_objects.
user_data.reserve(new_size);
2401 tria_objects.
user_data.resize(new_size);
2425 tria_object.
used.size()));
2466 template <
int dim,
int spacedim>
2488 std::vector<unsigned int> &line_cell_count,
2489 std::vector<unsigned int> &quad_cell_count) = 0;
2496 const bool check_for_distorted_cells) = 0;
2525 virtual std::unique_ptr<Policy<dim, spacedim>>
2536 template <
int dim,
int spacedim,
typename T>
2543 T::update_neighbors(tria);
2550 std::vector<unsigned int> &line_cell_count,
2551 std::vector<unsigned int> &quad_cell_count)
override
2553 T::delete_children(tria, cell, line_cell_count, quad_cell_count);
2558 const bool check_for_distorted_cells)
override
2560 return T::execute_refinement(triangulation, check_for_distorted_cells);
2567 T::prevent_distorted_boundary_cells(triangulation);
2574 T::prepare_refinement_dim_dependent(triangulation);
2585 std::unique_ptr<Policy<dim, spacedim>>
2588 return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
2702 template <
int dim,
int spacedim>
2706 const unsigned int level_objects,
2709 using line_iterator =
2712 number_cache.n_levels = 0;
2713 if (level_objects > 0)
2715 for (
unsigned int level = 0; level < level_objects; ++level)
2716 if (triangulation.
begin(level) != triangulation.
end(level))
2717 number_cache.n_levels = level + 1;
2725 number_cache.n_lines = 0;
2726 number_cache.n_active_lines = 0;
2732 number_cache.n_lines_level.resize(number_cache.n_levels);
2733 number_cache.n_active_lines_level.resize(number_cache.n_levels);
2735 for (
unsigned int level = 0; level < number_cache.n_levels; ++level)
2738 number_cache.n_lines_level[level] = 0;
2739 number_cache.n_active_lines_level[level] = 0;
2741 line_iterator line = triangulation.
begin_line(level),
2743 (level == number_cache.n_levels - 1 ?
2744 line_iterator(triangulation.
end_line()) :
2746 for (; line != endc; ++line)
2748 ++number_cache.n_lines_level[level];
2749 if (line->has_children() ==
false)
2750 ++number_cache.n_active_lines_level[level];
2754 number_cache.n_lines += number_cache.n_lines_level[level];
2755 number_cache.n_active_lines +=
2756 number_cache.n_active_lines_level[level];
2762 number_cache.n_lines_level.clear();
2763 number_cache.n_active_lines_level.clear();
2765 line_iterator line = triangulation.
begin_line(),
2767 for (; line != endc; ++line)
2769 ++number_cache.n_lines;
2770 if (line->has_children() ==
false)
2771 ++number_cache.n_active_lines;
2790 template <
int dim,
int spacedim>
2794 const unsigned int level_objects,
2811 using quad_iterator =
2817 number_cache.n_quads = 0;
2818 number_cache.n_active_quads = 0;
2828 unsigned int n_levels = 0;
2829 if (level_objects > 0)
2831 for (
unsigned int level = 0; level < level_objects; ++level)
2832 if (triangulation.
begin(level) != triangulation.
end(level))
2833 n_levels = level + 1;
2835 number_cache.n_quads_level.resize(n_levels);
2836 number_cache.n_active_quads_level.resize(n_levels);
2838 for (
unsigned int level = 0; level < n_levels; ++level)
2841 number_cache.n_quads_level[level] = 0;
2842 number_cache.n_active_quads_level[level] = 0;
2844 quad_iterator quad = triangulation.
begin_quad(level),
2846 (level == n_levels - 1 ?
2847 quad_iterator(triangulation.
end_quad()) :
2849 for (; quad != endc; ++quad)
2851 ++number_cache.n_quads_level[level];
2852 if (quad->has_children() ==
false)
2853 ++number_cache.n_active_quads_level[level];
2857 number_cache.n_quads += number_cache.n_quads_level[level];
2858 number_cache.n_active_quads +=
2859 number_cache.n_active_quads_level[level];
2865 number_cache.n_quads_level.clear();
2866 number_cache.n_active_quads_level.clear();
2868 quad_iterator quad = triangulation.
begin_quad(),
2870 for (; quad != endc; ++quad)
2872 ++number_cache.n_quads;
2873 if (quad->has_children() ==
false)
2874 ++number_cache.n_active_quads;
2879 update_lines.
join();
2897 template <
int dim,
int spacedim>
2901 const unsigned int level_objects,
2918 using hex_iterator =
2924 number_cache.n_hexes = 0;
2925 number_cache.n_active_hexes = 0;
2936 unsigned int n_levels = 0;
2937 if (level_objects > 0)
2939 for (
unsigned int level = 0; level < level_objects; ++level)
2940 if (triangulation.
begin(level) != triangulation.
end(level))
2941 n_levels = level + 1;
2943 number_cache.n_hexes_level.resize(n_levels);
2944 number_cache.n_active_hexes_level.resize(n_levels);
2946 for (
unsigned int level = 0; level < n_levels; ++level)
2949 number_cache.n_hexes_level[level] = 0;
2950 number_cache.n_active_hexes_level[level] = 0;
2952 hex_iterator hex = triangulation.
begin_hex(level),
2953 endc = (level == n_levels - 1 ?
2954 hex_iterator(triangulation.
end_hex()) :
2956 for (; hex != endc; ++hex)
2958 ++number_cache.n_hexes_level[level];
2959 if (hex->has_children() ==
false)
2960 ++number_cache.n_active_hexes_level[level];
2964 number_cache.n_hexes += number_cache.n_hexes_level[level];
2965 number_cache.n_active_hexes +=
2966 number_cache.n_active_hexes_level[level];
2972 number_cache.n_hexes_level.clear();
2973 number_cache.n_active_hexes_level.clear();
2975 hex_iterator hex = triangulation.
begin_hex(),
2976 endc = triangulation.
end_hex();
2977 for (; hex != endc; ++hex)
2979 ++number_cache.n_hexes;
2980 if (hex->has_children() ==
false)
2981 ++number_cache.n_active_hexes;
2986 update_quads_and_lines.
join();
2990 template <
int dim,
int spacedim>
2994 const unsigned int level_objects,
2999 number_cache.active_cell_index_partitioner =
3000 std::make_shared<const Utilities::MPI::Partitioner>(
3003 number_cache.level_cell_index_partitioners.resize(
3005 for (
unsigned int level = 0; level < triangulation.
n_levels(); ++level)
3006 number_cache.level_cell_index_partitioners[level] =
3007 std::make_shared<const Utilities::MPI::Partitioner>(
3008 triangulation.
n_cells(level));
3012 template <
int spacedim>
3018 template <
int dim,
int spacedim>
3073 static const unsigned int left_right_offset[2][6][2] = {
3082 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
3093 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
3094 adjacent_cells(2 * triangulation.
n_raw_faces(), dummy);
3102 const unsigned int offset =
3108 adjacent_cells[2 * face->index() + offset] = cell;
3118 if (cell->
is_active() && face->has_children())
3120 adjacent_cells[2 * face->child(0)->index() + offset] =
3122 adjacent_cells[2 * face->child(1)->index() + offset] =
3148 if (face->has_children() &&
3154 for (
unsigned int c = 0; c < face->n_children(); ++c)
3155 adjacent_cells[2 * face->child(c)->index() + offset] =
3157 if (face->child(0)->has_children())
3159 adjacent_cells[2 * face->child(0)->child(0)->index() +
3161 adjacent_cells[2 * face->child(0)->child(1)->index() +
3164 if (face->child(1)->has_children())
3166 adjacent_cells[2 * face->child(1)->child(0)->index() +
3168 adjacent_cells[2 * face->child(1)->child(1)->index() +
3183 const unsigned int offset =
3189 f, adjacent_cells[2 * cell->
face(f)->index() + 1 - offset]);
3197 template <
int dim,
int spacedim>
3211 if (dim == spacedim)
3212 for (
unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
3235 tria.
faces = std::make_unique<
3244 const unsigned int n_cell = cells.size();
3249 auto &lines_0 = tria.
faces->lines;
3252 const auto &crs = connectivity.entity_to_entities(1, 0);
3253 const unsigned int n_lines = crs.ptr.size() - 1;
3259 for (
unsigned int line = 0; line < n_lines; ++line)
3260 for (
unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
3269 auto &quads_0 = tria.
faces->quads;
3270 auto &faces = *tria.
faces;
3273 const auto &crs = connectivity.entity_to_entities(2, 1);
3274 const unsigned int n_quads = crs.ptr.size() - 1;
3281 for (
unsigned int q = 0, k = 0; q < n_quads; ++q)
3284 faces.set_quad_type(q, connectivity.entity_types(2)[q]);
3287 for (
unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
3295 const unsigned char combined_orientation =
3296 connectivity.entity_orientations(1)
3297 .get_combined_orientation(k);
3301 combined_orientation ==
3303 combined_orientation ==
3306 faces.quads_line_orientations
3308 combined_orientation ==
3316 auto &cells_0 = tria.
levels[0]->cells;
3317 auto &level = *tria.
levels[0];
3320 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
3321 const auto &nei = connectivity.entity_to_entities(dim, dim);
3325 bool orientation_needed =
false;
3327 orientation_needed =
true;
3330 const auto &orientations = connectivity.entity_orientations(1);
3331 for (
unsigned int i = 0; i < orientations.n_objects(); ++i)
3332 if (orientations.get_combined_orientation(i) !=
3335 orientation_needed =
true;
3345 for (
unsigned int cell = 0; cell < n_cell; ++cell)
3348 cells_0.boundary_or_material_id[cell].material_id =
3349 cells[cell].material_id;
3352 cells_0.manifold_id[cell] = cells[cell].manifold_id;
3355 level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
3358 for (
unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
3362 if (nei.col[i] !=
static_cast<unsigned int>(-1))
3364 j] = {0, nei.col[i]};
3371 if (orientation_needed)
3373 level.face_orientations.set_combined_orientation(
3375 connectivity.entity_orientations(dim - 1)
3376 .get_combined_orientation(i));
3385 auto &bids_face = dim == 3 ?
3386 tria.
faces->quads.boundary_or_material_id :
3387 tria.
faces->lines.boundary_or_material_id;
3390 std::vector<unsigned int> count(bids_face.size(), 0);
3393 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
3396 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
3397 for (
unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
3398 count[crs.col[i]]++;
3401 for (
unsigned int face = 0; face < count.size(); ++face)
3403 if (count[face] != 1)
3407 bids_face[face].boundary_id = 0;
3413 const auto &crs = connectivity.entity_to_entities(2, 1);
3414 for (
unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
3415 tria.
faces->lines.boundary_or_material_id[crs.col[i]]
3421 static const unsigned int t_tba =
static_cast<unsigned int>(-1);
3422 static const unsigned int t_inner =
static_cast<unsigned int>(-2);
3424 std::vector<unsigned int> type(vertices.size(), t_tba);
3426 const auto &crs = connectivity.entity_to_entities(1, 0);
3428 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
3429 for (
unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
3431 if (type[crs.col[i]] != t_inner)
3432 type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
3434 for (
unsigned int face = 0; face < type.size(); ++face)
3439 if (type[face] != t_inner && type[face] != t_tba)
3460 template <
int structdim,
int spacedim,
typename T>
3470 if (boundary_objects_in.empty())
3474 auto boundary_objects = boundary_objects_in;
3477 for (
auto &boundary_object : boundary_objects)
3478 std::sort(boundary_object.vertices.begin(),
3479 boundary_object.vertices.end());
3482 std::sort(boundary_objects.begin(),
3483 boundary_objects.end(),
3484 [](
const auto &a,
const auto &b) {
3485 return a.vertices < b.vertices;
3488 unsigned int counter = 0;
3490 std::vector<unsigned int> key;
3493 for (
unsigned int o = 0; o < obj.
n_objects(); ++o)
3507 key.assign(crs.
col.data() + crs.
ptr[o],
3508 crs.
col.data() + crs.
ptr[o + 1]);
3509 std::sort(key.begin(), key.end());
3512 const auto subcell_object =
3513 std::lower_bound(boundary_objects.begin(),
3514 boundary_objects.end(),
3516 [&](
const auto &cell,
const auto &key) {
3517 return cell.vertices < key;
3521 if (subcell_object == boundary_objects.end() ||
3522 subcell_object->vertices != key)
3528 manifold_id = subcell_object->manifold_id;
3531 if (subcell_object->boundary_id !=
3534 (void)vertex_locations;
3538 "The input arguments for creating a triangulation "
3539 "specified a boundary id for an internal face. This "
3542 "The object in question has vertex indices " +
3543 [subcell_object]() {
3545 for (
const auto v : subcell_object->vertices)
3546 s += std::to_string(v) +
',';
3549 " which are located at positions " +
3550 [vertex_locations, subcell_object]() {
3551 std::ostringstream s;
3552 for (
const auto v : subcell_object->vertices)
3553 s <<
'(' << vertex_locations[v] <<
')';
3557 boundary_id = subcell_object->boundary_id;
3571 const unsigned structdim,
3572 const unsigned int size)
3574 const unsigned int dim = faces.
dim;
3576 const unsigned int max_lines_per_face = 2 * structdim;
3578 if (dim == 3 && structdim == 2)
3593 const unsigned int spacedim,
3594 const unsigned int size,
3595 const bool orientation_needed)
3597 const unsigned int dim = level.
dim;
3599 const unsigned int max_faces_per_cell = 2 * dim;
3609 level.
parents.assign((size + 1) / 2, -1);
3611 if (dim == spacedim - 1)
3614 level.
neighbors.assign(size * max_faces_per_cell, {-1, -1});
3618 if (orientation_needed)
3633 const unsigned int structdim = obj.
structdim;
3636 const unsigned int max_faces_per_cell = 2 * structdim;
3638 obj.
used.assign(size,
true);
3642 BoundaryOrMaterialId());
3652 obj.
cells.assign(max_faces_per_cell * size, -1);
3682 template <
int spacedim>
3686 std::vector<unsigned int> &,
3687 std::vector<unsigned int> &)
3689 const unsigned int dim = 1;
3715 neighbor = neighbor->
child(1);
3727 neighbor = neighbor->
child(1);
3743 neighbor = neighbor->
child(0);
3755 neighbor = neighbor->
child(0);
3771 for (
unsigned int child = 0; child < cell->
n_children(); ++child)
3786 template <
int spacedim>
3790 std::vector<unsigned int> &line_cell_count,
3791 std::vector<unsigned int> &)
3793 const unsigned int dim = 2;
3801 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3804 lines_to_delete.reserve(4 * 2 + 4);
3809 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
3813 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3831 lines_to_delete.push_back(cell->
child(0)->
line(1));
3832 lines_to_delete.push_back(cell->
child(0)->
line(3));
3833 lines_to_delete.push_back(cell->
child(3)->
line(0));
3834 lines_to_delete.push_back(cell->
child(3)->
line(2));
3838 unsigned int inner_face_no =
3843 lines_to_delete.push_back(cell->
child(0)->
line(inner_face_no));
3847 for (
unsigned int child = 0; child < cell->
n_children(); ++child)
3864 for (
unsigned int line_no = 0;
3865 line_no < GeometryInfo<dim>::lines_per_cell;
3869 cell->
line(line_no);
3871 if (line->has_children())
3876 Assert((line_cell_count[line->child_index(0)] == 0 &&
3877 line_cell_count[line->child_index(1)] == 0) ||
3878 (line_cell_count[line->child_index(0)] > 0 &&
3879 line_cell_count[line->child_index(1)] > 0),
3882 if (line_cell_count[line->child_index(0)] == 0)
3884 for (
unsigned int c = 0; c < 2; ++c)
3885 Assert(!line->child(c)->has_children(),
3895 lines_to_delete.push_back(line->child(0));
3896 lines_to_delete.push_back(line->child(1));
3898 line->clear_children();
3910 typename std::vector<
3912 line = lines_to_delete.begin(),
3913 endline = lines_to_delete.end();
3914 for (; line != endline; ++line)
3916 (*line)->clear_user_data();
3917 (*line)->clear_user_flag();
3918 (*line)->clear_used_flag();
3924 template <
int spacedim>
3928 std::vector<unsigned int> &line_cell_count,
3929 std::vector<unsigned int> &quad_cell_count)
3931 const unsigned int dim = 3;
3943 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3945 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3948 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3949 quads_to_delete.reserve(6 * 4 + 12);
3953 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
3957 const auto line_indices = TriaAccessorImplementation::
3958 Implementation::get_line_indices_of_cell(*child);
3960 --line_cell_count[line_indices[l]];
3976 quads_to_delete.push_back(cell->
child(0)->
face(1));
3979 quads_to_delete.push_back(cell->
child(0)->
face(3));
3982 quads_to_delete.push_back(cell->
child(0)->
face(5));
3985 quads_to_delete.push_back(cell->
child(0)->
face(1));
3986 quads_to_delete.push_back(cell->
child(0)->
face(3));
3987 quads_to_delete.push_back(cell->
child(3)->
face(0));
3988 quads_to_delete.push_back(cell->
child(3)->
face(2));
3990 lines_to_delete.push_back(cell->
child(0)->
line(11));
3993 quads_to_delete.push_back(cell->
child(0)->
face(1));
3994 quads_to_delete.push_back(cell->
child(0)->
face(5));
3995 quads_to_delete.push_back(cell->
child(3)->
face(0));
3996 quads_to_delete.push_back(cell->
child(3)->
face(4));
3998 lines_to_delete.push_back(cell->
child(0)->
line(5));
4001 quads_to_delete.push_back(cell->
child(0)->
face(3));
4002 quads_to_delete.push_back(cell->
child(0)->
face(5));
4003 quads_to_delete.push_back(cell->
child(3)->
face(2));
4004 quads_to_delete.push_back(cell->
child(3)->
face(4));
4006 lines_to_delete.push_back(cell->
child(0)->
line(7));
4009 quads_to_delete.push_back(cell->
child(0)->
face(1));
4010 quads_to_delete.push_back(cell->
child(2)->
face(1));
4011 quads_to_delete.push_back(cell->
child(4)->
face(1));
4012 quads_to_delete.push_back(cell->
child(6)->
face(1));
4014 quads_to_delete.push_back(cell->
child(0)->
face(3));
4015 quads_to_delete.push_back(cell->
child(1)->
face(3));
4016 quads_to_delete.push_back(cell->
child(4)->
face(3));
4017 quads_to_delete.push_back(cell->
child(5)->
face(3));
4019 quads_to_delete.push_back(cell->
child(0)->
face(5));
4020 quads_to_delete.push_back(cell->
child(1)->
face(5));
4021 quads_to_delete.push_back(cell->
child(2)->
face(5));
4022 quads_to_delete.push_back(cell->
child(3)->
face(5));
4024 lines_to_delete.push_back(cell->
child(0)->
line(5));
4025 lines_to_delete.push_back(cell->
child(0)->
line(7));
4026 lines_to_delete.push_back(cell->
child(0)->
line(11));
4027 lines_to_delete.push_back(cell->
child(7)->
line(0));
4028 lines_to_delete.push_back(cell->
child(7)->
line(2));
4029 lines_to_delete.push_back(cell->
child(7)->
line(8));
4047 for (
unsigned int child = 0; child < cell->
n_children(); ++child)
4078 cell->
face(quad_no);
4082 quad->has_children()) ||
4087 switch (quad->refinement_case())
4099 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
4100 quad_cell_count[quad->child_index(1)] == 0) ||
4101 (quad_cell_count[quad->child_index(0)] > 0 &&
4102 quad_cell_count[quad->child_index(1)] > 0),
4108 unsigned int deleted_grandchildren = 0;
4109 unsigned int number_of_child_refinements = 0;
4111 for (
unsigned int c = 0; c < 2; ++c)
4112 if (quad->child(c)->has_children())
4114 ++number_of_child_refinements;
4119 (quad_cell_count[quad->child(c)->child_index(0)] ==
4121 quad_cell_count[quad->child(c)->child_index(1)] ==
4123 (quad_cell_count[quad->child(c)->child_index(0)] >
4125 quad_cell_count[quad->child(c)->child_index(1)] >
4128 if (quad_cell_count[quad->child(c)->child_index(0)] ==
4135 Assert(quad->refinement_case() +
4136 quad->child(c)->refinement_case() ==
4144 quads_to_delete.push_back(
4145 quad->child(c)->child(0));
4146 quads_to_delete.push_back(
4147 quad->child(c)->child(1));
4148 if (quad->child(c)->refinement_case() ==
4150 lines_to_delete.push_back(
4151 quad->child(c)->child(0)->line(1));
4153 lines_to_delete.push_back(
4154 quad->child(c)->child(0)->line(3));
4155 quad->child(c)->clear_children();
4156 quad->child(c)->clear_refinement_case();
4157 ++deleted_grandchildren;
4165 if (number_of_child_refinements > 0 &&
4166 deleted_grandchildren == number_of_child_refinements)
4171 middle_line = quad->child(0)->line(1);
4173 middle_line = quad->child(0)->line(3);
4175 lines_to_delete.push_back(middle_line->child(0));
4176 lines_to_delete.push_back(middle_line->child(1));
4179 middle_line)] =
false;
4180 middle_line->clear_children();
4185 if (quad_cell_count[quad->child_index(0)] == 0)
4191 quads_to_delete.push_back(quad->child(0));
4192 quads_to_delete.push_back(quad->child(1));
4194 lines_to_delete.push_back(quad->child(0)->line(1));
4196 lines_to_delete.push_back(quad->child(0)->line(3));
4216 if (quad->child(0)->has_children())
4217 if (quad->refinement_case() ==
4251 quad_iterator switch_1 =
4252 quad->child(0)->child(1),
4254 quad->child(1)->child(0);
4256 Assert(!switch_1->has_children(),
4258 Assert(!switch_2->has_children(),
4261 const int switch_1_index = switch_1->index();
4262 const int switch_2_index = switch_2->index();
4263 for (
unsigned int l = 0;
4264 l < triangulation.
levels.size();
4266 for (
unsigned int h = 0;
4268 triangulation.
levels[l]->cells.n_objects();
4270 for (
const unsigned int q :
4275 ->cells.get_bounding_object_indices(
4277 if (
index == switch_1_index)
4279 ->cells.get_bounding_object_indices(
4280 h)[q] = switch_2_index;
4281 else if (
index == switch_2_index)
4283 ->cells.get_bounding_object_indices(
4284 h)[q] = switch_1_index;
4289 const int switch_1_lines[4] = {
4290 static_cast<signed int>(
4291 switch_1->line_index(0)),
4292 static_cast<signed int>(
4293 switch_1->line_index(1)),
4294 static_cast<signed int>(
4295 switch_1->line_index(2)),
4296 static_cast<signed int>(
4297 switch_1->line_index(3))};
4298 const bool switch_1_line_orientations[4] = {
4299 switch_1->line_orientation(0),
4300 switch_1->line_orientation(1),
4301 switch_1->line_orientation(2),
4302 switch_1->line_orientation(3)};
4304 switch_1->boundary_id();
4305 const unsigned int switch_1_user_index =
4306 switch_1->user_index();
4307 const bool switch_1_user_flag =
4308 switch_1->user_flag_set();
4310 switch_1->set_bounding_object_indices(
4311 {switch_2->line_index(0),
4312 switch_2->line_index(1),
4313 switch_2->line_index(2),
4314 switch_2->line_index(3)});
4315 switch_1->set_line_orientation(
4316 0, switch_2->line_orientation(0));
4317 switch_1->set_line_orientation(
4318 1, switch_2->line_orientation(1));
4319 switch_1->set_line_orientation(
4320 2, switch_2->line_orientation(2));
4321 switch_1->set_line_orientation(
4322 3, switch_2->line_orientation(3));
4323 switch_1->set_boundary_id_internal(
4324 switch_2->boundary_id());
4325 switch_1->set_manifold_id(
4326 switch_2->manifold_id());
4327 switch_1->set_user_index(switch_2->user_index());
4328 if (switch_2->user_flag_set())
4329 switch_1->set_user_flag();
4331 switch_1->clear_user_flag();
4333 switch_2->set_bounding_object_indices(
4337 switch_1_lines[3]});
4338 switch_2->set_line_orientation(
4339 0, switch_1_line_orientations[0]);
4340 switch_2->set_line_orientation(
4341 1, switch_1_line_orientations[1]);
4342 switch_2->set_line_orientation(
4343 2, switch_1_line_orientations[2]);
4344 switch_2->set_line_orientation(
4345 3, switch_1_line_orientations[3]);
4346 switch_2->set_boundary_id_internal(
4347 switch_1_boundary_id);
4348 switch_2->set_manifold_id(
4349 switch_1->manifold_id());
4350 switch_2->set_user_index(switch_1_user_index);
4351 if (switch_1_user_flag)
4352 switch_2->set_user_flag();
4354 switch_2->clear_user_flag();
4356 const unsigned int child_0 =
4357 quad->child(0)->child_index(0);
4358 const unsigned int child_2 =
4359 quad->child(1)->child_index(0);
4360 quad->clear_children();
4361 quad->clear_refinement_case();
4362 quad->set_refinement_case(
4364 quad->set_children(0, child_0);
4365 quad->set_children(2, child_2);
4366 std::swap(quad_cell_count[child_0 + 1],
4367 quad_cell_count[child_2]);
4382 const unsigned int child_0 =
4383 quad->child(0)->child_index(0);
4384 const unsigned int child_2 =
4385 quad->child(1)->child_index(0);
4386 quad->clear_children();
4387 quad->clear_refinement_case();
4388 quad->set_refinement_case(
4390 quad->set_children(0, child_0);
4391 quad->set_children(2, child_2);
4395 quad->clear_children();
4396 quad->clear_refinement_case();
4407 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
4408 quad_cell_count[quad->child_index(1)] == 0 &&
4409 quad_cell_count[quad->child_index(2)] == 0 &&
4410 quad_cell_count[quad->child_index(3)] == 0) ||
4411 (quad_cell_count[quad->child_index(0)] > 0 &&
4412 quad_cell_count[quad->child_index(1)] > 0 &&
4413 quad_cell_count[quad->child_index(2)] > 0 &&
4414 quad_cell_count[quad->child_index(3)] > 0),
4417 if (quad_cell_count[quad->child_index(0)] == 0)
4423 lines_to_delete.push_back(quad->child(0)->line(1));
4424 lines_to_delete.push_back(quad->child(3)->line(0));
4425 lines_to_delete.push_back(quad->child(0)->line(3));
4426 lines_to_delete.push_back(quad->child(3)->line(2));
4428 for (
unsigned int child = 0; child < quad->n_children();
4430 quads_to_delete.push_back(quad->child(child));
4436 quad->clear_children();
4437 quad->clear_refinement_case();
4456 for (
unsigned int line_no = 0;
4457 line_no < GeometryInfo<dim>::lines_per_cell;
4461 cell->
line(line_no);
4465 line->has_children()) ||
4470 if (line->has_children())
4475 Assert((line_cell_count[line->child_index(0)] == 0 &&
4476 line_cell_count[line->child_index(1)] == 0) ||
4477 (line_cell_count[line->child_index(0)] > 0 &&
4478 line_cell_count[line->child_index(1)] > 0),
4481 if (line_cell_count[line->child_index(0)] == 0)
4483 for (
unsigned int c = 0; c < 2; ++c)
4484 Assert(!line->child(c)->has_children(),
4494 lines_to_delete.push_back(line->child(0));
4495 lines_to_delete.push_back(line->child(1));
4497 line->clear_children();
4509 typename std::vector<
4511 line = lines_to_delete.begin(),
4512 endline = lines_to_delete.end();
4513 for (; line != endline; ++line)
4515 (*line)->clear_user_data();
4516 (*line)->clear_user_flag();
4517 (*line)->clear_used_flag();
4520 typename std::vector<
4522 quad = quads_to_delete.begin(),
4523 endquad = quads_to_delete.end();
4524 for (; quad != endquad; ++quad)
4526 (*quad)->clear_user_data();
4527 (*quad)->clear_children();
4528 (*quad)->clear_refinement_case();
4529 (*quad)->clear_user_flag();
4530 (*quad)->clear_used_flag();
4552 template <
int spacedim>
4556 unsigned int &next_unused_vertex,
4563 const unsigned int dim = 2;
4659 int new_vertices[9];
4660 for (
unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
4661 new_vertices[vertex_no] = cell->
vertex_index(vertex_no);
4662 for (
unsigned int line_no = 0; line_no < 4; ++line_no)
4663 if (cell->
line(line_no)->has_children())
4664 new_vertices[4 + line_no] =
4665 cell->
line(line_no)->child(0)->vertex_index(1);
4674 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
4675 ++next_unused_vertex;
4678 "Internal error: During refinement, the triangulation "
4679 "wants to access an element of the 'vertices' array "
4680 "but it turns out that the array is not large enough."));
4683 new_vertices[8] = next_unused_vertex;
4688 triangulation.
vertices[next_unused_vertex] =
4689 cell->
center(
true,
true);
4695 unsigned int lmin = 8;
4696 unsigned int lmax = 12;
4703 for (
unsigned int l = lmin; l < lmax; ++l)
4705 while (next_unused_line->used() ==
true)
4707 new_lines[l] = next_unused_line;
4725 for (
unsigned int c = 0; c < 2; ++c, ++l)
4726 new_lines[l] = cell->
line(face_no)->child(c);
4729 new_lines[8]->set_bounding_object_indices(
4730 {new_vertices[6], new_vertices[8]});
4731 new_lines[9]->set_bounding_object_indices(
4732 {new_vertices[8], new_vertices[7]});
4733 new_lines[10]->set_bounding_object_indices(
4734 {new_vertices[4], new_vertices[8]});
4735 new_lines[11]->set_bounding_object_indices(
4736 {new_vertices[8], new_vertices[5]});
4745 new_lines[0] = cell->
line(0);
4746 new_lines[1] = cell->
line(1);
4747 new_lines[2] = cell->
line(2)->child(0);
4748 new_lines[3] = cell->
line(2)->child(1);
4749 new_lines[4] = cell->
line(3)->child(0);
4750 new_lines[5] = cell->
line(3)->child(1);
4751 new_lines[6]->set_bounding_object_indices(
4752 {new_vertices[6], new_vertices[7]});
4762 new_lines[0] = cell->
line(0)->child(0);
4763 new_lines[1] = cell->
line(0)->child(1);
4764 new_lines[2] = cell->
line(1)->child(0);
4765 new_lines[3] = cell->
line(1)->child(1);
4766 new_lines[4] = cell->
line(2);
4767 new_lines[5] = cell->
line(3);
4768 new_lines[6]->set_bounding_object_indices(
4769 {new_vertices[4], new_vertices[5]});
4772 for (
unsigned int l = lmin; l < lmax; ++l)
4774 new_lines[l]->set_used_flag();
4775 new_lines[l]->clear_user_flag();
4777 new_lines[l]->clear_children();
4779 new_lines[l]->set_boundary_id_internal(
4781 new_lines[l]->set_manifold_id(cell->
manifold_id());
4788 while (next_unused_cell->
used() ==
true)
4792 for (
unsigned int i = 0; i < n_children; ++i)
4795 subcells[i] = next_unused_cell;
4797 if (i % 2 == 1 && i < n_children - 1)
4798 while (next_unused_cell->
used() ==
true)
4817 new_lines[8]->index(),
4818 new_lines[4]->index(),
4819 new_lines[10]->index()});
4821 new_lines[2]->index(),
4822 new_lines[5]->index(),
4823 new_lines[11]->index()});
4825 new_lines[9]->index(),
4826 new_lines[10]->index(),
4827 new_lines[6]->index()});
4829 new_lines[3]->index(),
4830 new_lines[11]->index(),
4831 new_lines[7]->index()});
4848 new_lines[6]->index(),
4849 new_lines[2]->index(),
4850 new_lines[4]->index()});
4852 new_lines[1]->index(),
4853 new_lines[3]->index(),
4854 new_lines[5]->index()});
4872 new_lines[2]->index(),
4873 new_lines[4]->index(),
4874 new_lines[6]->index()});
4876 new_lines[3]->index(),
4877 new_lines[6]->index(),
4878 new_lines[5]->index()});
4883 for (
unsigned int i = 0; i < n_children; ++i)
4902 for (
unsigned int i = 0; i < n_children / 2; ++i)
4912 if (dim == spacedim - 1)
4913 for (
unsigned int c = 0; c < n_children; ++c)
4919 template <
int dim,
int spacedim>
4922 const bool check_for_distorted_cells)
4929 triangulation.
levels.size() - 1))
4932 triangulation.
levels.push_back(
4943 line->clear_user_flag();
4947 unsigned int n_single_lines = 0;
4948 unsigned int n_lines_in_pairs = 0;
4949 unsigned int needed_vertices = 0;
4951 for (
int level = triangulation.
levels.size() - 2; level >= 0; --level)
4955 unsigned int needed_cells = 0;
4957 for (
const auto &cell :
4964 needed_vertices += 0;
4965 n_single_lines += 3;
4971 needed_vertices += 1;
4972 n_single_lines += 4;
4981 auto line = cell->
line(line_no);
4982 if (line->has_children() ==
false)
4983 line->set_user_flag();
4988 const unsigned int used_cells =
4989 std::count(triangulation.
levels[level + 1]->cells.used.begin(),
4990 triangulation.
levels[level + 1]->cells.used.end(),
4995 used_cells + needed_cells,
5007 if (line->user_flag_set())
5010 n_lines_in_pairs += 2;
5011 needed_vertices += 1;
5016 needed_vertices += std::count(triangulation.
vertices_used.begin(),
5020 if (needed_vertices > triangulation.
vertices.size())
5026 unsigned int next_unused_vertex = 0;
5035 for (; line != endl; ++line)
5036 if (line->user_flag_set())
5040 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5041 ++next_unused_vertex;
5044 "Internal error: During refinement, the triangulation "
5045 "wants to access an element of the 'vertices' array "
5046 "but it turns out that the array is not large "
5050 triangulation.
vertices[next_unused_vertex] = line->center(
true);
5052 bool pair_found =
false;
5054 for (; next_unused_line != endl; ++next_unused_line)
5055 if (!next_unused_line->used() &&
5056 !(++next_unused_line)->used())
5064 line->set_children(0, next_unused_line->index());
5067 children[2] = {next_unused_line, ++next_unused_line};
5072 children[0]->set_bounding_object_indices(
5073 {line->vertex_index(0), next_unused_vertex});
5074 children[1]->set_bounding_object_indices(
5075 {next_unused_vertex, line->vertex_index(1)});
5077 for (
auto &child : children)
5079 child->set_used_flag();
5080 child->clear_children();
5081 child->clear_user_data();
5082 child->clear_user_flag();
5083 child->set_boundary_id_internal(line->boundary_id());
5084 child->set_manifold_id(line->manifold_id());
5089 line->clear_user_flag();
5096 cells_with_distorted_children;
5102 unsigned int &next_unused_vertex,
5103 auto &next_unused_line,
5104 auto &next_unused_cell,
5109 unsigned int n_new_vertices = 0;
5118 std::vector<unsigned int> new_vertices(n_new_vertices,
5120 for (
unsigned int vertex_no = 0; vertex_no < cell->
n_vertices();
5122 new_vertices[vertex_no] = cell->
vertex_index(vertex_no);
5123 for (
unsigned int line_no = 0; line_no < cell->
n_lines(); ++line_no)
5124 if (cell->
line(line_no)->has_children())
5126 cell->
line(line_no)->child(0)->vertex_index(1);
5130 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5131 ++next_unused_vertex;
5133 next_unused_vertex < triangulation.
vertices.size(),
5135 "Internal error: During refinement, the triangulation wants "
5136 "to access an element of the 'vertices' array but it turns "
5137 "out that the array is not large enough."));
5140 new_vertices[8] = next_unused_vertex;
5142 triangulation.
vertices[next_unused_vertex] =
5143 cell->
center(
true,
true);
5146 std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
5149 std::array<unsigned char, 12> inherited_orientations;
5150 inherited_orientations.fill(
5152 unsigned int lmin = 0;
5153 unsigned int lmax = 0;
5162 std::fill(inherited_orientations.begin() + lmin,
5163 inherited_orientations.begin() + lmax,
5176 for (
unsigned int l = lmin; l < lmax; ++l)
5178 while (next_unused_line->used() ==
true)
5180 new_lines[l] = next_unused_line;
5187 for (
const unsigned int face_no : cell->
face_indices())
5194 const unsigned char combined_orientation =
5196 Assert(combined_orientation ==
5198 combined_orientation ==
5201 for (
unsigned int c = 0; c < 2; ++c)
5203 new_lines[2 * face_no + c] = cell->
line(face_no)->child(c);
5204 inherited_orientations[2 * face_no + c] =
5207 if (combined_orientation ==
5209 std::swap(new_lines[2 * face_no], new_lines[2 * face_no + 1]);
5215 new_lines[6]->set_bounding_object_indices(
5216 {new_vertices[3], new_vertices[4]});
5217 new_lines[7]->set_bounding_object_indices(
5218 {new_vertices[4], new_vertices[5]});
5219 new_lines[8]->set_bounding_object_indices(
5220 {new_vertices[5], new_vertices[3]});
5224 new_lines[8]->set_bounding_object_indices(
5225 {new_vertices[6], new_vertices[8]});
5226 new_lines[9]->set_bounding_object_indices(
5227 {new_vertices[8], new_vertices[7]});
5228 new_lines[10]->set_bounding_object_indices(
5229 {new_vertices[4], new_vertices[8]});
5230 new_lines[11]->set_bounding_object_indices(
5231 {new_vertices[8], new_vertices[5]});
5238 for (
unsigned int l = lmin; l < lmax; ++l)
5240 new_lines[l]->set_used_flag();
5241 new_lines[l]->clear_user_flag();
5242 new_lines[l]->clear_user_data();
5243 new_lines[l]->clear_children();
5245 new_lines[l]->set_boundary_id_internal(
5247 new_lines[l]->set_manifold_id(cell->
manifold_id());
5252 while (next_unused_cell->used() ==
true)
5255 unsigned int n_children = 0;
5263 for (
unsigned int i = 0; i < n_children; ++i)
5266 subcells[i] = next_unused_cell;
5268 if (i % 2 == 1 && i < n_children - 1)
5269 while (next_unused_cell->used() ==
true)
5275 static constexpr ::ndarray<unsigned int, 4, 4> tri_child_lines =
5276 {{{{0, 8, 5, X}}, {{1, 2, 6, X}}, {{7, 3, 4, X}}, {{6, 7, 8, X}}}};
5277 static constexpr ::ndarray<unsigned int, 4, 4>
5278 quad_child_lines = {{{{0, 8, 4, 10}},
5284 const auto &child_lines =
5288 for (
unsigned int i = 0; i < n_children; ++i)
5292 {new_lines[child_lines[i][0]]->index(),
5293 new_lines[child_lines[i][1]]->index(),
5294 new_lines[child_lines[i][2]]->index()});
5297 {new_lines[child_lines[i][0]]->index(),
5298 new_lines[child_lines[i][1]]->index(),
5299 new_lines[child_lines[i][2]]->index(),
5300 new_lines[child_lines[i][3]]->index()});
5319 face_no, inherited_orientations[child_lines[i][face_no]]);
5332 for (
unsigned int i = 0; i < n_children / 2; ++i)
5337 if (dim == spacedim - 1)
5338 for (
unsigned int c = 0; c < n_children; ++c)
5343 level < static_cast<int>(triangulation.
levels.size()) - 1;
5347 next_unused_cell = triangulation.
begin_raw(level + 1);
5349 for (
const auto &cell :
5360 check_for_distorted_cells &&
5361 has_distorted_children<dim, spacedim>(cell))
5369 return cells_with_distorted_children;
5378 template <
int spacedim>
5383 const unsigned int dim = 1;
5388 triangulation.
levels.size() - 1))
5391 triangulation.
levels.push_back(
5402 unsigned int needed_vertices = 0;
5403 for (
int level = triangulation.
levels.size() - 2; level >= 0; --level)
5407 unsigned int flagged_cells = 0;
5409 for (
const auto &acell :
5411 if (acell->refine_flag_set())
5416 const unsigned int used_cells =
5417 std::count(triangulation.
levels[level + 1]->cells.used.begin(),
5418 triangulation.
levels[level + 1]->cells.used.end(),
5436 needed_vertices += flagged_cells;
5441 needed_vertices += std::count(triangulation.
vertices_used.begin(),
5447 if (needed_vertices > triangulation.
vertices.size())
5458 unsigned int next_unused_vertex = 0;
5460 for (
int level = triangulation.
levels.size() - 2; level >= 0; --level)
5463 next_unused_cell = triangulation.
begin_raw(level + 1);
5465 for (
const auto &cell :
5476 ++next_unused_vertex;
5478 next_unused_vertex < triangulation.
vertices.size(),
5480 "Internal error: During refinement, the triangulation "
5481 "wants to access an element of the 'vertices' array "
5482 "but it turns out that the array is not large enough."));
5487 triangulation.
vertices[next_unused_vertex] =
5497 while (next_unused_cell->
used() ==
true)
5499 first_child = next_unused_cell;
5504 second_child = next_unused_cell;
5518 if (dim == spacedim - 1)
5557 left_neighbor = left_neighbor->
child(nbnb);
5570 if (dim == spacedim - 1)
5589 right_neighbor = cell->
neighbor(1);
5592 right_neighbor = right_neighbor->
child(nbnb);
5613 template <
int spacedim>
5616 const bool check_for_distorted_cells)
5618 const unsigned int dim = 2;
5623 bool do_isotropic_refinement =
true;
5628 do_isotropic_refinement =
false;
5632 if (do_isotropic_refinement)
5634 check_for_distorted_cells);
5642 triangulation.
levels.size() - 1))
5645 triangulation.
levels.push_back(
5658 line->clear_user_flag();
5664 unsigned int n_single_lines = 0;
5668 unsigned int n_lines_in_pairs = 0;
5674 unsigned int needed_vertices = 0;
5675 for (
int level = triangulation.
levels.size() - 2; level >= 0; --level)
5679 unsigned int needed_cells = 0;
5681 for (
const auto &cell :
5694 n_single_lines += 4;
5706 n_single_lines += 1;
5714 for (
const unsigned int line_no :
5722 line = cell->
line(line_no);
5723 if (line->has_children() ==
false)
5724 line->set_user_flag();
5731 const unsigned int used_cells =
5732 std::count(triangulation.
levels[level + 1]->cells.used.begin(),
5733 triangulation.
levels[level + 1]->cells.used.end(),
5741 used_cells + needed_cells,
5757 if (line->user_flag_set())
5760 n_lines_in_pairs += 2;
5761 needed_vertices += 1;
5773 needed_vertices += std::count(triangulation.
vertices_used.begin(),
5779 if (needed_vertices > triangulation.
vertices.size())
5790 unsigned int next_unused_vertex = 0;
5802 for (; line != endl; ++line)
5803 if (line->user_flag_set())
5809 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5810 ++next_unused_vertex;
5812 next_unused_vertex < triangulation.
vertices.size(),
5814 "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."));
5817 triangulation.
vertices[next_unused_vertex] = line->center(
true);
5822 bool pair_found =
false;
5824 for (; next_unused_line != endl; ++next_unused_line)
5825 if (!next_unused_line->used() &&
5826 !(++next_unused_line)->used())
5839 line->set_children(0, next_unused_line->index());
5843 children[2] = {next_unused_line, ++next_unused_line};
5849 children[0]->set_bounding_object_indices(
5850 {line->vertex_index(0), next_unused_vertex});
5851 children[1]->set_bounding_object_indices(
5852 {next_unused_vertex, line->vertex_index(1)});
5854 children[0]->set_used_flag();
5855 children[1]->set_used_flag();
5856 children[0]->clear_children();
5857 children[1]->clear_children();
5860 children[0]->clear_user_flag();
5861 children[1]->clear_user_flag();
5864 children[0]->set_boundary_id_internal(line->boundary_id());
5865 children[1]->set_boundary_id_internal(line->boundary_id());
5867 children[0]->set_manifold_id(line->manifold_id());
5868 children[1]->set_manifold_id(line->manifold_id());
5872 line->clear_user_flag();
5884 cells_with_distorted_children;
5892 level < static_cast<int>(triangulation.
levels.size()) - 1;
5896 next_unused_cell = triangulation.
begin_raw(level + 1);
5898 for (
const auto &cell :
5910 if (check_for_distorted_cells &&
5911 has_distorted_children<dim, spacedim>(cell))
5919 return cells_with_distorted_children;
5923 template <
int spacedim>
5926 const bool check_for_distorted_cells)
5928 static const int dim = 3;
5930 using raw_line_iterator =
5932 using raw_quad_iterator =
5944 triangulation.
levels.size() - 1))
5947 triangulation.
levels.push_back(
5955 triangulation.
faces->quads.clear_user_data();
5956 triangulation.
faces->lines.clear_user_flags();
5957 triangulation.
faces->quads.clear_user_flags();
5971 unsigned int needed_vertices = 0;
5972 unsigned int needed_lines_single = 0;
5973 unsigned int needed_quads_single = 0;
5974 unsigned int needed_lines_pair = 0;
5975 unsigned int needed_quads_pair = 0;
5976 for (
int level = triangulation.
levels.size() - 2; level >= 0; --level)
5978 unsigned int new_cells = 0;
5980 for (
const auto &cell :
5996 needed_lines_single += 6;
5997 needed_quads_single += 12;
6002 needed_lines_single += 1;
6003 needed_quads_single += 8;
6015 if (cell->
face(face)->n_children() == 0)
6016 cell->
face(face)->set_user_flag();
6025 if (cell->
line(line)->has_children() ==
false)
6026 cell->
line(line)->set_user_flag();
6032 const unsigned int used_cells =
6033 std::count(triangulation.
levels[level + 1]->cells.used.begin(),
6034 triangulation.
levels[level + 1]->cells.used.end(),
6039 used_cells + new_cells,
6045 used_cells + new_cells,
6060 if (quad->user_flag_set() ==
false)
6065 needed_quads_pair += 4;
6066 needed_lines_pair += 4;
6067 needed_vertices += 1;
6071 needed_quads_pair += 4;
6072 needed_lines_single += 3;
6085 if (line->user_flag_set() ==
false)
6088 needed_lines_pair += 2;
6089 needed_vertices += 1;
6094 needed_lines_single);
6097 needed_quads_single);
6100 needed_quads_single);
6104 needed_vertices += std::count(triangulation.
vertices_used.begin(),
6108 if (needed_vertices > triangulation.
vertices.size())
6127 for (
unsigned int line_n = 0; line_n < cell->
n_lines(); ++line_n)
6128 if (cell->
line(line_n)->has_children())
6129 for (
unsigned int c = 0; c < 2; ++c)
6130 Assert(cell->
line(line_n)->child(c)->user_flag_set() ==
false,
6134 unsigned int current_vertex = 0;
6138 auto get_next_unused_vertex = [](
const unsigned int current_vertex,
6139 std::vector<bool> &vertices_used) {
6140 unsigned int next_vertex = current_vertex;
6141 while (next_vertex < vertices_used.size() &&
6142 vertices_used[next_vertex] ==
true)
6145 vertices_used[next_vertex] =
true;
6155 raw_line_iterator next_unused_line = triangulation.
begin_raw_line();
6157 for (; line != endl; ++line)
6159 if (line->user_flag_set() ==
false)
6163 triangulation.
faces->lines.template next_free_pair_object<1>(
6171 line->set_children(0, next_unused_line->index());
6173 const std::array<raw_line_iterator, 2> children{
6174 {next_unused_line, ++next_unused_line}};
6180 get_next_unused_vertex(current_vertex,
6182 triangulation.
vertices[current_vertex] = line->center(
true);
6184 children[0]->set_bounding_object_indices(
6185 {line->vertex_index(0), current_vertex});
6186 children[1]->set_bounding_object_indices(
6187 {current_vertex, line->vertex_index(1)});
6189 const auto manifold_id = line->manifold_id();
6190 const auto boundary_id = line->boundary_id();
6191 for (
const auto &child : children)
6193 child->set_used_flag();
6194 child->clear_children();
6195 child->clear_user_data();
6196 child->clear_user_flag();
6197 child->set_boundary_id_internal(boundary_id);
6198 child->set_manifold_id(manifold_id);
6201 line->clear_user_flag();
6211 for (; quad != endq; ++quad)
6213 if (quad->user_flag_set() ==
false)
6216 const auto reference_face_type = quad->reference_cell();
6220 std::array<raw_line_iterator, 4> new_lines;
6223 for (
unsigned int l = 0; l < 2; ++l)
6225 auto next_unused_line =
6226 triangulation.
faces->lines
6227 .template next_free_pair_object<1>(triangulation);
6228 new_lines[2 * l] = next_unused_line;
6229 new_lines[2 * l + 1] = ++next_unused_line;
6234 for (
unsigned int l = 0; l < 3; ++l)
6236 triangulation.
faces->lines
6237 .template next_free_single_object<1>(triangulation);
6244 for (
const unsigned int line : quad->line_indices())
6252 std::array<raw_quad_iterator, 4> new_quads;
6253 for (
unsigned int q = 0; q < 2; ++q)
6255 auto next_unused_quad =
6256 triangulation.
faces->quads
6257 .template next_free_pair_object<2>(triangulation);
6259 new_quads[2 * q] = next_unused_quad;
6260 new_quads[2 * q + 1] = ++next_unused_quad;
6262 quad->set_children(2 * q, new_quads[2 * q]->
index());
6266 for (
const auto &quad : new_quads)
6276 std::array<unsigned int, 9> vertex_indices = {};
6278 for (
const auto i : quad->vertex_indices())
6279 vertex_indices[k++] = quad->vertex_index(i);
6281 for (
const auto i : quad->line_indices())
6282 vertex_indices[k++] = quad->line(i)->child(0)->vertex_index(1);
6287 get_next_unused_vertex(current_vertex,
6289 vertex_indices[k++] = current_vertex;
6291 triangulation.
vertices[current_vertex] =
6292 quad->center(
true,
true);
6296 std::array<raw_line_iterator, 12> lines;
6297 unsigned int n_lines = 0;
6298 for (
unsigned int l = 0; l < quad->
n_lines(); ++l)
6299 for (
unsigned int c = 0; c < 2; ++c)
6301 static constexpr ::ndarray<unsigned int, 2, 2>
index =
6308 quad->line(l)->child(
index[c][quad->line_orientation(l)]);
6311 for (
unsigned int l = 0; l < quad->
n_lines(); ++l)
6312 lines[n_lines++] = new_lines[l];
6314 std::array<int, 12> line_indices;
6315 for (
unsigned int i = 0; i < n_lines; ++i)
6316 line_indices[i] = lines[i]->
index();
6318 static constexpr ::ndarray<unsigned int, 12, 2>
6319 line_vertices_quad{{{{0, 4}},
6332 static constexpr ::ndarray<unsigned int, 4, 4>
6333 quad_lines_quad{{{{0, 8, 4, 10}},
6338 static constexpr ::ndarray<unsigned int, 12, 2>
6339 line_vertices_tri{{{{0, 3}},
6352 static constexpr ::ndarray<unsigned int, 4, 4>
6353 quad_lines_tri{{{{0, 8, 5, X}},
6358 static constexpr ::ndarray<unsigned int, 4, 4, 2>
6359 quad_line_vertices_tri{
6360 {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
6361 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
6362 {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
6363 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
6365 const auto &line_vertices =
6367 line_vertices_quad :
6369 const auto &quad_lines =
6374 for (
unsigned int i = 0, j = 2 * quad->
n_lines();
6375 i < quad->n_lines();
6378 auto &new_line = new_lines[i];
6379 new_line->set_bounding_object_indices(
6380 {vertex_indices[line_vertices[j][0]],
6381 vertex_indices[line_vertices[j][1]]});
6382 new_line->set_used_flag();
6383 new_line->clear_user_flag();
6384 new_line->clear_user_data();
6385 new_line->clear_children();
6386 new_line->set_boundary_id_internal(quad->boundary_id());
6387 new_line->set_manifold_id(quad->manifold_id());
6391 for (
unsigned int i = 0; i < new_quads.size(); ++i)
6393 auto &new_quad = new_quads[i];
6397 triangulation.
faces->set_quad_type(new_quad->index(),
6398 reference_face_type);
6401 new_quad->set_bounding_object_indices(
6402 {line_indices[quad_lines[i][0]],
6403 line_indices[quad_lines[i][1]],
6404 line_indices[quad_lines[i][2]]});
6406 new_quad->set_bounding_object_indices(
6407 {line_indices[quad_lines[i][0]],
6408 line_indices[quad_lines[i][1]],
6409 line_indices[quad_lines[i][2]],
6410 line_indices[quad_lines[i][3]]});
6414 new_quad->set_used_flag();
6415 new_quad->clear_user_flag();
6416 new_quad->clear_user_data();
6417 new_quad->clear_children();
6418 new_quad->set_boundary_id_internal(quad->boundary_id());
6419 new_quad->set_manifold_id(quad->manifold_id());
6422 std::set<unsigned int> s;
6430 for (
const auto f : new_quad->line_indices())
6432 const std::array<unsigned int, 2> vertices_0 = {
6433 {lines[quad_lines[i][f]]->vertex_index(0),
6434 lines[quad_lines[i][f]]->vertex_index(1)}};
6436 const std::array<unsigned int, 2> vertices_1 = {
6437 {vertex_indices[quad_line_vertices_tri[i][f][0]],
6438 vertex_indices[quad_line_vertices_tri[i][f][1]]}};
6440 const auto orientation =
6446 for (
const auto i : vertices_0)
6448 for (
const auto i : vertices_1)
6452 new_quad->set_line_orientation(
6456 default_combined_face_orientation());
6468 static constexpr ::ndarray<unsigned int, 4, 2>
6469 quad_child_boundary_lines{
6470 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
6472 for (
unsigned int i = 0; i < 4; ++i)
6473 for (
unsigned int j = 0; j < 2; ++j)
6474 new_quads[quad_child_boundary_lines[i][j]]
6475 ->set_line_orientation(i, quad->line_orientation(i));
6478 quad->clear_user_flag();
6483 cells_with_distorted_children;
6487 for (
unsigned int level = 0; level != triangulation.
levels.size() - 1;
6493 hex->level() >=
static_cast<int>(level),
6496 for (; hex != triangulation.
end() &&
6497 hex->level() ==
static_cast<int>(level);
6500 if (hex->refine_flag_set() ==
6504 const auto &reference_cell_type = hex->reference_cell();
6507 hex->clear_refine_flag();
6508 hex->set_refinement_case(ref_case);
6510 unsigned int n_new_lines = 0;
6511 unsigned int n_new_quads = 0;
6512 unsigned int n_new_hexes = 0;
6529 std::array<raw_line_iterator, 6> new_lines;
6530 for (
unsigned int i = 0; i < n_new_lines; ++i)
6533 triangulation.
faces->lines
6534 .template next_free_single_object<1>(triangulation);
6537 new_lines[i]->set_used_flag();
6538 new_lines[i]->clear_user_flag();
6539 new_lines[i]->clear_user_data();
6540 new_lines[i]->clear_children();
6541 new_lines[i]->set_boundary_id_internal(
6543 new_lines[i]->set_manifold_id(hex->manifold_id());
6546 std::array<raw_quad_iterator, 12> new_quads;
6547 for (
unsigned int i = 0; i < n_new_quads; ++i)
6550 triangulation.
faces->quads
6551 .template next_free_single_object<2>(triangulation);
6553 auto &new_quad = new_quads[i];
6557 triangulation.
faces->set_quad_type(
6559 reference_cell_type.face_reference_cell(0));
6562 new_quad->set_used_flag();
6563 new_quad->clear_user_flag();
6564 new_quad->clear_user_data();
6565 new_quad->clear_children();
6566 new_quad->set_boundary_id_internal(
6568 new_quad->set_manifold_id(hex->manifold_id());
6569 for (
const auto j : new_quads[i]->line_indices())
6570 new_quad->set_line_orientation(j,
true);
6579 for (
unsigned int i = 0; i < n_new_hexes; ++i)
6583 triangulation.
levels[level + 1]->cells.next_free_hex(
6584 triangulation, level + 1);
6588 new_hexes[i] = next_unused_hex;
6590 auto &new_hex = new_hexes[i];
6593 triangulation.
levels[new_hex->level()]
6594 ->reference_cell[new_hex->index()] =
6595 reference_cell_type;
6598 new_hex->set_used_flag();
6599 new_hex->clear_user_flag();
6600 new_hex->clear_user_data();
6601 new_hex->clear_children();
6602 new_hex->set_material_id(hex->material_id());
6603 new_hex->set_manifold_id(hex->manifold_id());
6604 new_hex->set_subdomain_id(hex->subdomain_id());
6607 new_hex->set_parent(hex->index());
6612 for (
const auto f : new_hex->face_indices())
6613 new_hex->set_combined_face_orientation(
6617 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
6618 hex->set_children(2 * i, new_hexes[2 * i]->index());
6623 std::array<unsigned int, 27> vertex_indices = {};
6630 const unsigned int n_vertices =
6632 for (
unsigned int i = 0; i < n_vertices; ++i)
6633 vertex_indices[k++] = hex->vertex_index(i);
6635 const std::array<unsigned int, 12> line_indices =
6636 TriaAccessorImplementation::Implementation::
6637 get_line_indices_of_cell(*hex);
6649 for (
unsigned int l = 0; l < n_lines; ++l)
6651 raw_line_iterator line(&triangulation,
6654 vertex_indices[k++] = line->child(0)->vertex_index(1);
6659 for (
const unsigned int i : hex->face_indices())
6660 vertex_indices[k++] =
6661 hex->face(i)->child(0)->vertex_index(3);
6665 get_next_unused_vertex(current_vertex,
6667 vertex_indices[k++] = current_vertex;
6669 triangulation.
vertices[current_vertex] =
6670 hex->center(
true,
true);
6674 unsigned int chosen_line_tetrahedron = 0;
6678 static constexpr ::ndarray<unsigned int, 6, 2>
6679 new_line_vertices = {{{{22, 26}},
6685 for (
unsigned int i = 0; i < n_new_lines; ++i)
6686 new_lines[i]->set_bounding_object_indices(
6687 {vertex_indices[new_line_vertices[i][0]],
6688 vertex_indices[new_line_vertices[i][1]]});
6696 static constexpr ::ndarray<unsigned int, 3, 2>
6697 new_line_vertices = {{{{6, 8}}, {{5, 7}}, {{4, 9}}}};
6701 std::uint8_t refinement_choice = hex->refine_choice();
6702 if (refinement_choice ==
6707 double min_distance =
6708 std::numeric_limits<double>::infinity();
6709 for (
unsigned int i = 0; i < new_line_vertices.size();
6712 const double current_distance =
6714 [vertex_indices[new_line_vertices[i][0]]]
6716 vertices[vertex_indices
6717 [new_line_vertices[i][1]]]);
6718 if (current_distance < min_distance)
6720 chosen_line_tetrahedron = i;
6721 min_distance = current_distance;
6725 else if (refinement_choice ==
6728 chosen_line_tetrahedron = 0;
6729 else if (refinement_choice ==
6732 chosen_line_tetrahedron = 1;
6733 else if (refinement_choice ==
6736 chosen_line_tetrahedron = 2;
6740 hex->set_refinement_case(
6743 new_lines[0]->set_bounding_object_indices(
6745 [new_line_vertices[chosen_line_tetrahedron][0]],
6747 [new_line_vertices[chosen_line_tetrahedron][1]]});
6752 boost::container::small_vector<raw_line_iterator, 30>
6757 relevant_lines.resize(30);
6758 for (
unsigned int f = 0, k = 0; f < 6; ++f)
6759 for (
unsigned int c = 0; c < 4; ++c, ++k)
6762 ndarray<unsigned int, 4, 2>
6764 {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}};
6770 standard_to_real_face_vertex(
6772 hex->face_orientation(f),
6774 hex->face_rotation(f)))
6776 standard_to_real_face_line(
6778 hex->face_orientation(f),
6780 hex->face_rotation(f)));
6783 for (
unsigned int i = 0, k = 24; i < 6; ++i, ++k)
6784 relevant_lines[k] = new_lines[i];
6799 relevant_lines.resize(13);
6802 for (
unsigned int f = 0; f < 4; ++f)
6803 for (
unsigned int l = 0; l < 3; ++l, ++k)
6807 array<std::array<unsigned int, 3>, 6>
6808 table = {{{{1, 0, 2}},
6815 const unsigned char combined_orientation =
6816 hex->combined_face_orientation(f);
6820 ->line(table[combined_orientation][l]);
6823 relevant_lines[k++] = new_lines[0];
6829 boost::container::small_vector<unsigned int, 30>
6830 relevant_line_indices(relevant_lines.size());
6831 for (
unsigned int i = 0; i < relevant_line_indices.size();
6833 relevant_line_indices[i] = relevant_lines[i]->
index();
6842 const auto &new_quad_lines =
6843 hex->reference_cell().new_isotropic_child_face_lines(
6844 chosen_line_tetrahedron);
6857 hex->reference_cell()
6858 .new_isotropic_child_face_line_vertices(
6859 chosen_line_tetrahedron);
6861 static constexpr ::ndarray<unsigned int, 4, 2>
6862 representative_lines{
6863 {{{0, 2}}, {{2, 0}}, {{3, 3}}, {{1, 1}}}};
6865 for (
unsigned int q = 0; q < n_new_quads; ++q)
6867 auto &new_quad = new_quads[q];
6869 if (new_quad->n_lines() == 3)
6870 new_quad->set_bounding_object_indices(
6871 {relevant_line_indices[new_quad_lines[q][0]],
6872 relevant_line_indices[new_quad_lines[q][1]],
6873 relevant_line_indices[new_quad_lines[q][2]]});
6874 else if (new_quad->n_lines() == 4)
6875 new_quad->set_bounding_object_indices(
6876 {relevant_line_indices[new_quad_lines[q][0]],
6877 relevant_line_indices[new_quad_lines[q][1]],
6878 relevant_line_indices[new_quad_lines[q][2]],
6879 relevant_line_indices[new_quad_lines[q][3]]});
6887 const unsigned int n_compute_lines =
6890 new_quad->n_lines();
6891 for (
unsigned int line = 0; line < n_compute_lines;
6894 const unsigned int l =
6895 (reference_cell_type ==
6897 representative_lines[q % 4][0] :
6900 const std::array<unsigned int, 2> vertices_0 = {
6901 {relevant_lines[new_quad_lines[q][l]]
6903 relevant_lines[new_quad_lines[q][l]]
6904 ->vertex_index(1)}};
6906 const std::array<unsigned int, 2> vertices_1 = {
6907 {vertex_indices[table[q][l][0]],
6908 vertex_indices[table[q][l][1]]}};
6910 const auto orientation =
6915 new_quad->set_line_orientation(
6919 default_combined_face_orientation());
6924 if (reference_cell_type ==
6926 new_quads[representative_lines[q % 4][1] + q -
6928 ->set_line_orientation(
6932 default_combined_face_orientation());
6939 std::array<int, 36> quad_indices;
6943 for (
unsigned int i = 0; i < n_new_quads; ++i)
6944 quad_indices[i] = new_quads[i]->
index();
6946 for (
unsigned int f = 0, k = n_new_quads; f < 6; ++f)
6947 for (
unsigned int c = 0; c < 4; ++c, ++k)
6949 hex->face(f)->isotropic_child_index(
6952 hex->face_orientation(f),
6954 hex->face_rotation(f)));
6964 for (
unsigned int i = 0; i < n_new_quads; ++i)
6965 quad_indices[i] = new_quads[i]->
index();
6967 for (
unsigned int f = 0, k = n_new_quads; f < 4; ++f)
6968 for (
unsigned int c = 0; c < 4; ++c, ++k)
6970 const unsigned char combined_orientation =
6971 hex->combined_face_orientation(f);
6972 quad_indices[k] = hex->face(f)->child_index(
6975 .standard_to_real_face_vertex(
6976 c, f, combined_orientation));
6992 const auto &cell_quads =
6993 hex->reference_cell().new_isotropic_child_cell_faces(
6994 chosen_line_tetrahedron);
6996 for (
unsigned int c = 0;
6997 c < GeometryInfo<dim>::max_children_per_cell;
7000 auto &new_hex = new_hexes[c];
7002 if (new_hex->n_faces() == 4)
7004 new_hex->set_bounding_object_indices(
7005 {quad_indices[cell_quads[c][0]],
7006 quad_indices[cell_quads[c][1]],
7007 quad_indices[cell_quads[c][2]],
7008 quad_indices[cell_quads[c][3]]});
7013 for (
const auto f : new_hex->face_indices())
7015 const auto &face = new_hex->face(f);
7017 Assert(face->n_vertices() == 3,
7020 const std::array<unsigned int, 3> vertices_0 = {
7021 {face->vertex_index(0),
7022 face->vertex_index(1),
7023 face->vertex_index(2)}};
7032 const auto new_hex_vertices =
7033 hex->reference_cell()
7034 .new_isotropic_child_cell_vertices(
7035 chosen_line_tetrahedron)[c];
7039 const std::array<unsigned int, 3> vertices_1 = {
7044 .face_to_cell_vertices(f, 0, 1)]],
7048 .face_to_cell_vertices(f, 1, 1)]],
7052 .face_to_cell_vertices(f, 2, 1)]],
7055 new_hex->set_combined_face_orientation(
7057 face->reference_cell()
7058 .get_combined_orientation(
7063 else if (new_hex->n_faces() == 6)
7064 new_hex->set_bounding_object_indices(
7065 {quad_indices[cell_quads[c][0]],
7066 quad_indices[cell_quads[c][1]],
7067 quad_indices[cell_quads[c][2]],
7068 quad_indices[cell_quads[c][3]],
7069 quad_indices[cell_quads[c][4]],
7070 quad_indices[cell_quads[c][5]]});
7079 static constexpr ::ndarray<unsigned int, 6, 4>
7080 face_to_child_indices_hex{{{{0, 2, 4, 6}},
7087 for (
const auto f : hex->face_indices())
7089 const unsigned char combined_orientation =
7090 hex->combined_face_orientation(f);
7091 for (
unsigned int c = 0; c < 4; ++c)
7092 new_hexes[face_to_child_indices_hex[f][c]]
7093 ->set_combined_face_orientation(
7094 f, combined_orientation);
7099 if (check_for_distorted_cells &&
7100 has_distorted_children<dim, spacedim>(hex))
7107 triangulation.
faces->quads.clear_user_data();
7109 return cells_with_distorted_children;
7116 template <
int spacedim>
7119 const bool check_for_distorted_cells)
7121 const unsigned int dim = 3;
7124 bool flag_isotropic_mesh =
true;
7126 cell = triangulation.
begin(),
7127 endc = triangulation.
end();
7128 for (; cell != endc; ++cell)
7138 flag_isotropic_mesh =
false;
7142 if (flag_isotropic_mesh)
7144 check_for_distorted_cells);
7156 triangulation.
levels.size() - 1))
7159 triangulation.
levels.push_back(
7168 triangulation.
faces->quads.clear_user_data();
7174 line->clear_user_flag();
7179 quad->clear_user_flag();
7202 unsigned int needed_vertices = 0;
7203 unsigned int needed_lines_single = 0;
7204 unsigned int needed_quads_single = 0;
7205 unsigned int needed_lines_pair = 0;
7206 unsigned int needed_quads_pair = 0;
7207 for (
int level = triangulation.
levels.size() - 2; level >= 0; --level)
7211 unsigned int new_cells = 0;
7213 for (
const auto &acell :
7215 if (acell->refine_flag_set())
7225 ++needed_quads_single;
7233 ++needed_lines_single;
7234 needed_quads_single += 4;
7241 needed_lines_single += 6;
7242 needed_quads_single += 12;
7256 for (
const unsigned int face :
7260 aface = acell->face(face);
7267 acell->face_orientation(face),
7268 acell->face_flip(face),
7269 acell->face_rotation(face));
7274 if (face_ref_case ==
7277 if (aface->n_active_descendants() < 4)
7280 aface->set_user_flag();
7282 else if (aface->refinement_case() != face_ref_case)
7293 Assert(aface->refinement_case() ==
7296 aface->refinement_case() ==
7299 aface->set_user_index(face_ref_case);
7305 for (
unsigned int line = 0;
7306 line < GeometryInfo<dim>::lines_per_cell;
7310 !acell->line(line)->has_children())
7311 acell->line(line)->set_user_flag();
7317 const unsigned int used_cells =
7318 std::count(triangulation.
levels[level + 1]->cells.used.begin(),
7319 triangulation.
levels[level + 1]->cells.used.end(),
7327 used_cells + new_cells,
7341 if (quad->user_flag_set())
7347 needed_quads_pair += 4;
7348 needed_lines_pair += 4;
7349 needed_vertices += 1;
7351 if (quad->user_index())
7355 needed_quads_pair += 2;
7356 needed_lines_single += 1;
7369 if (quad->has_children())
7371 Assert(quad->refinement_case() ==
7374 if ((face_refinement_cases[quad->user_index()] ==
7376 (quad->child(0)->line_index(1) + 1 !=
7377 quad->child(2)->line_index(1))) ||
7378 (face_refinement_cases[quad->user_index()] ==
7380 (quad->child(0)->line_index(3) + 1 !=
7381 quad->child(1)->line_index(3))))
7382 needed_lines_pair += 2;
7391 if (line->user_flag_set())
7393 needed_lines_pair += 2;
7394 needed_vertices += 1;
7400 needed_lines_single);
7404 needed_quads_single);
7407 needed_quads_single);
7411 needed_vertices += std::count(triangulation.
vertices_used.begin(),
7417 if (needed_vertices > triangulation.
vertices.size())
7437 for (
unsigned int line = 0;
7438 line < GeometryInfo<dim>::lines_per_cell;
7440 if (cell->
line(line)->has_children())
7441 for (
unsigned int c = 0; c < 2; ++c)
7442 Assert(cell->
line(line)->child(c)->user_flag_set() ==
false,
7454 unsigned int next_unused_vertex = 0;
7465 for (; line != endl; ++line)
7466 if (line->user_flag_set())
7472 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
7473 ++next_unused_vertex;
7475 next_unused_vertex < triangulation.
vertices.size(),
7477 "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."));
7480 triangulation.
vertices[next_unused_vertex] = line->center(
true);
7486 triangulation.
faces->lines.template next_free_pair_object<1>(
7494 line->set_children(0, next_unused_line->index());
7498 children[2] = {next_unused_line, ++next_unused_line};
7505 children[0]->set_bounding_object_indices(
7506 {line->vertex_index(0), next_unused_vertex});
7507 children[1]->set_bounding_object_indices(
7508 {next_unused_vertex, line->vertex_index(1)});
7510 children[0]->set_used_flag();
7511 children[1]->set_used_flag();
7512 children[0]->clear_children();
7513 children[1]->clear_children();
7516 children[0]->clear_user_flag();
7517 children[1]->clear_user_flag();
7519 children[0]->set_boundary_id_internal(line->boundary_id());
7520 children[1]->set_boundary_id_internal(line->boundary_id());
7522 children[0]->set_manifold_id(line->manifold_id());
7523 children[1]->set_manifold_id(line->manifold_id());
7528 line->clear_user_flag();
7558 for (
unsigned int loop = 0; loop < 2; ++loop)
7571 for (; quad != endq; ++quad)
7573 if (quad->user_index())
7576 face_refinement_cases[quad->user_index()];
7585 if (aniso_quad_ref_case == quad->refinement_case())
7588 Assert(quad->refinement_case() ==
7590 quad->refinement_case() ==
7595 Assert(quad->user_index() ==
7597 quad->user_index() ==
7606 triangulation.
faces->lines
7607 .template next_free_single_object<1>(triangulation);
7623 unsigned int vertex_indices[2];
7627 quad->line(2)->child(0)->vertex_index(1);
7629 quad->line(3)->child(0)->vertex_index(1);
7634 quad->line(0)->child(0)->vertex_index(1);
7636 quad->line(1)->child(0)->vertex_index(1);
7639 new_line->set_bounding_object_indices(
7640 {vertex_indices[0], vertex_indices[1]});
7641 new_line->set_used_flag();
7642 new_line->clear_user_flag();
7644 new_line->clear_children();
7645 new_line->set_boundary_id_internal(quad->boundary_id());
7646 new_line->set_manifold_id(quad->manifold_id());
7654 const unsigned int index[2][2] = {
7664 triangulation.
faces->quads
7665 .template next_free_pair_object<2>(triangulation);
7666 new_quads[0] = next_unused_quad;
7670 new_quads[1] = next_unused_quad;
7675 new_quads[0]->set_bounding_object_indices(
7676 {
static_cast<int>(quad->line_index(0)),
7679 ->child(
index[0][quad->line_orientation(2)])
7682 ->child(
index[0][quad->line_orientation(3)])
7684 new_quads[1]->set_bounding_object_indices(
7686 static_cast<int>(quad->line_index(1)),
7688 ->child(
index[1][quad->line_orientation(2)])
7691 ->child(
index[1][quad->line_orientation(3)])
7696 new_quads[0]->set_bounding_object_indices(
7698 ->child(
index[0][quad->line_orientation(0)])
7701 ->child(
index[0][quad->line_orientation(1)])
7703 static_cast<int>(quad->line_index(2)),
7704 new_line->index()});
7705 new_quads[1]->set_bounding_object_indices(
7707 ->child(
index[1][quad->line_orientation(0)])
7710 ->child(
index[1][quad->line_orientation(1)])
7713 static_cast<int>(quad->line_index(3))});
7716 for (
const auto &new_quad : new_quads)
7718 new_quad->set_used_flag();
7719 new_quad->clear_user_flag();
7720 new_quad->clear_user_data();
7721 new_quad->clear_children();
7722 new_quad->set_boundary_id_internal(quad->boundary_id());
7723 new_quad->set_manifold_id(quad->manifold_id());
7727 for (
unsigned int j = 0;
7728 j < GeometryInfo<dim>::lines_per_face;
7730 new_quad->set_line_orientation(j,
true);
7736 new_quads[0]->set_line_orientation(
7737 0, quad->line_orientation(0));
7738 new_quads[0]->set_line_orientation(
7739 2, quad->line_orientation(2));
7740 new_quads[1]->set_line_orientation(
7741 1, quad->line_orientation(1));
7742 new_quads[1]->set_line_orientation(
7743 3, quad->line_orientation(3));
7746 new_quads[0]->set_line_orientation(
7747 3, quad->line_orientation(3));
7748 new_quads[1]->set_line_orientation(
7749 2, quad->line_orientation(2));
7753 new_quads[0]->set_line_orientation(
7754 1, quad->line_orientation(1));
7755 new_quads[1]->set_line_orientation(
7756 0, quad->line_orientation(0));
7762 if (quad->refinement_case() ==
7785 if (aniso_quad_ref_case ==
7788 old_child[0] = quad->child(0)->line(1);
7789 old_child[1] = quad->child(2)->line(1);
7793 Assert(aniso_quad_ref_case ==
7797 old_child[0] = quad->child(0)->line(3);
7798 old_child[1] = quad->child(1)->line(3);
7801 if (old_child[0]->
index() + 1 != old_child[1]->
index())
7807 spacedim>::raw_line_iterator
7810 new_child[0] = new_child[1] =
7811 triangulation.
faces->lines
7812 .template next_free_pair_object<1>(
7816 new_child[0]->set_used_flag();
7817 new_child[1]->set_used_flag();
7819 const int old_index_0 = old_child[0]->index(),
7820 old_index_1 = old_child[1]->index(),
7821 new_index_0 = new_child[0]->index(),
7822 new_index_1 = new_child[1]->index();
7826 for (
unsigned int q = 0;
7827 q < triangulation.
faces->quads.n_objects();
7829 for (
unsigned int l = 0;
7830 l < GeometryInfo<dim>::lines_per_face;
7833 const int this_index =
7834 triangulation.
faces->quads
7835 .get_bounding_object_indices(q)[l];
7836 if (this_index == old_index_0)
7837 triangulation.
faces->quads
7838 .get_bounding_object_indices(q)[l] =
7840 else if (this_index == old_index_1)
7841 triangulation.
faces->quads
7842 .get_bounding_object_indices(q)[l] =
7847 for (
unsigned int i = 0; i < 2; ++i)
7849 Assert(!old_child[i]->has_children(),
7852 new_child[i]->set_bounding_object_indices(
7853 {old_child[i]->vertex_index(0),
7854 old_child[i]->vertex_index(1)});
7855 new_child[i]->set_boundary_id_internal(
7856 old_child[i]->boundary_id());
7857 new_child[i]->set_manifold_id(
7858 old_child[i]->manifold_id());
7859 new_child[i]->set_user_index(
7860 old_child[i]->user_index());
7861 if (old_child[i]->user_flag_set())
7862 new_child[i]->set_user_flag();
7864 new_child[i]->clear_user_flag();
7866 new_child[i]->clear_children();
7868 old_child[i]->clear_user_flag();
7869 old_child[i]->clear_user_index();
7870 old_child[i]->clear_used_flag();
7876 if (aniso_quad_ref_case ==
7879 new_line->set_children(
7880 0, quad->child(0)->line_index(1));
7881 Assert(new_line->child(1) ==
7882 quad->child(2)->line(1),
7919 quad_iterator switch_1 = quad->child(1),
7920 switch_2 = quad->child(2);
7921 const int switch_1_index = switch_1->index();
7922 const int switch_2_index = switch_2->index();
7923 for (
unsigned int l = 0;
7924 l < triangulation.
levels.size();
7926 for (
unsigned int h = 0;
7928 triangulation.
levels[l]->cells.n_objects();
7930 for (
const unsigned int q :
7933 const int face_index =
7935 ->cells.get_bounding_object_indices(
7937 if (face_index == switch_1_index)
7939 ->cells.get_bounding_object_indices(
7940 h)[q] = switch_2_index;
7941 else if (face_index == switch_2_index)
7943 ->cells.get_bounding_object_indices(
7944 h)[q] = switch_1_index;
7948 const unsigned int switch_1_lines[4] = {
7949 switch_1->line_index(0),
7950 switch_1->line_index(1),
7951 switch_1->line_index(2),
7952 switch_1->line_index(3)};
7953 const bool switch_1_line_orientations[4] = {
7954 switch_1->line_orientation(0),
7955 switch_1->line_orientation(1),
7956 switch_1->line_orientation(2),
7957 switch_1->line_orientation(3)};
7959 switch_1->boundary_id();
7960 const unsigned int switch_1_user_index =
7961 switch_1->user_index();
7962 const bool switch_1_user_flag =
7963 switch_1->user_flag_set();
7965 switch_1_refinement_case =
7966 switch_1->refinement_case();
7967 const int switch_1_first_child_pair =
7968 (switch_1_refinement_case ?
7969 switch_1->child_index(0) :
7971 const int switch_1_second_child_pair =
7972 (switch_1_refinement_case ==
7974 switch_1->child_index(2) :
7977 switch_1->set_bounding_object_indices(
7978 {switch_2->line_index(0),
7979 switch_2->line_index(1),
7980 switch_2->line_index(2),
7981 switch_2->line_index(3)});
7982 switch_1->set_line_orientation(
7983 0, switch_2->line_orientation(0));
7984 switch_1->set_line_orientation(
7985 1, switch_2->line_orientation(1));
7986 switch_1->set_line_orientation(
7987 2, switch_2->line_orientation(2));
7988 switch_1->set_line_orientation(
7989 3, switch_2->line_orientation(3));
7990 switch_1->set_boundary_id_internal(
7991 switch_2->boundary_id());
7992 switch_1->set_manifold_id(switch_2->manifold_id());
7993 switch_1->set_user_index(switch_2->user_index());
7994 if (switch_2->user_flag_set())
7995 switch_1->set_user_flag();
7997 switch_1->clear_user_flag();
7998 switch_1->clear_refinement_case();
7999 switch_1->set_refinement_case(
8000 switch_2->refinement_case());
8001 switch_1->clear_children();
8002 if (switch_2->refinement_case())
8003 switch_1->set_children(0,
8004 switch_2->child_index(0));
8005 if (switch_2->refinement_case() ==
8007 switch_1->set_children(2,
8008 switch_2->child_index(2));
8010 switch_2->set_bounding_object_indices(
8014 switch_1_lines[3]});
8015 switch_2->set_line_orientation(
8016 0, switch_1_line_orientations[0]);
8017 switch_2->set_line_orientation(
8018 1, switch_1_line_orientations[1]);
8019 switch_2->set_line_orientation(
8020 2, switch_1_line_orientations[2]);
8021 switch_2->set_line_orientation(
8022 3, switch_1_line_orientations[3]);
8023 switch_2->set_boundary_id_internal(
8024 switch_1_boundary_id);
8025 switch_2->set_manifold_id(switch_1->manifold_id());
8026 switch_2->set_user_index(switch_1_user_index);
8027 if (switch_1_user_flag)
8028 switch_2->set_user_flag();
8030 switch_2->clear_user_flag();
8031 switch_2->clear_refinement_case();
8032 switch_2->set_refinement_case(
8033 switch_1_refinement_case);
8034 switch_2->clear_children();
8035 switch_2->set_children(0,
8036 switch_1_first_child_pair);
8037 switch_2->set_children(2,
8038 switch_1_second_child_pair);
8040 new_quads[0]->set_refinement_case(
8042 new_quads[0]->set_children(0, quad->child_index(0));
8043 new_quads[1]->set_refinement_case(
8045 new_quads[1]->set_children(0, quad->child_index(2));
8049 new_quads[0]->set_refinement_case(
8051 new_quads[0]->set_children(0, quad->child_index(0));
8052 new_quads[1]->set_refinement_case(
8054 new_quads[1]->set_children(0, quad->child_index(2));
8055 new_line->set_children(
8056 0, quad->child(0)->line_index(3));
8057 Assert(new_line->child(1) ==
8058 quad->child(1)->line(3),
8061 quad->clear_children();
8065 quad->set_children(0, new_quads[0]->
index());
8067 quad->set_refinement_case(aniso_quad_ref_case);
8074 if (quad->user_flag_set())
8086 ++next_unused_vertex;
8088 next_unused_vertex < triangulation.
vertices.size(),
8090 "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."));
8098 quad->refinement_case();
8104 quad->child(0)->set_user_index(
8106 quad->child(1)->set_user_index(
8112 middle_line = quad->child(0)->line(1);
8114 middle_line = quad->child(0)->line(3);
8121 if (!middle_line->has_children())
8128 triangulation.
vertices[next_unused_vertex] =
8129 middle_line->center(
true);
8136 triangulation.
faces->lines
8137 .template next_free_pair_object<1>(
8142 middle_line->set_children(
8143 0, next_unused_line->index());
8147 raw_line_iterator children[2] = {
8148 next_unused_line, ++next_unused_line};
8156 children[0]->set_bounding_object_indices(
8157 {middle_line->vertex_index(0),
8158 next_unused_vertex});
8159 children[1]->set_bounding_object_indices(
8160 {next_unused_vertex,
8161 middle_line->vertex_index(1)});
8163 children[0]->set_used_flag();
8164 children[1]->set_used_flag();
8165 children[0]->clear_children();
8166 children[1]->clear_children();
8169 children[0]->clear_user_flag();
8170 children[1]->clear_user_flag();
8172 children[0]->set_boundary_id_internal(
8173 middle_line->boundary_id());
8174 children[1]->set_boundary_id_internal(
8175 middle_line->boundary_id());
8177 children[0]->set_manifold_id(
8178 middle_line->manifold_id());
8179 children[1]->set_manifold_id(
8180 middle_line->manifold_id());
8186 quad->clear_user_flag();
8215 triangulation.
vertices[next_unused_vertex] =
8216 quad->center(
true,
true);
8225 for (
unsigned int i = 0; i < 4; ++i)
8235 triangulation.
faces->lines
8236 .template next_free_pair_object<1>(triangulation);
8238 new_lines[i] = next_unused_line;
8261 const unsigned int vertex_indices[5] = {
8262 quad->line(0)->child(0)->vertex_index(1),
8263 quad->line(1)->child(0)->vertex_index(1),
8264 quad->line(2)->child(0)->vertex_index(1),
8265 quad->line(3)->child(0)->vertex_index(1),
8266 next_unused_vertex};
8268 new_lines[0]->set_bounding_object_indices(
8269 {vertex_indices[2], vertex_indices[4]});
8270 new_lines[1]->set_bounding_object_indices(
8271 {vertex_indices[4], vertex_indices[3]});
8272 new_lines[2]->set_bounding_object_indices(
8273 {vertex_indices[0], vertex_indices[4]});
8274 new_lines[3]->set_bounding_object_indices(
8275 {vertex_indices[4], vertex_indices[1]});
8277 for (
const auto &new_line : new_lines)
8279 new_line->set_used_flag();
8280 new_line->clear_user_flag();
8281 new_line->clear_user_data();
8282 new_line->clear_children();
8283 new_line->set_boundary_id_internal(quad->boundary_id());
8284 new_line->set_manifold_id(quad->manifold_id());
8303 const unsigned int index[2][2] = {
8307 const int line_indices[12] = {
8309 ->child(
index[0][quad->line_orientation(0)])
8312 ->child(
index[1][quad->line_orientation(0)])
8315 ->child(
index[0][quad->line_orientation(1)])
8318 ->child(
index[1][quad->line_orientation(1)])
8321 ->child(
index[0][quad->line_orientation(2)])
8324 ->child(
index[1][quad->line_orientation(2)])
8327 ->child(
index[0][quad->line_orientation(3)])
8330 ->child(
index[1][quad->line_orientation(3)])
8332 new_lines[0]->index(),
8333 new_lines[1]->index(),
8334 new_lines[2]->index(),
8335 new_lines[3]->index()};
8344 triangulation.
faces->quads
8345 .template next_free_pair_object<2>(triangulation);
8347 new_quads[0] = next_unused_quad;
8351 new_quads[1] = next_unused_quad;
8355 triangulation.
faces->quads
8356 .template next_free_pair_object<2>(triangulation);
8357 new_quads[2] = next_unused_quad;
8361 new_quads[3] = next_unused_quad;
8365 quad->set_children(0, new_quads[0]->
index());
8366 quad->set_children(2, new_quads[2]->
index());
8369 new_quads[0]->set_bounding_object_indices(
8374 new_quads[1]->set_bounding_object_indices(
8379 new_quads[2]->set_bounding_object_indices(
8384 new_quads[3]->set_bounding_object_indices(
8389 for (
const auto &new_quad : new_quads)
8391 new_quad->set_used_flag();
8392 new_quad->clear_user_flag();
8393 new_quad->clear_user_data();
8394 new_quad->clear_children();
8395 new_quad->set_boundary_id_internal(quad->boundary_id());
8396 new_quad->set_manifold_id(quad->manifold_id());
8400 for (
unsigned int j = 0;
8401 j < GeometryInfo<dim>::lines_per_face;
8403 new_quad->set_line_orientation(j,
true);
8409 new_quads[0]->set_line_orientation(
8410 0, quad->line_orientation(0));
8411 new_quads[0]->set_line_orientation(
8412 2, quad->line_orientation(2));
8413 new_quads[1]->set_line_orientation(
8414 1, quad->line_orientation(1));
8415 new_quads[1]->set_line_orientation(
8416 2, quad->line_orientation(2));
8417 new_quads[2]->set_line_orientation(
8418 0, quad->line_orientation(0));
8419 new_quads[2]->set_line_orientation(
8420 3, quad->line_orientation(3));
8421 new_quads[3]->set_line_orientation(
8422 1, quad->line_orientation(1));
8423 new_quads[3]->set_line_orientation(
8424 3, quad->line_orientation(3));
8428 quad->clear_user_flag();
8439 cells_with_distorted_children;
8441 for (
unsigned int level = 0; level != triangulation.
levels.size() - 1;
8453 for (; hex != endh; ++hex)
8454 if (hex->refine_flag_set())
8462 hex->clear_refine_flag();
8463 hex->set_refinement_case(ref_case);
8474 unsigned int n_new_lines = 0;
8475 unsigned int n_new_quads = 0;
8476 unsigned int n_new_hexes = 0;
8507 new_lines(n_new_lines);
8508 for (
unsigned int i = 0; i < n_new_lines; ++i)
8511 triangulation.
faces->lines
8512 .template next_free_single_object<1>(triangulation);
8515 new_lines[i]->set_used_flag();
8516 new_lines[i]->clear_user_flag();
8517 new_lines[i]->clear_user_data();
8518 new_lines[i]->clear_children();
8520 new_lines[i]->set_boundary_id_internal(
8524 new_lines[i]->set_manifold_id(hex->manifold_id());
8531 new_quads(n_new_quads);
8532 for (
unsigned int i = 0; i < n_new_quads; ++i)
8535 triangulation.
faces->quads
8536 .template next_free_single_object<2>(triangulation);
8539 new_quads[i]->set_used_flag();
8540 new_quads[i]->clear_user_flag();
8541 new_quads[i]->clear_user_data();
8542 new_quads[i]->clear_children();
8544 new_quads[i]->set_boundary_id_internal(
8548 new_quads[i]->set_manifold_id(hex->manifold_id());
8551 for (
unsigned int j = 0;
8552 j < GeometryInfo<dim>::lines_per_face;
8554 new_quads[i]->set_line_orientation(j,
true);
8563 new_hexes(n_new_hexes);
8564 for (
unsigned int i = 0; i < n_new_hexes; ++i)
8568 triangulation.
levels[level + 1]->cells.next_free_hex(
8569 triangulation, level + 1);
8573 new_hexes[i] = next_unused_hex;
8576 new_hexes[i]->set_used_flag();
8577 new_hexes[i]->clear_user_flag();
8578 new_hexes[i]->clear_user_data();
8579 new_hexes[i]->clear_children();
8582 new_hexes[i]->set_material_id(hex->material_id());
8583 new_hexes[i]->set_manifold_id(hex->manifold_id());
8584 new_hexes[i]->set_subdomain_id(subdomainid);
8587 new_hexes[i]->set_parent(hex->index());
8599 for (
const unsigned int f :
8601 new_hexes[i]->set_combined_face_orientation(
8606 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
8607 hex->set_children(2 * i, new_hexes[2 * i]->index());
8614 const bool f_or[6] = {hex->face_orientation(0),
8615 hex->face_orientation(1),
8616 hex->face_orientation(2),
8617 hex->face_orientation(3),
8618 hex->face_orientation(4),
8619 hex->face_orientation(5)};
8622 const bool f_fl[6] = {hex->face_flip(0),
8630 const bool f_ro[6] = {hex->face_rotation(0),
8631 hex->face_rotation(1),
8632 hex->face_rotation(2),
8633 hex->face_rotation(3),
8634 hex->face_rotation(4),
8635 hex->face_rotation(5)};
8638 const unsigned char f_co[6] = {
8639 hex->combined_face_orientation(0),
8640 hex->combined_face_orientation(1),
8641 hex->combined_face_orientation(2),
8642 hex->combined_face_orientation(3),
8643 hex->combined_face_orientation(4),
8644 hex->combined_face_orientation(5)};
8654 const unsigned int child_at_origin[2][2][2] = {
8741 raw_line_iterator lines[4] = {
8742 hex->face(2)->child(0)->line(
8743 (hex->face(2)->refinement_case() ==
8747 hex->face(3)->child(0)->line(
8748 (hex->face(3)->refinement_case() ==
8752 hex->face(4)->child(0)->line(
8753 (hex->face(4)->refinement_case() ==
8757 hex->face(5)->child(0)->line(
8758 (hex->face(5)->refinement_case() ==
8764 unsigned int line_indices[4];
8765 for (
unsigned int i = 0; i < 4; ++i)
8766 line_indices[i] = lines[i]->
index();
8773 bool line_orientation[4];
8779 const unsigned int middle_vertices[2] = {
8780 hex->line(2)->child(0)->vertex_index(1),
8781 hex->line(7)->child(0)->vertex_index(1)};
8783 for (
unsigned int i = 0; i < 4; ++i)
8784 if (lines[i]->vertex_index(i % 2) ==
8785 middle_vertices[i % 2])
8786 line_orientation[i] =
true;
8791 Assert(lines[i]->vertex_index((i + 1) % 2) ==
8792 middle_vertices[i % 2],
8794 line_orientation[i] =
false;
8799 new_quads[0]->set_bounding_object_indices(
8805 new_quads[0]->set_line_orientation(
8806 0, line_orientation[0]);
8807 new_quads[0]->set_line_orientation(
8808 1, line_orientation[1]);
8809 new_quads[0]->set_line_orientation(
8810 2, line_orientation[2]);
8811 new_quads[0]->set_line_orientation(
8812 3, line_orientation[3]);
8845 const int quad_indices[11] = {
8846 new_quads[0]->index(),
8848 hex->face(0)->index(),
8850 hex->face(1)->index(),
8852 hex->face(2)->child_index(
8853 child_at_origin[hex->face(2)->refinement_case() -
8854 1][f_fl[2]][f_ro[2]]),
8855 hex->face(2)->child_index(
8857 child_at_origin[hex->face(2)->refinement_case() -
8858 1][f_fl[2]][f_ro[2]]),
8860 hex->face(3)->child_index(
8861 child_at_origin[hex->face(3)->refinement_case() -
8862 1][f_fl[3]][f_ro[3]]),
8863 hex->face(3)->child_index(
8865 child_at_origin[hex->face(3)->refinement_case() -
8866 1][f_fl[3]][f_ro[3]]),
8868 hex->face(4)->child_index(
8869 child_at_origin[hex->face(4)->refinement_case() -
8870 1][f_fl[4]][f_ro[4]]),
8871 hex->face(4)->child_index(
8873 child_at_origin[hex->face(4)->refinement_case() -
8874 1][f_fl[4]][f_ro[4]]),
8876 hex->face(5)->child_index(
8877 child_at_origin[hex->face(5)->refinement_case() -
8878 1][f_fl[5]][f_ro[5]]),
8879 hex->face(5)->child_index(
8881 child_at_origin[hex->face(5)->refinement_case() -
8882 1][f_fl[5]][f_ro[5]])
8886 new_hexes[0]->set_bounding_object_indices(
8893 new_hexes[1]->set_bounding_object_indices(
8968 raw_line_iterator lines[4] = {
8969 hex->face(0)->child(0)->line(
8970 (hex->face(0)->refinement_case() ==
8974 hex->face(1)->child(0)->line(
8975 (hex->face(1)->refinement_case() ==
8979 hex->face(4)->child(0)->line(
8980 (hex->face(4)->refinement_case() ==
8984 hex->face(5)->child(0)->line(
8985 (hex->face(5)->refinement_case() ==
8991 unsigned int line_indices[4];
8992 for (
unsigned int i = 0; i < 4; ++i)
8993 line_indices[i] = lines[i]->
index();
9000 bool line_orientation[4];
9006 const unsigned int middle_vertices[2] = {
9007 hex->line(0)->child(0)->vertex_index(1),
9008 hex->line(5)->child(0)->vertex_index(1)};
9010 for (
unsigned int i = 0; i < 4; ++i)
9011 if (lines[i]->vertex_index(i % 2) ==
9012 middle_vertices[i % 2])
9013 line_orientation[i] =
true;
9017 Assert(lines[i]->vertex_index((i + 1) % 2) ==
9018 middle_vertices[i % 2],
9020 line_orientation[i] =
false;
9025 new_quads[0]->set_bounding_object_indices(
9031 new_quads[0]->set_line_orientation(
9032 0, line_orientation[2]);
9033 new_quads[0]->set_line_orientation(
9034 1, line_orientation[3]);
9035 new_quads[0]->set_line_orientation(
9036 2, line_orientation[0]);
9037 new_quads[0]->set_line_orientation(
9038 3, line_orientation[1]);
9071 const int quad_indices[11] = {
9072 new_quads[0]->index(),
9074 hex->face(0)->child_index(
9075 child_at_origin[hex->face(0)->refinement_case() -
9076 1][f_fl[0]][f_ro[0]]),
9077 hex->face(0)->child_index(
9079 child_at_origin[hex->face(0)->refinement_case() -
9080 1][f_fl[0]][f_ro[0]]),
9082 hex->face(1)->child_index(
9083 child_at_origin[hex->face(1)->refinement_case() -
9084 1][f_fl[1]][f_ro[1]]),
9085 hex->face(1)->child_index(
9087 child_at_origin[hex->face(1)->refinement_case() -
9088 1][f_fl[1]][f_ro[1]]),
9090 hex->face(2)->index(),
9092 hex->face(3)->index(),
9094 hex->face(4)->child_index(
9095 child_at_origin[hex->face(4)->refinement_case() -
9096 1][f_fl[4]][f_ro[4]]),
9097 hex->face(4)->child_index(
9099 child_at_origin[hex->face(4)->refinement_case() -
9100 1][f_fl[4]][f_ro[4]]),
9102 hex->face(5)->child_index(
9103 child_at_origin[hex->face(5)->refinement_case() -
9104 1][f_fl[5]][f_ro[5]]),
9105 hex->face(5)->child_index(
9107 child_at_origin[hex->face(5)->refinement_case() -
9108 1][f_fl[5]][f_ro[5]])
9112 new_hexes[0]->set_bounding_object_indices(
9119 new_hexes[1]->set_bounding_object_indices(
9196 raw_line_iterator lines[4] = {
9197 hex->face(0)->child(0)->line(
9198 (hex->face(0)->refinement_case() ==
9202 hex->face(1)->child(0)->line(
9203 (hex->face(1)->refinement_case() ==
9207 hex->face(2)->child(0)->line(
9208 (hex->face(2)->refinement_case() ==
9212 hex->face(3)->child(0)->line(
9213 (hex->face(3)->refinement_case() ==
9219 unsigned int line_indices[4];
9220 for (
unsigned int i = 0; i < 4; ++i)
9221 line_indices[i] = lines[i]->
index();
9228 bool line_orientation[4];
9234 const unsigned int middle_vertices[2] = {
9235 middle_vertex_index<dim, spacedim>(hex->line(8)),
9236 middle_vertex_index<dim, spacedim>(hex->line(11))};
9238 for (
unsigned int i = 0; i < 4; ++i)
9239 if (lines[i]->vertex_index(i % 2) ==
9240 middle_vertices[i % 2])
9241 line_orientation[i] =
true;
9245 Assert(lines[i]->vertex_index((i + 1) % 2) ==
9246 middle_vertices[i % 2],
9248 line_orientation[i] =
false;
9253 new_quads[0]->set_bounding_object_indices(
9259 new_quads[0]->set_line_orientation(
9260 0, line_orientation[0]);
9261 new_quads[0]->set_line_orientation(
9262 1, line_orientation[1]);
9263 new_quads[0]->set_line_orientation(
9264 2, line_orientation[2]);
9265 new_quads[0]->set_line_orientation(
9266 3, line_orientation[3]);
9300 const int quad_indices[11] = {
9301 new_quads[0]->index(),
9303 hex->face(0)->child_index(
9304 child_at_origin[hex->face(0)->refinement_case() -
9305 1][f_fl[0]][f_ro[0]]),
9306 hex->face(0)->child_index(
9308 child_at_origin[hex->face(0)->refinement_case() -
9309 1][f_fl[0]][f_ro[0]]),
9311 hex->face(1)->child_index(
9312 child_at_origin[hex->face(1)->refinement_case() -
9313 1][f_fl[1]][f_ro[1]]),
9314 hex->face(1)->child_index(
9316 child_at_origin[hex->face(1)->refinement_case() -
9317 1][f_fl[1]][f_ro[1]]),
9319 hex->face(2)->child_index(
9320 child_at_origin[hex->face(2)->refinement_case() -
9321 1][f_fl[2]][f_ro[2]]),
9322 hex->face(2)->child_index(
9324 child_at_origin[hex->face(2)->refinement_case() -
9325 1][f_fl[2]][f_ro[2]]),
9327 hex->face(3)->child_index(
9328 child_at_origin[hex->face(3)->refinement_case() -
9329 1][f_fl[3]][f_ro[3]]),
9330 hex->face(3)->child_index(
9332 child_at_origin[hex->face(3)->refinement_case() -
9333 1][f_fl[3]][f_ro[3]]),
9335 hex->face(4)->index(),
9337 hex->face(5)->index()
9340 new_hexes[0]->set_bounding_object_indices(
9347 new_hexes[1]->set_bounding_object_indices(
9379 new_lines[0]->set_bounding_object_indices(
9380 {middle_vertex_index<dim, spacedim>(hex->face(4)),
9381 middle_vertex_index<dim, spacedim>(hex->face(5))});
9449 spacedim>::raw_line_iterator lines[13] = {
9450 hex->face(0)->child(0)->line(
9451 (hex->face(0)->refinement_case() ==
9455 hex->face(1)->child(0)->line(
9456 (hex->face(1)->refinement_case() ==
9460 hex->face(2)->child(0)->line(
9461 (hex->face(2)->refinement_case() ==
9465 hex->face(3)->child(0)->line(
9466 (hex->face(3)->refinement_case() ==
9474 0, f_or[4], f_fl[4], f_ro[4]))
9477 1, f_or[4], f_fl[4], f_ro[4])),
9481 3, f_or[4], f_fl[4], f_ro[4]))
9484 0, f_or[4], f_fl[4], f_ro[4])),
9488 0, f_or[4], f_fl[4], f_ro[4]))
9491 3, f_or[4], f_fl[4], f_ro[4])),
9495 3, f_or[4], f_fl[4], f_ro[4]))
9498 2, f_or[4], f_fl[4], f_ro[4])),
9503 0, f_or[5], f_fl[5], f_ro[5]))
9506 1, f_or[5], f_fl[5], f_ro[5])),
9510 3, f_or[5], f_fl[5], f_ro[5]))
9513 0, f_or[5], f_fl[5], f_ro[5])),
9517 0, f_or[5], f_fl[5], f_ro[5]))
9520 3, f_or[5], f_fl[5], f_ro[5])),
9524 3, f_or[5], f_fl[5], f_ro[5]))
9527 2, f_or[5], f_fl[5], f_ro[5])),
9532 unsigned int line_indices[13];
9533 for (
unsigned int i = 0; i < 13; ++i)
9534 line_indices[i] = lines[i]->
index();
9541 bool line_orientation[13];
9545 const unsigned int middle_vertices[4] = {
9546 hex->line(0)->child(0)->vertex_index(1),
9547 hex->line(1)->child(0)->vertex_index(1),
9548 hex->line(2)->child(0)->vertex_index(1),
9549 hex->line(3)->child(0)->vertex_index(1),
9555 for (
unsigned int i = 0; i < 4; ++i)
9556 if (lines[i]->vertex_index(0) == middle_vertices[i])
9557 line_orientation[i] =
true;
9561 Assert(lines[i]->vertex_index(1) ==
9564 line_orientation[i] =
false;
9573 for (
unsigned int i = 4; i < 12; ++i)
9574 if (lines[i]->vertex_index((i + 1) % 2) ==
9575 middle_vertex_index<dim, spacedim>(
9576 hex->face(3 + i / 4)))
9577 line_orientation[i] =
true;
9582 Assert(lines[i]->vertex_index(i % 2) ==
9583 (middle_vertex_index<dim, spacedim>(
9584 hex->face(3 + i / 4))),
9586 line_orientation[i] =
false;
9591 line_orientation[12] =
true;
9616 new_quads[0]->set_bounding_object_indices(
9621 new_quads[1]->set_bounding_object_indices(
9626 new_quads[2]->set_bounding_object_indices(
9631 new_quads[3]->set_bounding_object_indices(
9637 new_quads[0]->set_line_orientation(
9638 0, line_orientation[2]);
9639 new_quads[0]->set_line_orientation(
9640 2, line_orientation[4]);
9641 new_quads[0]->set_line_orientation(
9642 3, line_orientation[8]);
9644 new_quads[1]->set_line_orientation(
9645 1, line_orientation[3]);
9646 new_quads[1]->set_line_orientation(
9647 2, line_orientation[5]);
9648 new_quads[1]->set_line_orientation(
9649 3, line_orientation[9]);
9651 new_quads[2]->set_line_orientation(
9652 0, line_orientation[6]);
9653 new_quads[2]->set_line_orientation(
9654 1, line_orientation[10]);
9655 new_quads[2]->set_line_orientation(
9656 2, line_orientation[0]);
9658 new_quads[3]->set_line_orientation(
9659 0, line_orientation[7]);
9660 new_quads[3]->set_line_orientation(
9661 1, line_orientation[11]);
9662 new_quads[3]->set_line_orientation(
9663 3, line_orientation[1]);
9696 const int quad_indices[20] = {
9697 new_quads[0]->index(),
9698 new_quads[1]->index(),
9699 new_quads[2]->index(),
9700 new_quads[3]->index(),
9702 hex->face(0)->child_index(
9703 child_at_origin[hex->face(0)->refinement_case() -
9704 1][f_fl[0]][f_ro[0]]),
9705 hex->face(0)->child_index(
9707 child_at_origin[hex->face(0)->refinement_case() -
9708 1][f_fl[0]][f_ro[0]]),
9710 hex->face(1)->child_index(
9711 child_at_origin[hex->face(1)->refinement_case() -
9712 1][f_fl[1]][f_ro[1]]),
9713 hex->face(1)->child_index(
9715 child_at_origin[hex->face(1)->refinement_case() -
9716 1][f_fl[1]][f_ro[1]]),
9718 hex->face(2)->child_index(
9719 child_at_origin[hex->face(2)->refinement_case() -
9720 1][f_fl[2]][f_ro[2]]),
9721 hex->face(2)->child_index(
9723 child_at_origin[hex->face(2)->refinement_case() -
9724 1][f_fl[2]][f_ro[2]]),
9726 hex->face(3)->child_index(
9727 child_at_origin[hex->face(3)->refinement_case() -
9728 1][f_fl[3]][f_ro[3]]),
9729 hex->face(3)->child_index(
9731 child_at_origin[hex->face(3)->refinement_case() -
9732 1][f_fl[3]][f_ro[3]]),
9734 hex->face(4)->isotropic_child_index(
9736 0, f_or[4], f_fl[4], f_ro[4])),
9737 hex->face(4)->isotropic_child_index(
9739 1, f_or[4], f_fl[4], f_ro[4])),
9740 hex->face(4)->isotropic_child_index(
9742 2, f_or[4], f_fl[4], f_ro[4])),
9743 hex->face(4)->isotropic_child_index(
9745 3, f_or[4], f_fl[4], f_ro[4])),
9747 hex->face(5)->isotropic_child_index(
9749 0, f_or[5], f_fl[5], f_ro[5])),
9750 hex->face(5)->isotropic_child_index(
9752 1, f_or[5], f_fl[5], f_ro[5])),
9753 hex->face(5)->isotropic_child_index(
9755 2, f_or[5], f_fl[5], f_ro[5])),
9756 hex->face(5)->isotropic_child_index(
9758 3, f_or[5], f_fl[5], f_ro[5]))};
9760 new_hexes[0]->set_bounding_object_indices(
9767 new_hexes[1]->set_bounding_object_indices(
9774 new_hexes[2]->set_bounding_object_indices(
9781 new_hexes[3]->set_bounding_object_indices(
9813 new_lines[0]->set_bounding_object_indices(
9814 {middle_vertex_index<dim, spacedim>(hex->face(2)),
9815 middle_vertex_index<dim, spacedim>(hex->face(3))});
9883 spacedim>::raw_line_iterator lines[13] = {
9884 hex->face(0)->child(0)->line(
9885 (hex->face(0)->refinement_case() ==
9889 hex->face(1)->child(0)->line(
9890 (hex->face(1)->refinement_case() ==
9894 hex->face(4)->child(0)->line(
9895 (hex->face(4)->refinement_case() ==
9899 hex->face(5)->child(0)->line(
9900 (hex->face(5)->refinement_case() ==
9908 0, f_or[2], f_fl[2], f_ro[2]))
9911 3, f_or[2], f_fl[2], f_ro[2])),
9915 3, f_or[2], f_fl[2], f_ro[2]))
9918 2, f_or[2], f_fl[2], f_ro[2])),
9922 0, f_or[2], f_fl[2], f_ro[2]))
9925 1, f_or[2], f_fl[2], f_ro[2])),
9929 3, f_or[2], f_fl[2], f_ro[2]))
9932 0, f_or[2], f_fl[2], f_ro[2])),
9937 0, f_or[3], f_fl[3], f_ro[3]))
9940 3, f_or[3], f_fl[3], f_ro[3])),
9944 3, f_or[3], f_fl[3], f_ro[3]))
9947 2, f_or[3], f_fl[3], f_ro[3])),
9951 0, f_or[3], f_fl[3], f_ro[3]))
9954 1, f_or[3], f_fl[3], f_ro[3])),
9958 3, f_or[3], f_fl[3], f_ro[3]))
9961 0, f_or[3], f_fl[3], f_ro[3])),
9966 unsigned int line_indices[13];
9967 for (
unsigned int i = 0; i < 13; ++i)
9968 line_indices[i] = lines[i]->
index();
9975 bool line_orientation[13];
9979 const unsigned int middle_vertices[4] = {
9980 hex->line(8)->child(0)->vertex_index(1),
9981 hex->line(9)->child(0)->vertex_index(1),
9982 hex->line(2)->child(0)->vertex_index(1),
9983 hex->line(6)->child(0)->vertex_index(1),
9988 for (
unsigned int i = 0; i < 4; ++i)
9989 if (lines[i]->vertex_index(0) == middle_vertices[i])
9990 line_orientation[i] =
true;
9994 Assert(lines[i]->vertex_index(1) ==
9997 line_orientation[i] =
false;
10006 for (
unsigned int i = 4; i < 12; ++i)
10007 if (lines[i]->vertex_index((i + 1) % 2) ==
10008 middle_vertex_index<dim, spacedim>(
10009 hex->face(1 + i / 4)))
10010 line_orientation[i] =
true;
10015 Assert(lines[i]->vertex_index(i % 2) ==
10016 (middle_vertex_index<dim, spacedim>(
10017 hex->face(1 + i / 4))),
10019 line_orientation[i] =
false;
10024 line_orientation[12] =
true;
10050 new_quads[0]->set_bounding_object_indices(
10054 line_indices[10]});
10055 new_quads[1]->set_bounding_object_indices(
10059 line_indices[11]});
10060 new_quads[2]->set_bounding_object_indices(
10064 line_indices[12]});
10065 new_quads[3]->set_bounding_object_indices(
10071 new_quads[0]->set_line_orientation(
10072 0, line_orientation[0]);
10073 new_quads[0]->set_line_orientation(
10074 2, line_orientation[6]);
10075 new_quads[0]->set_line_orientation(
10076 3, line_orientation[10]);
10078 new_quads[1]->set_line_orientation(
10079 1, line_orientation[1]);
10080 new_quads[1]->set_line_orientation(
10081 2, line_orientation[7]);
10082 new_quads[1]->set_line_orientation(
10083 3, line_orientation[11]);
10085 new_quads[2]->set_line_orientation(
10086 0, line_orientation[4]);
10087 new_quads[2]->set_line_orientation(
10088 1, line_orientation[8]);
10089 new_quads[2]->set_line_orientation(
10090 2, line_orientation[2]);
10092 new_quads[3]->set_line_orientation(
10093 0, line_orientation[5]);
10094 new_quads[3]->set_line_orientation(
10095 1, line_orientation[9]);
10096 new_quads[3]->set_line_orientation(
10097 3, line_orientation[3]);
10131 const int quad_indices[20] = {
10132 new_quads[0]->index(),
10133 new_quads[1]->index(),
10134 new_quads[2]->index(),
10135 new_quads[3]->index(),
10137 hex->face(0)->child_index(
10138 child_at_origin[hex->face(0)->refinement_case() -
10139 1][f_fl[0]][f_ro[0]]),
10140 hex->face(0)->child_index(
10142 child_at_origin[hex->face(0)->refinement_case() -
10143 1][f_fl[0]][f_ro[0]]),
10145 hex->face(1)->child_index(
10146 child_at_origin[hex->face(1)->refinement_case() -
10147 1][f_fl[1]][f_ro[1]]),
10148 hex->face(1)->child_index(
10150 child_at_origin[hex->face(1)->refinement_case() -
10151 1][f_fl[1]][f_ro[1]]),
10153 hex->face(2)->isotropic_child_index(
10155 0, f_or[2], f_fl[2], f_ro[2])),
10156 hex->face(2)->isotropic_child_index(
10158 1, f_or[2], f_fl[2], f_ro[2])),
10159 hex->face(2)->isotropic_child_index(
10161 2, f_or[2], f_fl[2], f_ro[2])),
10162 hex->face(2)->isotropic_child_index(
10164 3, f_or[2], f_fl[2], f_ro[2])),
10166 hex->face(3)->isotropic_child_index(
10168 0, f_or[3], f_fl[3], f_ro[3])),
10169 hex->face(3)->isotropic_child_index(
10171 1, f_or[3], f_fl[3], f_ro[3])),
10172 hex->face(3)->isotropic_child_index(
10174 2, f_or[3], f_fl[3], f_ro[3])),
10175 hex->face(3)->isotropic_child_index(
10177 3, f_or[3], f_fl[3], f_ro[3])),
10179 hex->face(4)->child_index(
10180 child_at_origin[hex->face(4)->refinement_case() -
10181 1][f_fl[4]][f_ro[4]]),
10182 hex->face(4)->child_index(
10184 child_at_origin[hex->face(4)->refinement_case() -
10185 1][f_fl[4]][f_ro[4]]),
10187 hex->face(5)->child_index(
10188 child_at_origin[hex->face(5)->refinement_case() -
10189 1][f_fl[5]][f_ro[5]]),
10190 hex->face(5)->child_index(
10192 child_at_origin[hex->face(5)->refinement_case() -
10193 1][f_fl[5]][f_ro[5]])};
10204 new_hexes[0]->set_bounding_object_indices(
10211 new_hexes[1]->set_bounding_object_indices(
10217 quad_indices[18]});
10218 new_hexes[2]->set_bounding_object_indices(
10225 new_hexes[3]->set_bounding_object_indices(
10231 quad_indices[19]});
10258 new_lines[0]->set_bounding_object_indices(
10260 {middle_vertex_index<dim, spacedim>(hex->face(0)),
10261 middle_vertex_index<dim, spacedim>(hex->face(1))});
10330 spacedim>::raw_line_iterator lines[13] = {
10331 hex->face(2)->child(0)->line(
10332 (hex->face(2)->refinement_case() ==
10336 hex->face(3)->child(0)->line(
10337 (hex->face(3)->refinement_case() ==
10341 hex->face(4)->child(0)->line(
10342 (hex->face(4)->refinement_case() ==
10346 hex->face(5)->child(0)->line(
10347 (hex->face(5)->refinement_case() ==
10355 0, f_or[0], f_fl[0], f_ro[0]))
10358 1, f_or[0], f_fl[0], f_ro[0])),
10362 3, f_or[0], f_fl[0], f_ro[0]))
10365 0, f_or[0], f_fl[0], f_ro[0])),
10369 0, f_or[0], f_fl[0], f_ro[0]))
10372 3, f_or[0], f_fl[0], f_ro[0])),
10376 3, f_or[0], f_fl[0], f_ro[0]))
10379 2, f_or[0], f_fl[0], f_ro[0])),
10384 0, f_or[1], f_fl[1], f_ro[1]))
10387 1, f_or[1], f_fl[1], f_ro[1])),
10391 3, f_or[1], f_fl[1], f_ro[1]))
10394 0, f_or[1], f_fl[1], f_ro[1])),
10398 0, f_or[1], f_fl[1], f_ro[1]))
10401 3, f_or[1], f_fl[1], f_ro[1])),
10405 3, f_or[1], f_fl[1], f_ro[1]))
10408 2, f_or[1], f_fl[1], f_ro[1])),
10413 unsigned int line_indices[13];
10415 for (
unsigned int i = 0; i < 13; ++i)
10416 line_indices[i] = lines[i]->
index();
10423 bool line_orientation[13];
10427 const unsigned int middle_vertices[4] = {
10428 hex->line(8)->child(0)->vertex_index(1),
10429 hex->line(10)->child(0)->vertex_index(1),
10430 hex->line(0)->child(0)->vertex_index(1),
10431 hex->line(4)->child(0)->vertex_index(1),
10436 for (
unsigned int i = 0; i < 4; ++i)
10437 if (lines[i]->vertex_index(0) == middle_vertices[i])
10438 line_orientation[i] =
true;
10442 Assert(lines[i]->vertex_index(1) ==
10443 middle_vertices[i],
10445 line_orientation[i] =
false;
10454 for (
unsigned int i = 4; i < 12; ++i)
10455 if (lines[i]->vertex_index((i + 1) % 2) ==
10456 middle_vertex_index<dim, spacedim>(
10457 hex->face(i / 4 - 1)))
10458 line_orientation[i] =
true;
10463 Assert(lines[i]->vertex_index(i % 2) ==
10464 (middle_vertex_index<dim, spacedim>(
10465 hex->face(i / 4 - 1))),
10467 line_orientation[i] =
false;
10472 line_orientation[12] =
true;
10492 new_quads[0]->set_bounding_object_indices(
10496 line_indices[12]});
10497 new_quads[1]->set_bounding_object_indices(
10502 new_quads[2]->set_bounding_object_indices(
10507 new_quads[3]->set_bounding_object_indices(
10513 new_quads[0]->set_line_orientation(
10514 0, line_orientation[6]);
10515 new_quads[0]->set_line_orientation(
10516 1, line_orientation[10]);
10517 new_quads[0]->set_line_orientation(
10518 2, line_orientation[0]);
10520 new_quads[1]->set_line_orientation(
10521 0, line_orientation[7]);
10522 new_quads[1]->set_line_orientation(
10523 1, line_orientation[11]);
10524 new_quads[1]->set_line_orientation(
10525 3, line_orientation[1]);
10527 new_quads[2]->set_line_orientation(
10528 0, line_orientation[2]);
10529 new_quads[2]->set_line_orientation(
10530 2, line_orientation[4]);
10531 new_quads[2]->set_line_orientation(
10532 3, line_orientation[8]);
10534 new_quads[3]->set_line_orientation(
10535 1, line_orientation[3]);
10536 new_quads[3]->set_line_orientation(
10537 2, line_orientation[5]);
10538 new_quads[3]->set_line_orientation(
10539 3, line_orientation[9]);
10573 const int quad_indices[20] = {
10574 new_quads[0]->index(),
10575 new_quads[1]->index(),
10576 new_quads[2]->index(),
10577 new_quads[3]->index(),
10579 hex->face(0)->isotropic_child_index(
10581 0, f_or[0], f_fl[0], f_ro[0])),
10582 hex->face(0)->isotropic_child_index(
10584 1, f_or[0], f_fl[0], f_ro[0])),
10585 hex->face(0)->isotropic_child_index(
10587 2, f_or[0], f_fl[0], f_ro[0])),
10588 hex->face(0)->isotropic_child_index(
10590 3, f_or[0], f_fl[0], f_ro[0])),
10592 hex->face(1)->isotropic_child_index(
10594 0, f_or[1], f_fl[1], f_ro[1])),
10595 hex->face(1)->isotropic_child_index(
10597 1, f_or[1], f_fl[1], f_ro[1])),
10598 hex->face(1)->isotropic_child_index(
10600 2, f_or[1], f_fl[1], f_ro[1])),
10601 hex->face(1)->isotropic_child_index(
10603 3, f_or[1], f_fl[1], f_ro[1])),
10605 hex->face(2)->child_index(
10606 child_at_origin[hex->face(2)->refinement_case() -
10607 1][f_fl[2]][f_ro[2]]),
10608 hex->face(2)->child_index(
10610 child_at_origin[hex->face(2)->refinement_case() -
10611 1][f_fl[2]][f_ro[2]]),
10613 hex->face(3)->child_index(
10614 child_at_origin[hex->face(3)->refinement_case() -
10615 1][f_fl[3]][f_ro[3]]),
10616 hex->face(3)->child_index(
10618 child_at_origin[hex->face(3)->refinement_case() -
10619 1][f_fl[3]][f_ro[3]]),
10621 hex->face(4)->child_index(
10622 child_at_origin[hex->face(4)->refinement_case() -
10623 1][f_fl[4]][f_ro[4]]),
10624 hex->face(4)->child_index(
10626 child_at_origin[hex->face(4)->refinement_case() -
10627 1][f_fl[4]][f_ro[4]]),
10629 hex->face(5)->child_index(
10630 child_at_origin[hex->face(5)->refinement_case() -
10631 1][f_fl[5]][f_ro[5]]),
10632 hex->face(5)->child_index(
10634 child_at_origin[hex->face(5)->refinement_case() -
10635 1][f_fl[5]][f_ro[5]])};
10637 new_hexes[0]->set_bounding_object_indices(
10644 new_hexes[1]->set_bounding_object_indices(
10651 new_hexes[2]->set_bounding_object_indices(
10657 quad_indices[18]});
10658 new_hexes[3]->set_bounding_object_indices(
10664 quad_indices[19]});
10696 ++next_unused_vertex;
10698 next_unused_vertex < triangulation.
vertices.size(),
10700 "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."));
10710 triangulation.
vertices[next_unused_vertex] =
10711 hex->center(
true,
true);
10732 const unsigned int vertex_indices[7] = {
10733 middle_vertex_index<dim, spacedim>(hex->face(0)),
10734 middle_vertex_index<dim, spacedim>(hex->face(1)),
10735 middle_vertex_index<dim, spacedim>(hex->face(2)),
10736 middle_vertex_index<dim, spacedim>(hex->face(3)),
10737 middle_vertex_index<dim, spacedim>(hex->face(4)),
10738 middle_vertex_index<dim, spacedim>(hex->face(5)),
10739 next_unused_vertex};
10741 new_lines[0]->set_bounding_object_indices(
10742 {vertex_indices[2], vertex_indices[6]});
10743 new_lines[1]->set_bounding_object_indices(
10744 {vertex_indices[6], vertex_indices[3]});
10745 new_lines[2]->set_bounding_object_indices(
10746 {vertex_indices[0], vertex_indices[6]});
10747 new_lines[3]->set_bounding_object_indices(
10748 {vertex_indices[6], vertex_indices[1]});
10749 new_lines[4]->set_bounding_object_indices(
10750 {vertex_indices[4], vertex_indices[6]});
10751 new_lines[5]->set_bounding_object_indices(
10752 {vertex_indices[6], vertex_indices[5]});
10821 spacedim>::raw_line_iterator lines[30] = {
10825 0, f_or[0], f_fl[0], f_ro[0]))
10828 1, f_or[0], f_fl[0], f_ro[0])),
10832 3, f_or[0], f_fl[0], f_ro[0]))
10835 0, f_or[0], f_fl[0], f_ro[0])),
10839 0, f_or[0], f_fl[0], f_ro[0]))
10842 3, f_or[0], f_fl[0], f_ro[0])),
10846 3, f_or[0], f_fl[0], f_ro[0]))
10849 2, f_or[0], f_fl[0], f_ro[0])),
10854 0, f_or[1], f_fl[1], f_ro[1]))
10857 1, f_or[1], f_fl[1], f_ro[1])),
10861 3, f_or[1], f_fl[1], f_ro[1]))
10864 0, f_or[1], f_fl[1], f_ro[1])),
10868 0, f_or[1], f_fl[1], f_ro[1]))
10871 3, f_or[1], f_fl[1], f_ro[1])),
10875 3, f_or[1], f_fl[1], f_ro[1]))
10878 2, f_or[1], f_fl[1], f_ro[1])),
10883 0, f_or[2], f_fl[2], f_ro[2]))
10886 1, f_or[2], f_fl[2], f_ro[2])),
10890 3, f_or[2], f_fl[2], f_ro[2]))
10893 0, f_or[2], f_fl[2], f_ro[2])),
10897 0, f_or[2], f_fl[2], f_ro[2]))
10900 3, f_or[2], f_fl[2], f_ro[2])),
10904 3, f_or[2], f_fl[2], f_ro[2]))
10907 2, f_or[2], f_fl[2], f_ro[2])),
10912 0, f_or[3], f_fl[3], f_ro[3]))
10915 1, f_or[3], f_fl[3], f_ro[3])),
10919 3, f_or[3], f_fl[3], f_ro[3]))
10922 0, f_or[3], f_fl[3], f_ro[3])),
10926 0, f_or[3], f_fl[3], f_ro[3]))
10929 3, f_or[3], f_fl[3], f_ro[3])),
10933 3, f_or[3], f_fl[3], f_ro[3]))
10936 2, f_or[3], f_fl[3], f_ro[3])),
10941 0, f_or[4], f_fl[4], f_ro[4]))
10944 1, f_or[4], f_fl[4], f_ro[4])),
10948 3, f_or[4], f_fl[4], f_ro[4]))
10951 0, f_or[4], f_fl[4], f_ro[4])),
10955 0, f_or[4], f_fl[4], f_ro[4]))
10958 3, f_or[4], f_fl[4], f_ro[4])),
10962 3, f_or[4], f_fl[4], f_ro[4]))
10965 2, f_or[4], f_fl[4], f_ro[4])),
10970 0, f_or[5], f_fl[5], f_ro[5]))
10973 1, f_or[5], f_fl[5], f_ro[5])),
10977 3, f_or[5], f_fl[5], f_ro[5]))
10980 0, f_or[5], f_fl[5], f_ro[5])),
10984 0, f_or[5], f_fl[5], f_ro[5]))
10987 3, f_or[5], f_fl[5], f_ro[5])),
10991 3, f_or[5], f_fl[5], f_ro[5]))
10994 2, f_or[5], f_fl[5], f_ro[5])),
11004 unsigned int line_indices[30];
11005 for (
unsigned int i = 0; i < 30; ++i)
11006 line_indices[i] = lines[i]->
index();
11013 bool line_orientation[30];
11021 for (
unsigned int i = 0; i < 24; ++i)
11022 if (lines[i]->vertex_index((i + 1) % 2) ==
11023 vertex_indices[i / 4])
11024 line_orientation[i] =
true;
11029 Assert(lines[i]->vertex_index(i % 2) ==
11030 vertex_indices[i / 4],
11032 line_orientation[i] =
false;
11037 for (
unsigned int i = 24; i < 30; ++i)
11038 line_orientation[i] =
true;
11070 new_quads[0]->set_bounding_object_indices(
11074 line_indices[24]});
11075 new_quads[1]->set_bounding_object_indices(
11079 line_indices[25]});
11080 new_quads[2]->set_bounding_object_indices(
11084 line_indices[20]});
11085 new_quads[3]->set_bounding_object_indices(
11089 line_indices[21]});
11090 new_quads[4]->set_bounding_object_indices(
11094 line_indices[28]});
11095 new_quads[5]->set_bounding_object_indices(
11099 line_indices[29]});
11100 new_quads[6]->set_bounding_object_indices(
11105 new_quads[7]->set_bounding_object_indices(
11110 new_quads[8]->set_bounding_object_indices(
11114 line_indices[26]});
11115 new_quads[9]->set_bounding_object_indices(
11119 line_indices[27]});
11120 new_quads[10]->set_bounding_object_indices(
11124 line_indices[12]});
11125 new_quads[11]->set_bounding_object_indices(
11129 line_indices[13]});
11134 new_quads[0]->set_line_orientation(
11135 0, line_orientation[10]);
11136 new_quads[0]->set_line_orientation(
11137 2, line_orientation[16]);
11139 new_quads[1]->set_line_orientation(
11140 1, line_orientation[14]);
11141 new_quads[1]->set_line_orientation(
11142 2, line_orientation[17]);
11144 new_quads[2]->set_line_orientation(
11145 0, line_orientation[11]);
11146 new_quads[2]->set_line_orientation(
11147 3, line_orientation[20]);
11149 new_quads[3]->set_line_orientation(
11150 1, line_orientation[15]);
11151 new_quads[3]->set_line_orientation(
11152 3, line_orientation[21]);
11154 new_quads[4]->set_line_orientation(
11155 0, line_orientation[18]);
11156 new_quads[4]->set_line_orientation(
11157 2, line_orientation[0]);
11159 new_quads[5]->set_line_orientation(
11160 1, line_orientation[22]);
11161 new_quads[5]->set_line_orientation(
11162 2, line_orientation[1]);
11164 new_quads[6]->set_line_orientation(
11165 0, line_orientation[19]);
11166 new_quads[6]->set_line_orientation(
11167 3, line_orientation[4]);
11169 new_quads[7]->set_line_orientation(
11170 1, line_orientation[23]);
11171 new_quads[7]->set_line_orientation(
11172 3, line_orientation[5]);
11174 new_quads[8]->set_line_orientation(
11175 0, line_orientation[2]);
11176 new_quads[8]->set_line_orientation(
11177 2, line_orientation[8]);
11179 new_quads[9]->set_line_orientation(
11180 1, line_orientation[6]);
11181 new_quads[9]->set_line_orientation(
11182 2, line_orientation[9]);
11184 new_quads[10]->set_line_orientation(
11185 0, line_orientation[3]);
11186 new_quads[10]->set_line_orientation(
11187 3, line_orientation[12]);
11189 new_quads[11]->set_line_orientation(
11190 1, line_orientation[7]);
11191 new_quads[11]->set_line_orientation(
11192 3, line_orientation[13]);
11233 const int quad_indices[36] = {
11234 new_quads[0]->index(),
11235 new_quads[1]->index(),
11236 new_quads[2]->index(),
11237 new_quads[3]->index(),
11238 new_quads[4]->index(),
11239 new_quads[5]->index(),
11240 new_quads[6]->index(),
11241 new_quads[7]->index(),
11242 new_quads[8]->index(),
11243 new_quads[9]->index(),
11244 new_quads[10]->index(),
11245 new_quads[11]->index(),
11247 hex->face(0)->isotropic_child_index(
11249 0, f_or[0], f_fl[0], f_ro[0])),
11250 hex->face(0)->isotropic_child_index(
11252 1, f_or[0], f_fl[0], f_ro[0])),
11253 hex->face(0)->isotropic_child_index(
11255 2, f_or[0], f_fl[0], f_ro[0])),
11256 hex->face(0)->isotropic_child_index(
11258 3, f_or[0], f_fl[0], f_ro[0])),
11260 hex->face(1)->isotropic_child_index(
11262 0, f_or[1], f_fl[1], f_ro[1])),
11263 hex->face(1)->isotropic_child_index(
11265 1, f_or[1], f_fl[1], f_ro[1])),
11266 hex->face(1)->isotropic_child_index(
11268 2, f_or[1], f_fl[1], f_ro[1])),
11269 hex->face(1)->isotropic_child_index(
11271 3, f_or[1], f_fl[1], f_ro[1])),
11273 hex->face(2)->isotropic_child_index(
11275 0, f_or[2], f_fl[2], f_ro[2])),
11276 hex->face(2)->isotropic_child_index(
11278 1, f_or[2], f_fl[2], f_ro[2])),
11279 hex->face(2)->isotropic_child_index(
11281 2, f_or[2], f_fl[2], f_ro[2])),
11282 hex->face(2)->isotropic_child_index(
11284 3, f_or[2], f_fl[2], f_ro[2])),
11286 hex->face(3)->isotropic_child_index(
11288 0, f_or[3], f_fl[3], f_ro[3])),
11289 hex->face(3)->isotropic_child_index(
11291 1, f_or[3], f_fl[3], f_ro[3])),
11292 hex->face(3)->isotropic_child_index(
11294 2, f_or[3], f_fl[3], f_ro[3])),
11295 hex->face(3)->isotropic_child_index(
11297 3, f_or[3], f_fl[3], f_ro[3])),
11299 hex->face(4)->isotropic_child_index(
11301 0, f_or[4], f_fl[4], f_ro[4])),
11302 hex->face(4)->isotropic_child_index(
11304 1, f_or[4], f_fl[4], f_ro[4])),
11305 hex->face(4)->isotropic_child_index(
11307 2, f_or[4], f_fl[4], f_ro[4])),
11308 hex->face(4)->isotropic_child_index(
11310 3, f_or[4], f_fl[4], f_ro[4])),
11312 hex->face(5)->isotropic_child_index(
11314 0, f_or[5], f_fl[5], f_ro[5])),
11315 hex->face(5)->isotropic_child_index(
11317 1, f_or[5], f_fl[5], f_ro[5])),
11318 hex->face(5)->isotropic_child_index(
11320 2, f_or[5], f_fl[5], f_ro[5])),
11321 hex->face(5)->isotropic_child_index(
11323 3, f_or[5], f_fl[5], f_ro[5]))};
11326 new_hexes[0]->set_bounding_object_indices(
11333 new_hexes[1]->set_bounding_object_indices(
11340 new_hexes[2]->set_bounding_object_indices(
11346 quad_indices[10]});
11347 new_hexes[3]->set_bounding_object_indices(
11353 quad_indices[11]});
11356 new_hexes[4]->set_bounding_object_indices(
11362 quad_indices[32]});
11363 new_hexes[5]->set_bounding_object_indices(
11369 quad_indices[33]});
11370 new_hexes[6]->set_bounding_object_indices(
11376 quad_indices[34]});
11377 new_hexes[7]->set_bounding_object_indices(
11383 quad_indices[35]});
11417 for (
unsigned int s = 0;
11424 const unsigned int current_child =
11433 ref_case, f, f_or[f], f_fl[f], f_ro[f]));
11434 new_hexes[current_child]->set_combined_face_orientation(
11440 if (check_for_distorted_cells &&
11441 has_distorted_children<dim, spacedim>(hex))
11458 triangulation.
faces->quads.clear_user_data();
11461 return cells_with_distorted_children;
11477 template <
int spacedim>
11484 template <
int dim,
int spacedim>
11491 if (spacedim > dim)
11503 for (
const unsigned int face_no :
11505 if (cell->
face(face_no)->at_boundary())
11516 spacedim>::face_iterator
11517 face = cell->
face(face_no);
11527 .transform_real_to_unit_cell(cell, new_bound);
11537 const double dist =
11538 std::fabs(new_unit[face_no / 2] - face_no % 2);
11543 const double allowed = 0.25;
11545 if (dist > allowed)
11565 template <
int dim,
int spacedim>
11570 ExcMessage(
"Wrong function called -- there should "
11571 "be a specialization."));
11575 template <
int spacedim>
11580 const unsigned int dim = 3;
11581 using raw_line_iterator =
11586 bool mesh_changed =
false;
11590 mesh_changed =
false;
11606 const std::array<unsigned int, 12> line_indices =
11607 TriaAccessorImplementation::Implementation::
11608 get_line_indices_of_cell(*cell);
11609 for (
unsigned int l = 0; l < cell->
n_lines(); ++l)
11614 raw_line_iterator line(&triangulation,
11618 line->set_user_flag();
11624 const std::array<unsigned int, 12> line_indices =
11625 TriaAccessorImplementation::Implementation::
11626 get_line_indices_of_cell(*cell);
11627 for (
unsigned int l = 0; l < cell->
n_lines(); ++l)
11632 raw_line_iterator line(&triangulation,
11636 line->set_user_flag();
11650 cell != triangulation.
end();
11653 const std::array<unsigned int, 12> line_indices =
11654 TriaAccessorImplementation::Implementation::
11655 get_line_indices_of_cell(*cell);
11656 for (
unsigned int l = 0; l < cell->
n_lines(); ++l)
11658 raw_line_iterator line(&triangulation, 0, line_indices[l]);
11659 if (line->has_children())
11668 bool offending_line_found =
false;
11670 for (
unsigned int c = 0; c < 2; ++c)
11672 Assert(line->child(c)->has_children() ==
false,
11675 if (line->child(c)->user_flag_set() &&
11688 allow_anisotropic_smoothing)
11693 for (
unsigned int k = 0; k < cell->
n_lines();
11699 raw_line_iterator(&triangulation,
11705 offending_line_found =
true;
11711 for (
unsigned int k = 0; k < cell->
n_lines();
11715 raw_line_iterator(&triangulation,
11718 if (!line->has_children() &&
11720 line_refinement_case(
11723 line->set_user_flag();
11730 if (offending_line_found)
11732 mesh_changed =
true;
11751 triangulation.
last();
11752 cell != triangulation.
end();
11756 const std::array<unsigned int, 12> line_indices =
11757 TriaAccessorImplementation::Implementation::
11758 get_line_indices_of_cell(*cell);
11759 for (
unsigned int l = 0; l < cell->
n_lines(); ++l)
11761 raw_line_iterator line(&triangulation,
11764 if (line->has_children() &&
11765 (line->child(0)->user_flag_set() ||
11766 line->child(1)->user_flag_set()))
11768 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
11771 for (
unsigned int k = 0; k < cell->
n_lines(); ++k)
11777 raw_line_iterator(&triangulation,
11781 mesh_changed =
true;
11787 while (mesh_changed ==
true);
11798 template <
int dim,
int spacedim>
11818 const unsigned int n_subfaces =
11823 for (
unsigned int c = 0; c < n_subfaces; ++c)
11826 child = cell->
child(
11830 child_neighbor = child->
neighbor(n);
11870 template <
int spacedim>
11875 template <
int dim,
int spacedim>
11878 std::vector<std::pair<unsigned int, unsigned int>> adjacent_cells(
11880 {numbers::invalid_unsigned_int, numbers::invalid_unsigned_int});
11882 const auto set_entry = [&](
const auto &face_index,
const auto &cell) {
11883 const std::pair<unsigned int, unsigned int> cell_pair = {
11885 unsigned int index;
11887 if (adjacent_cells[2 * face_index].first ==
11889 adjacent_cells[2 * face_index].second ==
11892 index = 2 * face_index + 0;
11896 Assert(((adjacent_cells[2 * face_index + 1].first ==
11898 (adjacent_cells[2 * face_index + 1].second ==
11901 index = 2 * face_index + 1;
11904 adjacent_cells[
index] = cell_pair;
11907 const auto get_entry =
11908 [&](
const auto &face_index,
11910 auto test = adjacent_cells[2 * face_index];
11912 if (test == std::pair<unsigned int, unsigned int>(cell->
level(),
11914 test = adjacent_cells[2 * face_index + 1];
11928 set_entry(face->index(), cell);
11930 if (cell->
is_active() && face->has_children())
11931 for (
unsigned int c = 0; c < face->n_children(); ++c)
11932 set_entry(face->child(c)->index(), cell);
11940 template <
int dim,
int spacedim>
11945 std::vector<unsigned int> &line_cell_count,
11946 std::vector<unsigned int> &quad_cell_count)
11949 (void)triangulation;
11951 (void)line_cell_count;
11952 (void)quad_cell_count;
11955 template <
int dim,
int spacedim>
11958 const bool check_for_distorted_cells)
11961 triangulation, check_for_distorted_cells);
11964 template <
int dim,
int spacedim>
11970 (void)triangulation;
11973 template <
int dim,
int spacedim>
11981 template <
int dim,
int spacedim>
11993 template <
int dim,
int spacedim>
11998 return flat_manifold;
12005template <
int dim,
int spacedim>
12011template <
int dim,
int spacedim>
12014 const MeshSmoothing smooth_grid,
12015 const bool check_for_distorted_cells)
12016 : cell_attached_data({0, 0, {}, {}})
12017 , smooth_grid(smooth_grid)
12018 , anisotropic_refinement(
false)
12019 , check_for_distorted_cells(check_for_distorted_cells)
12023 vertex_to_boundary_id_map_1d =
12024 std::make_unique<std::map<unsigned int, types::boundary_id>>();
12025 vertex_to_manifold_id_map_1d =
12026 std::make_unique<std::map<unsigned int, types::manifold_id>>();
12030 signals.create.connect(signals.any_change);
12031 signals.post_refinement.connect(signals.any_change);
12032 signals.clear.connect(signals.any_change);
12033 signals.mesh_movement.connect(signals.any_change);
12038template <
int dim,
int spacedim>
12041 Triangulation<dim, spacedim> &&tria) noexcept
12042 : Subscriptor(std::move(tria))
12043 , smooth_grid(tria.smooth_grid)
12044 , reference_cells(std::move(tria.reference_cells))
12045 , periodic_face_pairs_level_0(std::move(tria.periodic_face_pairs_level_0))
12046 , periodic_face_map(std::move(tria.periodic_face_map))
12047 , levels(std::move(tria.levels))
12048 , faces(std::move(tria.faces))
12049 , vertices(std::move(tria.vertices))
12050 , vertices_used(std::move(tria.vertices_used))
12051 , manifolds(std::move(tria.manifolds))
12052 , anisotropic_refinement(tria.anisotropic_refinement)
12053 , check_for_distorted_cells(tria.check_for_distorted_cells)
12054 , number_cache(std::move(tria.number_cache))
12055 , vertex_to_boundary_id_map_1d(std::move(tria.vertex_to_boundary_id_map_1d))
12056 , vertex_to_manifold_id_map_1d(std::move(tria.vertex_to_manifold_id_map_1d))
12058 tria.number_cache = internal::TriangulationImplementation::NumberCache<dim>();
12061 this->policy = tria.policy->clone();
12065template <
int dim,
int spacedim>
12068 Triangulation<dim, spacedim> &&tria)
noexcept
12072 smooth_grid = tria.smooth_grid;
12073 reference_cells = std::move(tria.reference_cells);
12074 periodic_face_pairs_level_0 = std::move(tria.periodic_face_pairs_level_0);
12075 periodic_face_map = std::move(tria.periodic_face_map);
12076 levels = std::move(tria.levels);
12077 faces = std::move(tria.faces);
12078 vertices = std::move(tria.vertices);
12079 vertices_used = std::move(tria.vertices_used);
12080 manifolds = std::move(tria.manifolds);
12081 anisotropic_refinement = tria.anisotropic_refinement;
12082 number_cache = tria.number_cache;
12083 vertex_to_boundary_id_map_1d = std::move(tria.vertex_to_boundary_id_map_1d);
12084 vertex_to_manifold_id_map_1d = std::move(tria.vertex_to_manifold_id_map_1d);
12086 tria.number_cache = internal::TriangulationImplementation::NumberCache<dim>();
12089 this->policy = tria.policy->clone();
12096template <
int dim,
int spacedim>
12114 AssertNothrow((dim == 1) || (vertex_to_boundary_id_map_1d ==
nullptr),
12119 AssertNothrow((dim == 1) || (vertex_to_manifold_id_map_1d ==
nullptr),
12125template <
int dim,
int spacedim>
12133 clear_despite_subscriptions();
12134 periodic_face_pairs_level_0.clear();
12135 periodic_face_map.clear();
12136 reference_cells.clear();
12138 cell_attached_data = {0, 0, {}, {}};
12139 data_serializer.clear();
12142template <
int dim,
int spacedim>
12146 return MPI_COMM_SELF;
12151template <
int dim,
int spacedim>
12156 return number_cache.active_cell_index_partitioner;
12161template <
int dim,
int spacedim>
12168 return number_cache.level_cell_index_partitioners[level];
12173template <
int dim,
int spacedim>
12176 const MeshSmoothing mesh_smoothing)
12178 smooth_grid = mesh_smoothing;
12183template <
int dim,
int spacedim>
12188 return smooth_grid;
12193template <
int dim,
int spacedim>
12197 const Manifold<dim, spacedim> &manifold_object)
12201 manifolds[m_number] = manifold_object.
clone();
12206template <
int dim,
int spacedim>
12214 manifolds[m_number] =
12221template <
int dim,
int spacedim>
12225 for (
auto &m : manifolds)
12232template <
int dim,
int spacedim>
12240 "Error: set_all_manifold_ids() can not be called on an empty Triangulation."));
12242 for (
const auto &cell : this->active_cell_iterators())
12247template <
int dim,
int spacedim>
12255 "Error: set_all_manifold_ids_on_boundary() can not be called on an empty Triangulation."));
12257 for (
const auto &cell : this->active_cell_iterators())
12259 if (cell->
face(f)->at_boundary())
12260 cell->
face(f)->set_all_manifold_ids(m_number);
12264template <
int dim,
int spacedim>
12273 "Error: set_all_manifold_ids_on_boundary() can not be called on an empty Triangulation."));
12275 bool boundary_found =
false;
12277 for (
const auto &cell : this->active_cell_iterators())
12281 if (cell->
face(f)->at_boundary() &&
12282 cell->
face(f)->boundary_id() == b_id)
12284 boundary_found =
true;
12285 cell->
face(f)->set_manifold_id(m_number);
12290 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++
e)
12291 if (cell->
line(e)->at_boundary() &&
12292 cell->
line(e)->boundary_id() == b_id)
12294 boundary_found =
true;
12295 cell->
line(e)->set_manifold_id(m_number);
12299 (void)boundary_found;
12300 Assert(boundary_found, ExcBoundaryIdNotFound(b_id));
12305template <
int dim,
int spacedim>
12317 const auto it = manifolds.find(m_number);
12319 if (it != manifolds.end())
12322 return *(it->second);
12328 "No manifold of the manifold id " + std::to_string(m_number) +
12329 " has been attached to the triangulation. "
12330 "Please attach the right manifold with Triangulation::set_manifold()."));
12338template <
int dim,
int spacedim>
12347 std::vector<types::boundary_id> boundary_ids;
12348 for (std::map<unsigned int, types::boundary_id>::const_iterator p =
12349 vertex_to_boundary_id_map_1d->begin();
12350 p != vertex_to_boundary_id_map_1d->end();
12352 boundary_ids.push_back(p->second);
12354 return boundary_ids;
12358 std::set<types::boundary_id> b_ids;
12359 for (
auto cell : active_cell_iterators())
12363 b_ids.insert(cell->
face(face)->boundary_id());
12364 std::vector<types::boundary_id> boundary_ids(b_ids.begin(), b_ids.end());
12365 return boundary_ids;
12371template <
int dim,
int spacedim>
12376 std::set<types::manifold_id> m_ids;
12377 for (
const auto &cell : active_cell_iterators())
12382 m_ids.insert(face->manifold_id());
12385 const auto line_indices = internal::TriaAccessorImplementation::
12386 Implementation::get_line_indices_of_cell(*cell);
12387 for (
unsigned int l = 0;
l < cell->
n_lines(); ++
l)
12389 raw_line_iterator line(
this, 0, line_indices[l]);
12390 m_ids.insert(line->manifold_id());
12394 return {m_ids.begin(), m_ids.end()};
12402template <
int dim,
int spacedim>
12405 const Triangulation<dim, spacedim> &other_tria)
12407 Assert((vertices.empty()) && (levels.empty()) && (faces ==
nullptr),
12408 ExcTriangulationNotEmpty(vertices.size(), levels.size()));
12410 (dim == 1 || other_tria.
faces !=
nullptr),
12412 "When calling Triangulation::copy_triangulation(), "
12413 "the target triangulation must be empty but the source "
12414 "triangulation (the argument to this function) must contain "
12415 "something. Here, it seems like the source does not "
12416 "contain anything at all."));
12427 faces = std::make_unique<internal::TriangulationImplementation::TriaFaces>(
12428 *other_tria.
faces);
12430 for (
const auto &p : other_tria.
manifolds)
12431 set_manifold(p.first, *p.second);
12434 levels.reserve(other_tria.
levels.size());
12435 for (
unsigned int level = 0; level < other_tria.
levels.size(); ++level)
12437 std::make_unique<internal::TriangulationImplementation::TriaLevel>(
12438 *other_tria.
levels[level]));
12444 vertex_to_boundary_id_map_1d =
12445 std::make_unique<std::map<unsigned int, types::boundary_id>>(
12448 vertex_to_manifold_id_map_1d =
12449 std::make_unique<std::map<unsigned int, types::manifold_id>>(
12454 this->policy = other_tria.
policy->clone();
12457 this->periodic_face_pairs_level_0.reserve(
12462 auto entry = other_entry;
12464 cell_iterator(
this, entry.cell[0]->level(), entry.cell[0]->index());
12466 cell_iterator(
this, entry.cell[1]->level(), entry.cell[1]->index());
12467 periodic_face_pairs_level_0.emplace_back(entry);
12470 for (
auto [first_cell_, second_cell_and_orientation] :
12473 auto first_cell = first_cell_;
12474 first_cell.first = cell_iterator(
this,
12475 first_cell.first->level(),
12476 first_cell.first->index());
12477 second_cell_and_orientation.first.first =
12478 cell_iterator(
this,
12479 second_cell_and_orientation.first.first->level(),
12480 second_cell_and_orientation.first.first->index());
12482 this->periodic_face_map[first_cell] = second_cell_and_orientation;
12497template <
int dim,
int spacedim>
12501 this->update_reference_cells();
12503 if (this->all_reference_cells_are_hyper_cube())
12506 std::make_unique<internal::TriangulationImplementation::PolicyWrapper<
12509 internal::TriangulationImplementation::Implementation>>();
12514 std::make_unique<internal::TriangulationImplementation::PolicyWrapper<
12517 internal::TriangulationImplementation::ImplementationMixedMesh>>();
12523template <
int dim,
int spacedim>
12526 const std::vector<Point<spacedim>> &v,
12527 const std::vector<CellData<dim>> &cells,
12528 const SubCellData &subcelldata)
12530 Assert((vertices.empty()) && (levels.empty()) && (faces ==
nullptr),
12531 ExcTriangulationNotEmpty(vertices.size(), levels.size()));
12546 clear_despite_subscriptions();
12554 reset_cell_vertex_indices_cache();
12556 *
this, levels.size(), number_cache);
12557 reset_active_cell_indices();
12558 reset_global_cell_indices();
12563 if (check_for_distorted_cells)
12565 DistortedCellList distorted_cells = collect_distorted_coarse_cells(*
this);
12569 AssertThrow(distorted_cells.distorted_cells.empty(), distorted_cells);
12602 if ((dim == spacedim - 1) && all_reference_cells_are_hyper_cube())
12610 const bool values[][2] = {{
false,
true}, {
true,
false}};
12613 correct(i, j) =
values[i][j];
12618 const bool values[][4] = {{
false,
true,
true,
false},
12619 {
true,
false,
false,
true},
12620 {
true,
false,
false,
true},
12621 {
false,
true,
true,
false}};
12624 correct(i, j) = (
values[i][j]);
12632 std::list<active_cell_iterator> this_round, next_round;
12633 active_cell_iterator neighbor;
12637 this_round.push_back(begin_active());
12638 begin_active()->set_direction_flag(
true);
12639 begin_active()->set_user_flag();
12641 while (this_round.size() > 0)
12643 for (
const auto &cell : this_round)
12647 if (cell->
face(i)->at_boundary() ==
false)
12653 const unsigned int nb_of_nb =
12658 if (neighbor->user_flag_set())
12661 !(correct(i, nb_of_nb) ^
12662 (neighbor->direction_flag() ==
12665 "The triangulation you are trying to create is not orientable."));
12674 if (correct(i, nb_of_nb) ^
12675 (neighbor->direction_flag() ==
12677 neighbor->set_direction_flag(
12678 !neighbor->direction_flag());
12679 neighbor->set_user_flag();
12680 next_round.push_back(neighbor);
12691 if (next_round.empty())
12692 for (
const auto &cell : this->active_cell_iterators())
12695 next_round.push_back(cell);
12702 next_round.swap(this_round);
12703 next_round.clear();
12705 clear_user_flags();
12708 this->update_cell_relations();
12716template <
int dim,
int spacedim>
12719 const TriangulationDescription::Description<dim, spacedim> &construction_data)
12727 auto cell_infos = construction_data.
cell_infos;
12730 for (
auto &cell_info : cell_infos)
12734 [&](TriangulationDescription::CellData<dim> a,
12735 TriangulationDescription::CellData<dim> b) {
12736 const CellId a_id(a.id);
12737 const CellId b_id(b.id);
12739 const auto a_coarse_cell_index =
12740 this->coarse_cell_id_to_coarse_cell_index(a_id.get_coarse_cell_id());
12741 const auto b_coarse_cell_index =
12742 this->coarse_cell_id_to_coarse_cell_index(b_id.get_coarse_cell_id());
12749 if (a_coarse_cell_index != b_coarse_cell_index)
12750 return a_coarse_cell_index < b_coarse_cell_index;
12752 return a_id < b_id;
12758 for (
unsigned int level = 0;
12759 level < cell_infos.size() && !cell_infos[level].empty();
12765 auto cell = this->
begin(level);
12766 auto cell_info = cell_infos[level].begin();
12767 for (; cell_info != cell_infos[level].end(); ++cell_info)
12769 while (cell_info->id != cell->
id().template to_binary<dim>())
12773 cell->
face(face)->set_manifold_id(
12774 cell_info->manifold_line_ids[face]);
12778 cell->
face(face)->set_manifold_id(
12779 cell_info->manifold_quad_ids[face]);
12781 const auto line_indices = internal::TriaAccessorImplementation::
12782 Implementation::get_line_indices_of_cell(*cell);
12783 for (
unsigned int l = 0;
l < cell->
n_lines(); ++
l)
12785 raw_line_iterator line(
this, 0, line_indices[l]);
12786 line->set_manifold_id(cell_info->manifold_line_ids[l]);
12795 if (level + 1 != cell_infos.size())
12799 auto coarse_cell = this->
begin(level);
12800 auto fine_cell_info = cell_infos[level + 1].begin();
12803 for (; fine_cell_info != cell_infos[level + 1].end();
12808 !coarse_cell->id().is_parent_of(CellId(fine_cell_info->id)))
12812 coarse_cell->set_refine_flag();
12816 ::Triangulation<dim,
12817 spacedim>::execute_coarsening_and_refinement();
12822 for (
unsigned int level = 0;
12823 level < cell_infos.size() && !cell_infos[level].empty();
12826 auto cell = this->
begin(level);
12827 auto cell_info = cell_infos[level].begin();
12828 for (; cell_info != cell_infos[level].end(); ++cell_info)
12831 while (cell_info->id != cell->
id().template to_binary<dim>())
12835 for (
auto pair : cell_info->boundary_ids)
12836 if (cell->
face(pair.first)->at_boundary())
12837 cell->
face(pair.first)->set_boundary_id(pair.second);
12846template <
int dim,
int spacedim>
12852 "This function can only be called if dim == spacedim-1."));
12853 for (
const auto &cell : this->active_cell_iterators())
12859template <
int dim,
int spacedim>
12864 ExcMessage(
"Error: An empty Triangulation can not be refined."));
12866 for (
const auto &cell : this->active_cell_iterators())
12876template <
int dim,
int spacedim>
12880 for (
unsigned int i = 0; i < times; ++i)
12882 set_all_refine_flags();
12883 execute_coarsening_and_refinement();
12889template <
int dim,
int spacedim>
12893 for (
unsigned int i = 0; i < times; ++i)
12895 for (
const auto &cell : this->active_cell_iterators())
12900 execute_coarsening_and_refinement();
12910template <
int dim,
int spacedim>
12915 std::vector<bool>::iterator i = v.begin();
12917 for (
const auto &cell : this->active_cell_iterators())
12918 for (
unsigned int j = 0; j < dim; ++j, ++i)
12927template <
int dim,
int spacedim>
12931 std::vector<bool> v;
12932 save_refine_flags(v);
12941template <
int dim,
int spacedim>
12945 std::vector<bool> v;
12947 load_refine_flags(v);
12952template <
int dim,
int spacedim>
12958 std::vector<bool>::const_iterator i = v.begin();
12959 for (
const auto &cell : this->active_cell_iterators())
12961 unsigned int ref_case = 0;
12963 for (
unsigned int j = 0; j < dim; ++j, ++i)
12965 ref_case += 1 << j;
12967 ExcGridReadError());
12979template <
int dim,
int spacedim>
12982 std::vector<bool> &v)
const
12985 std::vector<bool>::iterator i = v.begin();
12986 for (
const auto &cell : this->active_cell_iterators())
12997template <
int dim,
int spacedim>
13001 std::vector<bool> v;
13002 save_coarsen_flags(v);
13011template <
int dim,
int spacedim>
13015 std::vector<bool> v;
13020 load_coarsen_flags(v);
13025template <
int dim,
int spacedim>
13028 const std::vector<bool> &v)
13032 std::vector<bool>::const_iterator i = v.begin();
13033 for (
const auto &cell : this->active_cell_iterators())
13046template <
int dim,
int spacedim>
13050 return anisotropic_refinement;
13060 std::vector<std::vector<bool>>
13061 extract_raw_coarsen_flags(
13062 const std::vector<std::unique_ptr<
13063 ::internal::TriangulationImplementation::TriaLevel>> &levels)
13065 std::vector<std::vector<bool>> coarsen_flags(levels.size());
13066 for (
unsigned int level = 0; level < levels.size(); ++level)
13067 coarsen_flags[level] = levels[level]->coarsen_flags;
13068 return coarsen_flags;
13071 std::vector<std::vector<std::uint8_t>>
13072 extract_raw_refine_flags(
13073 const std::vector<std::unique_ptr<
13074 ::internal::TriangulationImplementation::TriaLevel>> &levels)
13076 std::vector<std::vector<std::uint8_t>> refine_flags(levels.size());
13077 for (
unsigned int level = 0; level < levels.size(); ++level)
13078 refine_flags[level] = levels[level]->refine_flags;
13079 return refine_flags;
13092 clear_user_data(std::vector<std::unique_ptr<
13093 internal::TriangulationImplementation::TriaLevel>> &levels)
13095 for (
auto &level : levels)
13096 level->cells.clear_user_data();
13102 clear_user_data(internal::TriangulationImplementation::TriaFaces *faces)
13104 if (faces->
dim == 2)
13110 if (faces->
dim == 3)
13120template <
int dim,
int spacedim>
13125 ::clear_user_data(levels);
13127 ::clear_user_data(faces.get());
13135 clear_user_flags_line(
13138 std::unique_ptr<internal::TriangulationImplementation::TriaLevel>>
13140 internal::TriangulationImplementation::TriaFaces *faces)
13144 for (
const auto &level : levels)
13145 level->cells.clear_user_flags();
13147 else if (dim == 2 || dim == 3)
13159template <
int dim,
int spacedim>
13163 ::clear_user_flags_line(dim, levels, faces.get());
13171 clear_user_flags_quad(
13174 std::unique_ptr<internal::TriangulationImplementation::TriaLevel>>
13176 internal::TriangulationImplementation::TriaFaces *faces)
13184 for (
const auto &level : levels)
13185 level->cells.clear_user_flags();
13199template <
int dim,
int spacedim>
13203 ::clear_user_flags_quad(dim, levels, faces.get());
13211 clear_user_flags_hex(
13214 std::unique_ptr<internal::TriangulationImplementation::TriaLevel>>
13216 internal::TriangulationImplementation::TriaFaces *)
13228 for (
const auto &level : levels)
13229 level->cells.clear_user_flags();
13239template <
int dim,
int spacedim>
13243 ::clear_user_flags_hex(dim, levels, faces.get());
13248template <
int dim,
int spacedim>
13252 clear_user_flags_line();
13253 clear_user_flags_quad();
13254 clear_user_flags_hex();
13259template <
int dim,
int spacedim>
13263 save_user_flags_line(out);
13266 save_user_flags_quad(out);
13269 save_user_flags_hex(out);
13277template <
int dim,
int spacedim>
13285 std::vector<bool> tmp;
13287 save_user_flags_line(tmp);
13288 v.insert(v.end(), tmp.begin(), tmp.end());
13292 save_user_flags_quad(tmp);
13293 v.insert(v.end(), tmp.begin(), tmp.end());
13298 save_user_flags_hex(tmp);
13299 v.insert(v.end(), tmp.begin(), tmp.end());
13308template <
int dim,
int spacedim>
13312 load_user_flags_line(in);
13315 load_user_flags_quad(in);
13318 load_user_flags_hex(in);
13326template <
int dim,
int spacedim>
13331 std::vector<bool> tmp;
13335 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
13337 load_user_flags_line(tmp);
13342 tmp.insert(tmp.end(),
13343 v.begin() + n_lines(),
13344 v.begin() + n_lines() + n_quads());
13345 load_user_flags_quad(tmp);
13351 tmp.insert(tmp.end(),
13352 v.begin() + n_lines() + n_quads(),
13353 v.begin() + n_lines() + n_quads() + n_hexs());
13354 load_user_flags_hex(tmp);
13363template <
int dim,
int spacedim>
13366 std::vector<bool> &v)
const
13368 v.resize(n_lines(),
false);
13369 std::vector<bool>::iterator i = v.begin();
13370 line_iterator line = begin_line(), endl = end_line();
13371 for (; line != endl; ++line, ++i)
13372 *i = line->user_flag_set();
13379template <
int dim,
int spacedim>
13383 std::vector<bool> v;
13384 save_user_flags_line(v);
13393template <
int dim,
int spacedim>
13397 std::vector<bool> v;
13402 load_user_flags_line(v);
13407template <
int dim,
int spacedim>
13410 const std::vector<bool> &v)
13412 Assert(v.size() == n_lines(), ExcGridReadError());
13414 line_iterator line = begin_line(), endl = end_line();
13415 std::vector<bool>::const_iterator i = v.begin();
13416 for (; line != endl; ++line, ++i)
13418 line->set_user_flag();
13420 line->clear_user_flag();
13429 template <
typename Iterator>
13431 get_user_flag(
const Iterator &i)
13433 return i->user_flag_set();
13438 template <
int structdim,
int dim,
int spacedim>
13440 get_user_flag(
const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
13448 template <
typename Iterator>
13450 set_user_flag(
const Iterator &i)
13452 i->set_user_flag();
13457 template <
int structdim,
int dim,
int spacedim>
13459 set_user_flag(
const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
13466 template <
typename Iterator>
13468 clear_user_flag(
const Iterator &i)
13470 i->clear_user_flag();
13475 template <
int structdim,
int dim,
int spacedim>
13478 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
13486template <
int dim,
int spacedim>
13489 std::vector<bool> &v)
const
13491 v.resize(n_quads(),
false);
13495 std::vector<bool>::iterator i = v.begin();
13496 quad_iterator quad = begin_quad(), endq = end_quad();
13497 for (; quad != endq; ++quad, ++i)
13498 *i = get_user_flag(quad);
13506template <
int dim,
int spacedim>
13510 std::vector<bool> v;
13511 save_user_flags_quad(v);
13520template <
int dim,
int spacedim>
13524 std::vector<bool> v;
13529 load_user_flags_quad(v);
13534template <
int dim,
int spacedim>
13537 const std::vector<bool> &v)
13539 Assert(v.size() == n_quads(), ExcGridReadError());
13543 quad_iterator quad = begin_quad(), endq = end_quad();
13544 std::vector<bool>::const_iterator i = v.begin();
13545 for (; quad != endq; ++quad, ++i)
13547 set_user_flag(quad);
13549 clear_user_flag(quad);
13557template <
int dim,
int spacedim>
13560 std::vector<bool> &v)
const
13562 v.resize(n_hexs(),
false);
13566 std::vector<bool>::iterator i = v.begin();
13567 hex_iterator hex = begin_hex(), endh = end_hex();
13568 for (; hex != endh; ++hex, ++i)
13569 *i = get_user_flag(hex);
13577template <
int dim,
int spacedim>
13581 std::vector<bool> v;
13582 save_user_flags_hex(v);
13591template <
int dim,
int spacedim>
13595 std::vector<bool> v;
13600 load_user_flags_hex(v);
13605template <
int dim,
int spacedim>
13608 const std::vector<bool> &v)
13610 Assert(v.size() == n_hexs(), ExcGridReadError());
13614 hex_iterator hex = begin_hex(), endh = end_hex();
13615 std::vector<bool>::const_iterator i = v.begin();
13616 for (; hex != endh; ++hex, ++i)
13618 set_user_flag(hex);
13620 clear_user_flag(hex);
13628template <
int dim,
int spacedim>
13631 std::vector<unsigned int> &v)
const
13637 std::vector<unsigned int> tmp;
13639 save_user_indices_line(tmp);
13640 v.insert(v.end(), tmp.begin(), tmp.end());
13644 save_user_indices_quad(tmp);
13645 v.insert(v.end(), tmp.begin(), tmp.end());
13650 save_user_indices_hex(tmp);
13651 v.insert(v.end(), tmp.begin(), tmp.end());
13660template <
int dim,
int spacedim>
13663 const std::vector<unsigned int> &v)
13666 std::vector<unsigned int> tmp;
13670 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
13672 load_user_indices_line(tmp);
13677 tmp.insert(tmp.end(),
13678 v.begin() + n_lines(),
13679 v.begin() + n_lines() + n_quads());
13680 load_user_indices_quad(tmp);
13686 tmp.insert(tmp.end(),
13687 v.begin() + n_lines() + n_quads(),
13688 v.begin() + n_lines() + n_quads() + n_hexs());
13689 load_user_indices_hex(tmp);
13698template <
int dim,
int spacedim>
13703 std::ofstream ofs(file_basename +
"_triangulation.data");
13704 boost::archive::text_oarchive oa(ofs, boost::archive::no_header);
13709 std::ofstream ifs(file_basename +
".info");
13711 <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_active_cells"
13713 << internal::CellAttachedDataSerializer<dim, spacedim>::version_number
13714 <<
" " << 1 <<
" " << this->cell_attached_data.pack_callbacks_fixed.size()
13715 <<
" " << this->cell_attached_data.pack_callbacks_variable.size() <<
" "
13716 << this->n_global_active_cells() << std::endl;
13719 this->save_attached_data(0, this->n_global_active_cells(), file_basename);
13724template <
int dim,
int spacedim>
13729 std::ifstream ifs(file_basename +
"_triangulation.data");
13732 boost::archive::text_iarchive ia(ifs, boost::archive::no_header);
13736 unsigned int version, numcpus, attached_count_fixed, attached_count_variable,
13737 n_global_active_cells;
13739 std::ifstream ifs(std::string(file_basename) +
".info");
13741 std::string firstline;
13742 getline(ifs, firstline);
13743 ifs >> version >> numcpus >> attached_count_fixed >>
13744 attached_count_variable >> n_global_active_cells;
13748 ExcMessage(
"Incompatible number of CPUs found in .info file."));
13750 const auto expected_version =
13751 ::internal::CellAttachedDataSerializer<dim,
13752 spacedim>::version_number;
13755 "The information saved in the file you are trying "
13756 "to read the triangulation from was written with an "
13757 "incompatible file format version and cannot be read."));
13758 Assert(this->n_global_active_cells() == n_global_active_cells,
13759 ExcMessage(
"The number of cells of the triangulation differs "
13760 "from the number of cells written into the .info file."));
13764 this->cell_attached_data.n_attached_data_sets = 0;
13765 this->cell_attached_data.n_attached_deserialize =
13766 attached_count_fixed + attached_count_variable;
13768 this->load_attached_data(0,
13769 this->n_global_active_cells(),
13772 attached_count_fixed,
13773 attached_count_variable);
13775 this->update_cell_relations();
13781 template <
typename Iterator>
13783 get_user_index(
const Iterator &i)
13785 return i->user_index();
13790 template <
int structdim,
int dim,
int spacedim>
13793 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
13801 template <
typename Iterator>
13803 set_user_index(
const Iterator &i,
const unsigned int x)
13805 i->set_user_index(x);
13810 template <
int structdim,
int dim,
int spacedim>
13813 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &,
13814 const unsigned int)
13822template <
int dim,
int spacedim>
13825 std::vector<unsigned int> &v)
const
13827 v.resize(n_lines(), 0);
13828 std::vector<unsigned int>::iterator i = v.begin();
13829 line_iterator line = begin_line(), endl = end_line();
13830 for (; line != endl; ++line, ++i)
13831 *i = line->user_index();
13836template <
int dim,
int spacedim>
13839 const std::vector<unsigned int> &v)
13841 Assert(v.size() == n_lines(), ExcGridReadError());
13843 line_iterator line = begin_line(), endl = end_line();
13844 std::vector<unsigned int>::const_iterator i = v.begin();
13845 for (; line != endl; ++line, ++i)
13846 line->set_user_index(*i);
13850template <
int dim,
int spacedim>
13853 std::vector<unsigned int> &v)
const
13855 v.resize(n_quads(), 0);
13859 std::vector<unsigned int>::iterator i = v.begin();
13860 quad_iterator quad = begin_quad(), endq = end_quad();
13861 for (; quad != endq; ++quad, ++i)
13862 *i = get_user_index(quad);
13868template <
int dim,
int spacedim>
13871 const std::vector<unsigned int> &v)
13873 Assert(v.size() == n_quads(), ExcGridReadError());
13877 quad_iterator quad = begin_quad(), endq = end_quad();
13878 std::vector<unsigned int>::const_iterator i = v.begin();
13879 for (; quad != endq; ++quad, ++i)
13880 set_user_index(quad, *i);
13885template <
int dim,
int spacedim>
13888 std::vector<unsigned int> &v)
const
13890 v.resize(n_hexs(), 0);
13894 std::vector<unsigned int>::iterator i = v.begin();
13895 hex_iterator hex = begin_hex(), endh = end_hex();
13896 for (; hex != endh; ++hex, ++i)
13897 *i = get_user_index(hex);
13903template <
int dim,
int spacedim>
13906 const std::vector<unsigned int> &v)
13908 Assert(v.size() == n_hexs(), ExcGridReadError());
13912 hex_iterator hex = begin_hex(), endh = end_hex();
13913 std::vector<unsigned int>::const_iterator i = v.begin();
13914 for (; hex != endh; ++hex, ++i)
13915 set_user_index(hex, *i);
13927 template <
typename Iterator>
13929 get_user_pointer(
const Iterator &i)
13931 return i->user_pointer();
13936 template <
int structdim,
int dim,
int spacedim>
13939 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
13947 template <
typename Iterator>
13949 set_user_pointer(
const Iterator &i,
void *x)
13951 i->set_user_pointer(x);
13956 template <
int structdim,
int dim,
int spacedim>
13959 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &,
13968template <
int dim,
int spacedim>
13971 std::vector<void *> &v)
const
13977 std::vector<void *> tmp;
13979 save_user_pointers_line(tmp);
13980 v.insert(v.end(), tmp.begin(), tmp.end());
13984 save_user_pointers_quad(tmp);
13985 v.insert(v.end(), tmp.begin(), tmp.end());
13990 save_user_pointers_hex(tmp);
13991 v.insert(v.end(), tmp.begin(), tmp.end());
14000template <
int dim,
int spacedim>
14003 const std::vector<void *> &v)
14006 std::vector<void *> tmp;
14010 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
14012 load_user_pointers_line(tmp);
14017 tmp.insert(tmp.end(),
14018 v.begin() + n_lines(),
14019 v.begin() + n_lines() + n_quads());
14020 load_user_pointers_quad(tmp);
14026 tmp.insert(tmp.end(),
14027 v.begin() + n_lines() + n_quads(),
14028 v.begin() + n_lines() + n_quads() + n_hexs());
14029 load_user_pointers_hex(tmp);
14038template <
int dim,
int spacedim>
14041 std::vector<void *> &v)
const
14043 v.resize(n_lines(),
nullptr);
14044 std::vector<void *>::iterator i = v.begin();
14045 line_iterator line = begin_line(), endl = end_line();
14046 for (; line != endl; ++line, ++i)
14047 *i = line->user_pointer();
14052template <
int dim,
int spacedim>
14055 const std::vector<void *> &v)
14057 Assert(v.size() == n_lines(), ExcGridReadError());
14059 line_iterator line = begin_line(), endl = end_line();
14060 std::vector<void *>::const_iterator i = v.begin();
14061 for (; line != endl; ++line, ++i)
14062 line->set_user_pointer(*i);
14067template <
int dim,
int spacedim>
14070 std::vector<void *> &v)
const
14072 v.resize(n_quads(),
nullptr);
14076 std::vector<void *>::iterator i = v.begin();
14077 quad_iterator quad = begin_quad(), endq = end_quad();
14078 for (; quad != endq; ++quad, ++i)
14079 *i = get_user_pointer(quad);
14085template <
int dim,
int spacedim>
14088 const std::vector<void *> &v)
14090 Assert(v.size() == n_quads(), ExcGridReadError());
14094 quad_iterator quad = begin_quad(), endq = end_quad();
14095 std::vector<void *>::const_iterator i = v.begin();
14096 for (; quad != endq; ++quad, ++i)
14097 set_user_pointer(quad, *i);
14102template <
int dim,
int spacedim>
14105 std::vector<void *> &v)
const
14107 v.resize(n_hexs(),
nullptr);
14111 std::vector<void *>::iterator i = v.begin();
14112 hex_iterator hex = begin_hex(), endh = end_hex();
14113 for (; hex != endh; ++hex, ++i)
14114 *i = get_user_pointer(hex);
14120template <
int dim,
int spacedim>
14123 const std::vector<void *> &v)
14125 Assert(v.size() == n_hexs(), ExcGridReadError());
14129 hex_iterator hex = begin_hex(), endh = end_hex();
14130 std::vector<void *>::const_iterator i = v.begin();
14131 for (; hex != endh; ++hex, ++i)
14132 set_user_pointer(hex, *i);
14142template <
int dim,
int spacedim>
14150 return begin_raw_line(level);
14152 return begin_raw_quad(level);
14154 return begin_raw_hex(level);
14157 return raw_cell_iterator();
14163template <
int dim,
int spacedim>
14171 return begin_line(level);
14173 return begin_quad(level);
14175 return begin_hex(level);
14178 return cell_iterator();
14184template <
int dim,
int spacedim>
14192 return begin_active_line(level);
14194 return begin_active_quad(level);
14196 return begin_active_hex(level);
14199 return active_cell_iterator();
14205template <
int dim,
int spacedim>
14210 const unsigned int level = levels.size() - 1;
14211 if (levels[level]->cells.n_objects() == 0)
14216 raw_cell_iterator ri(
const_cast<Triangulation<dim, spacedim> *
>(
this),
14218 levels[level]->cells.n_objects() - 1);
14221 if (ri->used() ==
true)
14224 if (ri->used() ==
true)
14231template <
int dim,
int spacedim>
14237 cell_iterator cell = last();
14242 if (cell->is_active() ==
true)
14245 if (cell->is_active() ==
true)
14253template <
int dim,
int spacedim>
14257 const CellId &cell_id)
const
14260 this->contains_cell(cell_id),
14262 "CellId is invalid for this triangulation.\n"
14263 "Either the provided CellId does not correspond to a cell in this "
14264 "triangulation object, or, in case you are using a parallel "
14265 "triangulation, may correspond to an artificial cell that is less "
14266 "refined on this processor. In the case of "
14267 "parallel::fullydistributed::Triangulation, the corresponding coarse "
14268 "cell might not be accessible by the current process."));
14270 cell_iterator cell(
14274 cell = cell->child(
static_cast<unsigned int>(child_index));
14281template <
int dim,
int spacedim>
14285 const auto coarse_cell_index =
14291 cell_iterator cell(
this, 0, coarse_cell_index);
14295 if (cell->has_children() ==
false)
14297 cell = cell->child(
static_cast<unsigned int>(child_index));
14305template <
int dim,
int spacedim>
14310 return cell_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
14317template <
int dim,
int spacedim>
14330 if (level >= levels.size())
14332 Assert(level < n_global_levels(),
14333 ExcInvalidLevel(level, n_global_levels()));
14339 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
14340 if (level < levels.size() - 1)
14341 return begin_raw(level + 1);
14347template <
int dim,
int spacedim>
14360 if (level >= levels.size())
14362 Assert(level < n_global_levels(),
14363 ExcInvalidLevel(level, n_global_levels()));
14369 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
14370 if (level < levels.size() - 1)
14371 return begin(level + 1);
14377template <
int dim,
int spacedim>
14390 if (level >= levels.size())
14392 Assert(level < n_global_levels(),
14393 ExcInvalidLevel(level, n_global_levels()));
14399 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
14400 return (level >= levels.size() - 1 ? active_cell_iterator(
end()) :
14401 begin_active(level + 1));
14406template <
int dim,
int spacedim>
14408IteratorRange<typename Triangulation<dim, spacedim>::
14412 return IteratorRange<typename Triangulation<dim, spacedim>::cell_iterator>(
14417template <
int dim,
int spacedim>
14419IteratorRange<typename Triangulation<dim, spacedim>::
14423 return IteratorRange<
14430template <
int dim,
int spacedim>
14432IteratorRange<typename Triangulation<dim, spacedim>::
14436 return IteratorRange<typename Triangulation<dim, spacedim>::cell_iterator>(
14442template <
int dim,
int spacedim>
14444IteratorRange<typename Triangulation<dim, spacedim>::
14448 return IteratorRange<
14450 begin_active(level), end_active(level));
14458template <
int dim,
int spacedim>
14467 return raw_face_iterator();
14469 return begin_line();
14471 return begin_quad();
14474 return face_iterator();
14480template <
int dim,
int spacedim>
14489 return raw_face_iterator();
14491 return begin_active_line();
14493 return begin_active_quad();
14496 return active_face_iterator();
14502template <
int dim,
int spacedim>
14511 return raw_face_iterator();
14518 return raw_face_iterator();
14524template <
int dim,
int spacedim>
14526IteratorRange<typename Triangulation<dim, spacedim>::
14530 return IteratorRange<
14532 begin_active_face(), end_face());
14538template <
int dim,
int spacedim>
14543 vertex_iterator i =
14544 raw_vertex_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this), 0, 0);
14548 while (i->used() ==
false)
14556template <
int dim,
int spacedim>
14562 return begin_vertex();
14567template <
int dim,
int spacedim>
14572 return raw_vertex_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
14584template <
int dim,
int spacedim>
14597 if (level >= levels.size())
14599 Assert(level < n_global_levels(),
14600 ExcInvalidLevel(level, n_global_levels()));
14609 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
14611 if (level >= levels.size() || levels[level]->cells.n_objects() == 0)
14614 return raw_line_iterator(
14615 const_cast<Triangulation<dim, spacedim> *
>(
this), level, 0);
14618 Assert(level == 0, ExcFacesHaveNoLevel());
14619 return raw_line_iterator(
14620 const_cast<Triangulation<dim, spacedim> *
>(
this), 0, 0);
14625template <
int dim,
int spacedim>
14631 raw_line_iterator ri = begin_raw_line(level);
14634 while (ri->used() ==
false)
14642template <
int dim,
int spacedim>
14646 const unsigned int level)
const
14649 line_iterator i = begin_line(level);
14652 while (i->has_children())
14660template <
int dim,
int spacedim>
14665 return raw_line_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
14676template <
int dim,
int spacedim>
14689 if (level >= levels.size())
14691 Assert(level < n_global_levels(),
14692 ExcInvalidLevel(level, n_global_levels()));
14700 return raw_hex_iterator();
14705 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
14707 if (level >= levels.size() || levels[level]->cells.n_objects() == 0)
14710 return raw_quad_iterator(
14711 const_cast<Triangulation<dim, spacedim> *
>(
this), level, 0);
14716 Assert(level == 0, ExcFacesHaveNoLevel());
14718 return raw_quad_iterator(
14719 const_cast<Triangulation<dim, spacedim> *
>(
this), 0, 0);
14725 return raw_hex_iterator();
14731template <
int dim,
int spacedim>
14737 raw_quad_iterator ri = begin_raw_quad(level);
14740 while (ri->used() ==
false)
14748template <
int dim,
int spacedim>
14752 const unsigned int level)
const
14755 quad_iterator i = begin_quad(level);
14758 while (i->has_children())
14766template <
int dim,
int spacedim>
14771 return raw_quad_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
14782template <
int dim,
int spacedim>
14795 if (level >= levels.size())
14797 Assert(level < n_global_levels(),
14798 ExcInvalidLevel(level, n_global_levels()));
14807 return raw_hex_iterator();
14812 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
14814 if (level >= levels.size() || levels[level]->cells.n_objects() == 0)
14817 return raw_hex_iterator(
14818 const_cast<Triangulation<dim, spacedim> *
>(
this), level, 0);
14823 return raw_hex_iterator();
14829template <
int dim,
int spacedim>
14835 raw_hex_iterator ri = begin_raw_hex(level);
14838 while (ri->used() ==
false)
14846template <
int dim,
int spacedim>
14852 hex_iterator i = begin_hex(level);
14855 while (i->has_children())
14863template <
int dim,
int spacedim>
14868 return raw_hex_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
14882 inline unsigned int
14889 inline unsigned int
14893 return c.n_active_lines;
14897 inline unsigned int
14904 inline unsigned int
14908 return c.n_active_quads;
14912 inline unsigned int
14919 inline unsigned int
14923 return c.n_active_hexes;
14930template <
int dim,
int spacedim>
14938template <
int dim,
int spacedim>
14945template <
int dim,
int spacedim>
14953template <
int dim,
int spacedim>
14961template <
int dim,
int spacedim>
14968 return n_used_vertices();
14980template <
int dim,
int spacedim>
14987 return n_vertices();
14989 return n_raw_lines();
14991 return n_raw_quads();
14999template <
int dim,
int spacedim>
15006 return n_used_vertices();
15008 return n_active_lines();
15010 return n_active_quads();
15018template <
int dim,
int spacedim>
15021 const unsigned int level)
const
15026 return n_raw_lines(level);
15028 return n_raw_quads(level);
15030 return n_raw_hexs(level);
15039template <
int dim,
int spacedim>
15042 const unsigned int level)
const
15047 return n_lines(level);
15049 return n_quads(level);
15051 return n_hexs(level);
15060template <
int dim,
int spacedim>
15063 const unsigned int level)
const
15068 return n_active_lines(level);
15070 return n_active_quads(level);
15072 return n_active_hexs(level);
15080template <
int dim,
int spacedim>
15084 if (anisotropic_refinement ==
false)
15086 for (
unsigned int lvl = 0; lvl < n_global_levels() - 1; ++lvl)
15092 for (
const auto &cell : active_cell_iterators())
15093 for (
const auto &i : cell->face_indices())
15094 if (cell->face(i)->has_children())
15101template <
int dim,
int spacedim>
15105 return number_cache.n_lines;
15110template <
int dim,
int spacedim>
15113 const unsigned int level)
const
15118 return levels[level]->cells.n_objects();
15121 Assert(
false, ExcFacesHaveNoLevel());
15126template <
int dim,
int spacedim>
15140template <
int dim,
int spacedim>
15143 const unsigned int level)
const
15146 Assert(dim == 1, ExcFacesHaveNoLevel());
15147 return number_cache.n_lines_level[level];
15151template <
int dim,
int spacedim>
15155 return number_cache.n_active_lines;
15159template <
int dim,
int spacedim>
15162 const unsigned int level)
const
15165 Assert(dim == 1, ExcFacesHaveNoLevel());
15167 return number_cache.n_active_lines_level[level];
15317template <
int dim,
int spacedim>
15321 return number_cache.n_quads;
15325template <
int dim,
int spacedim>
15328 const unsigned int level)
const
15330 Assert(dim == 2, ExcFacesHaveNoLevel());
15332 return number_cache.n_quads_level[level];
15342 return levels[level]->cells.n_objects();
15352 return levels[level]->cells.n_objects();
15366template <
int dim,
int spacedim>
15380 return faces->quads.n_objects();
15385template <
int dim,
int spacedim>
15389 return number_cache.n_active_quads;
15393template <
int dim,
int spacedim>
15396 const unsigned int level)
const
15399 Assert(dim == 2, ExcFacesHaveNoLevel());
15401 return number_cache.n_active_quads_level[level];
15405template <
int dim,
int spacedim>
15414template <
int dim,
int spacedim>
15423template <
int dim,
int spacedim>
15431template <
int dim,
int spacedim>
15440template <
int dim,
int spacedim>
15443 const unsigned int)
const
15475 return levels[level]->cells.n_objects();
15499template <
int dim,
int spacedim>
15503 return std::count(vertices_used.begin(), vertices_used.end(),
true);
15508template <
int dim,
int spacedim>
15512 return vertices_used;
15543template <
int dim,
int spacedim>
15547 cell_iterator cell = begin(0),
15548 endc = (n_levels() > 1 ? begin(1) : cell_iterator(end()));
15551 unsigned int max_vertex_index = 0;
15552 for (; cell != endc; ++cell)
15554 if (cell->vertex_index(vertex) > max_vertex_index)
15555 max_vertex_index = cell->vertex_index(vertex);
15561 std::vector<unsigned short int> usage_count(max_vertex_index + 1, 0);
15565 for (cell = begin(); cell != endc; ++cell)
15567 ++usage_count[cell->vertex_index(vertex)];
15570 static_cast<unsigned int>(
15571 *std::max_element(usage_count.begin(), usage_count.end())));
15576template <
int dim,
int spacedim>
15586template <
int dim,
int spacedim>
15595template <
int dim,
int spacedim>
15597const Triangulation<dim, spacedim>
15605template <
int dim,
int spacedim>
15608 const std::vector<GridTools::PeriodicFacePair<cell_iterator>>
15609 &periodicity_vector)
15611 periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(),
15612 periodicity_vector.begin(),
15613 periodicity_vector.end());
15616 update_periodic_face_map();
15621template <
int dim,
int spacedim>
15623const typename std::map<
15624 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
unsigned int>,
15625 std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
15630 return periodic_face_map;
15634template <
int dim,
int spacedim>
15641 if (
dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *
>(
15645 this->local_cell_relations.clear();
15648 for (
const auto &cell : this->active_cell_iterators())
15649 this->local_cell_relations.emplace_back(
15655template <
int dim,
int spacedim>
15669 if (smooth_grid & limit_level_difference_at_vertices)
15673 signals.pre_refinement();
15675 execute_coarsening();
15677 const DistortedCellList cells_with_distorted_children = execute_refinement();
15679 reset_cell_vertex_indices_cache();
15684 if (smooth_grid & limit_level_difference_at_vertices)
15689 this->policy->update_neighbors(*
this);
15690 reset_active_cell_indices();
15692 reset_global_cell_indices();
15695 signals.post_refinement();
15697 AssertThrow(cells_with_distorted_children.distorted_cells.empty(),
15698 cells_with_distorted_children);
15700 update_periodic_face_map();
15702 if (this->cell_attached_data.n_attached_data_sets == 0)
15703 this->update_cell_relations();
15729 for (
const auto &cell : cell_iterators())
15731 if (cell->has_children() && cell->reference_cell().is_hyper_cube())
15732 for (
const unsigned int f : cell->face_indices())
15733 if (cell->at_boundary(f) ==
false)
15735 for (
const auto &child : cell->child_iterators())
15738 child->at_boundary(f) ==
false,
15740 "We ended up with a triangulation whose child cells "
15741 "are not connected to their neighbors as expected. "
15742 "When you created the triangulation, did you forget "
15743 "to call GridTools::consistently_order_cells() "
15744 "before calling Triangulation::create_triangulation()?"));
15753template <
int dim,
int spacedim>
15757 unsigned int active_cell_index = 0;
15758 for (raw_cell_iterator cell = begin_raw(); cell !=
end(); ++cell)
15759 if ((cell->used() ==
false) || cell->has_children())
15763 cell->set_active_cell_index(active_cell_index);
15764 ++active_cell_index;
15772template <
int dim,
int spacedim>
15778 for (
const auto &cell : active_cell_iterators())
15779 cell->set_global_active_cell_index(cell_index++);
15782 for (
unsigned int l = 0;
l < levels.size(); ++
l)
15785 for (
const auto &cell : cell_iterators_on_level(l))
15786 cell->set_global_level_cell_index(cell_index++);
15792template <
int dim,
int spacedim>
15796 for (
unsigned int l = 0;
l < levels.size(); ++
l)
15798 constexpr unsigned int max_vertices_per_cell = 1 << dim;
15799 std::vector<unsigned int> &cache = levels[
l]->cell_vertex_indices_cache;
15801 cache.resize(levels[l]->refine_flags.size() * max_vertices_per_cell,
15803 for (
const auto &cell : cell_iterators_on_level(l))
15805 const unsigned int my_index = cell->index() * max_vertices_per_cell;
15810 const ReferenceCell ref_cell = cell->reference_cell();
15812 for (
unsigned int face = 4; face < 6; ++face)
15814 const auto face_iter = cell->face(face);
15815 const std::array<bool, 2> line_orientations{
15816 {face_iter->line_orientation(0),
15817 face_iter->line_orientation(1)}};
15818 std::array<unsigned int, 4> raw_vertex_indices{
15819 {face_iter->line(0)->vertex_index(1 - line_orientations[0]),
15820 face_iter->line(1)->vertex_index(1 - line_orientations[1]),
15821 face_iter->line(0)->vertex_index(line_orientations[0]),
15822 face_iter->line(1)->vertex_index(line_orientations[1])}};
15824 const unsigned char combined_orientation =
15825 levels[
l]->face_orientations.get_combined_orientation(
15827 std::array<unsigned int, 4> vertex_order{
15830 combined_orientation),
15833 combined_orientation),
15836 combined_orientation),
15838 3, face, combined_orientation)}};
15840 const unsigned int index = my_index + 4 * (face - 4);
15841 for (
unsigned int i = 0; i < 4; ++i)
15842 cache[
index + i] = raw_vertex_indices[vertex_order[i]];
15846 const std::array<bool, 2> line_orientations{
15847 {cell->line_orientation(0), cell->line_orientation(1)}};
15848 std::array<unsigned int, 4> raw_vertex_indices{
15849 {cell->line(0)->vertex_index(1 - line_orientations[0]),
15850 cell->line(1)->vertex_index(1 - line_orientations[1]),
15851 cell->line(0)->vertex_index(line_orientations[0]),
15852 cell->line(1)->vertex_index(line_orientations[1])}};
15853 for (
unsigned int i = 0; i < 4; ++i)
15854 cache[my_index + i] = raw_vertex_indices[i];
15857 for (
const unsigned int i : cell->vertex_indices())
15858 cache[my_index + i] = internal::TriaAccessorImplementation::
15859 Implementation::vertex_index(*cell, i);
15866template <
int dim,
int spacedim>
15871 periodic_face_map.clear();
15873 typename std::vector<
15874 GridTools::PeriodicFacePair<cell_iterator>>::const_iterator it;
15875 for (it = periodic_face_pairs_level_0.begin();
15876 it != periodic_face_pairs_level_0.end();
15879 update_periodic_face_map_recursively<dim, spacedim>(it->cell[0],
15884 periodic_face_map);
15886 const auto face_reference_cell =
15887 it->cell[0]->reference_cell().face_reference_cell(it->face_idx[0]);
15889 update_periodic_face_map_recursively<dim, spacedim>(
15894 face_reference_cell.get_inverse_combined_orientation(it->orientation),
15895 periodic_face_map);
15899 typename std::map<std::pair<cell_iterator, unsigned int>,
15900 std::pair<std::pair<cell_iterator, unsigned int>,
15901 unsigned char>>::const_iterator it_test;
15902 for (it_test = periodic_face_map.begin(); it_test != periodic_face_map.end();
15906 it_test->first.first;
15908 it_test->second.first.first;
15913 Assert(periodic_face_map[it_test->second.first].first ==
15922template <
int dim,
int spacedim>
15926 std::set<ReferenceCell> reference_cells_set;
15927 for (
auto cell : active_cell_iterators())
15928 if (cell->is_locally_owned())
15929 reference_cells_set.insert(cell->reference_cell());
15931 this->reference_cells =
15932 std::vector<ReferenceCell>(reference_cells_set.begin(),
15933 reference_cells_set.end());
15938template <
int dim,
int spacedim>
15940const std::vector<ReferenceCell>
15943 return this->reference_cells;
15948template <
int dim,
int spacedim>
15952 Assert(this->reference_cells.size() > 0,
15953 ExcMessage(
"You can't ask about the kinds of reference "
15954 "cells used by this triangulation if the "
15955 "triangulation doesn't yet have any cells in it."));
15956 return (this->reference_cells.size() == 1 &&
15957 this->reference_cells[0].is_hyper_cube());
15962template <
int dim,
int spacedim>
15966 Assert(this->reference_cells.size() > 0,
15967 ExcMessage(
"You can't ask about the kinds of reference "
15968 "cells used by this triangulation if the "
15969 "triangulation doesn't yet have any cells in it."));
15970 return (this->reference_cells.size() == 1 &&
15971 this->reference_cells[0].is_simplex());
15976template <
int dim,
int spacedim>
15980 Assert(this->reference_cells.size() > 0,
15981 ExcMessage(
"You can't ask about the kinds of reference "
15982 "cells used by this triangulation if the "
15983 "triangulation doesn't yet have any cells in it."));
15984 return reference_cells.size() > 1 ||
15985 ((reference_cells[0].is_hyper_cube() ==
false) &&
15986 (reference_cells[0].is_simplex() ==
false));
15991template <
int dim,
int spacedim>
15994 const std::function<std::vector<char>(
const cell_iterator &,
15995 const ::CellStatus)>
15997 const bool returns_variable_size_data)
16003 if (returns_variable_size_data)
16005 handle = 2 * this->cell_attached_data.pack_callbacks_variable.size();
16006 this->cell_attached_data.pack_callbacks_variable.push_back(pack_callback);
16010 handle = 2 * this->cell_attached_data.pack_callbacks_fixed.size() + 1;
16011 this->cell_attached_data.pack_callbacks_fixed.push_back(pack_callback);
16015 ++this->cell_attached_data.n_attached_data_sets;
16022template <
int dim,
int spacedim>
16025 const unsigned int handle,
16026 const std::function<
16027 void(
const cell_iterator &,
16028 const ::CellStatus,
16029 const boost::iterator_range<std::vector<char>::const_iterator> &)>
16033 this->data_serializer.unpack_data(this->local_cell_relations,
16038 --this->cell_attached_data.n_attached_data_sets;
16039 if (this->cell_attached_data.n_attached_deserialize > 0)
16040 --this->cell_attached_data.n_attached_deserialize;
16049 if (this->cell_attached_data.n_attached_data_sets == 0 &&
16050 this->cell_attached_data.n_attached_deserialize == 0)
16053 this->cell_attached_data.pack_callbacks_fixed.clear();
16054 this->cell_attached_data.pack_callbacks_variable.clear();
16055 this->data_serializer.clear();
16058 for (
auto &cell_rel : this->local_cell_relations)
16065template <
int dim,
int spacedim>
16068 const unsigned int global_first_cell,
16069 const unsigned int global_num_cells,
16070 const std::string &file_basename)
const
16073 auto tria =
const_cast<Triangulation<dim, spacedim> *
>(
this);
16075 if (this->cell_attached_data.n_attached_data_sets > 0)
16078 tria->data_serializer.pack_data(
16079 tria->local_cell_relations,
16080 tria->cell_attached_data.pack_callbacks_fixed,
16081 tria->cell_attached_data.pack_callbacks_variable,
16082 this->get_communicator());
16085 tria->data_serializer.save(global_first_cell,
16088 this->get_communicator());
16091 tria->data_serializer.clear();
16097 tria->cell_attached_data.n_attached_data_sets = 0;
16098 tria->cell_attached_data.pack_callbacks_fixed.clear();
16099 tria->cell_attached_data.pack_callbacks_variable.clear();
16104template <
int dim,
int spacedim>
16107 const unsigned int global_first_cell,
16108 const unsigned int global_num_cells,
16109 const unsigned int local_num_cells,
16110 const std::string &file_basename,
16111 const unsigned int n_attached_deserialize_fixed,
16112 const unsigned int n_attached_deserialize_variable)
16115 if (this->cell_attached_data.n_attached_deserialize > 0)
16117 this->data_serializer.load(global_first_cell,
16121 n_attached_deserialize_fixed,
16122 n_attached_deserialize_variable,
16123 this->get_communicator());
16125 this->data_serializer.unpack_cell_status(this->local_cell_relations);
16129 for (
const auto &cell_rel : this->local_cell_relations)
16132 Assert((cell_rel.second ==
16140template <
int dim,
int spacedim>
16148 vertices_used.clear();
16156 vertex_to_boundary_id_map_1d->clear();
16159 vertex_to_manifold_id_map_1d->clear();
16169 number_cache = internal::TriangulationImplementation::NumberCache<dim>();
16174template <
int dim,
int spacedim>
16176typename Triangulation<dim, spacedim>::DistortedCellList
16179 const DistortedCellList cells_with_distorted_children =
16180 this->policy->execute_refinement(*
this, check_for_distorted_cells);
16186 *
this, levels.size(), number_cache);
16189 for (
const auto &level : levels)
16196 for (
const auto &cell : this->cell_iterators())
16200 return cells_with_distorted_children;
16205template <
int dim,
int spacedim>
16211 const cell_iterator endc =
end();
16212 bool do_coarsen =
false;
16213 if (levels.size() >= 2)
16214 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
16215 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
16227 std::vector<unsigned int> line_cell_count(dim > 1 ? this->n_raw_lines() : 0);
16228 std::vector<unsigned int> quad_cell_count(dim > 2 ? this->n_raw_quads() : 0);
16230 for (
const auto &cell : this->cell_iterators())
16234 const auto line_indices = internal::TriaAccessorImplementation::
16235 Implementation::get_line_indices_of_cell(*cell);
16238 const unsigned int n_lines =
std::min(cell->n_lines(), 12u);
16239 for (
unsigned int l = 0;
l < n_lines; ++
l)
16240 ++line_cell_count[line_indices[l]];
16241 for (
const unsigned int q : cell->face_indices())
16242 ++quad_cell_count[cell->face_index(q)];
16245 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
16246 ++line_cell_count[cell->line(l)->index()];
16263 if (levels.size() >= 2)
16264 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
16265 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
16267 for (
unsigned int child = 0; child < cell->n_children(); ++child)
16269 Assert(cell->child(child)->coarsen_flag_set(),
16271 cell->child(child)->clear_coarsen_flag();
16274 signals.pre_coarsening_on_cell(cell);
16276 this->policy->delete_children(*
this,
16284 *
this, levels.size(), number_cache);
16289template <
int dim,
int spacedim>
16306 auto previous_coarsen_flags = internal::extract_raw_coarsen_flags(levels);
16308 bool continue_iterating =
true;
16312 if (smooth_grid & limit_level_difference_at_vertices)
16314 Assert(!anisotropic_refinement,
16315 ExcMessage(
"In case of anisotropic refinement the "
16316 "limit_level_difference_at_vertices flag for "
16317 "mesh smoothing must not be set!"));
16321 std::vector<int> vertex_level(vertices.size(), 0);
16322 for (
const auto &cell : this->active_cell_iterators())
16324 if (cell->refine_flag_set())
16325 for (
const unsigned int vertex :
16327 vertex_level[cell->vertex_index(vertex)] =
16328 std::max(vertex_level[cell->vertex_index(vertex)],
16329 cell->level() + 1);
16330 else if (!cell->coarsen_flag_set())
16331 for (
const unsigned int vertex :
16333 vertex_level[cell->vertex_index(vertex)] =
16334 std::max(vertex_level[cell->vertex_index(vertex)],
16345 for (
const unsigned int vertex :
16347 vertex_level[cell->vertex_index(vertex)] =
16348 std::max(vertex_level[cell->vertex_index(vertex)],
16349 cell->level() - 1);
16363 active_cell_iterator endc =
end();
16364 for (active_cell_iterator cell = last_active(); cell != endc; --cell)
16365 if (cell->refine_flag_set() ==
false)
16367 for (
const unsigned int vertex :
16369 if (vertex_level[cell->vertex_index(vertex)] >=
16373 cell->clear_coarsen_flag();
16378 if (vertex_level[cell->vertex_index(vertex)] >
16381 cell->set_refine_flag();
16383 for (
const unsigned int v :
16385 vertex_level[cell->vertex_index(v)] =
16386 std::max(vertex_level[cell->vertex_index(v)],
16387 cell->level() + 1);
16403 for (
const auto &acell : this->active_cell_iterators_on_level(0))
16404 acell->clear_coarsen_flag();
16406 const cell_iterator endc =
end();
16407 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
16410 if (cell->is_active())
16413 const unsigned int n_children = cell->n_children();
16414 unsigned int flagged_children = 0;
16415 for (
unsigned int child = 0; child < n_children; ++child)
16417 const auto child_cell = cell->child(child);
16418 if (child_cell->is_active() && child_cell->coarsen_flag_set())
16420 ++flagged_children;
16422 child_cell->clear_coarsen_flag();
16428 if (flagged_children == n_children &&
16429 this->policy->coarsening_allowed(cell))
16430 for (
unsigned int c = 0; c < n_children; ++c)
16432 Assert(cell->child(c)->refine_flag_set() ==
false,
16435 cell->child(c)->set_coarsen_flag();
16441 auto current_coarsen_flags = internal::extract_raw_coarsen_flags(levels);
16443 continue_iterating = (current_coarsen_flags != previous_coarsen_flags);
16444 previous_coarsen_flags.swap(current_coarsen_flags);
16446 while (continue_iterating ==
true);
16458 const auto flags_before = internal::extract_raw_coarsen_flags(
levels);
16463 const auto flags_after = internal::extract_raw_coarsen_flags(
levels);
16465 return (flags_before != flags_after);
16476 const auto flags_before = internal::extract_raw_coarsen_flags(
levels);
16481 const auto flags_after = internal::extract_raw_coarsen_flags(
levels);
16483 return (flags_before != flags_after);
16494 const auto flags_before = internal::extract_raw_coarsen_flags(
levels);
16499 const auto flags_after = internal::extract_raw_coarsen_flags(
levels);
16501 return (flags_before != flags_after);
16512 template <
int dim,
int spacedim>
16514 possibly_do_not_produce_unrefined_islands(
16519 unsigned int n_neighbors = 0;
16522 unsigned int count = 0;
16530 if (face_will_be_refined_by_neighbor(cell, n))
16537 if (count == n_neighbors ||
16538 (count >= n_neighbors - 1 &&
16541 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
16547 (cell_will_be_coarsened(cell->
neighbor(face))))
16548 possibly_do_not_produce_unrefined_islands<dim, spacedim>(
16564 template <
int dim,
int spacedim>
16566 possibly_refine_unrefined_island(
16568 const bool allow_anisotropic_smoothing)
16575 if (
dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim>
16586 "The triangulation is trying to avoid unrefined islands "
16587 "during mesh refinement/coarsening, as you had requested "
16588 " by passing the appropriate 'smoothing flags' to the "
16589 "constructor of the triangulation. However, for objects "
16590 "of type parallel::distributed::Triangulation, control "
16591 "over which cells get refined rests with p4est, not the "
16592 "deal.II triangulation, and consequently it is not "
16593 "always possible to avoid unrefined islands in the mesh. "
16594 "Please remove the constructor argument to the triangulation "
16595 "object that requests mesh smoothing."));
16607 if (allow_anisotropic_smoothing ==
false)
16610 unsigned int refined_neighbors = 0, unrefined_neighbors = 0;
16614 if (face_will_be_refined_by_neighbor(cell, face))
16615 ++refined_neighbors;
16617 ++unrefined_neighbors;
16620 if (unrefined_neighbors < refined_neighbors)
16628 if (unrefined_neighbors > 0)
16631 (face_will_be_refined_by_neighbor(cell, face) ==
false) &&
16634 possibly_refine_unrefined_island<dim, spacedim>(
16635 cell->
neighbor(face), allow_anisotropic_smoothing);
16642 RefinementCase<dim> smoothing_cell_refinement_case =
16647 for (
unsigned int face_pair = 0;
16648 face_pair < GeometryInfo<dim>::faces_per_cell / 2;
16654 RefinementCase<dim> directional_cell_refinement_case =
16657 for (
unsigned int face_index = 0; face_index < 2; ++face_index)
16659 unsigned int face = 2 * face_pair + face_index;
16662 RefinementCase<dim - 1> expected_face_ref_case =
16663 RefinementCase<dim - 1>::no_refinement;
16666 face_will_be_refined_by_neighbor<dim, spacedim>(
16667 cell, face, expected_face_ref_case);
16681 directional_cell_refinement_case =
16682 (directional_cell_refinement_case &
16685 expected_face_ref_case,
16695 Assert(directional_cell_refinement_case <
16698 smoothing_cell_refinement_case =
16699 smoothing_cell_refinement_case | directional_cell_refinement_case;
16704 if (smoothing_cell_refinement_case)
16708 smoothing_cell_refinement_case);
16715template <
int dim,
int spacedim>
16721 const auto coarsen_flags_before = internal::extract_raw_coarsen_flags(levels);
16722 const auto refine_flags_before = internal::extract_raw_refine_flags(levels);
16740 auto coarsen_flags_before_loop = coarsen_flags_before;
16741 auto refine_flags_before_loop = refine_flags_before;
16818 if (((smooth_grid & coarsest_level_1) || (smooth_grid & patch_level_1)) &&
16821 for (
const auto &cell : active_cell_iterators_on_level(1))
16825 bool mesh_changed_in_this_loop =
false;
16834 if (smooth_grid & do_not_produce_unrefined_islands &&
16835 !(smooth_grid & patch_level_1))
16837 for (
const auto &cell : cell_iterators())
16841 if (!cell->
is_active() && cell_will_be_coarsened(cell))
16842 possibly_do_not_produce_unrefined_islands<dim, spacedim>(cell);
16866 if (smooth_grid & (eliminate_refined_inner_islands |
16867 eliminate_refined_boundary_islands) &&
16868 !(smooth_grid & patch_level_1))
16870 for (
const auto &cell : cell_iterators())
16880 bool all_children_active =
true;
16882 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
16887 all_children_active =
false;
16891 if (all_children_active)
16898 unsigned int unrefined_neighbors = 0, total_neighbors = 0;
16909 bool at_periodic_boundary =
false;
16911 for (
const unsigned int n :
16914 const cell_iterator neighbor = cell->
neighbor(n);
16919 if (!face_will_be_refined_by_neighbor(cell, n))
16920 ++unrefined_neighbors;
16925 at_periodic_boundary =
true;
16940 if ((unrefined_neighbors == total_neighbors) &&
16942 (smooth_grid & eliminate_refined_inner_islands)) ||
16945 eliminate_refined_boundary_islands))) &&
16946 (total_neighbors != 0))
16949 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
16970 if (smooth_grid & limit_level_difference_at_vertices)
16972 Assert(!anisotropic_refinement,
16973 ExcMessage(
"In case of anisotropic refinement the "
16974 "limit_level_difference_at_vertices flag for "
16975 "mesh smoothing must not be set!"));
16979 std::vector<int> vertex_level(vertices.size(), 0);
16980 for (
const auto &cell : active_cell_iterators())
16983 for (
const unsigned int vertex :
16987 cell->
level() + 1);
16989 for (
const unsigned int vertex :
17001 for (
const unsigned int vertex :
17005 cell->
level() - 1);
17019 for (active_cell_iterator cell = last_active(); cell !=
end(); --cell)
17022 for (
const unsigned int vertex :
17038 for (
const unsigned int v :
17042 cell->
level() + 1);
17063 if (smooth_grid & eliminate_unrefined_islands)
17065 for (active_cell_iterator cell = last_active(); cell !=
end(); --cell)
17070 possibly_refine_unrefined_island<dim, spacedim>(
17071 cell, (smooth_grid & allow_anisotropic_smoothing) != 0);
17099 if (smooth_grid & patch_level_1)
17109 for (
const auto &cell : cell_iterators())
17121 RefinementCase<dim> combined_ref_case =
17123 for (
unsigned int i = 0; i < cell->
n_children(); ++i)
17124 combined_ref_case =
17127 for (
unsigned int i = 0; i < cell->
n_children(); ++i)
17129 cell_iterator child = cell->
child(i);
17131 child->clear_coarsen_flag();
17132 child->set_refine_flag(combined_ref_case);
17145 for (
const auto &cell : cell_iterators())
17158 const unsigned int n_children = cell->
n_children();
17159 bool has_active_grandchildren =
false;
17161 for (
unsigned int i = 0; i < n_children; ++i)
17164 has_active_grandchildren =
true;
17168 if (has_active_grandchildren ==
false)
17174 unsigned int n_grandchildren = 0;
17177 unsigned int n_coarsen_flags = 0;
17182 for (
unsigned int c = 0; c < n_children; ++c)
17187 cell_iterator child = cell->
child(c);
17189 const unsigned int nn_children = child->n_children();
17190 n_grandchildren += nn_children;
17195 if (child->child(0)->is_active())
17196 for (
unsigned int cc = 0; cc < nn_children; ++cc)
17197 if (child->child(cc)->coarsen_flag_set())
17212 if ((n_coarsen_flags != n_grandchildren) && (n_coarsen_flags > 0))
17213 for (
unsigned int c = 0; c < n_children; ++c)
17215 const cell_iterator child = cell->
child(c);
17216 if (child->child(0)->is_active())
17217 for (
unsigned int cc = 0; cc < child->n_children(); ++cc)
17218 child->child(cc)->clear_coarsen_flag();
17229 this->policy->prevent_distorted_boundary_cells(*
this);
17245 bool changed =
true;
17249 active_cell_iterator cell = last_active(), endc =
end();
17251 for (; cell != endc; --cell)
17252 if (cell->refine_flag_set())
17255 for (
const auto i : cell->face_indices())
17260 const bool has_periodic_neighbor =
17261 cell->has_periodic_neighbor(i);
17262 const bool has_neighbor_or_periodic_neighbor =
17263 !cell->at_boundary(i) || has_periodic_neighbor;
17264 if (has_neighbor_or_periodic_neighbor &&
17266 cell->refine_flag_set(), i) !=
17279 if (cell->neighbor_or_periodic_neighbor(i)->is_active())
17281 if ((!has_periodic_neighbor &&
17282 cell->neighbor_is_coarser(i)) ||
17283 (has_periodic_neighbor &&
17284 cell->periodic_neighbor_is_coarser(i)))
17286 if (cell->neighbor_or_periodic_neighbor(i)
17287 ->coarsen_flag_set())
17288 cell->neighbor_or_periodic_neighbor(i)
17289 ->clear_coarsen_flag();
17306 allow_anisotropic_smoothing)
17308 has_periodic_neighbor ?
17309 cell->periodic_neighbor(i)
17310 ->flag_for_face_refinement(
17312 ->periodic_neighbor_of_coarser_periodic_neighbor(
17315 RefinementCase<dim - 1>::cut_x) :
17317 ->flag_for_face_refinement(
17319 ->neighbor_of_coarser_neighbor(
17322 RefinementCase<dim - 1>::cut_x);
17326 ->neighbor_or_periodic_neighbor(
17328 ->refine_flag_set())
17330 cell->neighbor_or_periodic_neighbor(i)
17331 ->set_refine_flag();
17463 std::pair<unsigned int, unsigned int>
17465 has_periodic_neighbor ?
17467 ->periodic_neighbor_of_coarser_periodic_neighbor(
17469 cell->neighbor_of_coarser_neighbor(i);
17470 unsigned int refined_along_x = 0,
17471 refined_along_y = 0,
17472 to_be_refined_along_x = 0,
17473 to_be_refined_along_y = 0;
17475 const int this_face_index =
17476 cell->face_index(i);
17485 const auto parent_face = [&]() {
17486 if (has_periodic_neighbor)
17488 const auto neighbor =
17489 cell->periodic_neighbor(i);
17490 const auto parent_face_no =
17492 ->periodic_neighbor_of_periodic_neighbor(
17495 neighbor->periodic_neighbor(
17497 return parent->face(parent_face_no);
17500 return cell->neighbor(i)->face(
17504 if ((this_face_index ==
17505 parent_face->child_index(0)) ||
17506 (this_face_index ==
17507 parent_face->child_index(1)))
17515 RefinementCase<dim - 1> frc =
17516 parent_face->refinement_case();
17517 if (frc & RefinementCase<dim>::cut_x)
17519 if (frc & RefinementCase<dim>::cut_y)
17532 RefinementCase<dim - 1> flagged_frc =
17534 cell->refine_flag_set(),
17536 cell->face_orientation(i),
17537 cell->face_flip(i),
17538 cell->face_rotation(i));
17540 RefinementCase<dim>::cut_x)
17541 ++to_be_refined_along_x;
17543 RefinementCase<dim>::cut_y)
17544 ++to_be_refined_along_y;
17549 allow_anisotropic_smoothing) ||
17550 cell->neighbor_or_periodic_neighbor(i)
17551 ->refine_flag_set())
17553 if (refined_along_x +
17554 to_be_refined_along_x >
17558 ->neighbor_or_periodic_neighbor(i)
17559 ->flag_for_face_refinement(
17561 RefinementCase<dim -
17563 if (refined_along_y +
17564 to_be_refined_along_y >
17568 ->neighbor_or_periodic_neighbor(i)
17569 ->flag_for_face_refinement(
17571 RefinementCase<dim -
17577 ->neighbor_or_periodic_neighbor(i)
17578 ->refine_flag_set() !=
17582 cell->neighbor_or_periodic_neighbor(i)
17583 ->set_refine_flag();
17589 cell->neighbor_or_periodic_neighbor(i);
17590 RefinementCase<dim - 1> nb_frc =
17592 nb->refine_flag_set(),
17594 nb->face_orientation(nb_indices.first),
17595 nb->face_flip(nb_indices.first),
17596 nb->face_rotation(nb_indices.first));
17597 if ((nb_frc & RefinementCase<dim>::cut_x) &&
17598 !((refined_along_x != 0u) ||
17599 (to_be_refined_along_x != 0u)))
17600 changed |= cell->flag_for_face_refinement(
17603 if ((nb_frc & RefinementCase<dim>::cut_y) &&
17604 !((refined_along_y != 0u) ||
17605 (to_be_refined_along_y != 0u)))
17606 changed |= cell->flag_for_face_refinement(
17613 cell->neighbor_or_periodic_neighbor(i)
17614 ->clear_coarsen_flag();
17615 const unsigned int nb_nb =
17616 has_periodic_neighbor ?
17618 ->periodic_neighbor_of_periodic_neighbor(
17620 cell->neighbor_of_neighbor(i);
17621 const cell_iterator neighbor =
17622 cell->neighbor_or_periodic_neighbor(i);
17623 RefinementCase<dim - 1> face_ref_case =
17625 neighbor->refine_flag_set(),
17627 neighbor->face_orientation(nb_nb),
17628 neighbor->face_flip(nb_nb),
17629 neighbor->face_rotation(nb_nb));
17630 RefinementCase<dim - 1> needed_face_ref_case =
17632 cell->refine_flag_set(),
17634 cell->face_orientation(i),
17635 cell->face_flip(i),
17636 cell->face_rotation(i));
17641 if ((face_ref_case ==
17642 RefinementCase<dim>::cut_x &&
17643 needed_face_ref_case ==
17644 RefinementCase<dim>::cut_y) ||
17646 RefinementCase<dim>::cut_y &&
17647 needed_face_ref_case ==
17648 RefinementCase<dim>::cut_x))
17650 changed = cell->flag_for_face_refinement(
17652 neighbor->flag_for_face_refinement(
17653 nb_nb, needed_face_ref_case);
17659 RefinementCase<dim - 1>
17660 face_ref_case = cell->face(i)->refinement_case(),
17661 needed_face_ref_case =
17663 cell->refine_flag_set(),
17665 cell->face_orientation(i),
17666 cell->face_flip(i),
17667 cell->face_rotation(i));
17671 if ((face_ref_case == RefinementCase<dim>::cut_x &&
17672 needed_face_ref_case ==
17673 RefinementCase<dim>::cut_y) ||
17674 (face_ref_case == RefinementCase<dim>::cut_y &&
17675 needed_face_ref_case ==
17676 RefinementCase<dim>::cut_x))
17678 cell->flag_for_face_refinement(i,
17690 this->policy->prepare_refinement_dim_dependent(*
this);
17696 fix_coarsen_flags();
17699 auto coarsen_flags_after_loop =
17700 internal::extract_raw_coarsen_flags(levels);
17701 auto refine_flags_after_loop = internal::extract_raw_refine_flags(levels);
17704 mesh_changed_in_this_loop =
17705 ((coarsen_flags_before_loop != coarsen_flags_after_loop) ||
17706 (refine_flags_before_loop != refine_flags_after_loop));
17709 coarsen_flags_before_loop.swap(coarsen_flags_after_loop);
17710 refine_flags_before_loop.swap(refine_flags_after_loop);
17712 while (mesh_changed_in_this_loop);
17718 return ((coarsen_flags_before != coarsen_flags_before_loop) ||
17719 (refine_flags_before != refine_flags_before_loop));
17724template <
int dim,
int spacedim>
17727 const unsigned int magic_number1,
17728 const std::vector<bool> &v,
17729 const unsigned int magic_number2,
17732 const unsigned int N = v.size();
17733 unsigned char *flags =
new unsigned char[N / 8 + 1];
17734 for (
unsigned int i = 0; i < N / 8 + 1; ++i)
17737 for (
unsigned int position = 0; position < N; ++position)
17738 flags[position / 8] |= (v[position] ? (1 << (position % 8)) : 0);
17747 out << magic_number1 <<
' ' << N << std::endl;
17748 for (
unsigned int i = 0; i < N / 8 + 1; ++i)
17749 out <<
static_cast<unsigned int>(flags[i]) <<
' ';
17751 out << std::endl << magic_number2 << std::endl;
17759template <
int dim,
int spacedim>
17762 const unsigned int magic_number1,
17763 std::vector<bool> &v,
17764 const unsigned int magic_number2,
17769 unsigned int magic_number;
17770 in >> magic_number;
17771 AssertThrow(magic_number == magic_number1, ExcGridReadError());
17777 unsigned char *flags =
new unsigned char[N / 8 + 1];
17778 unsigned short int tmp;
17779 for (
unsigned int i = 0; i < N / 8 + 1; ++i)
17785 for (
unsigned int position = 0; position != N; ++position)
17786 v[position] = ((flags[position / 8] & (1 << (position % 8))) != 0);
17788 in >> magic_number;
17789 AssertThrow(magic_number == magic_number2, ExcGridReadError());
17798template <
int dim,
int spacedim>
17802 std::size_t mem = 0;
17803 mem +=
sizeof(MeshSmoothing);
17807 for (
const auto &level : levels)
17811 mem +=
sizeof(manifolds);
17812 mem +=
sizeof(smooth_grid);
17814 mem +=
sizeof(faces);
17823template <
int dim,
int spacedim>
17831#include "tria.inst"