560 const
unsigned int handle,
564 const
boost::iterator_range<
std::vector<
char>::const_iterator> &)>
565 &unpack_callback)
const
571 const bool callback_variable_transfer = (handle % 2 == 0);
572 const unsigned int callback_index = handle / 2;
580 if (cell_relations.size() > 0)
586 std::vector<char>::const_iterator dest_data_it;
587 std::vector<char>::const_iterator dest_sizes_cell_it;
597 if (callback_variable_transfer)
610 const unsigned int offset_variable_data_sizes =
611 sizes_fixed_cumulative[sizes_fixed_cumulative.size() - 2];
618 dest_sizes_cell_it = dest_data_fixed.cbegin() +
619 offset_variable_data_sizes +
620 callback_index *
sizeof(
unsigned int);
623 dest_data_it = dest_data_variable.cbegin();
630 offset = sizes_fixed_cumulative[callback_index];
631 size = sizes_fixed_cumulative[callback_index + 1] - offset;
632 data_increment = sizes_fixed_cumulative.back();
638 if (cell_relations.begin() != cell_relations.end())
639 dest_data_it = dest_data_fixed.cbegin() + offset;
643 auto cell_rel_it = cell_relations.begin();
644 auto dest_sizes_it = dest_sizes_variable.cbegin();
645 for (; cell_rel_it != cell_relations.end(); ++cell_rel_it)
647 const auto &dealii_cell = cell_rel_it->first;
648 const auto &cell_status = cell_rel_it->second;
650 if (callback_variable_transfer)
654 data_increment = *dest_sizes_it;
661 if (callback_index == 0)
665 &(*(dest_sizes_cell_it -
sizeof(
unsigned int))),
666 sizeof(
unsigned int));
669 &(*dest_sizes_cell_it),
670 sizeof(
unsigned int));
677 dest_data_it += offset;
678 data_increment -= offset;
683 if (cell_rel_it != cell_relations.end() - 1)
684 dest_sizes_cell_it += sizes_fixed_cumulative.back();
692 unpack_callback(dealii_cell,
694 boost::make_iterator_range(dest_data_it,
695 dest_data_it + size));
699 unpack_callback(dealii_cell->parent(),
701 boost::make_iterator_range(dest_data_it,
702 dest_data_it + size));
714 if (cell_rel_it != cell_relations.end() - 1)
715 dest_data_it += data_increment;
721 template <
int dim,
int spacedim>
724 const
unsigned int global_first_cell,
725 const
unsigned int global_num_cells,
726 const
std::
string &file_basename,
727 const
MPI_Comm &mpi_communicator)
const
732#ifdef DEAL_II_WITH_MPI
737 const unsigned int myrank =
739 const unsigned int mpisize =
750 const std::string fname_fixed =
751 std::string(file_basename) +
"_fixed.data";
754 int ierr = MPI_Info_create(&info);
758 ierr = MPI_File_open(mpi_communicator,
760 MPI_MODE_CREATE | MPI_MODE_WRONLY,
765 ierr = MPI_File_set_size(fh, 0);
769 ierr = MPI_Barrier(mpi_communicator);
771 ierr = MPI_Info_free(&info);
792 const MPI_Offset size_header =
797 const MPI_Offset my_global_file_position =
799 static_cast<MPI_Offset
>(global_first_cell) * bytes_per_cell;
803 my_global_file_position,
810 ierr = MPI_File_close(&fh);
821 const std::string fname_variable =
822 std::string(file_basename) +
"_variable.data";
825 int ierr = MPI_Info_create(&info);
829 ierr = MPI_File_open(mpi_communicator,
830 fname_variable.c_str(),
831 MPI_MODE_CREATE | MPI_MODE_WRONLY,
836 ierr = MPI_File_set_size(fh, 0);
840 ierr = MPI_Barrier(mpi_communicator);
842 ierr = MPI_Info_free(&info);
847 const MPI_Offset my_global_file_position =
848 static_cast<MPI_Offset
>(global_first_cell) *
849 sizeof(
unsigned int);
854 static_cast<std::size_t
>(
855 std::numeric_limits<int>::max()),
860 my_global_file_position,
872 std::uint64_t prefix_sum = 0;
873 ierr = MPI_Exscan(&size_on_proc,
881 const MPI_Offset my_global_file_position =
882 static_cast<MPI_Offset
>(global_num_cells) *
sizeof(
unsigned int) +
888 my_global_file_position,
896 ierr = MPI_File_close(&fh);
903 (void)global_first_cell;
904 (void)global_num_cells;
905 (void)mpi_communicator;
911 const std::string fname_fixed =
912 std::string(file_basename) +
"_fixed.data";
914 std::ofstream file(fname_fixed, std::ios::binary | std::ios::out);
918 file.write(
reinterpret_cast<const char *
>(
932 const std::string fname_variable =
933 std::string(file_basename) +
"_variable.data";
935 std::ofstream file(fname_variable,
936 std::ios::binary | std::ios::out);
940 file.write(
reinterpret_cast<const char *
>(
952 template <
int dim,
int spacedim>
955 const
unsigned int global_first_cell,
956 const
unsigned int global_num_cells,
957 const
unsigned int local_num_cells,
958 const
std::
string &file_basename,
959 const
unsigned int n_attached_deserialize_fixed,
960 const
unsigned int n_attached_deserialize_variable,
964 ExcMessage(
"Previously loaded data has not been released yet!"));
968#ifdef DEAL_II_WITH_MPI
973 const unsigned int mpisize =
982 const std::string fname_fixed =
983 std::string(file_basename) +
"_fixed.data";
986 int ierr = MPI_Info_create(&info);
990 ierr = MPI_File_open(
991 mpi_communicator, fname_fixed.c_str(), MPI_MODE_RDONLY, info, &fh);
994 ierr = MPI_Info_free(&info);
1018 const MPI_Offset size_header =
1023 const MPI_Offset my_global_file_position =
1025 static_cast<MPI_Offset
>(global_first_cell) * bytes_per_cell;
1029 my_global_file_position,
1037 ierr = MPI_File_close(&fh);
1046 const std::string fname_variable =
1047 std::string(file_basename) +
"_variable.data";
1050 int ierr = MPI_Info_create(&info);
1054 ierr = MPI_File_open(mpi_communicator,
1055 fname_variable.c_str(),
1061 ierr = MPI_Info_free(&info);
1067 const MPI_Offset my_global_file_position_sizes =
1068 static_cast<MPI_Offset
>(global_first_cell) *
sizeof(
unsigned int);
1072 my_global_file_position_sizes,
1082 const std::uint64_t size_on_proc =
1087 std::uint64_t prefix_sum = 0;
1088 ierr = MPI_Exscan(&size_on_proc,
1096 const MPI_Offset my_global_file_position =
1097 static_cast<MPI_Offset
>(global_num_cells) *
sizeof(
unsigned int) +
1104 my_global_file_position,
1111 ierr = MPI_File_close(&fh);
1118 (void)mpi_communicator;
1119 (void)global_first_cell;
1120 (void)global_num_cells;
1126 const std::string fname_fixed =
1127 std::string(file_basename) +
"_fixed.data";
1129 std::ifstream file(fname_fixed, std::ios::binary | std::ios::in);
1152 const std::string fname_variable =
1153 std::string(file_basename) +
"_variable.data";
1155 std::ifstream file(fname_variable, std::ios::binary | std::ios::in);
1164 const std::uint64_t size =
1176 template <
int dim,
int spacedim>
1218 template <
int dim,
int spacedim>
1220 cell_is_patch_level_1(
1225 unsigned int n_active_children = 0;
1226 for (
unsigned int i = 0; i < cell->n_children(); ++i)
1227 if (cell->child(i)->is_active())
1228 ++n_active_children;
1230 return (n_active_children == 0) ||
1231 (n_active_children == cell->n_children());
1241 template <
int dim,
int spacedim>
1243 cell_will_be_coarsened(
1249 if (cell->has_children())
1251 unsigned int children_to_coarsen = 0;
1252 const unsigned int n_children = cell->n_children();
1254 for (
unsigned int c = 0; c < n_children; ++c)
1255 if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
1256 ++children_to_coarsen;
1257 if (children_to_coarsen == n_children)
1260 for (
unsigned int c = 0; c < n_children; ++c)
1261 if (cell->child(c)->is_active())
1262 cell->child(c)->clear_coarsen_flag();
1298 template <
int dim,
int spacedim>
1300 face_will_be_refined_by_neighbor_internal(
1302 const unsigned int face_no,
1311 cell->neighbor(face_no);
1323 if (cell_will_be_coarsened(neighbor))
1332 expected_face_ref_case = cell->face(face_no)->refinement_case();
1344 const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
1358 neighbor_face = neighbor->
face(neighbor_neighbor);
1359 const int this_face_index = cell->face_index(face_no);
1365 if (neighbor_face->index() == this_face_index)
1371 expected_face_ref_case = face_ref_case;
1389 for (
unsigned int c = 0; c < neighbor_face->n_children(); ++c)
1390 if (neighbor_face->child_index(c) == this_face_index)
1399 if ((neighbor_face->refinement_case() | face_ref_case) ==
1400 neighbor_face->refinement_case())
1419 expected_face_ref_case =
1420 ~neighbor_face->refinement_case();
1447 template <
int dim,
int spacedim>
1449 face_will_be_refined_by_neighbor(
1451 const unsigned int face_no)
1454 return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
1461 template <
int dim,
int spacedim>
1463 face_will_be_refined_by_neighbor(
1465 const unsigned int face_no,
1468 return face_will_be_refined_by_neighbor_internal(cell,
1470 expected_face_ref_case);
1475 template <
int dim,
int spacedim>
1477 satisfies_level1_at_vertex_rule(
1480 std::vector<unsigned int> min_adjacent_cell_level(
1482 std::vector<unsigned int> max_adjacent_cell_level(
1486 for (
const unsigned int v : cell->vertex_indices())
1488 min_adjacent_cell_level[cell->vertex_index(v)] =
1490 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
1491 max_adjacent_cell_level[cell->vertex_index(v)] =
1493 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
1496 for (
unsigned int k = 0; k < triangulation.
n_vertices(); ++k)
1498 if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
1522 template <
int dim,
int spacedim>
1524 middle_vertex_index(
1527 if (line->has_children())
1528 return line->child(0)->vertex_index(1);
1533 template <
int dim,
int spacedim>
1535 middle_vertex_index(
1538 switch (
static_cast<unsigned char>(quad->refinement_case()))
1541 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
1544 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
1547 return quad->child(0)->vertex_index(3);
1556 template <
int dim,
int spacedim>
1558 middle_vertex_index(
1561 switch (
static_cast<unsigned char>(hex->refinement_case()))
1564 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
1567 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
1570 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
1573 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
1576 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
1579 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
1582 return hex->child(0)->vertex_index(7);
1603 template <
class TRIANGULATION>
1604 typename TRIANGULATION::DistortedCellList
1605 collect_distorted_coarse_cells(
const TRIANGULATION &)
1607 return typename TRIANGULATION::DistortedCellList();
1629 vertices[i] = cell->vertex(i);
1635 if (determinants[i] <=
1643 return distorted_cells;
1655 has_distorted_children(
1660 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
1670 if (determinants[i] <=
1686 template <
int dim,
int spacedim>
1688 has_distorted_children(
1695 template <
int dim,
int spacedim>
1697 update_periodic_face_map_recursively(
1700 unsigned int n_face_1,
1701 unsigned int n_face_2,
1711 const FaceIterator face_1 = cell_1->
face(n_face_1);
1712 const FaceIterator face_2 = cell_2->
face(n_face_2);
1714 const auto inverse_orientation =
1715 face_1->reference_cell().get_inverse_combined_orientation(orientation);
1719 const auto [face_orientation, face_rotation, face_flip] =
1722 Assert((dim != 1) || (face_orientation ==
true && face_flip ==
false &&
1723 face_rotation ==
false),
1725 "(face_orientation, face_flip, face_rotation) "
1726 "is invalid for 1d"));
1728 Assert((dim != 2) || (face_flip ==
false && face_rotation ==
false),
1730 "(face_orientation, face_flip, face_rotation) "
1731 "is invalid for 2d"));
1736 Assert(face_1->at_boundary() && face_2->at_boundary(),
1737 ExcMessage(
"Periodic faces must be on the boundary"));
1748 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
1750 const CellFace cell_face_1(cell_1, n_face_1);
1751 const CellFace cell_face_2(cell_2, n_face_2);
1752 const std::pair<CellFace, types::geometric_orientation>
1753 cell_face_orientation_2(cell_face_2, orientation);
1755 const std::pair<CellFace, std::pair<CellFace, types::geometric_orientation>>
1756 periodic_faces(cell_face_1, cell_face_orientation_2);
1760 periodic_face_map.insert(periodic_faces);
1768 update_periodic_face_map_recursively<dim, spacedim>(
1769 cell_1->
child(n_face_1),
1770 cell_2->
child(n_face_2),
1778 update_periodic_face_map_recursively<dim, spacedim>(
1779 cell_1->
child(n_face_1),
1798 Assert(face_1->n_children() ==
1800 face_2->n_children() ==
1806 for (
unsigned int i = 0;
1807 i < GeometryInfo<dim>::max_children_per_face;
1811 const unsigned int j =
1813 i, n_face_1, inverse_orientation);
1816 unsigned int child_cell_1 =
1824 face_1->refinement_case());
1825 unsigned int child_cell_2 =
1833 face_2->refinement_case());
1844 update_periodic_face_map_recursively<dim, spacedim>(
1845 cell_1->
child(child_cell_1),
1846 cell_2->
child(child_cell_2),
1855 for (
unsigned int i = 0;
1856 i < GeometryInfo<dim>::max_children_per_face;
1860 unsigned int child_cell_1 =
1868 face_1->refinement_case());
1871 update_periodic_face_map_recursively<dim, spacedim>(
1872 cell_1->
child(child_cell_1),
1887 child_line_index(
const unsigned int child_no,
1896 return line_orientation == D ? 0 : 1;
1898 return line_orientation == D ? 1 : 0;
1921 <<
"Something went wrong upon construction of cell "
1935 <<
" has negative measure. This typically "
1936 <<
"indicates some distortion in the cell, or a mistakenly "
1937 <<
"swapped pair of vertices in the input to "
1938 <<
"Triangulation::create_triangulation().");
1950 <<
"Error while creating cell " << arg1
1951 <<
": the vertex index " << arg2 <<
" must be between 0 and "
1962 <<
"The input data for creating a triangulation contained "
1963 <<
"information about a line with indices " << arg1 <<
" and " << arg2
1964 <<
" that is described to have boundary indicator "
1965 <<
static_cast<int>(arg3)
1966 <<
". However, this is an internal line not located on the "
1967 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
1969 <<
"If this happened at a place where you call "
1970 <<
"Triangulation::create_triangulation() yourself, you need "
1971 <<
"to check the SubCellData object you pass to this function."
1974 <<
"If this happened in a place where you are reading a mesh "
1975 <<
"from a file, then you need to investigate why such a line "
1976 <<
"ended up in the input file. A typical case is a geometry "
1977 <<
"that consisted of multiple parts and for which the mesh "
1978 <<
"generator program assumes that the interface between "
1979 <<
"two parts is a boundary when that isn't supposed to be "
1980 <<
"the case, or where the mesh generator simply assigns "
1981 <<
"'geometry indicators' to lines at the perimeter of "
1982 <<
"a part that are not supposed to be interpreted as "
1983 <<
"'boundary indicators'.");
1995 <<
"The input data for creating a triangulation contained "
1996 <<
"information about a quad with indices " << arg1 <<
", " << arg2
1997 <<
", " << arg3 <<
", and " << arg4
1998 <<
" that is described to have boundary indicator "
1999 <<
static_cast<int>(arg5)
2000 <<
". However, this is an internal quad not located on the "
2001 <<
"boundary. You cannot assign a boundary indicator to it." << std::endl
2003 <<
"If this happened at a place where you call "
2004 <<
"Triangulation::create_triangulation() yourself, you need "
2005 <<
"to check the SubCellData object you pass to this function."
2008 <<
"If this happened in a place where you are reading a mesh "
2009 <<
"from a file, then you need to investigate why such a quad "
2010 <<
"ended up in the input file. A typical case is a geometry "
2011 <<
"that consisted of multiple parts and for which the mesh "
2012 <<
"generator program assumes that the interface between "
2013 <<
"two parts is a boundary when that isn't supposed to be "
2014 <<
"the case, or where the mesh generator simply assigns "
2015 <<
"'geometry indicators' to quads at the surface of "
2016 <<
"a part that are not supposed to be interpreted as "
2017 <<
"'boundary indicators'.");
2026 <<
"In SubCellData the line info of the line with vertex indices " << arg1
2027 <<
" and " << arg2 <<
" appears more than once. "
2028 <<
"This is not allowed.");
2038 <<
"In SubCellData the line info of the line with vertex indices " << arg1
2039 <<
" and " << arg2 <<
" appears multiple times with different (valid) "
2040 << arg3 <<
". This is not allowed.");
2052 <<
"In SubCellData the quad info of the quad with line indices " << arg1
2053 <<
", " << arg2 <<
", " << arg3 <<
" and " << arg4
2054 <<
" appears multiple times with different (valid) " << arg5
2055 <<
". This is not allowed.");
2063 <<
"The containers have sizes " << arg1 <<
" and " << arg2
2064 <<
", which is not as expected.");
2071 template <
int dim,
int spacedim>
2074 const unsigned int true_dimension)
2077 tria_level.neighbors.
size(),
2079 tria_level.neighbors.
size()));
2081 tria_level.neighbors.
size(),
2083 tria_level.neighbors.
size()));
2098 tria_object.
used.size()));
2129 template <
int dim,
int spacedim>
2151 std::vector<unsigned int> &line_cell_count,
2152 std::vector<unsigned int> &quad_cell_count) = 0;
2159 const bool check_for_distorted_cells) = 0;
2188 virtual std::unique_ptr<Policy<dim, spacedim>>
2199 template <
int dim,
int spacedim,
typename T>
2206 T::update_neighbors(tria);
2213 std::vector<unsigned int> &line_cell_count,
2214 std::vector<unsigned int> &quad_cell_count)
override
2216 T::delete_children(tria, cell, line_cell_count, quad_cell_count);
2221 const bool check_for_distorted_cells)
override
2223 return T::execute_refinement(triangulation, check_for_distorted_cells);
2230 T::prevent_distorted_boundary_cells(triangulation);
2237 T::prepare_refinement_dim_dependent(triangulation);
2248 std::unique_ptr<Policy<dim, spacedim>>
2251 return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
2282 const std::vector<unsigned int> &
columns)
2292 std::vector<unsigned int> &&
columns)
2381 for (
const auto &reference_cell :
cell_types)
2386 const std::set<ReferenceCell<dim>> &
2392 std::vector<types::geometric_orientation> &
2403 const std::vector<types::geometric_orientation> &
2414 template <
int structdim>
2415 std::vector<ReferenceCell<structdim>> &
2418 if constexpr (structdim == dim)
2420 else if constexpr ((structdim == 2) && (dim == 3))
2430 template <
int structdim>
2431 const std::vector<ReferenceCell<structdim>> &
2434 if constexpr (structdim == dim)
2436 else if constexpr ((structdim == 2) && (dim == 3))
2449 if (from == dim && to == dim)
2451 else if (from == dim && to == dim - 1)
2453 else if (dim == 3 && from == 2 && to == 0)
2455 else if (dim == 3 && from == 2 && to == 1)
2457 else if (from == 1 && to == 0)
2468 if (from == dim && to == dim)
2470 else if (from == dim && to == dim - 1)
2472 else if (dim == 3 && from == 2 && to == 0)
2474 else if (dim == 3 && from == 2 && to == 1)
2476 else if (from == 1 && to == 0)
2518 unsigned int n_faces = 0;
2519 for (
unsigned int cell = 0; cell < cells_to_faces.
size(); ++cell)
2522 cells_to_faces[cell];
2524 *std::max_element(face_indices.
begin(),
2525 face_indices.
end()) +
2529 std::vector<std::size_t> cells_to_cells_offsets(cells_to_faces.
size() + 1,
2531 for (
unsigned int cell = 0; cell < cells_to_faces.
size(); ++cell)
2534 cells_to_faces[cell];
2535 cells_to_cells_offsets[cell + 1] =
2536 cells_to_cells_offsets[cell] + face_indices.
size();
2538 std::vector<unsigned int> cells_to_cells_columns(
2539 cells_to_cells_offsets.back(),
2542 std::vector<std::pair<unsigned int, unsigned int>> neighbors(
2548 unsigned int global_face_index = 0;
2549 for (
unsigned int cell = 0; cell < cells_to_faces.
size(); ++cell)
2552 const auto faces = cells_to_faces[cell];
2553 for (
unsigned int f = 0; f < faces.size(); ++f, ++global_face_index)
2559 neighbors[faces[f]] = std::make_pair(cell, global_face_index);
2566 cells_to_cells_columns[global_face_index] =
2567 neighbors[faces[f]].first;
2568 cells_to_cells_columns[neighbors[faces[f]].second] = cell;
2574 std::move(cells_to_cells_columns));
2587 template <
int max_n_vertices,
int face_dim,
int dim,
typename FU>
2594 std::vector<types::geometric_orientation> &orientations,
2595 const FU &second_key_function)
2597 const bool compatibility_mode =
true;
2601 std::vector<std::size_t> offsets_0;
2602 std::vector<unsigned int> columns_0;
2604 unsigned int n_entities = 0;
2606 for (
const auto &c : cell_types)
2610 Assert((face_dim == dim - 1) || ((dim == 3) && (face_dim == 1)),
2612 n_entities += (face_dim == dim - 1 ? c.n_faces() : c.n_lines());
2621 std::tuple<std::array<unsigned int, max_n_vertices>,
unsigned int>>
2624 std::vector<std::array<unsigned int, max_n_vertices>> ad_entity_vertices;
2625 std::vector<ReferenceCell<face_dim>> ad_entity_types;
2626 std::vector<std::array<unsigned int, max_n_vertices>> ad_compatibility;
2628 keys.reserve(n_entities);
2629 ad_entity_vertices.reserve(n_entities);
2630 ad_entity_types.reserve(n_entities);
2631 ad_compatibility.reserve(n_entities);
2633 std::vector<std::size_t> cells_to_faces_offsets(cell_types.size() + 1);
2634 cells_to_faces_offsets[0] = 0;
2636 constexpr unsigned int offset = 1;
2639 for (
unsigned int c = 0, counter = 0; c < cell_types.size(); ++c)
2641 const auto &cell_type = cell_types[c];
2645 Assert((face_dim == dim - 1) || ((dim == 3) && (face_dim == 1)),
2647 const unsigned int n_face_entities =
2648 (face_dim == dim - 1 ? cell_type.n_faces() : cell_type.n_lines());
2649 cells_to_faces_offsets[c + 1] =
2650 cells_to_faces_offsets[c] + n_face_entities;
2654 cells_to_vertices[c];
2655 for (
unsigned int e = 0; e < n_face_entities; ++e)
2658 std::array<unsigned int, max_n_vertices> entity_vertices;
2659 std::fill(entity_vertices.begin(), entity_vertices.end(), 0);
2664 Assert((face_dim == dim - 1) || ((dim == 3) && (face_dim == 1)),
2666 for (
unsigned int i = 0;
2667 i < (face_dim == dim - 1 ?
2668 cell_type.face_reference_cell(e).n_vertices() :
2671 entity_vertices[i] =
2672 local_vertices[face_dim == dim - 1 ?
2673 cell_type.face_to_cell_vertices(
2677 cell_type.line_to_cell_vertices(e, i)] +
2681 std::array<unsigned int, max_n_vertices> key = entity_vertices;
2682 std::sort(key.begin(), key.end());
2683 keys.emplace_back(key, counter++);
2685 ad_entity_vertices.emplace_back(entity_vertices);
2687 if constexpr (face_dim == dim - 1)
2688 ad_entity_types.emplace_back(cell_type.face_reference_cell(e));
2689 else if constexpr (face_dim == dim - 2)
2695 if constexpr (face_dim == 0)
2697 else if constexpr (face_dim == 1)
2707 if (compatibility_mode)
2708 ad_compatibility.emplace_back(
2709 second_key_function(entity_vertices, cell_type, c, e));
2713 std::vector<unsigned int> cells_to_faces_indices(keys.size());
2718 std::sort(keys.begin(), keys.end());
2720 if (compatibility_mode)
2722 unsigned int n_unique_entities = 0;
2723 unsigned int n_unique_entity_vertices = 0;
2725 std::array<unsigned int, max_n_vertices> ref_key, new_key;
2726 std::fill(ref_key.begin(), ref_key.end(), 0);
2727 for (
unsigned int i = 0; i < keys.size(); ++i)
2729 const auto offset_i = std::get<1>(keys[i]);
2731 if (ref_key != std::get<0>(keys[i]))
2733 ref_key = std::get<0>(keys[i]);
2735 ++n_unique_entities;
2736 n_unique_entity_vertices +=
2737 ad_entity_types[offset_i].n_vertices();
2739 new_key = ad_compatibility[offset_i];
2742 std::get<0>(keys[i]) = new_key;
2745 std::sort(keys.begin(), keys.end());
2747 offsets_0.reserve(n_unique_entities + 1);
2748 columns_0.reserve(n_unique_entity_vertices);
2752 std::array<unsigned int, max_n_vertices> ref_key;
2753 std::array<unsigned int, max_n_vertices> ref_indices;
2754 std::fill(ref_key.begin(), ref_key.end(), 0);
2757 for (
unsigned int i = 0; i < keys.size(); i++)
2759 const auto offset_i = std::get<1>(keys[i]);
2761 if (ref_key != std::get<0>(keys[i]))
2765 ref_key = std::get<0>(keys[i]);
2766 ref_indices = ad_entity_vertices[offset_i];
2768 offsets_0.push_back(columns_0.size());
2769 for (
const auto j : ad_entity_vertices[offset_i])
2771 columns_0.push_back(j - offset);
2777 orientations[offset_i] =
2778 ad_entity_types[offset_i]
2779 .template get_combined_orientation<unsigned int>(
2781 ad_entity_vertices[offset_i].
begin() +
2782 ad_entity_types[offset_i].n_vertices()),
2784 ref_indices.begin() +
2785 ad_entity_types[offset_i].n_vertices()));
2787 cells_to_faces_indices[offset_i] = counter;
2789 offsets_0.push_back(columns_0.size());
2792 crs_0 =
ArrayOfArrays(std::move(offsets_0), std::move(columns_0));
2793 cells_to_faces =
ArrayOfArrays(std::move(cells_to_faces_offsets),
2794 std::move(cells_to_faces_indices));
2803 template <
int face_dim,
int dim,
typename FU>
2809 std::vector<types::geometric_orientation> &orientations,
2810 const FU &second_key_function)
2812 unsigned int max_n_vertices = 0;
2817 for (
const auto &c : cell_types)
2818 if (face_dim == dim - 1)
2820 for (
unsigned int f = 0; f < c.n_faces(); ++f)
2822 std::max(max_n_vertices, c.face_reference_cell(f).n_vertices());
2824 else if (face_dim == 1)
2825 max_n_vertices =
std::max(max_n_vertices, 2u);
2829 if (max_n_vertices == 2)
2835 second_key_function);
2836 else if (max_n_vertices == 3)
2842 second_key_function);
2843 else if (max_n_vertices == 4)
2849 second_key_function);
2872 const std::vector<types::geometric_orientation> &ori_cq,
2874 std::vector<types::geometric_orientation> &ori_ql,
2878 quad_t_id.resize(quads_to_vertices.
size());
2887 std::vector<std::size_t> quads_to_lines_offsets(quads_to_vertices.
size() +
2889 for (
unsigned int c = 0; c < cells_to_quads.
size(); ++c)
2893 for (
unsigned int f = 0; f < faces.
size(); ++f)
2895 quads_to_lines_offsets[faces[f] + 1] =
2896 cell_types[c].face_reference_cell(f).n_lines();
2899 quads_to_lines_offsets[0] = 0;
2900 for (
unsigned int i = 0; i < quads_to_lines_offsets.size() - 1; ++i)
2901 quads_to_lines_offsets[i + 1] += quads_to_lines_offsets[i];
2904 std::vector<unsigned int> quads_to_lines_indices(
2905 quads_to_lines_offsets.back());
2906 ori_ql.assign(quads_to_lines_offsets.back(),
2910 unsigned int global_face_index = 0;
2911 for (
unsigned int c = 0; c < cells_to_quads.
size(); ++c)
2917 for (
unsigned int f_index = 0; f_index < faces.
size();
2918 ++f_index, ++global_face_index)
2920 const unsigned int f = faces[f_index];
2923 if (ori_cq[global_face_index] !=
2929 for (
unsigned int l = 0; l < quad_t_id[f].n_lines(); ++l)
2932 const unsigned int local_line_index =
2935 const unsigned int global_line_index =
2936 cells_to_lines[c][local_line_index];
2937 quads_to_lines_indices[quads_to_lines_offsets[f] + l] =
2942 for (
unsigned int v = 0; v < 2; ++v)
2943 if (cells_to_vertices
2949 lines_to_vertices[global_line_index][v])
2956 ori_ql[quads_to_lines_offsets[f] + l] =
2964 quads_to_lines =
ArrayOfArrays(std::move(quads_to_lines_offsets),
2965 std::move(quads_to_lines_indices));
2988 if constexpr (dim == 1)
2991 if constexpr (dim == 2 || dim == 3)
2993 std::vector<types::geometric_orientation> dummy;
2996 connectivity.template entity_types<dim>(),
3001 [](
auto key,
const auto &,
const auto &,
const auto &) {
3007 if constexpr (dim == 3)
3010 connectivity.template entity_types<3>(),
3018 const unsigned int &c,
3019 const unsigned int &f) {
3021 AssertIndexRange(cell_type.face_reference_cell(f).n_lines(),
3026 for (; l < cell_type.face_reference_cell(f).n_lines(); ++l)
3028 AssertIndexRange(l, key.size());
3029 key[l] = temp1[c][cell_type.face_to_cell_lines(
3030 f, l, numbers::default_geometric_orientation)] +
3034 for (; l < key.size(); ++l)
3050 connectivity.template entity_types<2>());
3054 connectivity.entity_to_entities(dim, dim) =
3057 return connectivity;
3072 std::vector<unsigned int> cell_vertices;
3073 cell_vertices.reserve(std::accumulate(cells.begin(),
3076 [](
const unsigned int accumulator,
3078 return accumulator +
3079 cell.vertices.size();
3082 std::vector<std::size_t> cell_vertices_offsets;
3083 cell_vertices_offsets.reserve(cells.size() + 1);
3084 cell_vertices_offsets.push_back(0);
3086 std::vector<ReferenceCell<dim>> cell_types;
3087 cell_types.reserve(cells.size());
3091 std::set<ReferenceCell<dim>> all_reference_cells;
3093 for (
const auto &cell : cells)
3097 auto vertices_unique = cell.
vertices;
3098 std::sort(vertices_unique.begin(), vertices_unique.end());
3099 vertices_unique.erase(std::unique(vertices_unique.begin(),
3100 vertices_unique.end()),
3101 vertices_unique.end());
3106 "The definition of a cell refers to the same vertex several "
3107 "times. This is not possible. A common reason is that "
3108 "CellData::vertices has a size that does not match the "
3109 "size expected from the reference cell. Please resize "
3110 "CellData::vertices or use the appropriate constructor of "
3114 cell_types.push_back(
3117 all_reference_cells.insert(cell_types.back());
3120 cell_vertices.insert(cell_vertices.end(),
3123 cell_vertices_offsets.push_back(cell_vertices.size());
3128 all_reference_cells,
3129 {std::move(cell_vertices_offsets),
3130 std::move(cell_vertices)});
3242 template <
int dim,
int spacedim>
3246 const unsigned int level_objects,
3249 using line_iterator =
3252 number_cache.n_levels = 0;
3253 if (level_objects > 0)
3255 for (
unsigned int level = 0; level < level_objects; ++level)
3256 if (triangulation.
begin(level) != triangulation.
end(level))
3257 number_cache.n_levels = level + 1;
3265 number_cache.n_lines = 0;
3266 number_cache.n_active_lines = 0;
3272 number_cache.n_lines_level.resize(number_cache.n_levels);
3273 number_cache.n_active_lines_level.resize(number_cache.n_levels);
3275 for (
unsigned int level = 0; level < number_cache.n_levels; ++level)
3278 number_cache.n_lines_level[level] = 0;
3279 number_cache.n_active_lines_level[level] = 0;
3281 line_iterator line = triangulation.
begin_line(level),
3283 (level == number_cache.n_levels - 1 ?
3284 line_iterator(triangulation.
end_line()) :
3286 for (; line != endc; ++line)
3288 ++number_cache.n_lines_level[level];
3289 if (line->has_children() ==
false)
3290 ++number_cache.n_active_lines_level[level];
3294 number_cache.n_lines += number_cache.n_lines_level[level];
3295 number_cache.n_active_lines +=
3296 number_cache.n_active_lines_level[level];
3302 number_cache.n_lines_level.clear();
3303 number_cache.n_active_lines_level.clear();
3305 line_iterator line = triangulation.
begin_line(),
3307 for (; line != endc; ++line)
3309 ++number_cache.n_lines;
3310 if (line->has_children() ==
false)
3311 ++number_cache.n_active_lines;
3330 template <
int dim,
int spacedim>
3334 const unsigned int level_objects,
3351 using quad_iterator =
3357 number_cache.n_quads = 0;
3358 number_cache.n_active_quads = 0;
3368 unsigned int n_levels = 0;
3369 if (level_objects > 0)
3371 for (
unsigned int level = 0; level < level_objects; ++level)
3372 if (triangulation.
begin(level) != triangulation.
end(level))
3373 n_levels = level + 1;
3375 number_cache.n_quads_level.resize(n_levels);
3376 number_cache.n_active_quads_level.resize(n_levels);
3378 for (
unsigned int level = 0; level < n_levels; ++level)
3381 number_cache.n_quads_level[level] = 0;
3382 number_cache.n_active_quads_level[level] = 0;
3384 quad_iterator quad = triangulation.
begin_quad(level),
3386 (level == n_levels - 1 ?
3387 quad_iterator(triangulation.
end_quad()) :
3389 for (; quad != endc; ++quad)
3391 ++number_cache.n_quads_level[level];
3392 if (quad->has_children() ==
false)
3393 ++number_cache.n_active_quads_level[level];
3397 number_cache.n_quads += number_cache.n_quads_level[level];
3398 number_cache.n_active_quads +=
3399 number_cache.n_active_quads_level[level];
3405 number_cache.n_quads_level.clear();
3406 number_cache.n_active_quads_level.clear();
3408 quad_iterator quad = triangulation.
begin_quad(),
3410 for (; quad != endc; ++quad)
3412 ++number_cache.n_quads;
3413 if (quad->has_children() ==
false)
3414 ++number_cache.n_active_quads;
3419 update_lines.
join();
3437 template <
int dim,
int spacedim>
3441 const unsigned int level_objects,
3458 using hex_iterator =
3464 number_cache.n_hexes = 0;
3465 number_cache.n_active_hexes = 0;
3476 unsigned int n_levels = 0;
3477 if (level_objects > 0)
3479 for (
unsigned int level = 0; level < level_objects; ++level)
3480 if (triangulation.
begin(level) != triangulation.
end(level))
3481 n_levels = level + 1;
3483 number_cache.n_hexes_level.resize(n_levels);
3484 number_cache.n_active_hexes_level.resize(n_levels);
3486 for (
unsigned int level = 0; level < n_levels; ++level)
3489 number_cache.n_hexes_level[level] = 0;
3490 number_cache.n_active_hexes_level[level] = 0;
3492 hex_iterator hex = triangulation.
begin_hex(level),
3493 endc = (level == n_levels - 1 ?
3494 hex_iterator(triangulation.
end_hex()) :
3496 for (; hex != endc; ++hex)
3498 ++number_cache.n_hexes_level[level];
3499 if (hex->has_children() ==
false)
3500 ++number_cache.n_active_hexes_level[level];
3504 number_cache.n_hexes += number_cache.n_hexes_level[level];
3505 number_cache.n_active_hexes +=
3506 number_cache.n_active_hexes_level[level];
3512 number_cache.n_hexes_level.clear();
3513 number_cache.n_active_hexes_level.clear();
3515 hex_iterator hex = triangulation.
begin_hex(),
3516 endc = triangulation.
end_hex();
3517 for (; hex != endc; ++hex)
3519 ++number_cache.n_hexes;
3520 if (hex->has_children() ==
false)
3521 ++number_cache.n_active_hexes;
3526 update_quads_and_lines.
join();
3530 template <
int dim,
int spacedim>
3534 const unsigned int level_objects,
3539 number_cache.active_cell_index_partitioner =
3540 std::make_shared<const Utilities::MPI::Partitioner>(
3543 number_cache.level_cell_index_partitioners.resize(
3545 for (
unsigned int level = 0; level < triangulation.
n_levels(); ++level)
3546 number_cache.level_cell_index_partitioners[level] =
3547 std::make_shared<const Utilities::MPI::Partitioner>(
3548 triangulation.
n_cells(level));
3552 template <
int spacedim>
3558 template <
int dim,
int spacedim>
3591 constexpr unsigned int left_right_offset[2][6][2] = {
3600 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
3608 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
3609 adjacent_cells(2 * triangulation.
n_raw_faces(), dummy);
3617 const unsigned int offset =
3623 adjacent_cells[2 * face->index() + offset] = cell;
3631 if (cell->
is_active() && face->has_children())
3633 adjacent_cells[2 * face->child(0)->index() + offset] =
3635 adjacent_cells[2 * face->child(1)->index() + offset] =
3653 if (face->has_children() &&
3659 for (
unsigned int c = 0; c < face->n_children(); ++c)
3660 adjacent_cells[2 * face->child(c)->index() + offset] =
3662 if (face->child(0)->has_children())
3664 adjacent_cells[2 * face->child(0)->child(0)->index() +
3666 adjacent_cells[2 * face->child(0)->child(1)->index() +
3669 if (face->child(1)->has_children())
3671 adjacent_cells[2 * face->child(1)->child(0)->index() +
3673 adjacent_cells[2 * face->child(1)->child(1)->index() +
3687 const unsigned int offset =
3693 f, adjacent_cells[2 * cell->
face(f)->index() + 1 - offset]);
3701 template <
int dim,
int spacedim>
3715 if (dim == spacedim)
3716 for (
unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
3738 const auto n_cells = cells.size();
3741 const auto all_reference_cells =
3745 all_reference_cells.end());
3754 tria.
strides.max_children_per_cell,
3755 tria.
strides.max_faces_per_cell,
3756 tria.
strides.max_vertices_per_cell));
3759 tria.
faces = std::make_unique<
3761 tria.
strides.max_children_per_face,
3762 tria.
strides.max_lines_per_face);
3766 dim > 1 ? connectivity.entity_to_entities(1, 0) : empty;
3768 dim == 3 ? connectivity.entity_to_entities(2, 1) : empty;
3769 const auto n_lines = lines_to_vertices.
size();
3771 n_lines <= std::size_t(std::numeric_limits<int>::max()),
3772 ExcMessage(
"The number of lines must be less than the largest int."));
3773 const auto n_quads = quads_to_lines.
size();
3775 n_quads <= std::size_t(std::numeric_limits<int>::max()),
3776 ExcMessage(
"The number of quads must be less than the largest int."));
3779 tria.
faces->allocate(n_lines, n_quads);
3784 auto &faces = *tria.
faces;
3786 for (
int line = 0; line <
int(n_lines); ++line)
3788 const auto vertices = lines_to_vertices[line];
3789 for (
unsigned int v = 0; v < vertices.size(); ++v)
3790 faces.lines.set_bounding_object(line, v, vertices[v]);
3795 if constexpr (dim == 3)
3797 auto &faces = *tria.
faces;
3801 for (
int q = 0, k = 0; q <
int(n_quads); ++q)
3804 const auto face_reference_cell =
3805 connectivity.template entity_types<2>()[q];
3806 faces.set_quad_type(q, face_reference_cell);
3809 const auto lines = quads_to_lines[q];
3810 for (
unsigned int l = 0; l < lines.size(); ++l, ++k)
3814 faces.quads.set_bounding_object(q, l, lines[l]);
3817 faces.set_line_orientation(
3818 q, l, connectivity.entity_orientations(1)[k]);
3825 auto &cells_0 = tria.
levels[0]->cells;
3826 auto &level = *tria.
levels[0];
3829 const auto &cells_to_faces =
3830 connectivity.entity_to_entities(dim, dim - 1);
3831 const auto &cells_to_neighbor_cells =
3832 connectivity.entity_to_entities(dim, dim);
3836 bool orientation_needed =
false;
3838 orientation_needed =
true;
3841 const auto &orientations = connectivity.entity_orientations(1);
3842 orientation_needed = std::any_of(
3843 orientations.begin(), orientations.end(), [](
const auto &a) {
3844 return a != numbers::default_geometric_orientation;
3848 level.allocate(
n_cells, orientation_needed);
3850 unsigned int global_face_index = 0;
3851 for (
unsigned int cell = 0; cell <
n_cells; ++cell)
3854 cells_0.boundary_or_material_id[cell].material_id =
3855 cells[cell].material_id;
3858 cells_0.manifold_id[cell] = cells[cell].manifold_id;
3861 level.reference_cell[cell] =
3862 connectivity.template entity_types<dim>()[cell];
3863 Assert(all_reference_cells.count(level.reference_cell[cell]) ==
3868 const auto faces = cells_to_faces[cell];
3869 const auto neighbors = cells_to_neighbor_cells[cell];
3870 for (
unsigned int f = 0; f < faces.size();
3871 ++f, ++global_face_index)
3875 level.set_neighbor(cell, f, 0, neighbors[f]);
3878 cells_0.set_bounding_object(cell, f, faces[f]);
3881 if (orientation_needed)
3883 level.face_orientations.set_combined_orientation(
3886 connectivity.entity_orientations(dim -
3887 1)[global_face_index]);
3896 auto &bids_face = dim == 3 ?
3897 tria.
faces->quads.boundary_or_material_id :
3898 tria.
faces->lines.boundary_or_material_id;
3901 std::vector<unsigned int> count(bids_face.size(), 0);
3904 const auto &cells_to_faces =
3905 connectivity.entity_to_entities(dim, dim - 1);
3908 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
3909 for (
const unsigned int face : cells_to_faces[cell])
3913 for (
unsigned int face = 0; face < count.size(); ++face)
3915 if (count[face] != 1)
3919 bids_face[face].boundary_id = 0;
3925 const auto &quads_to_lines =
3926 connectivity.entity_to_entities(2, 1);
3927 for (
const unsigned int line : quads_to_lines[face])
3928 tria.
faces->lines.boundary_or_material_id[line].boundary_id =
3934 constexpr unsigned int t_tba =
static_cast<unsigned int>(-1);
3935 constexpr unsigned int t_inner =
static_cast<unsigned int>(-2);
3937 std::vector<unsigned int> type(vertices.size(), t_tba);
3939 const auto &cells_to_vertices =
3940 connectivity.entity_to_entities(1, 0);
3942 for (
unsigned int cell = 0; cell < cells.size(); ++cell)
3944 const auto vertices = cells_to_vertices[cell];
3945 for (
unsigned int v = 0; v < vertices.size(); ++v)
3946 if (type[vertices[v]] != t_inner)
3948 (type[vertices[v]] == t_tba ? v : t_inner);
3951 for (
unsigned int face = 0; face < type.size(); ++face)
3956 if (type[face] != t_inner && type[face] != t_tba)
3977 template <
int structdim,
int spacedim>
3987 if (boundary_objects_in.empty())
3991 auto boundary_objects = boundary_objects_in;
3994 for (
auto &boundary_object : boundary_objects)
3995 std::sort(boundary_object.vertices.begin(),
3996 boundary_object.vertices.end());
3999 std::sort(boundary_objects.begin(),
4000 boundary_objects.end(),
4001 [](
const auto &a,
const auto &b) {
4002 return a.vertices < b.vertices;
4005 [[maybe_unused]]
unsigned int counter = 0;
4011 for (
unsigned int o = 0; o < obj.
n_objects(); ++o)
4026 key.assign(indices.
begin(), indices.
end());
4027 std::sort(key.begin(), key.end());
4030 const auto subcell_object =
4031 std::lower_bound(boundary_objects.begin(),
4032 boundary_objects.end(),
4034 [&](
const auto &cell,
const auto &key) {
4035 return cell.vertices < key;
4039 if (subcell_object == boundary_objects.end() ||
4040 subcell_object->vertices != key)
4046 manifold_id = subcell_object->manifold_id;
4049 if (subcell_object->boundary_id !=
4055 "The input arguments for creating a triangulation "
4056 "specified a boundary id for an internal face. This "
4059 "The object in question has vertex indices " +
4060 [subcell_object]() {
4062 for (
const auto v : subcell_object->vertices)
4063 s += std::to_string(v) +
',';
4066 " which are located at coordinates " +
4067 [vertex_locations, subcell_object]() {
4068 std::ostringstream s;
4069 for (
unsigned int i = 0;
4070 i < subcell_object->vertices.size();
4073 << vertex_locations[subcell_object->vertices[i]]
4074 << (i != subcell_object->vertices.size() - 1 ?
"), " :
4079 boundary_id = subcell_object->boundary_id;
4105 template <
int spacedim>
4109 std::vector<unsigned int> &,
4110 std::vector<unsigned int> &)
4112 const unsigned int dim = 1;
4138 neighbor = neighbor->
child(1);
4150 neighbor = neighbor->
child(1);
4166 neighbor = neighbor->
child(0);
4178 neighbor = neighbor->
child(0);
4194 for (
unsigned int child = 0; child < cell->
n_children(); ++child)
4209 template <
int spacedim>
4213 std::vector<unsigned int> &line_cell_count,
4214 std::vector<unsigned int> &)
4216 const unsigned int dim = 2;
4224 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
4227 lines_to_delete.reserve(4 * 2 + 4);
4232 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
4236 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
4254 lines_to_delete.push_back(cell->
child(0)->
line(1));
4255 lines_to_delete.push_back(cell->
child(0)->
line(3));
4256 lines_to_delete.push_back(cell->
child(3)->
line(0));
4257 lines_to_delete.push_back(cell->
child(3)->
line(2));
4261 unsigned int inner_face_no =
4266 lines_to_delete.push_back(cell->
child(0)->
line(inner_face_no));
4270 for (
unsigned int child = 0; child < cell->
n_children(); ++child)
4287 for (
unsigned int line_no = 0;
4288 line_no < GeometryInfo<dim>::lines_per_cell;
4292 cell->
line(line_no);
4294 if (line->has_children())
4299 Assert((line_cell_count[line->child_index(0)] == 0 &&
4300 line_cell_count[line->child_index(1)] == 0) ||
4301 (line_cell_count[line->child_index(0)] > 0 &&
4302 line_cell_count[line->child_index(1)] > 0),
4305 if (line_cell_count[line->child_index(0)] == 0)
4307 for (
unsigned int c = 0; c < 2; ++c)
4308 Assert(!line->child(c)->has_children(),
4318 lines_to_delete.push_back(line->child(0));
4319 lines_to_delete.push_back(line->child(1));
4321 line->clear_children();
4333 typename std::vector<
4335 line = lines_to_delete.begin(),
4336 endline = lines_to_delete.end();
4337 for (; line != endline; ++line)
4339 (*line)->clear_user_data();
4340 (*line)->clear_user_flag();
4341 (*line)->clear_used_flag();
4347 template <
int spacedim>
4351 std::vector<unsigned int> &line_cell_count,
4352 std::vector<unsigned int> &quad_cell_count)
4354 const unsigned int dim = 3;
4366 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
4368 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
4371 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
4372 quads_to_delete.reserve(6 * 4 + 12);
4376 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
4380 const auto line_indices = TriaAccessorImplementation::
4381 Implementation::get_line_indices_of_cell(*child);
4383 --line_cell_count[line_indices[l]];
4399 quads_to_delete.push_back(cell->
child(0)->
face(1));
4402 quads_to_delete.push_back(cell->
child(0)->
face(3));
4405 quads_to_delete.push_back(cell->
child(0)->
face(5));
4408 quads_to_delete.push_back(cell->
child(0)->
face(1));
4409 quads_to_delete.push_back(cell->
child(0)->
face(3));
4410 quads_to_delete.push_back(cell->
child(3)->
face(0));
4411 quads_to_delete.push_back(cell->
child(3)->
face(2));
4413 lines_to_delete.push_back(cell->
child(0)->
line(11));
4416 quads_to_delete.push_back(cell->
child(0)->
face(1));
4417 quads_to_delete.push_back(cell->
child(0)->
face(5));
4418 quads_to_delete.push_back(cell->
child(3)->
face(0));
4419 quads_to_delete.push_back(cell->
child(3)->
face(4));
4421 lines_to_delete.push_back(cell->
child(0)->
line(5));
4424 quads_to_delete.push_back(cell->
child(0)->
face(3));
4425 quads_to_delete.push_back(cell->
child(0)->
face(5));
4426 quads_to_delete.push_back(cell->
child(3)->
face(2));
4427 quads_to_delete.push_back(cell->
child(3)->
face(4));
4429 lines_to_delete.push_back(cell->
child(0)->
line(7));
4432 quads_to_delete.push_back(cell->
child(0)->
face(1));
4433 quads_to_delete.push_back(cell->
child(2)->
face(1));
4434 quads_to_delete.push_back(cell->
child(4)->
face(1));
4435 quads_to_delete.push_back(cell->
child(6)->
face(1));
4437 quads_to_delete.push_back(cell->
child(0)->
face(3));
4438 quads_to_delete.push_back(cell->
child(1)->
face(3));
4439 quads_to_delete.push_back(cell->
child(4)->
face(3));
4440 quads_to_delete.push_back(cell->
child(5)->
face(3));
4442 quads_to_delete.push_back(cell->
child(0)->
face(5));
4443 quads_to_delete.push_back(cell->
child(1)->
face(5));
4444 quads_to_delete.push_back(cell->
child(2)->
face(5));
4445 quads_to_delete.push_back(cell->
child(3)->
face(5));
4447 lines_to_delete.push_back(cell->
child(0)->
line(5));
4448 lines_to_delete.push_back(cell->
child(0)->
line(7));
4449 lines_to_delete.push_back(cell->
child(0)->
line(11));
4450 lines_to_delete.push_back(cell->
child(7)->
line(0));
4451 lines_to_delete.push_back(cell->
child(7)->
line(2));
4452 lines_to_delete.push_back(cell->
child(7)->
line(8));
4470 for (
unsigned int child = 0; child < cell->
n_children(); ++child)
4501 cell->
face(quad_no);
4505 quad->has_children()) ||
4510 switch (quad->refinement_case())
4522 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
4523 quad_cell_count[quad->child_index(1)] == 0) ||
4524 (quad_cell_count[quad->child_index(0)] > 0 &&
4525 quad_cell_count[quad->child_index(1)] > 0),
4531 unsigned int deleted_grandchildren = 0;
4532 unsigned int number_of_child_refinements = 0;
4534 for (
unsigned int c = 0; c < 2; ++c)
4535 if (quad->child(c)->has_children())
4537 ++number_of_child_refinements;
4542 (quad_cell_count[quad->child(c)->child_index(0)] ==
4544 quad_cell_count[quad->child(c)->child_index(1)] ==
4546 (quad_cell_count[quad->child(c)->child_index(0)] >
4548 quad_cell_count[quad->child(c)->child_index(1)] >
4551 if (quad_cell_count[quad->child(c)->child_index(0)] ==
4558 Assert(quad->refinement_case() +
4559 quad->child(c)->refinement_case() ==
4567 quads_to_delete.push_back(
4568 quad->child(c)->child(0));
4569 quads_to_delete.push_back(
4570 quad->child(c)->child(1));
4571 if (quad->child(c)->refinement_case() ==
4573 lines_to_delete.push_back(
4574 quad->child(c)->child(0)->line(1));
4576 lines_to_delete.push_back(
4577 quad->child(c)->child(0)->line(3));
4578 quad->child(c)->clear_children();
4579 quad->child(c)->clear_refinement_case();
4580 ++deleted_grandchildren;
4588 if (number_of_child_refinements > 0 &&
4589 deleted_grandchildren == number_of_child_refinements)
4594 middle_line = quad->child(0)->line(1);
4596 middle_line = quad->child(0)->line(3);
4598 lines_to_delete.push_back(middle_line->child(0));
4599 lines_to_delete.push_back(middle_line->child(1));
4602 middle_line)] =
false;
4603 middle_line->clear_children();
4608 if (quad_cell_count[quad->child_index(0)] == 0)
4614 quads_to_delete.push_back(quad->child(0));
4615 quads_to_delete.push_back(quad->child(1));
4617 lines_to_delete.push_back(quad->child(0)->line(1));
4619 lines_to_delete.push_back(quad->child(0)->line(3));
4639 if (quad->child(0)->has_children())
4640 if (quad->refinement_case() ==
4674 quad_iterator switch_1 =
4675 quad->child(0)->child(1),
4677 quad->child(1)->child(0);
4679 Assert(!switch_1->has_children(),
4681 Assert(!switch_2->has_children(),
4684 const int switch_1_index = switch_1->index();
4685 const int switch_2_index = switch_2->index();
4686 for (
unsigned int l = 0;
4687 l < triangulation.
levels.size();
4689 for (
unsigned int h = 0;
4691 triangulation.
levels[l]->cells.n_objects();
4693 for (
const unsigned int q :
4698 ->cells.get_bounding_object_indices(
4700 if (
index == switch_1_index)
4702 ->cells.get_bounding_object_indices(
4703 h)[q] = switch_2_index;
4704 else if (
index == switch_2_index)
4706 ->cells.get_bounding_object_indices(
4707 h)[q] = switch_1_index;
4712 const int switch_1_lines[4] = {
4713 static_cast<signed int>(
4714 switch_1->line_index(0)),
4715 static_cast<signed int>(
4716 switch_1->line_index(1)),
4717 static_cast<signed int>(
4718 switch_1->line_index(2)),
4719 static_cast<signed int>(
4720 switch_1->line_index(3))};
4722 switch_1_line_orientations[4] = {
4723 switch_1->line_orientation(0),
4724 switch_1->line_orientation(1),
4725 switch_1->line_orientation(2),
4726 switch_1->line_orientation(3)};
4728 switch_1->boundary_id();
4729 const unsigned int switch_1_user_index =
4730 switch_1->user_index();
4731 const bool switch_1_user_flag =
4732 switch_1->user_flag_set();
4734 switch_1->set_bounding_object_indices(
4735 {switch_2->line_index(0),
4736 switch_2->line_index(1),
4737 switch_2->line_index(2),
4738 switch_2->line_index(3)});
4739 switch_1->set_line_orientation(
4740 0, switch_2->line_orientation(0));
4741 switch_1->set_line_orientation(
4742 1, switch_2->line_orientation(1));
4743 switch_1->set_line_orientation(
4744 2, switch_2->line_orientation(2));
4745 switch_1->set_line_orientation(
4746 3, switch_2->line_orientation(3));
4747 switch_1->set_boundary_id_internal(
4748 switch_2->boundary_id());
4749 switch_1->set_manifold_id(
4750 switch_2->manifold_id());
4751 switch_1->set_user_index(switch_2->user_index());
4752 if (switch_2->user_flag_set())
4753 switch_1->set_user_flag();
4755 switch_1->clear_user_flag();
4757 switch_2->set_bounding_object_indices(
4761 switch_1_lines[3]});
4762 switch_2->set_line_orientation(
4763 0, switch_1_line_orientations[0]);
4764 switch_2->set_line_orientation(
4765 1, switch_1_line_orientations[1]);
4766 switch_2->set_line_orientation(
4767 2, switch_1_line_orientations[2]);
4768 switch_2->set_line_orientation(
4769 3, switch_1_line_orientations[3]);
4770 switch_2->set_boundary_id_internal(
4771 switch_1_boundary_id);
4772 switch_2->set_manifold_id(
4773 switch_1->manifold_id());
4774 switch_2->set_user_index(switch_1_user_index);
4775 if (switch_1_user_flag)
4776 switch_2->set_user_flag();
4778 switch_2->clear_user_flag();
4780 const unsigned int child_0 =
4781 quad->child(0)->child_index(0);
4782 const unsigned int child_2 =
4783 quad->child(1)->child_index(0);
4784 quad->clear_children();
4785 quad->clear_refinement_case();
4786 quad->set_refinement_case(
4788 quad->set_children(0, child_0);
4789 quad->set_children(2, child_2);
4790 std::swap(quad_cell_count[child_0 + 1],
4791 quad_cell_count[child_2]);
4806 const unsigned int child_0 =
4807 quad->child(0)->child_index(0);
4808 const unsigned int child_2 =
4809 quad->child(1)->child_index(0);
4810 quad->clear_children();
4811 quad->clear_refinement_case();
4812 quad->set_refinement_case(
4814 quad->set_children(0, child_0);
4815 quad->set_children(2, child_2);
4819 quad->clear_children();
4820 quad->clear_refinement_case();
4831 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
4832 quad_cell_count[quad->child_index(1)] == 0 &&
4833 quad_cell_count[quad->child_index(2)] == 0 &&
4834 quad_cell_count[quad->child_index(3)] == 0) ||
4835 (quad_cell_count[quad->child_index(0)] > 0 &&
4836 quad_cell_count[quad->child_index(1)] > 0 &&
4837 quad_cell_count[quad->child_index(2)] > 0 &&
4838 quad_cell_count[quad->child_index(3)] > 0),
4841 if (quad_cell_count[quad->child_index(0)] == 0)
4847 lines_to_delete.push_back(quad->child(0)->line(1));
4848 lines_to_delete.push_back(quad->child(3)->line(0));
4849 lines_to_delete.push_back(quad->child(0)->line(3));
4850 lines_to_delete.push_back(quad->child(3)->line(2));
4852 for (
unsigned int child = 0; child < quad->n_children();
4854 quads_to_delete.push_back(quad->child(child));
4860 quad->clear_children();
4861 quad->clear_refinement_case();
4880 for (
unsigned int line_no = 0;
4881 line_no < GeometryInfo<dim>::lines_per_cell;
4885 cell->
line(line_no);
4889 line->has_children()) ||
4894 if (line->has_children())
4899 Assert((line_cell_count[line->child_index(0)] == 0 &&
4900 line_cell_count[line->child_index(1)] == 0) ||
4901 (line_cell_count[line->child_index(0)] > 0 &&
4902 line_cell_count[line->child_index(1)] > 0),
4905 if (line_cell_count[line->child_index(0)] == 0)
4907 for (
unsigned int c = 0; c < 2; ++c)
4908 Assert(!line->child(c)->has_children(),
4918 lines_to_delete.push_back(line->child(0));
4919 lines_to_delete.push_back(line->child(1));
4921 line->clear_children();
4933 typename std::vector<
4935 line = lines_to_delete.begin(),
4936 endline = lines_to_delete.end();
4937 for (; line != endline; ++line)
4939 (*line)->clear_user_data();
4940 (*line)->clear_user_flag();
4941 (*line)->clear_used_flag();
4944 typename std::vector<
4946 quad = quads_to_delete.begin(),
4947 endquad = quads_to_delete.end();
4948 for (; quad != endquad; ++quad)
4950 (*quad)->clear_user_data();
4951 (*quad)->clear_children();
4952 (*quad)->clear_refinement_case();
4953 (*quad)->clear_user_flag();
4954 (*quad)->clear_used_flag();
4976 template <
int spacedim>
4980 unsigned int &next_unused_vertex,
4987 const unsigned int dim = 2;
5083 int new_vertices[9];
5084 for (
unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
5085 new_vertices[vertex_no] = cell->
vertex_index(vertex_no);
5086 for (
unsigned int line_no = 0; line_no < 4; ++line_no)
5087 if (cell->
line(line_no)->has_children())
5088 new_vertices[4 + line_no] =
5089 cell->
line(line_no)->child(0)->vertex_index(1);
5098 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5099 ++next_unused_vertex;
5102 "Internal error: During refinement, the triangulation "
5103 "wants to access an element of the 'vertices' array "
5104 "but it turns out that the array is not large enough."));
5107 new_vertices[8] = next_unused_vertex;
5112 triangulation.
vertices[next_unused_vertex] =
5113 cell->
center(
true,
true);
5119 unsigned int lmin = 8;
5120 unsigned int lmax = 12;
5127 for (
unsigned int l = lmin; l < lmax; ++l)
5129 while (next_unused_line->used() ==
true)
5131 new_lines[l] = next_unused_line;
5149 for (
unsigned int c = 0; c < 2; ++c, ++l)
5150 new_lines[l] = cell->
line(face_no)->child(c);
5153 new_lines[8]->set_bounding_object_indices(
5154 {new_vertices[6], new_vertices[8]});
5155 new_lines[9]->set_bounding_object_indices(
5156 {new_vertices[8], new_vertices[7]});
5157 new_lines[10]->set_bounding_object_indices(
5158 {new_vertices[4], new_vertices[8]});
5159 new_lines[11]->set_bounding_object_indices(
5160 {new_vertices[8], new_vertices[5]});
5169 new_lines[0] = cell->
line(0);
5170 new_lines[1] = cell->
line(1);
5171 new_lines[2] = cell->
line(2)->child(0);
5172 new_lines[3] = cell->
line(2)->child(1);
5173 new_lines[4] = cell->
line(3)->child(0);
5174 new_lines[5] = cell->
line(3)->child(1);
5175 new_lines[6]->set_bounding_object_indices(
5176 {new_vertices[6], new_vertices[7]});
5186 new_lines[0] = cell->
line(0)->child(0);
5187 new_lines[1] = cell->
line(0)->child(1);
5188 new_lines[2] = cell->
line(1)->child(0);
5189 new_lines[3] = cell->
line(1)->child(1);
5190 new_lines[4] = cell->
line(2);
5191 new_lines[5] = cell->
line(3);
5192 new_lines[6]->set_bounding_object_indices(
5193 {new_vertices[4], new_vertices[5]});
5196 for (
unsigned int l = lmin; l < lmax; ++l)
5198 new_lines[l]->set_used_flag();
5199 new_lines[l]->clear_user_flag();
5201 new_lines[l]->clear_children();
5203 new_lines[l]->set_boundary_id_internal(
5205 new_lines[l]->set_manifold_id(cell->
manifold_id());
5212 while (next_unused_cell->
used() ==
true)
5216 for (
unsigned int i = 0; i < n_children; ++i)
5219 subcells[i] = next_unused_cell;
5221 if (i % 2 == 1 && i < n_children - 1)
5222 while (next_unused_cell->
used() ==
true)
5241 new_lines[8]->index(),
5242 new_lines[4]->index(),
5243 new_lines[10]->index()});
5245 new_lines[2]->index(),
5246 new_lines[5]->index(),
5247 new_lines[11]->index()});
5249 new_lines[9]->index(),
5250 new_lines[10]->index(),
5251 new_lines[6]->index()});
5253 new_lines[3]->index(),
5254 new_lines[11]->index(),
5255 new_lines[7]->index()});
5272 new_lines[6]->index(),
5273 new_lines[2]->index(),
5274 new_lines[4]->index()});
5276 new_lines[1]->index(),
5277 new_lines[3]->index(),
5278 new_lines[5]->index()});
5296 new_lines[2]->index(),
5297 new_lines[4]->index(),
5298 new_lines[6]->index()});
5300 new_lines[3]->index(),
5301 new_lines[6]->index(),
5302 new_lines[5]->index()});
5307 for (
unsigned int i = 0; i < n_children; ++i)
5331 for (
unsigned int i = 0; i < n_children / 2; ++i)
5341 if (dim == spacedim - 1)
5342 for (
unsigned int c = 0; c < n_children; ++c)
5348 template <
int spacedim>
5359 template <
int spacedim>
5362 const bool check_for_distorted_cells)
5364 constexpr int dim = 2;
5369 triangulation.
levels.size() - 1))
5372 triangulation.
levels.push_back(
5375 triangulation.
strides.max_children_per_cell,
5376 triangulation.
strides.max_faces_per_cell,
5377 triangulation.
strides.max_vertices_per_cell));
5386 line->clear_user_flag();
5390 unsigned int n_single_lines = 0;
5391 unsigned int n_lines_in_pairs = 0;
5392 unsigned int needed_vertices = 0;
5394 const bool orientation_needed =
5397 for (
int level_no = triangulation.
levels.size() - 2; level_no >= 0;
5402 std::size_t needed_cells = 0;
5404 for (
const auto &cell :
5411 needed_vertices += 0;
5412 n_single_lines += 3;
5418 needed_vertices += 1;
5419 n_single_lines += 4;
5428 auto line = cell->
line(line_no);
5429 if (line->has_children() ==
false)
5430 line->set_user_flag();
5435 auto &next_level = *triangulation.
levels[level_no + 1];
5436 const auto used_cells = std::count(next_level.cells.used.begin(),
5437 next_level.cells.used.end(),
5440 if (used_cells + needed_cells > next_level.size())
5441 next_level.allocate_end((used_cells + needed_cells) -
5448 next_level.cells.allocate_end(needed_cells, 0);
5454 if (line->user_flag_set())
5457 n_lines_in_pairs += 2;
5458 needed_vertices += 1;
5461 triangulation.
faces->lines.allocate_end(n_lines_in_pairs, 0);
5463 needed_vertices += std::count(triangulation.
vertices_used.begin(),
5467 if (needed_vertices > triangulation.
vertices.size())
5473 unsigned int next_unused_vertex = 0;
5482 for (; line != endl; ++line)
5483 if (line->user_flag_set())
5487 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5488 ++next_unused_vertex;
5491 "Internal error: During refinement, the triangulation "
5492 "wants to access an element of the 'vertices' array "
5493 "but it turns out that the array is not large "
5497 triangulation.
vertices[next_unused_vertex] = line->center(
true);
5499 bool pair_found =
false;
5500 for (; next_unused_line != endl; ++next_unused_line)
5501 if (!next_unused_line->used() &&
5502 !(++next_unused_line)->used())
5510 line->set_children(0, next_unused_line->index());
5513 children[2] = {next_unused_line, ++next_unused_line};
5518 children[0]->set_bounding_object_indices(
5519 {line->vertex_index(0), next_unused_vertex});
5520 children[1]->set_bounding_object_indices(
5521 {next_unused_vertex, line->vertex_index(1)});
5523 for (
auto &child : children)
5525 child->set_used_flag();
5526 child->clear_children();
5527 child->clear_user_data();
5528 child->clear_user_flag();
5529 child->set_boundary_id_internal(line->boundary_id());
5530 child->set_manifold_id(line->manifold_id());
5535 line->clear_user_flag();
5539 triangulation.
faces->lines.allocate_end(0, n_single_lines);
5543 cells_with_distorted_children;
5549 unsigned int &next_unused_vertex,
5550 auto &next_unused_line,
5551 auto &next_unused_cell,
5556 unsigned int n_new_vertices = 0;
5566 for (
unsigned int vertex_no = 0; vertex_no < cell->
n_vertices();
5568 new_vertices[vertex_no] = cell->
vertex_index(vertex_no);
5569 for (
unsigned int line_no = 0; line_no < cell->
n_lines(); ++line_no)
5570 if (cell->
line(line_no)->has_children())
5572 cell->
line(line_no)->child(0)->vertex_index(1);
5576 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5577 ++next_unused_vertex;
5579 next_unused_vertex < triangulation.
vertices.size(),
5581 "Internal error: During refinement, the triangulation wants "
5582 "to access an element of the 'vertices' array but it turns "
5583 "out that the array is not large enough."));
5586 new_vertices[8] = next_unused_vertex;
5588 triangulation.
vertices[next_unused_vertex] =
5589 cell->
center(
true,
true);
5592 std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
5595 std::array<types::geometric_orientation, 12> inherited_orientations;
5597 unsigned int lmin = 0;
5598 unsigned int lmax = 0;
5607 std::fill(inherited_orientations.begin() + lmin,
5608 inherited_orientations.begin() + lmax,
5621 for (
unsigned int l = lmin; l < lmax; ++l)
5623 while (next_unused_line->used() ==
true)
5625 new_lines[l] = next_unused_line;
5632 for (
const unsigned int face_no : cell->
face_indices())
5639 const auto combined_orientation =
5641 Assert(combined_orientation ==
5643 combined_orientation ==
5646 for (
unsigned int c = 0; c < 2; ++c)
5648 new_lines[2 * face_no + c] = cell->
line(face_no)->child(c);
5649 inherited_orientations[2 * face_no + c] =
5653 std::swap(new_lines[2 * face_no], new_lines[2 * face_no + 1]);
5659 new_lines[6]->set_bounding_object_indices(
5660 {new_vertices[3], new_vertices[4]});
5661 new_lines[7]->set_bounding_object_indices(
5662 {new_vertices[4], new_vertices[5]});
5663 new_lines[8]->set_bounding_object_indices(
5664 {new_vertices[5], new_vertices[3]});
5668 new_lines[8]->set_bounding_object_indices(
5669 {new_vertices[6], new_vertices[8]});
5670 new_lines[9]->set_bounding_object_indices(
5671 {new_vertices[8], new_vertices[7]});
5672 new_lines[10]->set_bounding_object_indices(
5673 {new_vertices[4], new_vertices[8]});
5674 new_lines[11]->set_bounding_object_indices(
5675 {new_vertices[8], new_vertices[5]});
5682 for (
unsigned int l = lmin; l < lmax; ++l)
5684 new_lines[l]->set_used_flag();
5685 new_lines[l]->clear_user_flag();
5686 new_lines[l]->clear_user_data();
5687 new_lines[l]->clear_children();
5689 new_lines[l]->set_boundary_id_internal(
5691 new_lines[l]->set_manifold_id(cell->
manifold_id());
5696 while (next_unused_cell->used() ==
true)
5699 unsigned int n_children = 0;
5707 for (
unsigned int i = 0; i < n_children; ++i)
5710 subcells[i] = next_unused_cell;
5712 if (i % 2 == 1 && i < n_children - 1)
5713 while (next_unused_cell->used() ==
true)
5719 constexpr ::ndarray<unsigned int, 4, 4> tri_child_lines = {
5724 constexpr ::ndarray<unsigned int, 4, 4> quad_child_lines = {
5731 const auto &child_lines =
5735 for (
unsigned int i = 0; i < n_children; ++i)
5739 {new_lines[child_lines[i][0]]->index(),
5740 new_lines[child_lines[i][1]]->index(),
5741 new_lines[child_lines[i][2]]->index()});
5744 {new_lines[child_lines[i][0]]->index(),
5745 new_lines[child_lines[i][1]]->index(),
5746 new_lines[child_lines[i][2]]->index(),
5747 new_lines[child_lines[i][3]]->index()});
5765 ->face_orientations.n_objects() > 0)
5768 face_no, inherited_orientations[child_lines[i][face_no]]);
5773 Assert(inherited_orientations[child_lines[i][face_no]] ==
5793 for (
unsigned int i = 0; i < n_children / 2; ++i)
5798 if (dim == spacedim - 1)
5799 for (
unsigned int c = 0; c < n_children; ++c)
5804 level < static_cast<int>(triangulation.
levels.size()) - 1;
5808 next_unused_cell = triangulation.
begin_raw(level + 1);
5810 for (
const auto &cell :
5821 check_for_distorted_cells &&
5822 has_distorted_children<dim, spacedim>(cell))
5830 return cells_with_distorted_children;
5839 template <
int spacedim>
5844 const unsigned int dim = 1;
5849 triangulation.
levels.size() - 1))
5852 triangulation.
levels.push_back(
5855 triangulation.
strides.max_children_per_cell,
5856 triangulation.
strides.max_faces_per_cell,
5857 triangulation.
strides.max_vertices_per_cell));
5865 std::size_t needed_vertices = 0;
5866 for (
int level_no = triangulation.
levels.size() - 2; level_no >= 0;
5869 std::size_t needed_cells = 0;
5871 for (
const auto &cell :
5877 auto &next_level = *triangulation.
levels[level_no + 1];
5878 const auto used_cells = std::count(next_level.cells.used.begin(),
5879 next_level.cells.used.end(),
5882 if (used_cells + needed_cells > next_level.size())
5883 next_level.allocate_end((used_cells + needed_cells) -
5891 next_level.cells.allocate_end(needed_cells, 0);
5892 needed_vertices += needed_cells / 2;
5897 needed_vertices += std::count(triangulation.
vertices_used.begin(),
5903 if (needed_vertices > triangulation.
vertices.size())
5914 unsigned int next_unused_vertex = 0;
5916 for (
int level = triangulation.
levels.size() - 2; level >= 0; --level)
5919 next_unused_cell = triangulation.
begin_raw(level + 1);
5921 for (
const auto &cell :
5932 ++next_unused_vertex;
5934 next_unused_vertex < triangulation.
vertices.size(),
5936 "Internal error: During refinement, the triangulation "
5937 "wants to access an element of the 'vertices' array "
5938 "but it turns out that the array is not large enough."));
5943 triangulation.
vertices[next_unused_vertex] =
5953 while (next_unused_cell->
used() ==
true)
5955 first_child = next_unused_cell;
5960 second_child = next_unused_cell;
5974 if (dim == spacedim - 1)
6018 left_neighbor = left_neighbor->
child(nbnb);
6031 if (dim == spacedim - 1)
6050 right_neighbor = cell->
neighbor(1);
6053 right_neighbor = right_neighbor->
child(nbnb);
6074 template <
int spacedim>
6077 const bool check_for_distorted_cells)
6079 const unsigned int dim = 2;
6084 bool do_isotropic_refinement =
true;
6089 do_isotropic_refinement =
false;
6093 if (do_isotropic_refinement)
6095 check_for_distorted_cells);
6103 triangulation.
levels.size() - 1))
6106 triangulation.
levels.push_back(
6109 triangulation.
strides.max_children_per_cell,
6110 triangulation.
strides.max_faces_per_cell,
6111 triangulation.
strides.max_vertices_per_cell));
6122 line->clear_user_flag();
6129 unsigned int n_single_lines = 0;
6133 unsigned int n_lines_in_pairs = 0;
6139 unsigned int needed_vertices = 0;
6141 const bool orientation_needed =
6144 for (
int level_no = triangulation.
levels.size() - 2; level_no >= 0;
6149 std::size_t needed_cells = 0;
6151 for (
const auto &cell :
6166 n_single_lines += 4;
6178 n_single_lines += 1;
6186 for (
const unsigned int line_no :
6194 line = cell->
line(line_no);
6195 if (line->has_children() ==
false)
6196 line->set_user_flag();
6201 auto &next_level = *triangulation.
levels[level_no + 1];
6202 const auto used_cells = std::count(next_level.cells.used.begin(),
6203 next_level.cells.used.end(),
6205 if (used_cells + needed_cells > next_level.size())
6206 next_level.allocate_end((used_cells + needed_cells) -
6213 next_level.cells.allocate_end(needed_cells, 0);
6221 if (line->user_flag_set())
6224 n_lines_in_pairs += 2;
6225 needed_vertices += 1;
6234 triangulation.
faces->lines.allocate_end(n_lines_in_pairs, 0);
6237 needed_vertices += std::count(triangulation.
vertices_used.begin(),
6243 if (needed_vertices > triangulation.
vertices.size())
6254 unsigned int next_unused_vertex = 0;
6266 for (; line != endl; ++line)
6267 if (line->user_flag_set())
6273 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
6274 ++next_unused_vertex;
6276 next_unused_vertex < triangulation.
vertices.size(),
6278 "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."));
6281 triangulation.
vertices[next_unused_vertex] = line->center(
true);
6286 bool pair_found =
false;
6287 for (; next_unused_line != endl; ++next_unused_line)
6288 if (!next_unused_line->used() &&
6289 !(++next_unused_line)->used())
6302 line->set_children(0, next_unused_line->index());
6306 children[2] = {next_unused_line, ++next_unused_line};
6312 children[0]->set_bounding_object_indices(
6313 {line->vertex_index(0), next_unused_vertex});
6314 children[1]->set_bounding_object_indices(
6315 {next_unused_vertex, line->vertex_index(1)});
6317 children[0]->set_used_flag();
6318 children[1]->set_used_flag();
6319 children[0]->clear_children();
6320 children[1]->clear_children();
6323 children[0]->clear_user_flag();
6324 children[1]->clear_user_flag();
6327 children[0]->set_boundary_id_internal(line->boundary_id());
6328 children[1]->set_boundary_id_internal(line->boundary_id());
6330 children[0]->set_manifold_id(line->manifold_id());
6331 children[1]->set_manifold_id(line->manifold_id());
6335 line->clear_user_flag();
6344 triangulation.
faces->lines.allocate_end(0, n_single_lines);
6347 cells_with_distorted_children;
6355 level < static_cast<int>(triangulation.
levels.size()) - 1;
6359 next_unused_cell = triangulation.
begin_raw(level + 1);
6361 for (
const auto &cell :
6373 if (check_for_distorted_cells &&
6374 has_distorted_children<dim, spacedim>(cell))
6382 return cells_with_distorted_children;
6386 template <
int spacedim>
6389 const bool check_for_distorted_cells)
6391 constexpr int dim = 3;
6393 using raw_line_iterator =
6395 using raw_quad_iterator =
6407 triangulation.
levels.size() - 1))
6413 triangulation.
levels.push_back(
6416 triangulation.
strides.max_children_per_cell,
6417 triangulation.
strides.max_faces_per_cell,
6418 triangulation.
strides.max_vertices_per_cell));
6424 triangulation.
faces->quads.clear_user_data();
6425 triangulation.
faces->lines.clear_user_flags();
6426 triangulation.
faces->quads.clear_user_flags();
6441 unsigned int needed_vertices = 0;
6442 unsigned int needed_lines_single = 0;
6443 unsigned int needed_faces_single = 0;
6444 unsigned int needed_lines_pair = 0;
6445 unsigned int needed_faces_pair = 0;
6446 for (
int level_no = triangulation.
levels.size() - 2; level_no >= 0;
6449 std::size_t needed_cells = 0;
6451 for (
const auto &cell :
6465 needed_cells += cell_reference_cell.n_isotropic_children();
6467 switch (cell_reference_cell)
6470 needed_lines_single += 1;
6471 needed_faces_single += 8;
6479 needed_lines_single += 4;
6480 needed_faces_single += 13;
6493 needed_lines_single += 3;
6494 needed_faces_single += 10;
6500 needed_lines_single += 6;
6501 needed_faces_single += 12;
6515 if (cell->
face(face)->n_children() == 0)
6516 cell->
face(face)->set_user_flag();
6520 cell_reference_cell.face_reference_cell(face)
6521 .n_isotropic_children(),
6526 if (cell->
line(line)->has_children() ==
false)
6527 cell->
line(line)->set_user_flag();
6533 auto &next_level = *triangulation.
levels[level_no + 1];
6534 const auto used_cells = std::count(next_level.cells.used.begin(),
6535 next_level.cells.used.end(),
6537 if (used_cells + needed_cells > next_level.size())
6538 next_level.allocate_end(
6539 (used_cells + needed_cells) - next_level.size(),
6549 next_level.cells.allocate_end(needed_cells, 0);
6559 if (quad->user_flag_set() ==
false)
6564 needed_faces_pair += 4;
6565 needed_lines_pair += 4;
6566 needed_vertices += 1;
6570 needed_faces_pair += 4;
6571 needed_lines_single += 3;
6584 if (line->user_flag_set() ==
false)
6587 needed_lines_pair += 2;
6588 needed_vertices += 1;
6591 triangulation.
faces->allocate_end(needed_faces_pair,
6592 needed_faces_single);
6593 triangulation.
faces->lines.allocate_end(needed_lines_pair,
6594 needed_lines_single);
6597 needed_vertices += std::count(triangulation.
vertices_used.begin(),
6601 if (needed_vertices > triangulation.
vertices.size())
6622 for (
unsigned int line_n = 0; line_n < cell->
n_lines();
6624 if (cell->
line(line_n)->has_children())
6625 for (
unsigned int c = 0; c < 2; ++c)
6626 Assert(cell->
line(line_n)->child(c)->user_flag_set() ==
6631 unsigned int current_vertex = 0;
6635 auto get_next_unused_vertex = [](
const unsigned int current_vertex,
6636 std::vector<bool> &vertices_used) {
6637 unsigned int next_vertex = current_vertex;
6638 while (next_vertex < vertices_used.size() &&
6639 vertices_used[next_vertex] ==
true)
6642 vertices_used[next_vertex] =
true;
6652 raw_line_iterator next_unused_line = triangulation.
begin_raw_line();
6654 for (; line != endl; ++line)
6656 if (line->user_flag_set() ==
false)
6660 triangulation.
faces->lines.template next_free_pair_object<1>(
6668 line->set_children(0, next_unused_line->index());
6670 const std::array<raw_line_iterator, 2> children{
6671 {next_unused_line, ++next_unused_line}};
6677 get_next_unused_vertex(current_vertex,
6679 triangulation.
vertices[current_vertex] = line->center(
true);
6681 children[0]->set_bounding_object_indices(
6682 {line->vertex_index(0), current_vertex});
6683 children[1]->set_bounding_object_indices(
6684 {current_vertex, line->vertex_index(1)});
6686 const auto manifold_id = line->manifold_id();
6687 const auto boundary_id = line->boundary_id();
6688 for (
const auto &child : children)
6690 child->set_used_flag();
6691 child->clear_children();
6692 child->clear_user_data();
6693 child->clear_user_flag();
6694 child->set_boundary_id_internal(boundary_id);
6695 child->set_manifold_id(manifold_id);
6698 line->clear_user_flag();
6708 for (; face != endf; ++face)
6710 if (face->user_flag_set() ==
false)
6713 const auto reference_face_type = face->reference_cell();
6717 std::array<raw_line_iterator, 4> new_lines;
6718 switch (reference_face_type)
6721 for (
unsigned int l = 0; l < 3; ++l)
6723 triangulation.
faces->lines
6724 .template next_free_single_object<1>(triangulation);
6728 for (
unsigned int l = 0; l < 2; ++l)
6730 auto next_unused_line =
6731 triangulation.
faces->lines
6732 .template next_free_pair_object<1>(triangulation);
6733 new_lines[2 * l] = next_unused_line;
6734 new_lines[2 * l + 1] = ++next_unused_line;
6744 for (
const unsigned int line : face->line_indices())
6753 std::array<raw_quad_iterator, 4> new_faces;
6754 for (
unsigned int f = 0; f < 2; ++f)
6756 auto next_unused_quad =
6757 triangulation.
faces->quads
6758 .template next_free_pair_object<2>(triangulation);
6760 new_faces[2 * f] = next_unused_quad;
6761 new_faces[2 * f + 1] = ++next_unused_quad;
6763 face->set_children(2 * f, new_faces[2 * f]->
index());
6770 for (
const auto &quad : new_faces)
6783 std::array<unsigned int, 9> vertex_indices = {};
6785 for (
const auto i : face->vertex_indices())
6786 vertex_indices[k++] = face->vertex_index(i);
6788 for (
const auto i : face->line_indices())
6789 vertex_indices[k++] = face->line(i)->child(0)->vertex_index(1);
6794 get_next_unused_vertex(current_vertex,
6796 vertex_indices[k++] = current_vertex;
6798 triangulation.
vertices[current_vertex] =
6799 face->center(
true,
true);
6805 std::array<raw_line_iterator, 12> lines;
6806 unsigned int n_lines = 0;
6807 for (
unsigned int l = 0; l < face->n_lines(); ++l)
6808 for (
unsigned int c = 0; c < 2; ++c)
6809 lines[n_lines++] = face->line(l)->child(
6810 child_line_index(c, face->line_orientation(l)));
6812 for (
unsigned int l = 0; l < face->n_lines(); ++l)
6813 lines[n_lines++] = new_lines[l];
6815 std::array<int, 12> line_indices;
6816 for (
unsigned int i = 0; i < n_lines; ++i)
6817 line_indices[i] = lines[i]->
index();
6824 constexpr ::ndarray<unsigned int, 12, 2> line_vertices_quad{
6838 constexpr ::ndarray<unsigned int, 4, 4> quad_lines_quad{
6853 constexpr ::ndarray<unsigned int, 12, 2> line_vertices_tri{
6867 constexpr ::ndarray<unsigned int, 4, 4> tri_lines_tri{
6877 constexpr ::ndarray<unsigned int, 4, 4, 2>
6878 tri_line_vertices_tri{
6879 {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
6880 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
6881 {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
6882 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
6885 const auto &line_vertices =
6887 line_vertices_quad :
6889 const auto &face_lines =
6897 for (
unsigned int i = 0, j = 2 * face->n_lines();
6898 i < face->n_lines();
6901 auto &new_line = new_lines[i];
6902 new_line->set_bounding_object_indices(
6903 {vertex_indices[line_vertices[j][0]],
6904 vertex_indices[line_vertices[j][1]]});
6905 new_line->set_used_flag();
6906 new_line->clear_user_flag();
6907 new_line->clear_user_data();
6908 new_line->clear_children();
6909 new_line->set_boundary_id_internal(face->boundary_id());
6910 new_line->set_manifold_id(face->manifold_id());
6914 for (
unsigned int i = 0; i < new_faces.size(); ++i)
6917 auto &new_face = new_faces[i];
6921 triangulation.
faces->set_quad_type(new_face->index(),
6922 reference_face_type);
6924 switch (reference_face_type)
6927 new_face->set_bounding_object_indices(
6928 {line_indices[face_lines[i][0]],
6929 line_indices[face_lines[i][1]],
6930 line_indices[face_lines[i][2]]});
6934 new_face->set_bounding_object_indices(
6935 {line_indices[face_lines[i][0]],
6936 line_indices[face_lines[i][1]],
6937 line_indices[face_lines[i][2]],
6938 line_indices[face_lines[i][3]]});
6945 new_face->set_used_flag();
6946 new_face->clear_user_flag();
6947 new_face->clear_user_data();
6948 new_face->clear_children();
6949 new_face->set_boundary_id_internal(face->boundary_id());
6950 new_face->set_manifold_id(face->manifold_id());
6952 [[maybe_unused]] std::set<unsigned int> s;
6960 for (
const auto f : new_face->line_indices())
6962 const std::array<unsigned int, 2> vertices_0 = {
6963 {lines[face_lines[i][f]]->vertex_index(0),
6964 lines[face_lines[i][f]]->vertex_index(1)}};
6966 const std::array<unsigned int, 2> vertices_1 = {
6967 {vertex_indices[tri_line_vertices_tri[i][f][0]],
6968 vertex_indices[tri_line_vertices_tri[i][f][1]]}};
6970 const auto orientation =
6978 for (
const auto i : vertices_0)
6980 for (
const auto i : vertices_1)
6984 new_face->set_line_orientation(f, orientation);
6998 constexpr ::ndarray<unsigned int, 4, 2>
6999 quad_child_line_vertices{{{{0, 2}},
7004 for (
unsigned int i = 0; i < 4; ++i)
7005 for (
unsigned int j = 0; j < 2; ++j)
7006 new_faces[quad_child_line_vertices[i][j]]
7007 ->set_line_orientation(i, face->line_orientation(i));
7010 face->clear_user_flag();
7016 cells_with_distorted_children;
7020 for (
unsigned int level = 0; level != triangulation.
levels.size() - 1;
7024 next_unused_cell = triangulation.
begin_raw(level + 1);
7026 cell->
level() >=
static_cast<int>(level),
7030 for (; cell != triangulation.
end() &&
7031 cell->
level() ==
static_cast<int>(level);
7046 unsigned int n_new_lines = 0;
7047 unsigned int n_new_faces = 0;
7048 const unsigned int n_new_cells =
7049 cell_reference_cell.n_isotropic_children();
7053 switch (cell_reference_cell)
7081 constexpr unsigned int max_n_new_lines = 6;
7082 constexpr unsigned int max_n_new_faces = 13;
7083 constexpr unsigned int max_n_new_children = 10;
7084 constexpr unsigned int max_relevant_vertices = 27;
7085 constexpr unsigned int max_relevant_lines = 30;
7086 constexpr unsigned int max_face_indices = 36;
7088 std::array<raw_line_iterator, max_n_new_lines> new_lines;
7089 for (
unsigned int i = 0; i < n_new_lines; ++i)
7092 triangulation.
faces->lines
7093 .template next_free_single_object<1>(triangulation);
7096 new_lines[i]->set_used_flag();
7097 new_lines[i]->clear_user_flag();
7098 new_lines[i]->clear_user_data();
7099 new_lines[i]->clear_children();
7100 new_lines[i]->set_boundary_id_internal(
7102 new_lines[i]->set_manifold_id(cell->
manifold_id());
7105 std::array<raw_quad_iterator, max_n_new_faces> new_faces;
7106 for (
unsigned int i = 0; i < n_new_faces; ++i)
7109 triangulation.
faces->quads
7110 .template next_free_single_object<2>(triangulation);
7112 auto &new_face = new_faces[i];
7115 new_face->set_used_flag();
7116 new_face->clear_user_flag();
7117 new_face->clear_user_data();
7118 new_face->clear_children();
7119 new_face->set_boundary_id_internal(
7126 for (
unsigned int j = 0;
7127 j < triangulation.
faces->lines_per_quad;
7129 new_face->set_line_orientation(
7182 std::array<unsigned int, max_relevant_vertices>
7183 vertex_indices = {};
7191 const unsigned int n_vertices =
7195 for (
unsigned int i = 0; i < n_vertices; ++i)
7204 const std::array<unsigned int, 12> line_indices =
7205 TriaAccessorImplementation::Implementation::
7206 get_line_indices_of_cell(*cell);
7208 for (
unsigned int l = 0; l < n_lines; ++l)
7210 raw_line_iterator line(&triangulation,
7213 vertex_indices[k++] = line->child(0)->vertex_index(1);
7218 for (
auto quad_face_index :
7219 cell_reference_cell.face_indices_by_type(
7222 vertex_indices[k++] = cell->
face(quad_face_index)
7232 get_next_unused_vertex(current_vertex,
7234 vertex_indices[k++] = current_vertex;
7236 triangulation.
vertices[current_vertex] =
7237 cell->
center(
true,
true);
7242 unsigned int chosen_line_tetrahedron = 0;
7245 switch (cell_reference_cell)
7253 constexpr ::ndarray<unsigned int, 3, 2>
7254 new_line_vertices = {{{{6, 8}},
7260 std::uint8_t refinement_choice =
7262 if (refinement_choice ==
7266 const auto &vertices =
7268 double min_distance =
7269 std::numeric_limits<double>::infinity();
7270 for (
unsigned int i = 0;
7271 i < new_line_vertices.size();
7274 const double current_distance =
7275 vertices[vertex_indices
7276 [new_line_vertices[i][0]]]
7278 vertices[vertex_indices
7279 [new_line_vertices[i][1]]]);
7280 if (current_distance < min_distance)
7282 chosen_line_tetrahedron = i;
7283 min_distance = current_distance;
7287 else if (refinement_choice ==
7290 chosen_line_tetrahedron = 0;
7291 else if (refinement_choice ==
7294 chosen_line_tetrahedron = 1;
7295 else if (refinement_choice ==
7298 chosen_line_tetrahedron = 2;
7305 new_lines[0]->set_bounding_object_indices(
7307 [new_line_vertices[chosen_line_tetrahedron][0]],
7308 vertex_indices[new_line_vertices
7309 [chosen_line_tetrahedron][1]]});
7318 static constexpr ::ndarray<unsigned int, 4, 2>
7319 new_line_vertices = {{{{13, 9}},
7324 for (
unsigned int i = 0; i < n_new_lines; ++i)
7325 new_lines[i]->set_bounding_object_indices(
7326 {vertex_indices[new_line_vertices[i][0]],
7327 vertex_indices[new_line_vertices[i][1]]});
7335 constexpr ::ndarray<unsigned int, 3, 2>
7336 new_line_vertices = {{{{15, 16}},
7340 for (
unsigned int i = 0; i < n_new_lines; ++i)
7341 new_lines[i]->set_bounding_object_indices(
7342 {vertex_indices[new_line_vertices[i][0]],
7343 vertex_indices[new_line_vertices[i][1]]});
7349 constexpr ::ndarray<unsigned int, 6, 2>
7350 new_line_vertices = {{{{22, 26}},
7356 for (
unsigned int i = 0; i < n_new_lines; ++i)
7357 new_lines[i]->set_bounding_object_indices(
7358 {vertex_indices[new_line_vertices[i][0]],
7359 vertex_indices[new_line_vertices[i][1]]});
7370 auto standard_to_real_quad_child_index =
7371 [](
const unsigned int child,
7382 switch (face_ref_cell_type)
7389 .standard_to_real_face_vertex(child,
7398 .standard_to_real_face_vertex(child,
7440 std::array<raw_line_iterator, max_relevant_lines>
7442 std::array<unsigned int, max_relevant_lines>
7443 relevant_line_indices;
7445 unsigned int relevant_lines_counter = 0;
7448 for (
auto f : cell_reference_cell.face_indices_by_type(
7468 constexpr ::ndarray<unsigned int, 4, 2>
7469 quad_relevant_children_and_lines = {{{{0, 1}},
7473 for (
unsigned int l = 0; l < 4;
7474 ++l, ++relevant_lines_counter)
7478 const unsigned int relevant_child =
7479 standard_to_real_quad_child_index(
7480 quad_relevant_children_and_lines[l][0],
7487 const unsigned int relevant_child_line =
7489 .standard_to_real_face_line(
7490 quad_relevant_children_and_lines[l][1],
7497 relevant_lines[relevant_lines_counter] =
7499 ->isotropic_child(relevant_child)
7500 ->line(relevant_child_line);
7505 for (
auto f : cell_reference_cell.face_indices_by_type(
7508 for (
unsigned int l = 0; l < 3;
7509 ++l, ++relevant_lines_counter)
7514 const std::array<std::array<unsigned int, 3>, 6>
7525 tri_line_perm = {{{{0, 1, 2}},
7532 const auto combined_orientation =
7534 relevant_lines[relevant_lines_counter] =
7537 ->line(tri_line_perm[combined_orientation][l]);
7542 for (
unsigned int i = 0; i < n_new_lines;
7543 ++i, ++relevant_lines_counter)
7544 relevant_lines[relevant_lines_counter] = new_lines[i];
7548 for (
unsigned int i = 0; i < relevant_lines_counter; ++i)
7549 relevant_line_indices[i] = relevant_lines[i]->
index();
7558 const auto &new_face_lines =
7559 cell_reference_cell.new_isotropic_child_face_lines(
7560 chosen_line_tetrahedron);
7568 const auto &new_face_lines_vert =
7570 .new_isotropic_child_face_line_vertices(
7571 chosen_line_tetrahedron);
7577 constexpr ::ndarray<unsigned int, 4, 2>
7578 representative_lines{{{{0, 2}},
7583 for (
unsigned int q = 0; q < n_new_faces; ++q)
7585 auto &new_face = new_faces[q];
7587 if (new_face_lines[q][3] == X)
7591 triangulation.
faces->set_quad_type(
7594 new_face->set_bounding_object_indices(
7595 {relevant_line_indices[new_face_lines[q][0]],
7596 relevant_line_indices[new_face_lines[q][1]],
7597 relevant_line_indices[new_face_lines[q][2]]});
7601 triangulation.
faces->set_quad_type(
7604 new_face->set_bounding_object_indices(
7605 {relevant_line_indices[new_face_lines[q][0]],
7606 relevant_line_indices[new_face_lines[q][1]],
7607 relevant_line_indices[new_face_lines[q][2]],
7608 relevant_line_indices[new_face_lines[q][3]]});
7615 const unsigned int n_compute_lines =
7618 new_face->n_lines();
7619 for (
unsigned int line = 0; line < n_compute_lines;
7622 const unsigned int l =
7623 (cell_reference_cell ==
7625 representative_lines[q % 4][0] :
7628 const std::array<unsigned int, 2> vertices_0 = {
7629 {relevant_lines[new_face_lines[q][l]]
7631 relevant_lines[new_face_lines[q][l]]
7632 ->vertex_index(1)}};
7634 const std::array<unsigned int, 2> vertices_1 = {
7635 {vertex_indices[new_face_lines_vert[q][l][0]],
7636 vertex_indices[new_face_lines_vert[q][l][1]]}};
7638 const auto orientation =
7643 new_face->set_line_orientation(l, orientation);
7648 if (cell_reference_cell ==
7650 new_faces[representative_lines[q % 4][1] + q -
7652 ->set_line_orientation(l, orientation);
7674 std::array<int, max_face_indices> face_indices;
7676 unsigned int face_indices_counter = 0;
7679 for (; face_indices_counter < n_new_faces;
7680 ++face_indices_counter)
7681 face_indices[face_indices_counter] =
7682 new_faces[face_indices_counter]->
index();
7685 for (
auto f_q : cell_reference_cell.face_indices_by_type(
7688 for (
unsigned int c = 0; c < 4;
7689 ++c, ++face_indices_counter)
7691 const auto combined_orientation =
7695 face_indices[face_indices_counter] =
7696 cell->
face(f_q)->isotropic_child_index(
7697 standard_to_real_quad_child_index(
7699 combined_orientation,
7705 for (
auto f_t : cell_reference_cell.face_indices_by_type(
7708 for (
unsigned int c = 0; c < 4;
7709 ++c, ++face_indices_counter)
7711 const auto combined_orientation =
7715 face_indices[face_indices_counter] =
7716 cell->
face(f_t)->child_index(
7717 standard_to_real_quad_child_index(
7719 combined_orientation,
7728 const auto cell_faces =
7729 cell_reference_cell.new_isotropic_child_cell_faces(
7730 chosen_line_tetrahedron);
7736 const auto new_cells_vertices =
7737 cell_reference_cell.new_isotropic_child_cell_vertices(
7738 chosen_line_tetrahedron);
7743 const ::types::subdomain_id cell_subdomain_id =
7751 for (
unsigned int c = 0; c < n_new_cells; ++c)
7758 triangulation.
levels[level + 1]
7759 ->cells.next_free_hex(triangulation, level + 1);
7763 new_cells[c] = next_unused_cell;
7765 auto &new_cell = new_cells[c];
7770 new_cell->set_used_flag();
7771 new_cell->clear_user_flag();
7772 new_cell->clear_user_data();
7773 new_cell->clear_children();
7776 new_cell->set_subdomain_id(cell_subdomain_id);
7782 new_cell->set_parent(cell->
index());
7796 triangulation.
levels[new_cell->level()]
7797 ->reference_cell[new_cell->index()] =
7798 child_reference_cell;
7804 switch (child_reference_cell)
7807 new_cell->set_bounding_object_indices(
7808 {face_indices[cell_faces[c][0]],
7809 face_indices[cell_faces[c][1]],
7810 face_indices[cell_faces[c][2]],
7811 face_indices[cell_faces[c][3]]});
7817 new_cell->set_bounding_object_indices(
7818 {face_indices[cell_faces[c][0]],
7819 face_indices[cell_faces[c][1]],
7820 face_indices[cell_faces[c][2]],
7821 face_indices[cell_faces[c][3]],
7822 face_indices[cell_faces[c][4]]});
7826 new_cell->set_bounding_object_indices(
7827 {face_indices[cell_faces[c][0]],
7828 face_indices[cell_faces[c][1]],
7829 face_indices[cell_faces[c][2]],
7830 face_indices[cell_faces[c][3]],
7831 face_indices[cell_faces[c][4]],
7832 face_indices[cell_faces[c][5]]});
7851 for (
const auto f : new_cell->face_indices())
7852 new_cell->set_combined_face_orientation(
7857 for (
const auto f : new_cell->face_indices())
7859 const auto &face = new_cell->face(f);
7862 auto new_cell_vertices = new_cells_vertices[c];
7865 std::array<unsigned int, 4> vertices_0,
7870 for (
unsigned int face_vertex_no :
7871 face->vertex_indices())
7876 const auto cell_vertex_no =
7877 child_reference_cell
7882 default_geometric_orientation);
7895 vertices_0[face_vertex_no] = vertex_indices
7896 [new_cell_vertices[cell_vertex_no]];
7900 for (
const auto i : face->vertex_indices())
7901 vertices_1[i] = face->vertex_index(i);
7906 new_cell->set_combined_face_orientation(
7908 face->reference_cell()
7909 .get_combined_orientation(
7911 vertices_0.cbegin() +
7912 face->n_vertices()),
7914 vertices_1.cbegin() +
7915 face->n_vertices())));
7924 constexpr ::ndarray<unsigned int, 6, 4>
7925 face_to_child_indices_hex{{{{0, 2, 4, 6}},
7934 const auto combined_orientation =
7936 for (
unsigned int c = 0; c < 4; ++c)
7937 new_cells[face_to_child_indices_hex[f][c]]
7938 ->set_combined_face_orientation(
7939 f, combined_orientation);
7945 if (check_for_distorted_cells &&
7946 has_distorted_children<dim, spacedim>(cell))
7953 triangulation.
faces->quads.clear_user_data();
7955 return cells_with_distorted_children;
7962 template <
int spacedim>
7965 const bool check_for_distorted_cells)
7967 const unsigned int dim = 3;
7970 bool flag_isotropic_mesh =
true;
7972 cell = triangulation.
begin(),
7973 endc = triangulation.
end();
7974 for (; cell != endc; ++cell)
7984 flag_isotropic_mesh =
false;
7988 if (flag_isotropic_mesh)
7990 check_for_distorted_cells);
8002 triangulation.
levels.size() - 1))
8008 triangulation.
levels.push_back(
8011 triangulation.
strides.max_children_per_cell,
8012 triangulation.
strides.max_faces_per_cell,
8013 triangulation.
strides.max_vertices_per_cell));
8020 triangulation.
faces->quads.clear_user_data();
8026 line->clear_user_flag();
8031 quad->clear_user_flag();
8054 unsigned int needed_vertices = 0;
8055 unsigned int needed_lines_single = 0;
8056 unsigned int needed_quads_single = 0;
8057 unsigned int needed_lines_pair = 0;
8058 unsigned int needed_quads_pair = 0;
8059 for (
int level_no = triangulation.
levels.size() - 2; level_no >= 0;
8064 std::size_t needed_cells = 0;
8066 for (
const auto &acell :
8068 if (acell->refine_flag_set())
8078 ++needed_quads_single;
8086 ++needed_lines_single;
8087 needed_quads_single += 4;
8094 needed_lines_single += 6;
8095 needed_quads_single += 12;
8109 for (
const unsigned int face :
8113 aface = acell->face(face);
8120 acell->face_orientation(face),
8121 acell->face_flip(face),
8122 acell->face_rotation(face));
8127 if (face_ref_case ==
8130 if (aface->n_active_descendants() < 4)
8133 aface->set_user_flag();
8135 else if (aface->refinement_case() != face_ref_case)
8146 Assert(aface->refinement_case() ==
8149 aface->refinement_case() ==
8152 aface->set_user_index(face_ref_case);
8158 for (
unsigned int line = 0;
8159 line < GeometryInfo<dim>::lines_per_cell;
8163 !acell->line(line)->has_children())
8164 acell->line(line)->set_user_flag();
8169 auto &next_level = *triangulation.
levels[level_no + 1];
8170 const auto used_cells = std::count(next_level.cells.used.begin(),
8171 next_level.cells.used.end(),
8173 if (used_cells + needed_cells > next_level.size())
8174 next_level.allocate_end(
8175 (used_cells + needed_cells) - next_level.size(),
8185 next_level.cells.allocate_end(needed_cells, 0);
8194 if (quad->user_flag_set())
8200 needed_quads_pair += 4;
8201 needed_lines_pair += 4;
8202 needed_vertices += 1;
8204 if (quad->user_index())
8208 needed_quads_pair += 2;
8209 needed_lines_single += 1;
8222 if (quad->has_children())
8224 Assert(quad->refinement_case() ==
8227 if ((face_refinement_cases[quad->user_index()] ==
8229 (quad->child(0)->line_index(1) + 1 !=
8230 quad->child(2)->line_index(1))) ||
8231 (face_refinement_cases[quad->user_index()] ==
8233 (quad->child(0)->line_index(3) + 1 !=
8234 quad->child(1)->line_index(3))))
8235 needed_lines_pair += 2;
8244 if (line->user_flag_set())
8246 needed_lines_pair += 2;
8247 needed_vertices += 1;
8251 triangulation.
faces->lines.allocate_end(needed_lines_pair,
8252 needed_lines_single);
8255 triangulation.
faces->allocate_end(needed_quads_pair,
8256 needed_quads_single);
8259 needed_vertices += std::count(triangulation.
vertices_used.begin(),
8265 if (needed_vertices > triangulation.
vertices.size())
8286 for (
unsigned int line = 0;
8287 line < GeometryInfo<dim>::lines_per_cell;
8289 if (cell->
line(line)->has_children())
8290 for (
unsigned int c = 0; c < 2; ++c)
8291 Assert(cell->
line(line)->child(c)->user_flag_set() ==
8304 unsigned int next_unused_vertex = 0;
8315 for (; line != endl; ++line)
8316 if (line->user_flag_set())
8322 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
8323 ++next_unused_vertex;
8325 next_unused_vertex < triangulation.
vertices.size(),
8327 "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."));
8330 triangulation.
vertices[next_unused_vertex] = line->center(
true);
8336 triangulation.
faces->lines.template next_free_pair_object<1>(
8344 line->set_children(0, next_unused_line->index());
8348 children[2] = {next_unused_line, ++next_unused_line};
8355 children[0]->set_bounding_object_indices(
8356 {line->vertex_index(0), next_unused_vertex});
8357 children[1]->set_bounding_object_indices(
8358 {next_unused_vertex, line->vertex_index(1)});
8360 children[0]->set_used_flag();
8361 children[1]->set_used_flag();
8362 children[0]->clear_children();
8363 children[1]->clear_children();
8366 children[0]->clear_user_flag();
8367 children[1]->clear_user_flag();
8369 children[0]->set_boundary_id_internal(line->boundary_id());
8370 children[1]->set_boundary_id_internal(line->boundary_id());
8372 children[0]->set_manifold_id(line->manifold_id());
8373 children[1]->set_manifold_id(line->manifold_id());
8378 line->clear_user_flag();
8408 for (
unsigned int loop = 0; loop < 2; ++loop)
8421 for (; quad != endq; ++quad)
8423 if (quad->user_index())
8426 face_refinement_cases[quad->user_index()];
8435 if (aniso_quad_ref_case == quad->refinement_case())
8438 Assert(quad->refinement_case() ==
8440 quad->refinement_case() ==
8445 Assert(quad->user_index() ==
8447 quad->user_index() ==
8456 triangulation.
faces->lines
8457 .template next_free_single_object<1>(triangulation);
8473 unsigned int vertex_indices[2];
8477 quad->line(2)->child(0)->vertex_index(1);
8479 quad->line(3)->child(0)->vertex_index(1);
8484 quad->line(0)->child(0)->vertex_index(1);
8486 quad->line(1)->child(0)->vertex_index(1);
8489 new_line->set_bounding_object_indices(
8490 {vertex_indices[0], vertex_indices[1]});
8491 new_line->set_used_flag();
8492 new_line->clear_user_flag();
8494 new_line->clear_children();
8495 new_line->set_boundary_id_internal(quad->boundary_id());
8496 new_line->set_manifold_id(quad->manifold_id());
8504 triangulation.
faces->quads
8505 .template next_free_pair_object<2>(triangulation);
8506 new_quads[0] = next_unused_quad;
8510 new_quads[1] = next_unused_quad;
8515 new_quads[0]->set_bounding_object_indices(
8516 {
static_cast<int>(quad->line_index(0)),
8520 child_line_index(0, quad->line_orientation(2)))
8524 child_line_index(0, quad->line_orientation(3)))
8526 new_quads[1]->set_bounding_object_indices(
8528 static_cast<int>(quad->line_index(1)),
8531 child_line_index(1, quad->line_orientation(2)))
8535 child_line_index(1, quad->line_orientation(3)))
8540 new_quads[0]->set_bounding_object_indices(
8543 child_line_index(0, quad->line_orientation(0)))
8547 child_line_index(0, quad->line_orientation(1)))
8549 static_cast<int>(quad->line_index(2)),
8550 new_line->index()});
8551 new_quads[1]->set_bounding_object_indices(
8554 child_line_index(1, quad->line_orientation(0)))
8558 child_line_index(1, quad->line_orientation(1)))
8561 static_cast<int>(quad->line_index(3))});
8564 for (
const auto &new_quad : new_quads)
8566 new_quad->set_used_flag();
8567 new_quad->clear_user_flag();
8568 new_quad->clear_user_data();
8569 new_quad->clear_children();
8570 new_quad->set_boundary_id_internal(quad->boundary_id());
8571 new_quad->set_manifold_id(quad->manifold_id());
8575 for (
unsigned int j = 0;
8576 j < GeometryInfo<dim>::lines_per_face;
8578 new_quad->set_line_orientation(
8585 new_quads[0]->set_line_orientation(
8586 0, quad->line_orientation(0));
8587 new_quads[0]->set_line_orientation(
8588 2, quad->line_orientation(2));
8589 new_quads[1]->set_line_orientation(
8590 1, quad->line_orientation(1));
8591 new_quads[1]->set_line_orientation(
8592 3, quad->line_orientation(3));
8595 new_quads[0]->set_line_orientation(
8596 3, quad->line_orientation(3));
8597 new_quads[1]->set_line_orientation(
8598 2, quad->line_orientation(2));
8602 new_quads[0]->set_line_orientation(
8603 1, quad->line_orientation(1));
8604 new_quads[1]->set_line_orientation(
8605 0, quad->line_orientation(0));
8611 if (quad->refinement_case() ==
8634 if (aniso_quad_ref_case ==
8637 old_child[0] = quad->child(0)->line(1);
8638 old_child[1] = quad->child(2)->line(1);
8642 Assert(aniso_quad_ref_case ==
8646 old_child[0] = quad->child(0)->line(3);
8647 old_child[1] = quad->child(1)->line(3);
8650 if (old_child[0]->
index() + 1 != old_child[1]->
index())
8656 spacedim>::raw_line_iterator
8659 new_child[0] = new_child[1] =
8660 triangulation.
faces->lines
8661 .template next_free_pair_object<1>(
8665 new_child[0]->set_used_flag();
8666 new_child[1]->set_used_flag();
8668 const int old_index_0 = old_child[0]->index(),
8669 old_index_1 = old_child[1]->index(),
8670 new_index_0 = new_child[0]->index(),
8671 new_index_1 = new_child[1]->index();
8675 for (
unsigned int q = 0;
8676 q < triangulation.
faces->quads.n_objects();
8678 for (
unsigned int l = 0;
8679 l < GeometryInfo<dim>::lines_per_face;
8682 const int this_index =
8683 triangulation.
faces->quads
8684 .get_bounding_object_indices(q)[l];
8685 if (this_index == old_index_0)
8686 triangulation.
faces->quads
8687 .get_bounding_object_indices(q)[l] =
8689 else if (this_index == old_index_1)
8690 triangulation.
faces->quads
8691 .get_bounding_object_indices(q)[l] =
8696 for (
unsigned int i = 0; i < 2; ++i)
8698 Assert(!old_child[i]->has_children(),
8701 new_child[i]->set_bounding_object_indices(
8702 {old_child[i]->vertex_index(0),
8703 old_child[i]->vertex_index(1)});
8704 new_child[i]->set_boundary_id_internal(
8705 old_child[i]->boundary_id());
8706 new_child[i]->set_manifold_id(
8707 old_child[i]->manifold_id());
8708 new_child[i]->set_user_index(
8709 old_child[i]->user_index());
8710 if (old_child[i]->user_flag_set())
8711 new_child[i]->set_user_flag();
8713 new_child[i]->clear_user_flag();
8715 new_child[i]->clear_children();
8717 old_child[i]->clear_user_flag();
8718 old_child[i]->clear_user_index();
8719 old_child[i]->clear_used_flag();
8725 if (aniso_quad_ref_case ==
8728 new_line->set_children(
8729 0, quad->child(0)->line_index(1));
8730 Assert(new_line->child(1) ==
8731 quad->child(2)->line(1),
8768 quad_iterator switch_1 = quad->child(1),
8769 switch_2 = quad->child(2);
8770 const int switch_1_index = switch_1->index();
8771 const int switch_2_index = switch_2->index();
8772 for (
unsigned int l = 0;
8773 l < triangulation.
levels.size();
8775 for (
unsigned int h = 0;
8777 triangulation.
levels[l]->cells.n_objects();
8779 for (
const unsigned int q :
8782 const int face_index =
8784 ->cells.get_bounding_object_indices(
8786 if (face_index == switch_1_index)
8788 ->cells.get_bounding_object_indices(
8789 h)[q] = switch_2_index;
8790 else if (face_index == switch_2_index)
8792 ->cells.get_bounding_object_indices(
8793 h)[q] = switch_1_index;
8797 const unsigned int switch_1_lines[4] = {
8798 switch_1->line_index(0),
8799 switch_1->line_index(1),
8800 switch_1->line_index(2),
8801 switch_1->line_index(3)};
8803 switch_1_line_orientations[4] = {
8804 switch_1->line_orientation(0),
8805 switch_1->line_orientation(1),
8806 switch_1->line_orientation(2),
8807 switch_1->line_orientation(3)};
8809 switch_1->boundary_id();
8810 const unsigned int switch_1_user_index =
8811 switch_1->user_index();
8812 const bool switch_1_user_flag =
8813 switch_1->user_flag_set();
8815 switch_1_refinement_case =
8816 switch_1->refinement_case();
8817 const int switch_1_first_child_pair =
8818 (switch_1_refinement_case ?
8819 switch_1->child_index(0) :
8821 const int switch_1_second_child_pair =
8822 (switch_1_refinement_case ==
8824 switch_1->child_index(2) :
8827 switch_1->set_bounding_object_indices(
8828 {switch_2->line_index(0),
8829 switch_2->line_index(1),
8830 switch_2->line_index(2),
8831 switch_2->line_index(3)});
8832 switch_1->set_line_orientation(
8833 0, switch_2->line_orientation(0));
8834 switch_1->set_line_orientation(
8835 1, switch_2->line_orientation(1));
8836 switch_1->set_line_orientation(
8837 2, switch_2->line_orientation(2));
8838 switch_1->set_line_orientation(
8839 3, switch_2->line_orientation(3));
8840 switch_1->set_boundary_id_internal(
8841 switch_2->boundary_id());
8842 switch_1->set_manifold_id(switch_2->manifold_id());
8843 switch_1->set_user_index(switch_2->user_index());
8844 if (switch_2->user_flag_set())
8845 switch_1->set_user_flag();
8847 switch_1->clear_user_flag();
8848 switch_1->clear_refinement_case();
8849 switch_1->set_refinement_case(
8850 switch_2->refinement_case());
8851 switch_1->clear_children();
8852 if (switch_2->refinement_case())
8853 switch_1->set_children(0,
8854 switch_2->child_index(0));
8855 if (switch_2->refinement_case() ==
8857 switch_1->set_children(2,
8858 switch_2->child_index(2));
8860 switch_2->set_bounding_object_indices(
8864 switch_1_lines[3]});
8865 switch_2->set_line_orientation(
8866 0, switch_1_line_orientations[0]);
8867 switch_2->set_line_orientation(
8868 1, switch_1_line_orientations[1]);
8869 switch_2->set_line_orientation(
8870 2, switch_1_line_orientations[2]);
8871 switch_2->set_line_orientation(
8872 3, switch_1_line_orientations[3]);
8873 switch_2->set_boundary_id_internal(
8874 switch_1_boundary_id);
8875 switch_2->set_manifold_id(switch_1->manifold_id());
8876 switch_2->set_user_index(switch_1_user_index);
8877 if (switch_1_user_flag)
8878 switch_2->set_user_flag();
8880 switch_2->clear_user_flag();
8881 switch_2->clear_refinement_case();
8882 switch_2->set_refinement_case(
8883 switch_1_refinement_case);
8884 switch_2->clear_children();
8885 switch_2->set_children(0,
8886 switch_1_first_child_pair);
8887 switch_2->set_children(2,
8888 switch_1_second_child_pair);
8890 new_quads[0]->set_refinement_case(
8892 new_quads[0]->set_children(0, quad->child_index(0));
8893 new_quads[1]->set_refinement_case(
8895 new_quads[1]->set_children(0, quad->child_index(2));
8899 new_quads[0]->set_refinement_case(
8901 new_quads[0]->set_children(0, quad->child_index(0));
8902 new_quads[1]->set_refinement_case(
8904 new_quads[1]->set_children(0, quad->child_index(2));
8905 new_line->set_children(
8906 0, quad->child(0)->line_index(3));
8907 Assert(new_line->child(1) ==
8908 quad->child(1)->line(3),
8911 quad->clear_children();
8915 quad->set_children(0, new_quads[0]->
index());
8917 quad->set_refinement_case(aniso_quad_ref_case);
8924 if (quad->user_flag_set())
8936 ++next_unused_vertex;
8938 next_unused_vertex < triangulation.
vertices.size(),
8940 "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."));
8948 quad->refinement_case();
8954 quad->child(0)->set_user_index(
8956 quad->child(1)->set_user_index(
8962 middle_line = quad->child(0)->line(1);
8964 middle_line = quad->child(0)->line(3);
8971 if (!middle_line->has_children())
8978 triangulation.
vertices[next_unused_vertex] =
8979 middle_line->center(
true);
8986 triangulation.
faces->lines
8987 .template next_free_pair_object<1>(
8992 middle_line->set_children(
8993 0, next_unused_line->index());
8997 raw_line_iterator children[2] = {
8998 next_unused_line, ++next_unused_line};
9006 children[0]->set_bounding_object_indices(
9007 {middle_line->vertex_index(0),
9008 next_unused_vertex});
9009 children[1]->set_bounding_object_indices(
9010 {next_unused_vertex,
9011 middle_line->vertex_index(1)});
9013 children[0]->set_used_flag();
9014 children[1]->set_used_flag();
9015 children[0]->clear_children();
9016 children[1]->clear_children();
9019 children[0]->clear_user_flag();
9020 children[1]->clear_user_flag();
9022 children[0]->set_boundary_id_internal(
9023 middle_line->boundary_id());
9024 children[1]->set_boundary_id_internal(
9025 middle_line->boundary_id());
9027 children[0]->set_manifold_id(
9028 middle_line->manifold_id());
9029 children[1]->set_manifold_id(
9030 middle_line->manifold_id());
9036 quad->clear_user_flag();
9065 triangulation.
vertices[next_unused_vertex] =
9066 quad->center(
true,
true);
9075 for (
unsigned int i = 0; i < 4; ++i)
9085 triangulation.
faces->lines
9086 .template next_free_pair_object<1>(triangulation);
9088 new_lines[i] = next_unused_line;
9111 const unsigned int vertex_indices[5] = {
9112 quad->line(0)->child(0)->vertex_index(1),
9113 quad->line(1)->child(0)->vertex_index(1),
9114 quad->line(2)->child(0)->vertex_index(1),
9115 quad->line(3)->child(0)->vertex_index(1),
9116 next_unused_vertex};
9118 new_lines[0]->set_bounding_object_indices(
9119 {vertex_indices[2], vertex_indices[4]});
9120 new_lines[1]->set_bounding_object_indices(
9121 {vertex_indices[4], vertex_indices[3]});
9122 new_lines[2]->set_bounding_object_indices(
9123 {vertex_indices[0], vertex_indices[4]});
9124 new_lines[3]->set_bounding_object_indices(
9125 {vertex_indices[4], vertex_indices[1]});
9127 for (
const auto &new_line : new_lines)
9129 new_line->set_used_flag();
9130 new_line->clear_user_flag();
9131 new_line->clear_user_data();
9132 new_line->clear_children();
9133 new_line->set_boundary_id_internal(quad->boundary_id());
9134 new_line->set_manifold_id(quad->manifold_id());
9147 const int line_indices[12] = {
9149 ->child(child_line_index(0, quad->line_orientation(0)))
9152 ->child(child_line_index(1, quad->line_orientation(0)))
9155 ->child(child_line_index(0, quad->line_orientation(1)))
9158 ->child(child_line_index(1, quad->line_orientation(1)))
9161 ->child(child_line_index(0, quad->line_orientation(2)))
9164 ->child(child_line_index(1, quad->line_orientation(2)))
9167 ->child(child_line_index(0, quad->line_orientation(3)))
9170 ->child(child_line_index(1, quad->line_orientation(3)))
9172 new_lines[0]->index(),
9173 new_lines[1]->index(),
9174 new_lines[2]->index(),
9175 new_lines[3]->index()};
9184 triangulation.
faces->quads
9185 .template next_free_pair_object<2>(triangulation);
9187 new_quads[0] = next_unused_quad;
9191 new_quads[1] = next_unused_quad;
9195 triangulation.
faces->quads
9196 .template next_free_pair_object<2>(triangulation);
9197 new_quads[2] = next_unused_quad;
9201 new_quads[3] = next_unused_quad;
9205 quad->set_children(0, new_quads[0]->
index());
9206 quad->set_children(2, new_quads[2]->
index());
9209 new_quads[0]->set_bounding_object_indices(
9214 new_quads[1]->set_bounding_object_indices(
9219 new_quads[2]->set_bounding_object_indices(
9224 new_quads[3]->set_bounding_object_indices(
9229 for (
const auto &new_quad : new_quads)
9231 new_quad->set_used_flag();
9232 new_quad->clear_user_flag();
9233 new_quad->clear_user_data();
9234 new_quad->clear_children();
9235 new_quad->set_boundary_id_internal(quad->boundary_id());
9236 new_quad->set_manifold_id(quad->manifold_id());
9240 for (
unsigned int j = 0;
9241 j < GeometryInfo<dim>::lines_per_face;
9243 new_quad->set_line_orientation(
9250 new_quads[0]->set_line_orientation(
9251 0, quad->line_orientation(0));
9252 new_quads[0]->set_line_orientation(
9253 2, quad->line_orientation(2));
9254 new_quads[1]->set_line_orientation(
9255 1, quad->line_orientation(1));
9256 new_quads[1]->set_line_orientation(
9257 2, quad->line_orientation(2));
9258 new_quads[2]->set_line_orientation(
9259 0, quad->line_orientation(0));
9260 new_quads[2]->set_line_orientation(
9261 3, quad->line_orientation(3));
9262 new_quads[3]->set_line_orientation(
9263 1, quad->line_orientation(1));
9264 new_quads[3]->set_line_orientation(
9265 3, quad->line_orientation(3));
9269 quad->clear_user_flag();
9280 cells_with_distorted_children;
9282 for (
unsigned int level = 0; level != triangulation.
levels.size() - 1;
9294 for (; hex != endh; ++hex)
9295 if (hex->refine_flag_set())
9303 hex->clear_refine_flag();
9304 hex->set_refinement_case(ref_case);
9315 unsigned int n_new_lines = 0;
9316 unsigned int n_new_quads = 0;
9317 unsigned int n_new_hexes = 0;
9348 new_lines(n_new_lines);
9349 for (
unsigned int i = 0; i < n_new_lines; ++i)
9352 triangulation.
faces->lines
9353 .template next_free_single_object<1>(triangulation);
9356 new_lines[i]->set_used_flag();
9357 new_lines[i]->clear_user_flag();
9358 new_lines[i]->clear_user_data();
9359 new_lines[i]->clear_children();
9361 new_lines[i]->set_boundary_id_internal(
9365 new_lines[i]->set_manifold_id(hex->manifold_id());
9372 new_quads(n_new_quads);
9373 for (
unsigned int i = 0; i < n_new_quads; ++i)
9376 triangulation.
faces->quads
9377 .template next_free_single_object<2>(triangulation);
9380 new_quads[i]->set_used_flag();
9381 new_quads[i]->clear_user_flag();
9382 new_quads[i]->clear_user_data();
9383 new_quads[i]->clear_children();
9385 new_quads[i]->set_boundary_id_internal(
9389 new_quads[i]->set_manifold_id(hex->manifold_id());
9392 for (
unsigned int j = 0;
9393 j < GeometryInfo<dim>::lines_per_face;
9395 new_quads[i]->set_line_orientation(
9405 new_hexes(n_new_hexes);
9406 for (
unsigned int i = 0; i < n_new_hexes; ++i)
9410 triangulation.
levels[level + 1]->cells.next_free_hex(
9411 triangulation, level + 1);
9415 new_hexes[i] = next_unused_hex;
9418 new_hexes[i]->set_used_flag();
9419 new_hexes[i]->clear_user_flag();
9420 new_hexes[i]->clear_user_data();
9421 new_hexes[i]->clear_children();
9424 new_hexes[i]->set_material_id(hex->material_id());
9425 new_hexes[i]->set_manifold_id(hex->manifold_id());
9426 new_hexes[i]->set_subdomain_id(subdomainid);
9433 new_hexes[i]->set_parent(hex->index());
9446 for (
const unsigned int f :
9448 new_hexes[i]->set_combined_face_orientation(
9452 for (
unsigned int i = 0; i < n_new_hexes / 2; ++i)
9453 hex->set_children(2 * i, new_hexes[2 * i]->index());
9460 const bool f_or[6] = {hex->face_orientation(0),
9461 hex->face_orientation(1),
9462 hex->face_orientation(2),
9463 hex->face_orientation(3),
9464 hex->face_orientation(4),
9465 hex->face_orientation(5)};
9468 const bool f_fl[6] = {hex->face_flip(0),
9476 const bool f_ro[6] = {hex->face_rotation(0),
9477 hex->face_rotation(1),
9478 hex->face_rotation(2),
9479 hex->face_rotation(3),
9480 hex->face_rotation(4),
9481 hex->face_rotation(5)};
9485 hex->combined_face_orientation(0),
9486 hex->combined_face_orientation(1),
9487 hex->combined_face_orientation(2),
9488 hex->combined_face_orientation(3),
9489 hex->combined_face_orientation(4),
9490 hex->combined_face_orientation(5)};
9500 const unsigned int child_at_origin[2][2][2] = {
9587 raw_line_iterator lines[4] = {
9588 hex->face(2)->child(0)->line(
9589 (hex->face(2)->refinement_case() ==
9593 hex->face(3)->child(0)->line(
9594 (hex->face(3)->refinement_case() ==
9598 hex->face(4)->child(0)->line(
9599 (hex->face(4)->refinement_case() ==
9603 hex->face(5)->child(0)->line(
9604 (hex->face(5)->refinement_case() ==
9610 unsigned int line_indices[4];
9611 for (
unsigned int i = 0; i < 4; ++i)
9612 line_indices[i] = lines[i]->
index();
9625 const unsigned int middle_vertices[2] = {
9626 hex->line(2)->child(0)->vertex_index(1),
9627 hex->line(7)->child(0)->vertex_index(1)};
9629 for (
unsigned int i = 0; i < 4; ++i)
9630 if (lines[i]->vertex_index(i % 2) ==
9631 middle_vertices[i % 2])
9632 line_orientation[i] =
9637 Assert(lines[i]->vertex_index((i + 1) % 2) ==
9638 middle_vertices[i % 2],
9640 line_orientation[i] =
9646 new_quads[0]->set_bounding_object_indices(
9652 new_quads[0]->set_line_orientation(
9653 0, line_orientation[0]);
9654 new_quads[0]->set_line_orientation(
9655 1, line_orientation[1]);
9656 new_quads[0]->set_line_orientation(
9657 2, line_orientation[2]);
9658 new_quads[0]->set_line_orientation(
9659 3, line_orientation[3]);
9692 const int quad_indices[11] = {
9693 new_quads[0]->index(),
9695 hex->face(0)->index(),
9697 hex->face(1)->index(),
9699 hex->face(2)->child_index(
9700 child_at_origin[hex->face(2)->refinement_case() -
9701 1][f_fl[2]][f_ro[2]]),
9702 hex->face(2)->child_index(
9704 child_at_origin[hex->face(2)->refinement_case() -
9705 1][f_fl[2]][f_ro[2]]),
9707 hex->face(3)->child_index(
9708 child_at_origin[hex->face(3)->refinement_case() -
9709 1][f_fl[3]][f_ro[3]]),
9710 hex->face(3)->child_index(
9712 child_at_origin[hex->face(3)->refinement_case() -
9713 1][f_fl[3]][f_ro[3]]),
9715 hex->face(4)->child_index(
9716 child_at_origin[hex->face(4)->refinement_case() -
9717 1][f_fl[4]][f_ro[4]]),
9718 hex->face(4)->child_index(
9720 child_at_origin[hex->face(4)->refinement_case() -
9721 1][f_fl[4]][f_ro[4]]),
9723 hex->face(5)->child_index(
9724 child_at_origin[hex->face(5)->refinement_case() -
9725 1][f_fl[5]][f_ro[5]]),
9726 hex->face(5)->child_index(
9728 child_at_origin[hex->face(5)->refinement_case() -
9729 1][f_fl[5]][f_ro[5]])
9733 new_hexes[0]->set_bounding_object_indices(
9740 new_hexes[1]->set_bounding_object_indices(
9815 raw_line_iterator lines[4] = {
9816 hex->face(0)->child(0)->line(
9817 (hex->face(0)->refinement_case() ==
9821 hex->face(1)->child(0)->line(
9822 (hex->face(1)->refinement_case() ==
9826 hex->face(4)->child(0)->line(
9827 (hex->face(4)->refinement_case() ==
9831 hex->face(5)->child(0)->line(
9832 (hex->face(5)->refinement_case() ==
9838 unsigned int line_indices[4];
9839 for (
unsigned int i = 0; i < 4; ++i)
9840 line_indices[i] = lines[i]->
index();
9853 const unsigned int middle_vertices[2] = {
9854 hex->line(0)->child(0)->vertex_index(1),
9855 hex->line(5)->child(0)->vertex_index(1)};
9857 for (
unsigned int i = 0; i < 4; ++i)
9858 if (lines[i]->vertex_index(i % 2) ==
9859 middle_vertices[i % 2])
9860 line_orientation[i] =
9865 Assert(lines[i]->vertex_index((i + 1) % 2) ==
9866 middle_vertices[i % 2],
9868 line_orientation[i] =
9874 new_quads[0]->set_bounding_object_indices(
9880 new_quads[0]->set_line_orientation(
9881 0, line_orientation[2]);
9882 new_quads[0]->set_line_orientation(
9883 1, line_orientation[3]);
9884 new_quads[0]->set_line_orientation(
9885 2, line_orientation[0]);
9886 new_quads[0]->set_line_orientation(
9887 3, line_orientation[1]);
9920 const int quad_indices[11] = {
9921 new_quads[0]->index(),
9923 hex->face(0)->child_index(
9924 child_at_origin[hex->face(0)->refinement_case() -
9925 1][f_fl[0]][f_ro[0]]),
9926 hex->face(0)->child_index(
9928 child_at_origin[hex->face(0)->refinement_case() -
9929 1][f_fl[0]][f_ro[0]]),
9931 hex->face(1)->child_index(
9932 child_at_origin[hex->face(1)->refinement_case() -
9933 1][f_fl[1]][f_ro[1]]),
9934 hex->face(1)->child_index(
9936 child_at_origin[hex->face(1)->refinement_case() -
9937 1][f_fl[1]][f_ro[1]]),
9939 hex->face(2)->index(),
9941 hex->face(3)->index(),
9943 hex->face(4)->child_index(
9944 child_at_origin[hex->face(4)->refinement_case() -
9945 1][f_fl[4]][f_ro[4]]),
9946 hex->face(4)->child_index(
9948 child_at_origin[hex->face(4)->refinement_case() -
9949 1][f_fl[4]][f_ro[4]]),
9951 hex->face(5)->child_index(
9952 child_at_origin[hex->face(5)->refinement_case() -
9953 1][f_fl[5]][f_ro[5]]),
9954 hex->face(5)->child_index(
9956 child_at_origin[hex->face(5)->refinement_case() -
9957 1][f_fl[5]][f_ro[5]])
9961 new_hexes[0]->set_bounding_object_indices(
9968 new_hexes[1]->set_bounding_object_indices(
10045 raw_line_iterator lines[4] = {
10046 hex->face(0)->child(0)->line(
10047 (hex->face(0)->refinement_case() ==
10051 hex->face(1)->child(0)->line(
10052 (hex->face(1)->refinement_case() ==
10056 hex->face(2)->child(0)->line(
10057 (hex->face(2)->refinement_case() ==
10061 hex->face(3)->child(0)->line(
10062 (hex->face(3)->refinement_case() ==
10068 unsigned int line_indices[4];
10069 for (
unsigned int i = 0; i < 4; ++i)
10070 line_indices[i] = lines[i]->
index();
10083 const unsigned int middle_vertices[2] = {
10084 middle_vertex_index<dim, spacedim>(hex->line(8)),
10085 middle_vertex_index<dim, spacedim>(hex->line(11))};
10087 for (
unsigned int i = 0; i < 4; ++i)
10088 if (lines[i]->vertex_index(i % 2) ==
10089 middle_vertices[i % 2])
10090 line_orientation[i] =
10095 Assert(lines[i]->vertex_index((i + 1) % 2) ==
10096 middle_vertices[i % 2],
10098 line_orientation[i] =
10104 new_quads[0]->set_bounding_object_indices(
10110 new_quads[0]->set_line_orientation(
10111 0, line_orientation[0]);
10112 new_quads[0]->set_line_orientation(
10113 1, line_orientation[1]);
10114 new_quads[0]->set_line_orientation(
10115 2, line_orientation[2]);
10116 new_quads[0]->set_line_orientation(
10117 3, line_orientation[3]);
10151 const int quad_indices[11] = {
10152 new_quads[0]->index(),
10154 hex->face(0)->child_index(
10155 child_at_origin[hex->face(0)->refinement_case() -
10156 1][f_fl[0]][f_ro[0]]),
10157 hex->face(0)->child_index(
10159 child_at_origin[hex->face(0)->refinement_case() -
10160 1][f_fl[0]][f_ro[0]]),
10162 hex->face(1)->child_index(
10163 child_at_origin[hex->face(1)->refinement_case() -
10164 1][f_fl[1]][f_ro[1]]),
10165 hex->face(1)->child_index(
10167 child_at_origin[hex->face(1)->refinement_case() -
10168 1][f_fl[1]][f_ro[1]]),
10170 hex->face(2)->child_index(
10171 child_at_origin[hex->face(2)->refinement_case() -
10172 1][f_fl[2]][f_ro[2]]),
10173 hex->face(2)->child_index(
10175 child_at_origin[hex->face(2)->refinement_case() -
10176 1][f_fl[2]][f_ro[2]]),
10178 hex->face(3)->child_index(
10179 child_at_origin[hex->face(3)->refinement_case() -
10180 1][f_fl[3]][f_ro[3]]),
10181 hex->face(3)->child_index(
10183 child_at_origin[hex->face(3)->refinement_case() -
10184 1][f_fl[3]][f_ro[3]]),
10186 hex->face(4)->index(),
10188 hex->face(5)->index()
10191 new_hexes[0]->set_bounding_object_indices(
10198 new_hexes[1]->set_bounding_object_indices(
10204 quad_indices[10]});
10230 new_lines[0]->set_bounding_object_indices(
10231 {middle_vertex_index<dim, spacedim>(hex->face(4)),
10232 middle_vertex_index<dim, spacedim>(hex->face(5))});
10300 spacedim>::raw_line_iterator lines[13] = {
10301 hex->face(0)->child(0)->line(
10302 (hex->face(0)->refinement_case() ==
10306 hex->face(1)->child(0)->line(
10307 (hex->face(1)->refinement_case() ==
10311 hex->face(2)->child(0)->line(
10312 (hex->face(2)->refinement_case() ==
10316 hex->face(3)->child(0)->line(
10317 (hex->face(3)->refinement_case() ==
10325 0, f_or[4], f_fl[4], f_ro[4]))
10328 1, f_or[4], f_fl[4], f_ro[4])),
10332 3, f_or[4], f_fl[4], f_ro[4]))
10335 0, f_or[4], f_fl[4], f_ro[4])),
10339 0, f_or[4], f_fl[4], f_ro[4]))
10342 3, f_or[4], f_fl[4], f_ro[4])),
10346 3, f_or[4], f_fl[4], f_ro[4]))
10349 2, f_or[4], f_fl[4], f_ro[4])),
10354 0, f_or[5], f_fl[5], f_ro[5]))
10357 1, f_or[5], f_fl[5], f_ro[5])),
10361 3, f_or[5], f_fl[5], f_ro[5]))
10364 0, f_or[5], f_fl[5], f_ro[5])),
10368 0, f_or[5], f_fl[5], f_ro[5]))
10371 3, f_or[5], f_fl[5], f_ro[5])),
10375 3, f_or[5], f_fl[5], f_ro[5]))
10378 2, f_or[5], f_fl[5], f_ro[5])),
10383 unsigned int line_indices[13];
10384 for (
unsigned int i = 0; i < 13; ++i)
10385 line_indices[i] = lines[i]->
index();
10396 const unsigned int middle_vertices[4] = {
10397 hex->line(0)->child(0)->vertex_index(1),
10398 hex->line(1)->child(0)->vertex_index(1),
10399 hex->line(2)->child(0)->vertex_index(1),
10400 hex->line(3)->child(0)->vertex_index(1),
10406 for (
unsigned int i = 0; i < 4; ++i)
10407 if (lines[i]->vertex_index(0) == middle_vertices[i])
10408 line_orientation[i] =
10413 Assert(lines[i]->vertex_index(1) ==
10414 middle_vertices[i],
10416 line_orientation[i] =
10426 for (
unsigned int i = 4; i < 12; ++i)
10427 if (lines[i]->vertex_index((i + 1) % 2) ==
10428 middle_vertex_index<dim, spacedim>(
10429 hex->face(3 + i / 4)))
10430 line_orientation[i] =
10435 Assert(lines[i]->vertex_index(i % 2) ==
10436 (middle_vertex_index<dim, spacedim>(
10437 hex->face(3 + i / 4))),
10439 line_orientation[i] =
10445 line_orientation[12] =
10471 new_quads[0]->set_bounding_object_indices(
10476 new_quads[1]->set_bounding_object_indices(
10481 new_quads[2]->set_bounding_object_indices(
10485 line_indices[12]});
10486 new_quads[3]->set_bounding_object_indices(
10492 new_quads[0]->set_line_orientation(
10493 0, line_orientation[2]);
10494 new_quads[0]->set_line_orientation(
10495 2, line_orientation[4]);
10496 new_quads[0]->set_line_orientation(
10497 3, line_orientation[8]);
10499 new_quads[1]->set_line_orientation(
10500 1, line_orientation[3]);
10501 new_quads[1]->set_line_orientation(
10502 2, line_orientation[5]);
10503 new_quads[1]->set_line_orientation(
10504 3, line_orientation[9]);
10506 new_quads[2]->set_line_orientation(
10507 0, line_orientation[6]);
10508 new_quads[2]->set_line_orientation(
10509 1, line_orientation[10]);
10510 new_quads[2]->set_line_orientation(
10511 2, line_orientation[0]);
10513 new_quads[3]->set_line_orientation(
10514 0, line_orientation[7]);
10515 new_quads[3]->set_line_orientation(
10516 1, line_orientation[11]);
10517 new_quads[3]->set_line_orientation(
10518 3, line_orientation[1]);
10551 const int quad_indices[20] = {
10552 new_quads[0]->index(),
10553 new_quads[1]->index(),
10554 new_quads[2]->index(),
10555 new_quads[3]->index(),
10557 hex->face(0)->child_index(
10558 child_at_origin[hex->face(0)->refinement_case() -
10559 1][f_fl[0]][f_ro[0]]),
10560 hex->face(0)->child_index(
10562 child_at_origin[hex->face(0)->refinement_case() -
10563 1][f_fl[0]][f_ro[0]]),
10565 hex->face(1)->child_index(
10566 child_at_origin[hex->face(1)->refinement_case() -
10567 1][f_fl[1]][f_ro[1]]),
10568 hex->face(1)->child_index(
10570 child_at_origin[hex->face(1)->refinement_case() -
10571 1][f_fl[1]][f_ro[1]]),
10573 hex->face(2)->child_index(
10574 child_at_origin[hex->face(2)->refinement_case() -
10575 1][f_fl[2]][f_ro[2]]),
10576 hex->face(2)->child_index(
10578 child_at_origin[hex->face(2)->refinement_case() -
10579 1][f_fl[2]][f_ro[2]]),
10581 hex->face(3)->child_index(
10582 child_at_origin[hex->face(3)->refinement_case() -
10583 1][f_fl[3]][f_ro[3]]),
10584 hex->face(3)->child_index(
10586 child_at_origin[hex->face(3)->refinement_case() -
10587 1][f_fl[3]][f_ro[3]]),
10589 hex->face(4)->isotropic_child_index(
10591 0, f_or[4], f_fl[4], f_ro[4])),
10592 hex->face(4)->isotropic_child_index(
10594 1, f_or[4], f_fl[4], f_ro[4])),
10595 hex->face(4)->isotropic_child_index(
10597 2, f_or[4], f_fl[4], f_ro[4])),
10598 hex->face(4)->isotropic_child_index(
10600 3, f_or[4], f_fl[4], f_ro[4])),
10602 hex->face(5)->isotropic_child_index(
10604 0, f_or[5], f_fl[5], f_ro[5])),
10605 hex->face(5)->isotropic_child_index(
10607 1, f_or[5], f_fl[5], f_ro[5])),
10608 hex->face(5)->isotropic_child_index(
10610 2, f_or[5], f_fl[5], f_ro[5])),
10611 hex->face(5)->isotropic_child_index(
10613 3, f_or[5], f_fl[5], f_ro[5]))};
10615 new_hexes[0]->set_bounding_object_indices(
10621 quad_indices[16]});
10622 new_hexes[1]->set_bounding_object_indices(
10628 quad_indices[17]});
10629 new_hexes[2]->set_bounding_object_indices(
10635 quad_indices[18]});
10636 new_hexes[3]->set_bounding_object_indices(
10642 quad_indices[19]});
10668 new_lines[0]->set_bounding_object_indices(
10669 {middle_vertex_index<dim, spacedim>(hex->face(2)),
10670 middle_vertex_index<dim, spacedim>(hex->face(3))});
10738 spacedim>::raw_line_iterator lines[13] = {
10739 hex->face(0)->child(0)->line(
10740 (hex->face(0)->refinement_case() ==
10744 hex->face(1)->child(0)->line(
10745 (hex->face(1)->refinement_case() ==
10749 hex->face(4)->child(0)->line(
10750 (hex->face(4)->refinement_case() ==
10754 hex->face(5)->child(0)->line(
10755 (hex->face(5)->refinement_case() ==
10763 0, f_or[2], f_fl[2], f_ro[2]))
10766 3, f_or[2], f_fl[2], f_ro[2])),
10770 3, f_or[2], f_fl[2], f_ro[2]))
10773 2, f_or[2], f_fl[2], f_ro[2])),
10777 0, f_or[2], f_fl[2], f_ro[2]))
10780 1, f_or[2], f_fl[2], f_ro[2])),
10784 3, f_or[2], f_fl[2], f_ro[2]))
10787 0, f_or[2], f_fl[2], f_ro[2])),
10792 0, f_or[3], f_fl[3], f_ro[3]))
10795 3, f_or[3], f_fl[3], f_ro[3])),
10799 3, f_or[3], f_fl[3], f_ro[3]))
10802 2, f_or[3], f_fl[3], f_ro[3])),
10806 0, f_or[3], f_fl[3], f_ro[3]))
10809 1, f_or[3], f_fl[3], f_ro[3])),
10813 3, f_or[3], f_fl[3], f_ro[3]))
10816 0, f_or[3], f_fl[3], f_ro[3])),
10821 unsigned int line_indices[13];
10822 for (
unsigned int i = 0; i < 13; ++i)
10823 line_indices[i] = lines[i]->
index();
10834 const unsigned int middle_vertices[4] = {
10835 hex->line(8)->child(0)->vertex_index(1),
10836 hex->line(9)->child(0)->vertex_index(1),
10837 hex->line(2)->child(0)->vertex_index(1),
10838 hex->line(6)->child(0)->vertex_index(1),
10843 for (
unsigned int i = 0; i < 4; ++i)
10844 if (lines[i]->vertex_index(0) == middle_vertices[i])
10845 line_orientation[i] =
10850 Assert(lines[i]->vertex_index(1) ==
10851 middle_vertices[i],
10853 line_orientation[i] =
10863 for (
unsigned int i = 4; i < 12; ++i)
10864 if (lines[i]->vertex_index((i + 1) % 2) ==
10865 middle_vertex_index<dim, spacedim>(
10866 hex->face(1 + i / 4)))
10867 line_orientation[i] =
10873 Assert(lines[i]->vertex_index(i % 2) ==
10874 (middle_vertex_index<dim, spacedim>(
10875 hex->face(1 + i / 4))),
10877 line_orientation[i] =
10883 line_orientation[12] =
10910 new_quads[0]->set_bounding_object_indices(
10914 line_indices[10]});
10915 new_quads[1]->set_bounding_object_indices(
10919 line_indices[11]});
10920 new_quads[2]->set_bounding_object_indices(
10924 line_indices[12]});
10925 new_quads[3]->set_bounding_object_indices(
10931 new_quads[0]->set_line_orientation(
10932 0, line_orientation[0]);
10933 new_quads[0]->set_line_orientation(
10934 2, line_orientation[6]);
10935 new_quads[0]->set_line_orientation(
10936 3, line_orientation[10]);
10938 new_quads[1]->set_line_orientation(
10939 1, line_orientation[1]);
10940 new_quads[1]->set_line_orientation(
10941 2, line_orientation[7]);
10942 new_quads[1]->set_line_orientation(
10943 3, line_orientation[11]);
10945 new_quads[2]->set_line_orientation(
10946 0, line_orientation[4]);
10947 new_quads[2]->set_line_orientation(
10948 1, line_orientation[8]);
10949 new_quads[2]->set_line_orientation(
10950 2, line_orientation[2]);
10952 new_quads[3]->set_line_orientation(
10953 0, line_orientation[5]);
10954 new_quads[3]->set_line_orientation(
10955 1, line_orientation[9]);
10956 new_quads[3]->set_line_orientation(
10957 3, line_orientation[3]);
10991 const int quad_indices[20] = {
10992 new_quads[0]->index(),
10993 new_quads[1]->index(),
10994 new_quads[2]->index(),
10995 new_quads[3]->index(),
10997 hex->face(0)->child_index(
10998 child_at_origin[hex->face(0)->refinement_case() -
10999 1][f_fl[0]][f_ro[0]]),
11000 hex->face(0)->child_index(
11002 child_at_origin[hex->face(0)->refinement_case() -
11003 1][f_fl[0]][f_ro[0]]),
11005 hex->face(1)->child_index(
11006 child_at_origin[hex->face(1)->refinement_case() -
11007 1][f_fl[1]][f_ro[1]]),
11008 hex->face(1)->child_index(
11010 child_at_origin[hex->face(1)->refinement_case() -
11011 1][f_fl[1]][f_ro[1]]),
11013 hex->face(2)->isotropic_child_index(
11015 0, f_or[2], f_fl[2], f_ro[2])),
11016 hex->face(2)->isotropic_child_index(
11018 1, f_or[2], f_fl[2], f_ro[2])),
11019 hex->face(2)->isotropic_child_index(
11021 2, f_or[2], f_fl[2], f_ro[2])),
11022 hex->face(2)->isotropic_child_index(
11024 3, f_or[2], f_fl[2], f_ro[2])),
11026 hex->face(3)->isotropic_child_index(
11028 0, f_or[3], f_fl[3], f_ro[3])),
11029 hex->face(3)->isotropic_child_index(
11031 1, f_or[3], f_fl[3], f_ro[3])),
11032 hex->face(3)->isotropic_child_index(
11034 2, f_or[3], f_fl[3], f_ro[3])),
11035 hex->face(3)->isotropic_child_index(
11037 3, f_or[3], f_fl[3], f_ro[3])),
11039 hex->face(4)->child_index(
11040 child_at_origin[hex->face(4)->refinement_case() -
11041 1][f_fl[4]][f_ro[4]]),
11042 hex->face(4)->child_index(
11044 child_at_origin[hex->face(4)->refinement_case() -
11045 1][f_fl[4]][f_ro[4]]),
11047 hex->face(5)->child_index(
11048 child_at_origin[hex->face(5)->refinement_case() -
11049 1][f_fl[5]][f_ro[5]]),
11050 hex->face(5)->child_index(
11052 child_at_origin[hex->face(5)->refinement_case() -
11053 1][f_fl[5]][f_ro[5]])};
11064 new_hexes[0]->set_bounding_object_indices(
11071 new_hexes[1]->set_bounding_object_indices(
11077 quad_indices[18]});
11078 new_hexes[2]->set_bounding_object_indices(
11085 new_hexes[3]->set_bounding_object_indices(
11091 quad_indices[19]});
11118 new_lines[0]->set_bounding_object_indices(
11120 {middle_vertex_index<dim, spacedim>(hex->face(0)),
11121 middle_vertex_index<dim, spacedim>(hex->face(1))});
11190 spacedim>::raw_line_iterator lines[13] = {
11191 hex->face(2)->child(0)->line(
11192 (hex->face(2)->refinement_case() ==
11196 hex->face(3)->child(0)->line(
11197 (hex->face(3)->refinement_case() ==
11201 hex->face(4)->child(0)->line(
11202 (hex->face(4)->refinement_case() ==
11206 hex->face(5)->child(0)->line(
11207 (hex->face(5)->refinement_case() ==
11215 0, f_or[0], f_fl[0], f_ro[0]))
11218 1, f_or[0], f_fl[0], f_ro[0])),
11222 3, f_or[0], f_fl[0], f_ro[0]))
11225 0, f_or[0], f_fl[0], f_ro[0])),
11229 0, f_or[0], f_fl[0], f_ro[0]))
11232 3, f_or[0], f_fl[0], f_ro[0])),
11236 3, f_or[0], f_fl[0], f_ro[0]))
11239 2, f_or[0], f_fl[0], f_ro[0])),
11244 0, f_or[1], f_fl[1], f_ro[1]))
11247 1, f_or[1], f_fl[1], f_ro[1])),
11251 3, f_or[1], f_fl[1], f_ro[1]))
11254 0, f_or[1], f_fl[1], f_ro[1])),
11258 0, f_or[1], f_fl[1], f_ro[1]))
11261 3, f_or[1], f_fl[1], f_ro[1])),
11265 3, f_or[1], f_fl[1], f_ro[1]))
11268 2, f_or[1], f_fl[1], f_ro[1])),
11273 unsigned int line_indices[13];
11275 for (
unsigned int i = 0; i < 13; ++i)
11276 line_indices[i] = lines[i]->
index();
11287 const unsigned int middle_vertices[4] = {
11288 hex->line(8)->child(0)->vertex_index(1),
11289 hex->line(10)->child(0)->vertex_index(1),
11290 hex->line(0)->child(0)->vertex_index(1),
11291 hex->line(4)->child(0)->vertex_index(1),
11296 for (
unsigned int i = 0; i < 4; ++i)
11297 if (lines[i]->vertex_index(0) == middle_vertices[i])
11298 line_orientation[i] =
11303 Assert(lines[i]->vertex_index(1) ==
11304 middle_vertices[i],
11306 line_orientation[i] =
11316 for (
unsigned int i = 4; i < 12; ++i)
11317 if (lines[i]->vertex_index((i + 1) % 2) ==
11318 middle_vertex_index<dim, spacedim>(
11319 hex->face(i / 4 - 1)))
11320 line_orientation[i] =
11325 Assert(lines[i]->vertex_index(i % 2) ==
11326 (middle_vertex_index<dim, spacedim>(
11327 hex->face(i / 4 - 1))),
11329 line_orientation[i] =
11334 line_orientation[12] =
11355 new_quads[0]->set_bounding_object_indices(
11359 line_indices[12]});
11360 new_quads[1]->set_bounding_object_indices(
11365 new_quads[2]->set_bounding_object_indices(
11370 new_quads[3]->set_bounding_object_indices(
11376 new_quads[0]->set_line_orientation(
11377 0, line_orientation[6]);
11378 new_quads[0]->set_line_orientation(
11379 1, line_orientation[10]);
11380 new_quads[0]->set_line_orientation(
11381 2, line_orientation[0]);
11383 new_quads[1]->set_line_orientation(
11384 0, line_orientation[7]);
11385 new_quads[1]->set_line_orientation(
11386 1, line_orientation[11]);
11387 new_quads[1]->set_line_orientation(
11388 3, line_orientation[1]);
11390 new_quads[2]->set_line_orientation(
11391 0, line_orientation[2]);
11392 new_quads[2]->set_line_orientation(
11393 2, line_orientation[4]);
11394 new_quads[2]->set_line_orientation(
11395 3, line_orientation[8]);
11397 new_quads[3]->set_line_orientation(
11398 1, line_orientation[3]);
11399 new_quads[3]->set_line_orientation(
11400 2, line_orientation[5]);
11401 new_quads[3]->set_line_orientation(
11402 3, line_orientation[9]);
11436 const int quad_indices[20] = {
11437 new_quads[0]->index(),
11438 new_quads[1]->index(),
11439 new_quads[2]->index(),
11440 new_quads[3]->index(),
11442 hex->face(0)->isotropic_child_index(
11444 0, f_or[0], f_fl[0], f_ro[0])),
11445 hex->face(0)->isotropic_child_index(
11447 1, f_or[0], f_fl[0], f_ro[0])),
11448 hex->face(0)->isotropic_child_index(
11450 2, f_or[0], f_fl[0], f_ro[0])),
11451 hex->face(0)->isotropic_child_index(
11453 3, f_or[0], f_fl[0], f_ro[0])),
11455 hex->face(1)->isotropic_child_index(
11457 0, f_or[1], f_fl[1], f_ro[1])),
11458 hex->face(1)->isotropic_child_index(
11460 1, f_or[1], f_fl[1], f_ro[1])),
11461 hex->face(1)->isotropic_child_index(
11463 2, f_or[1], f_fl[1], f_ro[1])),
11464 hex->face(1)->isotropic_child_index(
11466 3, f_or[1], f_fl[1], f_ro[1])),
11468 hex->face(2)->child_index(
11469 child_at_origin[hex->face(2)->refinement_case() -
11470 1][f_fl[2]][f_ro[2]]),
11471 hex->face(2)->child_index(
11473 child_at_origin[hex->face(2)->refinement_case() -
11474 1][f_fl[2]][f_ro[2]]),
11476 hex->face(3)->child_index(
11477 child_at_origin[hex->face(3)->refinement_case() -
11478 1][f_fl[3]][f_ro[3]]),
11479 hex->face(3)->child_index(
11481 child_at_origin[hex->face(3)->refinement_case() -
11482 1][f_fl[3]][f_ro[3]]),
11484 hex->face(4)->child_index(
11485 child_at_origin[hex->face(4)->refinement_case() -
11486 1][f_fl[4]][f_ro[4]]),
11487 hex->face(4)->child_index(
11489 child_at_origin[hex->face(4)->refinement_case() -
11490 1][f_fl[4]][f_ro[4]]),
11492 hex->face(5)->child_index(
11493 child_at_origin[hex->face(5)->refinement_case() -
11494 1][f_fl[5]][f_ro[5]]),
11495 hex->face(5)->child_index(
11497 child_at_origin[hex->face(5)->refinement_case() -
11498 1][f_fl[5]][f_ro[5]])};
11500 new_hexes[0]->set_bounding_object_indices(
11507 new_hexes[1]->set_bounding_object_indices(
11514 new_hexes[2]->set_bounding_object_indices(
11520 quad_indices[18]});
11521 new_hexes[3]->set_bounding_object_indices(
11527 quad_indices[19]});
11559 ++next_unused_vertex;
11561 next_unused_vertex < triangulation.
vertices.size(),
11563 "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."));
11573 triangulation.
vertices[next_unused_vertex] =
11574 hex->center(
true,
true);
11595 const unsigned int vertex_indices[7] = {
11596 middle_vertex_index<dim, spacedim>(hex->face(0)),
11597 middle_vertex_index<dim, spacedim>(hex->face(1)),
11598 middle_vertex_index<dim, spacedim>(hex->face(2)),
11599 middle_vertex_index<dim, spacedim>(hex->face(3)),
11600 middle_vertex_index<dim, spacedim>(hex->face(4)),
11601 middle_vertex_index<dim, spacedim>(hex->face(5)),
11602 next_unused_vertex};
11604 new_lines[0]->set_bounding_object_indices(
11605 {vertex_indices[2], vertex_indices[6]});
11606 new_lines[1]->set_bounding_object_indices(
11607 {vertex_indices[6], vertex_indices[3]});
11608 new_lines[2]->set_bounding_object_indices(
11609 {vertex_indices[0], vertex_indices[6]});
11610 new_lines[3]->set_bounding_object_indices(
11611 {vertex_indices[6], vertex_indices[1]});
11612 new_lines[4]->set_bounding_object_indices(
11613 {vertex_indices[4], vertex_indices[6]});
11614 new_lines[5]->set_bounding_object_indices(
11615 {vertex_indices[6], vertex_indices[5]});
11684 spacedim>::raw_line_iterator lines[30] = {
11688 0, f_or[0], f_fl[0], f_ro[0]))
11691 1, f_or[0], f_fl[0], f_ro[0])),
11695 3, f_or[0], f_fl[0], f_ro[0]))
11698 0, f_or[0], f_fl[0], f_ro[0])),
11702 0, f_or[0], f_fl[0], f_ro[0]))
11705 3, f_or[0], f_fl[0], f_ro[0])),
11709 3, f_or[0], f_fl[0], f_ro[0]))
11712 2, f_or[0], f_fl[0], f_ro[0])),
11717 0, f_or[1], f_fl[1], f_ro[1]))
11720 1, f_or[1], f_fl[1], f_ro[1])),
11724 3, f_or[1], f_fl[1], f_ro[1]))
11727 0, f_or[1], f_fl[1], f_ro[1])),
11731 0, f_or[1], f_fl[1], f_ro[1]))
11734 3, f_or[1], f_fl[1], f_ro[1])),
11738 3, f_or[1], f_fl[1], f_ro[1]))
11741 2, f_or[1], f_fl[1], f_ro[1])),
11746 0, f_or[2], f_fl[2], f_ro[2]))
11749 1, f_or[2], f_fl[2], f_ro[2])),
11753 3, f_or[2], f_fl[2], f_ro[2]))
11756 0, f_or[2], f_fl[2], f_ro[2])),
11760 0, f_or[2], f_fl[2], f_ro[2]))
11763 3, f_or[2], f_fl[2], f_ro[2])),
11767 3, f_or[2], f_fl[2], f_ro[2]))
11770 2, f_or[2], f_fl[2], f_ro[2])),
11775 0, f_or[3], f_fl[3], f_ro[3]))
11778 1, f_or[3], f_fl[3], f_ro[3])),
11782 3, f_or[3], f_fl[3], f_ro[3]))
11785 0, f_or[3], f_fl[3], f_ro[3])),
11789 0, f_or[3], f_fl[3], f_ro[3]))
11792 3, f_or[3], f_fl[3], f_ro[3])),
11796 3, f_or[3], f_fl[3], f_ro[3]))
11799 2, f_or[3], f_fl[3], f_ro[3])),
11804 0, f_or[4], f_fl[4], f_ro[4]))
11807 1, f_or[4], f_fl[4], f_ro[4])),
11811 3, f_or[4], f_fl[4], f_ro[4]))
11814 0, f_or[4], f_fl[4], f_ro[4])),
11818 0, f_or[4], f_fl[4], f_ro[4]))
11821 3, f_or[4], f_fl[4], f_ro[4])),
11825 3, f_or[4], f_fl[4], f_ro[4]))
11828 2, f_or[4], f_fl[4], f_ro[4])),
11833 0, f_or[5], f_fl[5], f_ro[5]))
11836 1, f_or[5], f_fl[5], f_ro[5])),
11840 3, f_or[5], f_fl[5], f_ro[5]))
11843 0, f_or[5], f_fl[5], f_ro[5])),
11847 0, f_or[5], f_fl[5], f_ro[5]))
11850 3, f_or[5], f_fl[5], f_ro[5])),
11854 3, f_or[5], f_fl[5], f_ro[5]))
11857 2, f_or[5], f_fl[5], f_ro[5])),
11867 unsigned int line_indices[30];
11868 for (
unsigned int i = 0; i < 30; ++i)
11869 line_indices[i] = lines[i]->
index();
11884 for (
unsigned int i = 0; i < 24; ++i)
11885 if (lines[i]->vertex_index((i + 1) % 2) ==
11886 vertex_indices[i / 4])
11887 line_orientation[i] =
11893 Assert(lines[i]->vertex_index(i % 2) ==
11894 vertex_indices[i / 4],
11896 line_orientation[i] =
11902 for (
unsigned int i = 24; i < 30; ++i)
11903 line_orientation[i] =
11936 new_quads[0]->set_bounding_object_indices(
11940 line_indices[24]});
11941 new_quads[1]->set_bounding_object_indices(
11945 line_indices[25]});
11946 new_quads[2]->set_bounding_object_indices(
11950 line_indices[20]});
11951 new_quads[3]->set_bounding_object_indices(
11955 line_indices[21]});
11956 new_quads[4]->set_bounding_object_indices(
11960 line_indices[28]});
11961 new_quads[5]->set_bounding_object_indices(
11965 line_indices[29]});
11966 new_quads[6]->set_bounding_object_indices(
11971 new_quads[7]->set_bounding_object_indices(
11976 new_quads[8]->set_bounding_object_indices(
11980 line_indices[26]});
11981 new_quads[9]->set_bounding_object_indices(
11985 line_indices[27]});
11986 new_quads[10]->set_bounding_object_indices(
11990 line_indices[12]});
11991 new_quads[11]->set_bounding_object_indices(
11995 line_indices[13]});
12000 new_quads[0]->set_line_orientation(
12001 0, line_orientation[10]);
12002 new_quads[0]->set_line_orientation(
12003 2, line_orientation[16]);
12005 new_quads[1]->set_line_orientation(
12006 1, line_orientation[14]);
12007 new_quads[1]->set_line_orientation(
12008 2, line_orientation[17]);
12010 new_quads[2]->set_line_orientation(
12011 0, line_orientation[11]);
12012 new_quads[2]->set_line_orientation(
12013 3, line_orientation[20]);
12015 new_quads[3]->set_line_orientation(
12016 1, line_orientation[15]);
12017 new_quads[3]->set_line_orientation(
12018 3, line_orientation[21]);
12020 new_quads[4]->set_line_orientation(
12021 0, line_orientation[18]);
12022 new_quads[4]->set_line_orientation(
12023 2, line_orientation[0]);
12025 new_quads[5]->set_line_orientation(
12026 1, line_orientation[22]);
12027 new_quads[5]->set_line_orientation(
12028 2, line_orientation[1]);
12030 new_quads[6]->set_line_orientation(
12031 0, line_orientation[19]);
12032 new_quads[6]->set_line_orientation(
12033 3, line_orientation[4]);
12035 new_quads[7]->set_line_orientation(
12036 1, line_orientation[23]);
12037 new_quads[7]->set_line_orientation(
12038 3, line_orientation[5]);
12040 new_quads[8]->set_line_orientation(
12041 0, line_orientation[2]);
12042 new_quads[8]->set_line_orientation(
12043 2, line_orientation[8]);
12045 new_quads[9]->set_line_orientation(
12046 1, line_orientation[6]);
12047 new_quads[9]->set_line_orientation(
12048 2, line_orientation[9]);
12050 new_quads[10]->set_line_orientation(
12051 0, line_orientation[3]);
12052 new_quads[10]->set_line_orientation(
12053 3, line_orientation[12]);
12055 new_quads[11]->set_line_orientation(
12056 1, line_orientation[7]);
12057 new_quads[11]->set_line_orientation(
12058 3, line_orientation[13]);
12099 const int quad_indices[36] = {
12100 new_quads[0]->index(),
12101 new_quads[1]->index(),
12102 new_quads[2]->index(),
12103 new_quads[3]->index(),
12104 new_quads[4]->index(),
12105 new_quads[5]->index(),
12106 new_quads[6]->index(),
12107 new_quads[7]->index(),
12108 new_quads[8]->index(),
12109 new_quads[9]->index(),
12110 new_quads[10]->index(),
12111 new_quads[11]->index(),
12113 hex->face(0)->isotropic_child_index(
12115 0, f_or[0], f_fl[0], f_ro[0])),
12116 hex->face(0)->isotropic_child_index(
12118 1, f_or[0], f_fl[0], f_ro[0])),
12119 hex->face(0)->isotropic_child_index(
12121 2, f_or[0], f_fl[0], f_ro[0])),
12122 hex->face(0)->isotropic_child_index(
12124 3, f_or[0], f_fl[0], f_ro[0])),
12126 hex->face(1)->isotropic_child_index(
12128 0, f_or[1], f_fl[1], f_ro[1])),
12129 hex->face(1)->isotropic_child_index(
12131 1, f_or[1], f_fl[1], f_ro[1])),
12132 hex->face(1)->isotropic_child_index(
12134 2, f_or[1], f_fl[1], f_ro[1])),
12135 hex->face(1)->isotropic_child_index(
12137 3, f_or[1], f_fl[1], f_ro[1])),
12139 hex->face(2)->isotropic_child_index(
12141 0, f_or[2], f_fl[2], f_ro[2])),
12142 hex->face(2)->isotropic_child_index(
12144 1, f_or[2], f_fl[2], f_ro[2])),
12145 hex->face(2)->isotropic_child_index(
12147 2, f_or[2], f_fl[2], f_ro[2])),
12148 hex->face(2)->isotropic_child_index(
12150 3, f_or[2], f_fl[2], f_ro[2])),
12152 hex->face(3)->isotropic_child_index(
12154 0, f_or[3], f_fl[3], f_ro[3])),
12155 hex->face(3)->isotropic_child_index(
12157 1, f_or[3], f_fl[3], f_ro[3])),
12158 hex->face(3)->isotropic_child_index(
12160 2, f_or[3], f_fl[3], f_ro[3])),
12161 hex->face(3)->isotropic_child_index(
12163 3, f_or[3], f_fl[3], f_ro[3])),
12165 hex->face(4)->isotropic_child_index(
12167 0, f_or[4], f_fl[4], f_ro[4])),
12168 hex->face(4)->isotropic_child_index(
12170 1, f_or[4], f_fl[4], f_ro[4])),
12171 hex->face(4)->isotropic_child_index(
12173 2, f_or[4], f_fl[4], f_ro[4])),
12174 hex->face(4)->isotropic_child_index(
12176 3, f_or[4], f_fl[4], f_ro[4])),
12178 hex->face(5)->isotropic_child_index(
12180 0, f_or[5], f_fl[5], f_ro[5])),
12181 hex->face(5)->isotropic_child_index(
12183 1, f_or[5], f_fl[5], f_ro[5])),
12184 hex->face(5)->isotropic_child_index(
12186 2, f_or[5], f_fl[5], f_ro[5])),
12187 hex->face(5)->isotropic_child_index(
12189 3, f_or[5], f_fl[5], f_ro[5]))};
12192 new_hexes[0]->set_bounding_object_indices(
12199 new_hexes[1]->set_bounding_object_indices(
12206 new_hexes[2]->set_bounding_object_indices(
12212 quad_indices[10]});
12213 new_hexes[3]->set_bounding_object_indices(
12219 quad_indices[11]});
12222 new_hexes[4]->set_bounding_object_indices(
12228 quad_indices[32]});
12229 new_hexes[5]->set_bounding_object_indices(
12235 quad_indices[33]});
12236 new_hexes[6]->set_bounding_object_indices(
12242 quad_indices[34]});
12243 new_hexes[7]->set_bounding_object_indices(
12249 quad_indices[35]});
12283 for (
unsigned int s = 0;
12290 const unsigned int current_child =
12299 ref_case, f, f_or[f], f_fl[f], f_ro[f]));
12300 new_hexes[current_child]->set_combined_face_orientation(
12306 if (check_for_distorted_cells &&
12307 has_distorted_children<dim, spacedim>(hex))
12324 triangulation.
faces->quads.clear_user_data();
12327 return cells_with_distorted_children;
12343 template <
int spacedim>
12350 template <
int dim,
int spacedim>
12357 if (spacedim > dim)
12369 for (
const unsigned int face_no :
12371 if (cell->
face(face_no)->at_boundary())
12382 spacedim>::face_iterator
12383 face = cell->
face(face_no);
12402 const double dist =
12403 std::fabs(new_unit[face_no / 2] - face_no % 2);
12408 const double allowed = 0.25;
12410 if (dist > allowed)
12430 template <
int dim,
int spacedim>
12435 ExcMessage(
"Wrong function called -- there should "
12436 "be a specialization."));
12440 template <
int spacedim>
12445 const unsigned int dim = 3;
12446 using raw_line_iterator =
12451 bool mesh_changed =
false;
12455 mesh_changed =
false;
12471 const std::array<unsigned int, 12> line_indices =
12472 TriaAccessorImplementation::Implementation::
12473 get_line_indices_of_cell(*cell);
12474 for (
unsigned int l = 0; l < cell->
n_lines(); ++l)
12479 raw_line_iterator line(&triangulation,
12483 line->set_user_flag();
12489 const std::array<unsigned int, 12> line_indices =
12490 TriaAccessorImplementation::Implementation::
12491 get_line_indices_of_cell(*cell);
12492 for (
unsigned int l = 0; l < cell->
n_lines(); ++l)
12497 raw_line_iterator line(&triangulation,
12501 line->set_user_flag();
12515 cell != triangulation.
end();
12518 const std::array<unsigned int, 12> line_indices =
12519 TriaAccessorImplementation::Implementation::
12520 get_line_indices_of_cell(*cell);
12521 for (
unsigned int l = 0; l < cell->
n_lines(); ++l)
12523 raw_line_iterator line(&triangulation, 0, line_indices[l]);
12524 if (line->has_children())
12533 bool offending_line_found =
false;
12535 for (
unsigned int c = 0; c < 2; ++c)
12537 Assert(line->child(c)->has_children() ==
false,
12540 if (line->child(c)->user_flag_set() &&
12553 allow_anisotropic_smoothing)
12558 for (
unsigned int k = 0; k < cell->
n_lines();
12564 raw_line_iterator(&triangulation,
12570 offending_line_found =
true;
12576 for (
unsigned int k = 0; k < cell->
n_lines();
12580 raw_line_iterator(&triangulation,
12583 if (!line->has_children() &&
12585 line_refinement_case(
12588 line->set_user_flag();
12595 if (offending_line_found)
12597 mesh_changed =
true;
12616 triangulation.
last();
12617 cell != triangulation.
end();
12621 const std::array<unsigned int, 12> line_indices =
12622 TriaAccessorImplementation::Implementation::
12623 get_line_indices_of_cell(*cell);
12624 for (
unsigned int l = 0; l < cell->
n_lines(); ++l)
12626 raw_line_iterator line(&triangulation,
12629 if (line->has_children() &&
12630 (line->child(0)->user_flag_set() ||
12631 line->child(1)->user_flag_set()))
12633 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
12636 for (
unsigned int k = 0; k < cell->
n_lines(); ++k)
12642 raw_line_iterator(&triangulation,
12646 mesh_changed =
true;
12652 while (mesh_changed ==
true);
12663 template <
int dim,
int spacedim>
12683 const unsigned int n_subfaces =
12688 for (
unsigned int c = 0; c < n_subfaces; ++c)
12691 child = cell->
child(
12695 child_neighbor = child->
neighbor(n);
12735 template <
int spacedim>
12740 template <
int dim,
int spacedim>
12743 std::vector<std::pair<unsigned int, unsigned int>> adjacent_cells(
12745 {numbers::invalid_unsigned_int, numbers::invalid_unsigned_int});
12747 const auto set_entry = [&](
const auto &face_index,
const auto &cell) {
12748 const std::pair<unsigned int, unsigned int> cell_pair = {
12750 unsigned int index;
12752 if (adjacent_cells[2 * face_index].first ==
12754 adjacent_cells[2 * face_index].second ==
12757 index = 2 * face_index + 0;
12761 Assert(((adjacent_cells[2 * face_index + 1].first ==
12763 (adjacent_cells[2 * face_index + 1].second ==
12766 index = 2 * face_index + 1;
12769 adjacent_cells[
index] = cell_pair;
12772 const auto get_entry =
12773 [&](
const auto &face_index,
12775 auto test = adjacent_cells[2 * face_index];
12777 if (test == std::pair<unsigned int, unsigned int>(cell->
level(),
12779 test = adjacent_cells[2 * face_index + 1];
12793 set_entry(face->index(), cell);
12795 if (cell->
is_active() && face->has_children())
12796 for (
unsigned int c = 0; c < face->n_children(); ++c)
12797 set_entry(face->child(c)->index(), cell);
12805 template <
int dim,
int spacedim>
12810 std::vector<unsigned int> & ,
12811 std::vector<unsigned int> & )
12816 template <
int dim,
int spacedim>
12819 const bool check_for_distorted_cells)
12822 triangulation, check_for_distorted_cells);
12825 template <
int dim,
int spacedim>
12833 template <
int dim,
int spacedim>
12841 template <
int dim,
int spacedim>
12856template <
int dim,
int spacedim>
12862template <
int dim,
int spacedim>
12865 const MeshSmoothing smooth_grid,
12866 const bool check_for_distorted_cells)
12867 : cell_attached_data({0, 0, {}, {}})
12868 , smooth_grid(smooth_grid)
12869 , default_flat_manifold(std::make_unique<
const FlatManifold<dim, spacedim>>())
12870 , anisotropic_refinement(
false)
12871 , check_for_distorted_cells(check_for_distorted_cells)
12875 vertex_to_boundary_id_map_1d =
12876 std::make_unique<std::map<unsigned int, types::boundary_id>>();
12877 vertex_to_manifold_id_map_1d =
12878 std::make_unique<std::map<unsigned int, types::manifold_id>>();
12882 signals.create.connect(signals.any_change);
12883 signals.post_refinement.connect(signals.any_change);
12884 signals.clear.connect(signals.any_change);
12885 signals.mesh_movement.connect(signals.any_change);
12890template <
int dim,
int spacedim>
12894 : EnableObserverPointer(std::move(tria))
12895 , cell_attached_data(std::move(tria.cell_attached_data))
12896 , smooth_grid(tria.smooth_grid)
12897 , reference_cells(std::move(tria.reference_cells))
12898 , strides(tria.strides)
12899 , periodic_face_pairs_level_0(std::move(tria.periodic_face_pairs_level_0))
12900 , periodic_face_map(std::move(tria.periodic_face_map))
12901 , levels(std::move(tria.levels))
12902 , faces(std::move(tria.faces))
12903 , vertices(std::move(tria.vertices))
12904 , vertices_used(std::move(tria.vertices_used))
12905 , manifolds(std::move(tria.manifolds))
12906 , default_flat_manifold(std::make_unique<
const FlatManifold<dim, spacedim>>())
12907 , anisotropic_refinement(tria.anisotropic_refinement)
12908 , check_for_distorted_cells(tria.check_for_distorted_cells)
12909 , number_cache(std::move(tria.number_cache))
12910 , vertex_to_boundary_id_map_1d(std::move(tria.vertex_to_boundary_id_map_1d))
12911 , vertex_to_manifold_id_map_1d(std::move(tria.vertex_to_manifold_id_map_1d))
12913 tria.number_cache = internal::TriangulationImplementation::NumberCache<dim>();
12917 this->policy = tria.policy->clone();
12921template <
int dim,
int spacedim>
12928 cell_attached_data = std::move(tria.cell_attached_data);
12929 smooth_grid = tria.smooth_grid;
12930 reference_cells = std::move(tria.reference_cells);
12931 strides = tria.strides;
12932 periodic_face_pairs_level_0 = std::move(tria.periodic_face_pairs_level_0);
12933 periodic_face_map = std::move(tria.periodic_face_map);
12934 levels = std::move(tria.levels);
12935 faces = std::move(tria.faces);
12936 vertices = std::move(tria.vertices);
12937 vertices_used = std::move(tria.vertices_used);
12938 manifolds = std::move(tria.manifolds);
12939 anisotropic_refinement = tria.anisotropic_refinement;
12940 number_cache = tria.number_cache;
12941 vertex_to_boundary_id_map_1d = std::move(tria.vertex_to_boundary_id_map_1d);
12942 vertex_to_manifold_id_map_1d = std::move(tria.vertex_to_manifold_id_map_1d);
12944 tria.number_cache = internal::TriangulationImplementation::NumberCache<dim>();
12948 this->policy = tria.policy->clone();
12955template <
int dim,
int spacedim>
12973 AssertNothrow((dim == 1) || (vertex_to_boundary_id_map_1d ==
nullptr),
12978 AssertNothrow((dim == 1) || (vertex_to_manifold_id_map_1d ==
nullptr),
12984template <
int dim,
int spacedim>
12992 clear_despite_subscriptions();
12993 periodic_face_pairs_level_0.clear();
12994 periodic_face_map.clear();
12995 reference_cells.clear();
12997 cell_attached_data = {0, 0, {}, {}};
12998 data_serializer.clear();
13001template <
int dim,
int spacedim>
13005 return MPI_COMM_SELF;
13010template <
int dim,
int spacedim>
13014 return get_mpi_communicator();
13019template <
int dim,
int spacedim>
13024 return number_cache.active_cell_index_partitioner;
13029template <
int dim,
int spacedim>
13036 return number_cache.level_cell_index_partitioners[level];
13041template <
int dim,
int spacedim>
13044 const MeshSmoothing mesh_smoothing)
13046 smooth_grid = mesh_smoothing;
13051template <
int dim,
int spacedim>
13056 return smooth_grid;
13061template <
int dim,
int spacedim>
13065 const Manifold<dim, spacedim> &manifold_object)
13069 manifolds[m_number] = manifold_object.
clone();
13074template <
int dim,
int spacedim>
13082 manifolds[m_number] = default_flat_manifold->clone();
13086template <
int dim,
int spacedim>
13090 for (
auto &m : manifolds)
13091 m.second = default_flat_manifold->clone();
13095template <
int dim,
int spacedim>
13103 "Error: set_all_manifold_ids() can not be called on an empty Triangulation."));
13105 for (
const auto &cell : this->active_cell_iterators())
13110template <
int dim,
int spacedim>
13118 "Error: set_all_manifold_ids_on_boundary() can not be called on an empty Triangulation."));
13120 for (
const auto &cell : this->active_cell_iterators())
13122 if (cell->
face(f)->at_boundary())
13123 cell->
face(f)->set_all_manifold_ids(m_number);
13127template <
int dim,
int spacedim>
13136 "Error: set_all_manifold_ids_on_boundary() can not be called on an empty Triangulation."));
13138 bool boundary_found =
false;
13139 for (
const auto &cell : this->active_cell_iterators())
13143 if (cell->
face(f)->at_boundary() &&
13144 cell->
face(f)->boundary_id() == b_id)
13146 boundary_found =
true;
13147 cell->
face(f)->set_manifold_id(m_number);
13152 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++
e)
13153 if (cell->
line(e)->at_boundary() &&
13154 cell->
line(e)->boundary_id() == b_id)
13156 boundary_found =
true;
13157 cell->
line(e)->set_manifold_id(m_number);
13161 Assert(boundary_found, ExcBoundaryIdNotFound(b_id));
13166template <
int dim,
int spacedim>
13173 return *default_flat_manifold;
13177 const auto it = manifolds.find(m_number);
13179 if (it != manifolds.end())
13182 return *(it->second);
13188 "No manifold of the manifold id " + std::to_string(m_number) +
13189 " has been attached to the triangulation. "
13190 "Please attach the right manifold with Triangulation::set_manifold()."));
13192 return *default_flat_manifold;
13197template <
int dim,
int spacedim>
13202 std::set<types::boundary_id> boundary_ids;
13203 for (
const auto &cell : active_cell_iterators())
13207 boundary_ids.insert(cell->
face(face)->boundary_id());
13209 return {boundary_ids.begin(), boundary_ids.end()};
13214template <
int dim,
int spacedim>
13219 std::set<types::manifold_id> m_ids;
13220 for (
const auto &cell : active_cell_iterators())
13225 m_ids.insert(face->manifold_id());
13228 const auto line_indices = internal::TriaAccessorImplementation::
13229 Implementation::get_line_indices_of_cell(*cell);
13230 for (
unsigned int l = 0;
l < cell->
n_lines(); ++
l)
13232 raw_line_iterator line(
this, 0, line_indices[l]);
13233 m_ids.insert(line->manifold_id());
13237 return {m_ids.begin(), m_ids.end()};
13245template <
int dim,
int spacedim>
13250 Assert((vertices.empty()) && (levels.empty()) && (faces ==
nullptr),
13251 ExcTriangulationNotEmpty(vertices.size(), levels.size()));
13253 (dim == 1 || other_tria.
faces !=
nullptr),
13255 "When calling Triangulation::copy_triangulation(), "
13256 "the target triangulation must be empty but the source "
13257 "triangulation (the argument to this function) must contain "
13258 "something. Here, it seems like the source does not "
13259 "contain anything at all."));
13268 strides = other_tria.
strides;
13272 std::make_unique<internal::TriangulationImplementation::TriaFaces<dim>>(
13273 *other_tria.
faces);
13275 for (
const auto &p : other_tria.
manifolds)
13276 set_manifold(p.first, *p.second);
13280 levels.reserve(other_tria.
levels.size());
13281 for (
const auto &level : other_tria.
levels)
13284 internal::TriangulationImplementation::TriaLevel<dim, spacedim>>(
13291 vertex_to_boundary_id_map_1d =
13292 std::make_unique<std::map<unsigned int, types::boundary_id>>(
13295 vertex_to_manifold_id_map_1d =
13296 std::make_unique<std::map<unsigned int, types::manifold_id>>(
13301 this->policy = other_tria.
policy->clone();
13304 this->periodic_face_pairs_level_0.reserve(
13309 auto entry = other_entry;
13311 cell_iterator(
this, entry.cell[0]->level(), entry.cell[0]->index());
13313 cell_iterator(
this, entry.cell[1]->level(), entry.cell[1]->index());
13314 periodic_face_pairs_level_0.emplace_back(entry);
13317 for (
auto [first_cell_, second_cell_and_orientation] :
13320 auto first_cell = first_cell_;
13321 first_cell.first = cell_iterator(
this,
13322 first_cell.first->level(),
13323 first_cell.first->index());
13324 second_cell_and_orientation.first.first =
13325 cell_iterator(
this,
13326 second_cell_and_orientation.first.first->level(),
13327 second_cell_and_orientation.first.first->index());
13329 this->periodic_face_map[first_cell] = second_cell_and_orientation;
13344template <
int dim,
int spacedim>
13348 if (this->all_reference_cells_are_hyper_cube())
13351 std::make_unique<internal::TriangulationImplementation::PolicyWrapper<
13354 internal::TriangulationImplementation::Implementation>>();
13359 std::make_unique<internal::TriangulationImplementation::PolicyWrapper<
13362 internal::TriangulationImplementation::ImplementationMixedMesh>>();
13368template <
int dim,
int spacedim>
13371 const std::vector<Point<spacedim>> &v,
13372 const std::vector<CellData<dim>> &cells,
13373 const SubCellData &subcelldata)
13375 Assert((vertices.empty()) && (levels.empty()) && (faces ==
nullptr),
13376 ExcTriangulationNotEmpty(vertices.size(), levels.size()));
13391 clear_despite_subscriptions();
13399 reset_cell_vertex_indices_cache();
13401 *
this, levels.size(), number_cache);
13402 reset_active_cell_indices();
13403 reset_global_cell_indices();
13408 if (check_for_distorted_cells)
13410 DistortedCellList distorted_cells = collect_distorted_coarse_cells(*
this);
13414 AssertThrow(distorted_cells.distorted_cells.empty(), distorted_cells);
13447 if ((dim == spacedim - 1) && all_reference_cells_are_hyper_cube())
13455 const bool values[][2] = {{
false,
true}, {
true,
false}};
13458 correct(i, j) =
values[i][j];
13463 const bool values[][4] = {{
false,
true,
true,
false},
13464 {
true,
false,
false,
true},
13465 {
true,
false,
false,
true},
13466 {
false,
true,
true,
false}};
13469 correct(i, j) = (
values[i][j]);
13477 std::list<active_cell_iterator> this_round, next_round;
13478 active_cell_iterator neighbor;
13482 this_round.push_back(begin_active());
13483 begin_active()->set_direction_flag(
true);
13484 begin_active()->set_user_flag();
13486 while (this_round.size() > 0)
13488 for (
const auto &cell : this_round)
13492 if (cell->
face(i)->at_boundary() ==
false)
13498 const unsigned int nb_of_nb =
13503 if (neighbor->user_flag_set())
13506 !(correct(i, nb_of_nb) ^
13507 (neighbor->direction_flag() ==
13510 "The triangulation you are trying to create is not orientable."));
13519 if (correct(i, nb_of_nb) ^
13520 (neighbor->direction_flag() ==
13522 neighbor->set_direction_flag(
13523 !neighbor->direction_flag());
13524 neighbor->set_user_flag();
13525 next_round.push_back(neighbor);
13536 if (next_round.empty())
13537 for (
const auto &cell : this->active_cell_iterators())
13540 next_round.push_back(cell);
13547 next_round.swap(this_round);
13548 next_round.clear();
13550 clear_user_flags();
13553 this->update_cell_relations();
13561template <
int dim,
int spacedim>
13564 const TriangulationDescription::Description<dim, spacedim> &construction_data)
13572 auto cell_infos = construction_data.
cell_infos;
13575 for (
auto &cell_info : cell_infos)
13579 [&](
const TriangulationDescription::CellData<dim> &a,
13580 const TriangulationDescription::CellData<dim> &b) {
13581 const CellId a_id(a.id);
13582 const CellId b_id(b.id);
13584 const auto a_coarse_cell_index =
13585 this->coarse_cell_id_to_coarse_cell_index(a_id.get_coarse_cell_id());
13586 const auto b_coarse_cell_index =
13587 this->coarse_cell_id_to_coarse_cell_index(b_id.get_coarse_cell_id());
13594 if (a_coarse_cell_index != b_coarse_cell_index)
13595 return a_coarse_cell_index < b_coarse_cell_index;
13597 return a_id < b_id;
13603 for (
unsigned int level = 0;
13604 level < cell_infos.size() && !cell_infos[level].empty();
13610 auto cell = this->
begin(level);
13611 auto cell_info = cell_infos[level].begin();
13612 for (; cell_info != cell_infos[level].end(); ++cell_info)
13614 while (cell_info->id != cell->
id().template to_binary<dim>())
13620 cell->
face(face)->set_manifold_id(
13621 cell_info->manifold_line_ids[face]);
13625 cell->
face(face)->set_manifold_id(
13626 cell_info->manifold_quad_ids[face]);
13628 const auto line_indices = internal::TriaAccessorImplementation::
13629 Implementation::get_line_indices_of_cell(*cell);
13630 for (
unsigned int l = 0;
l < cell->
n_lines(); ++
l)
13632 raw_line_iterator line(
this, 0, line_indices[l]);
13633 line->set_manifold_id(cell_info->manifold_line_ids[l]);
13642 if (level + 1 != cell_infos.size())
13646 auto coarse_cell = this->
begin(level);
13647 auto fine_cell_info = cell_infos[level + 1].begin();
13650 for (; fine_cell_info != cell_infos[level + 1].end();
13655 !coarse_cell->id().is_parent_of(CellId(fine_cell_info->id)))
13659 coarse_cell->set_refine_flag();
13663 ::Triangulation<dim,
13664 spacedim>::execute_coarsening_and_refinement();
13669 for (
unsigned int level = 0;
13670 level < cell_infos.size() && !cell_infos[level].empty();
13673 auto cell = this->
begin(level);
13674 auto cell_info = cell_infos[level].begin();
13675 for (; cell_info != cell_infos[level].end(); ++cell_info)
13678 while (cell_info->id != cell->
id().template to_binary<dim>())
13682 for (
auto pair : cell_info->boundary_ids)
13683 if (cell->
face(pair.first)->at_boundary())
13684 cell->
face(pair.first)->set_boundary_id(pair.second);
13693template <
int dim,
int spacedim>
13699 "This function can only be called if dim == spacedim-1."));
13700 for (
const auto &cell : this->active_cell_iterators())
13706template <
int dim,
int spacedim>
13710 if constexpr (dim == 3)
13713 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
13717 if (this->line_to_adjacent_cells_map)
13720 this->line_to_adjacent_cells_map->reinit(
13726 get_reference_cells()[0]
13733 this->line_to_adjacent_cells_map =
13734 std::make_optional<Table<2, std::set<active_cell_iterator>>>(
13735 Table<2, std::set<active_cell_iterator>>(
13738 get_reference_cells()[0].n_lines())));
13742 for (
const auto &cell : this->active_cell_iterators())
13751 &adjacent_cells_to_vertex_0 = vertex_to_cell[vertex_0];
13754 &adjacent_cells_to_vertex_1 = vertex_to_cell[vertex_1];
13757 std::set_intersection(
13758 adjacent_cells_to_vertex_0.begin(),
13759 adjacent_cells_to_vertex_0.end(),
13760 adjacent_cells_to_vertex_1.begin(),
13761 adjacent_cells_to_vertex_1.end(),
13762 std::inserter(line_to_adjacent_cells_map
13764 line_to_adjacent_cells_map
13773template <
int dim,
int spacedim>
13778 ExcMessage(
"Error: An empty Triangulation can not be refined."));
13780 for (
const auto &cell : this->active_cell_iterators())
13790template <
int dim,
int spacedim>
13795 ExcMessage(
"Error: An empty Triangulation can not be refined."));
13797 for (
unsigned int i = 0; i < times; ++i)
13799 set_all_refine_flags();
13800 execute_coarsening_and_refinement();
13806template <
int dim,
int spacedim>
13810 for (
unsigned int i = 0; i < times; ++i)
13812 for (
const auto &cell : this->active_cell_iterators())
13817 execute_coarsening_and_refinement();
13827template <
int dim,
int spacedim>
13832 std::vector<bool>::iterator i = v.begin();
13834 for (
const auto &cell : this->active_cell_iterators())
13835 for (
unsigned int j = 0; j < dim; ++j, ++i)
13844template <
int dim,
int spacedim>
13848 std::vector<bool> v;
13849 save_refine_flags(v);
13858template <
int dim,
int spacedim>
13862 std::vector<bool> v;
13864 load_refine_flags(v);
13869template <
int dim,
int spacedim>
13875 std::vector<bool>::const_iterator i = v.begin();
13876 for (
const auto &cell : this->active_cell_iterators())
13878 unsigned int ref_case = 0;
13880 for (
unsigned int j = 0; j < dim; ++j, ++i)
13882 ref_case += 1 << j;
13884 ExcGridReadError());
13896template <
int dim,
int spacedim>
13899 std::vector<bool> &v)
const
13902 std::vector<bool>::iterator i = v.begin();
13903 for (
const auto &cell : this->active_cell_iterators())
13914template <
int dim,
int spacedim>
13918 std::vector<bool> v;
13919 save_coarsen_flags(v);
13928template <
int dim,
int spacedim>
13932 std::vector<bool> v;
13937 load_coarsen_flags(v);
13942template <
int dim,
int spacedim>
13945 const std::vector<bool> &v)
13949 std::vector<bool>::const_iterator i = v.begin();
13950 for (
const auto &cell : this->active_cell_iterators())
13963template <
int dim,
int spacedim>
13967 return anisotropic_refinement;
13977 template <
int dim,
int spacedim>
13978 std::vector<std::vector<bool>>
13979 extract_raw_coarsen_flags(
13981 std::unique_ptr<::internal::TriangulationImplementation::
13982 TriaLevel<dim, spacedim>>> &levels)
13984 std::vector<std::vector<bool>> coarsen_flags(levels.size());
13985 for (
unsigned int level = 0; level < levels.size(); ++level)
13986 coarsen_flags[level] = levels[level]->coarsen_flags;
13987 return coarsen_flags;
13990 template <
int dim,
int spacedim>
13991 std::vector<std::vector<std::uint8_t>>
13992 extract_raw_refine_flags(
13994 std::unique_ptr<::internal::TriangulationImplementation::
13995 TriaLevel<dim, spacedim>>> &levels)
13997 std::vector<std::vector<std::uint8_t>> refine_flags(levels.size());
13998 for (
unsigned int level = 0; level < levels.size(); ++level)
13999 refine_flags[level] = levels[level]->refine_flags;
14000 return refine_flags;
14012 template <
int dim,
int spacedim>
14015 std::vector<std::unique_ptr<
14016 internal::TriangulationImplementation::TriaLevel<dim, spacedim>>> &levels)
14018 for (
auto &level : levels)
14019 level->cells.clear_user_data();
14026 clear_user_data(internal::TriangulationImplementation::TriaFaces<dim> *faces)
14044template <
int dim,
int spacedim>
14049 ::clear_user_data(levels);
14051 ::clear_user_data(faces.get());
14058 template <
int dim,
int spacedim>
14060 clear_user_flags_line(
14061 std::vector<std::unique_ptr<
14062 internal::TriangulationImplementation::TriaLevel<dim, spacedim>>> &levels,
14063 internal::TriangulationImplementation::TriaFaces<dim> *faces)
14067 for (
const auto &level : levels)
14068 level->cells.clear_user_flags();
14070 else if (dim == 2 || dim == 3)
14082template <
int dim,
int spacedim>
14086 ::clear_user_flags_line(levels, faces.get());
14093 template <
int dim,
int spacedim>
14095 clear_user_flags_quad(
14096 std::vector<std::unique_ptr<
14097 internal::TriangulationImplementation::TriaLevel<dim, spacedim>>> &levels,
14098 internal::TriangulationImplementation::TriaFaces<dim> *faces)
14106 for (
const auto &level : levels)
14107 level->cells.clear_user_flags();
14121template <
int dim,
int spacedim>
14125 ::clear_user_flags_quad(levels, faces.get());
14132 template <
int dim,
int spacedim>
14134 clear_user_flags_hex(
14135 std::vector<std::unique_ptr<
14136 internal::TriangulationImplementation::TriaLevel<dim, spacedim>>> &levels,
14137 internal::TriangulationImplementation::TriaFaces<dim> *)
14149 for (
const auto &level : levels)
14150 level->cells.clear_user_flags();
14160template <
int dim,
int spacedim>
14164 ::clear_user_flags_hex(levels, faces.get());
14169template <
int dim,
int spacedim>
14173 clear_user_flags_line();
14174 clear_user_flags_quad();
14175 clear_user_flags_hex();
14180template <
int dim,
int spacedim>
14184 save_user_flags_line(out);
14187 save_user_flags_quad(out);
14190 save_user_flags_hex(out);
14198template <
int dim,
int spacedim>
14206 std::vector<bool> tmp;
14208 save_user_flags_line(tmp);
14209 v.insert(v.end(), tmp.begin(), tmp.end());
14213 save_user_flags_quad(tmp);
14214 v.insert(v.end(), tmp.begin(), tmp.end());
14219 save_user_flags_hex(tmp);
14220 v.insert(v.end(), tmp.begin(), tmp.end());
14229template <
int dim,
int spacedim>
14233 load_user_flags_line(in);
14236 load_user_flags_quad(in);
14239 load_user_flags_hex(in);
14247template <
int dim,
int spacedim>
14252 std::vector<bool> tmp;
14256 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
14258 load_user_flags_line(tmp);
14263 tmp.insert(tmp.end(),
14264 v.begin() + n_lines(),
14265 v.begin() + n_lines() + n_quads());
14266 load_user_flags_quad(tmp);
14272 tmp.insert(tmp.end(),
14273 v.begin() + n_lines() + n_quads(),
14274 v.begin() + n_lines() + n_quads() + n_hexs());
14275 load_user_flags_hex(tmp);
14284template <
int dim,
int spacedim>
14287 std::vector<bool> &v)
const
14289 v.resize(n_lines(),
false);
14290 std::vector<bool>::iterator i = v.begin();
14291 line_iterator line = begin_line(), endl = end_line();
14292 for (; line != endl; ++line, ++i)
14293 *i = line->user_flag_set();
14300template <
int dim,
int spacedim>
14304 std::vector<bool> v;
14305 save_user_flags_line(v);
14314template <
int dim,
int spacedim>
14318 std::vector<bool> v;
14323 load_user_flags_line(v);
14328template <
int dim,
int spacedim>
14331 const std::vector<bool> &v)
14333 Assert(v.size() == n_lines(), ExcGridReadError());
14335 line_iterator line = begin_line(), endl = end_line();
14336 std::vector<bool>::const_iterator i = v.begin();
14337 for (; line != endl; ++line, ++i)
14339 line->set_user_flag();
14341 line->clear_user_flag();
14350 template <
typename Iterator>
14352 get_user_flag(
const Iterator &i)
14354 return i->user_flag_set();
14359 template <
int structdim,
int dim,
int spacedim>
14361 get_user_flag(
const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
14369 template <
typename Iterator>
14371 set_user_flag(
const Iterator &i)
14373 i->set_user_flag();
14378 template <
int structdim,
int dim,
int spacedim>
14380 set_user_flag(
const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
14387 template <
typename Iterator>
14389 clear_user_flag(
const Iterator &i)
14391 i->clear_user_flag();
14396 template <
int structdim,
int dim,
int spacedim>
14399 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
14407template <
int dim,
int spacedim>
14410 std::vector<bool> &v)
const
14412 v.resize(n_quads(),
false);
14416 std::vector<bool>::iterator i = v.begin();
14417 quad_iterator quad = begin_quad(), endq = end_quad();
14418 for (; quad != endq; ++quad, ++i)
14419 *i = get_user_flag(quad);
14427template <
int dim,
int spacedim>
14431 std::vector<bool> v;
14432 save_user_flags_quad(v);
14441template <
int dim,
int spacedim>
14445 std::vector<bool> v;
14450 load_user_flags_quad(v);
14455template <
int dim,
int spacedim>
14458 const std::vector<bool> &v)
14460 Assert(v.size() == n_quads(), ExcGridReadError());
14464 quad_iterator quad = begin_quad(), endq = end_quad();
14465 std::vector<bool>::const_iterator i = v.begin();
14466 for (; quad != endq; ++quad, ++i)
14468 set_user_flag(quad);
14470 clear_user_flag(quad);
14478template <
int dim,
int spacedim>
14481 std::vector<bool> &v)
const
14483 v.resize(n_hexs(),
false);
14487 std::vector<bool>::iterator i = v.begin();
14488 hex_iterator hex = begin_hex(), endh = end_hex();
14489 for (; hex != endh; ++hex, ++i)
14490 *i = get_user_flag(hex);
14498template <
int dim,
int spacedim>
14502 std::vector<bool> v;
14503 save_user_flags_hex(v);
14512template <
int dim,
int spacedim>
14516 std::vector<bool> v;
14521 load_user_flags_hex(v);
14526template <
int dim,
int spacedim>
14529 const std::vector<bool> &v)
14531 Assert(v.size() == n_hexs(), ExcGridReadError());
14535 hex_iterator hex = begin_hex(), endh = end_hex();
14536 std::vector<bool>::const_iterator i = v.begin();
14537 for (; hex != endh; ++hex, ++i)
14539 set_user_flag(hex);
14541 clear_user_flag(hex);
14549template <
int dim,
int spacedim>
14552 std::vector<unsigned int> &v)
const
14558 std::vector<unsigned int> tmp;
14560 save_user_indices_line(tmp);
14561 v.insert(v.end(), tmp.begin(), tmp.end());
14565 save_user_indices_quad(tmp);
14566 v.insert(v.end(), tmp.begin(), tmp.end());
14571 save_user_indices_hex(tmp);
14572 v.insert(v.end(), tmp.begin(), tmp.end());
14581template <
int dim,
int spacedim>
14584 const std::vector<unsigned int> &v)
14587 std::vector<unsigned int> tmp;
14591 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
14593 load_user_indices_line(tmp);
14598 tmp.insert(tmp.end(),
14599 v.begin() + n_lines(),
14600 v.begin() + n_lines() + n_quads());
14601 load_user_indices_quad(tmp);
14607 tmp.insert(tmp.end(),
14608 v.begin() + n_lines() + n_quads(),
14609 v.begin() + n_lines() + n_quads() + n_hexs());
14610 load_user_indices_hex(tmp);
14619template <
int dim,
int spacedim>
14625 std::ofstream ofs_tria(file_basename +
"_triangulation.data");
14628 boost::archive::text_oarchive oa(ofs_tria, boost::archive::no_header);
14630 internal::CellAttachedDataSerializer<dim, spacedim>::version_number);
14635 std::ofstream ofs_info(file_basename +
".info");
14637 <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_active_cells"
14639 << internal::CellAttachedDataSerializer<dim, spacedim>::version_number
14640 <<
" " << 1 <<
" " << this->cell_attached_data.pack_callbacks_fixed.size()
14641 <<
" " << this->cell_attached_data.pack_callbacks_variable.size() <<
" "
14642 << this->n_global_active_cells() << std::endl;
14645 this->save_attached_data(0, this->n_global_active_cells(), file_basename);
14650template <
int dim,
int spacedim>
14661 std::ifstream ifs_tria(file_basename +
"_triangulation.data");
14664 boost::archive::text_iarchive ia(ifs_tria, boost::archive::no_header);
14666 internal::CellAttachedDataSerializer<dim, spacedim>::version_number);
14670 unsigned int version, numcpus, attached_count_fixed, attached_count_variable,
14671 n_global_active_cells;
14673 std::ifstream ifs_info(std::string(file_basename) +
".info");
14675 std::string firstline;
14676 std::getline(ifs_info, firstline);
14677 ifs_info >> version >> numcpus >> attached_count_fixed >>
14678 attached_count_variable >> n_global_active_cells;
14682 ExcMessage(
"Incompatible number of CPUs found in .info file."));
14684 const auto expected_version =
14685 ::internal::CellAttachedDataSerializer<dim,
14686 spacedim>::version_number;
14689 "The information saved in the file you are trying "
14690 "to read the triangulation from was written with an "
14691 "incompatible file format version and cannot be read."));
14692 Assert(this->n_global_active_cells() == n_global_active_cells,
14693 ExcMessage(
"The number of cells of the triangulation differs "
14694 "from the number of cells written into the .info file."));
14698 this->cell_attached_data.n_attached_data_sets = 0;
14699 this->cell_attached_data.n_attached_deserialize =
14700 attached_count_fixed + attached_count_variable;
14702 this->load_attached_data(0,
14703 this->n_global_active_cells(),
14706 attached_count_fixed,
14707 attached_count_variable);
14709 this->update_cell_relations();
14715 template <
typename Iterator>
14717 get_user_index(
const Iterator &i)
14719 return i->user_index();
14724 template <
int structdim,
int dim,
int spacedim>
14727 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
14735 template <
typename Iterator>
14737 set_user_index(
const Iterator &i,
const unsigned int x)
14739 i->set_user_index(x);
14744 template <
int structdim,
int dim,
int spacedim>
14747 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &,
14748 const unsigned int)
14756template <
int dim,
int spacedim>
14759 std::vector<unsigned int> &v)
const
14761 v.resize(n_lines(), 0);
14762 std::vector<unsigned int>::iterator i = v.begin();
14763 line_iterator line = begin_line(), endl = end_line();
14764 for (; line != endl; ++line, ++i)
14765 *i = line->user_index();
14770template <
int dim,
int spacedim>
14773 const std::vector<unsigned int> &v)
14775 Assert(v.size() == n_lines(), ExcGridReadError());
14777 line_iterator line = begin_line(), endl = end_line();
14778 std::vector<unsigned int>::const_iterator i = v.begin();
14779 for (; line != endl; ++line, ++i)
14780 line->set_user_index(*i);
14784template <
int dim,
int spacedim>
14787 std::vector<unsigned int> &v)
const
14789 v.resize(n_quads(), 0);
14793 std::vector<unsigned int>::iterator i = v.begin();
14794 quad_iterator quad = begin_quad(), endq = end_quad();
14795 for (; quad != endq; ++quad, ++i)
14796 *i = get_user_index(quad);
14802template <
int dim,
int spacedim>
14805 const std::vector<unsigned int> &v)
14807 Assert(v.size() == n_quads(), ExcGridReadError());
14811 quad_iterator quad = begin_quad(), endq = end_quad();
14812 std::vector<unsigned int>::const_iterator i = v.begin();
14813 for (; quad != endq; ++quad, ++i)
14814 set_user_index(quad, *i);
14819template <
int dim,
int spacedim>
14822 std::vector<unsigned int> &v)
const
14824 v.resize(n_hexs(), 0);
14828 std::vector<unsigned int>::iterator i = v.begin();
14829 hex_iterator hex = begin_hex(), endh = end_hex();
14830 for (; hex != endh; ++hex, ++i)
14831 *i = get_user_index(hex);
14837template <
int dim,
int spacedim>
14840 const std::vector<unsigned int> &v)
14842 Assert(v.size() == n_hexs(), ExcGridReadError());
14846 hex_iterator hex = begin_hex(), endh = end_hex();
14847 std::vector<unsigned int>::const_iterator i = v.begin();
14848 for (; hex != endh; ++hex, ++i)
14849 set_user_index(hex, *i);
14861 template <
typename Iterator>
14863 get_user_pointer(
const Iterator &i)
14865 return i->user_pointer();
14870 template <
int structdim,
int dim,
int spacedim>
14873 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &)
14881 template <
typename Iterator>
14883 set_user_pointer(
const Iterator &i,
void *x)
14885 i->set_user_pointer(x);
14890 template <
int structdim,
int dim,
int spacedim>
14893 const TriaIterator<InvalidAccessor<structdim, dim, spacedim>> &,
14902template <
int dim,
int spacedim>
14905 std::vector<void *> &v)
const
14911 std::vector<void *> tmp;
14913 save_user_pointers_line(tmp);
14914 v.insert(v.end(), tmp.begin(), tmp.end());
14918 save_user_pointers_quad(tmp);
14919 v.insert(v.end(), tmp.begin(), tmp.end());
14924 save_user_pointers_hex(tmp);
14925 v.insert(v.end(), tmp.begin(), tmp.end());
14934template <
int dim,
int spacedim>
14937 const std::vector<void *> &v)
14940 std::vector<void *> tmp;
14944 tmp.insert(tmp.end(), v.begin(), v.begin() + n_lines());
14946 load_user_pointers_line(tmp);
14951 tmp.insert(tmp.end(),
14952 v.begin() + n_lines(),
14953 v.begin() + n_lines() + n_quads());
14954 load_user_pointers_quad(tmp);
14960 tmp.insert(tmp.end(),
14961 v.begin() + n_lines() + n_quads(),
14962 v.begin() + n_lines() + n_quads() + n_hexs());
14963 load_user_pointers_hex(tmp);
14972template <
int dim,
int spacedim>
14975 std::vector<void *> &v)
const
14977 v.resize(n_lines(),
nullptr);
14978 std::vector<void *>::iterator i = v.begin();
14979 line_iterator line = begin_line(), endl = end_line();
14980 for (; line != endl; ++line, ++i)
14981 *i = line->user_pointer();
14986template <
int dim,
int spacedim>
14989 const std::vector<void *> &v)
14991 Assert(v.size() == n_lines(), ExcGridReadError());
14993 line_iterator line = begin_line(), endl = end_line();
14994 std::vector<void *>::const_iterator i = v.begin();
14995 for (; line != endl; ++line, ++i)
14996 line->set_user_pointer(*i);
15001template <
int dim,
int spacedim>
15004 std::vector<void *> &v)
const
15006 v.resize(n_quads(),
nullptr);
15010 std::vector<void *>::iterator i = v.begin();
15011 quad_iterator quad = begin_quad(), endq = end_quad();
15012 for (; quad != endq; ++quad, ++i)
15013 *i = get_user_pointer(quad);
15019template <
int dim,
int spacedim>
15022 const std::vector<void *> &v)
15024 Assert(v.size() == n_quads(), ExcGridReadError());
15028 quad_iterator quad = begin_quad(), endq = end_quad();
15029 std::vector<void *>::const_iterator i = v.begin();
15030 for (; quad != endq; ++quad, ++i)
15031 set_user_pointer(quad, *i);
15036template <
int dim,
int spacedim>
15039 std::vector<void *> &v)
const
15041 v.resize(n_hexs(),
nullptr);
15045 std::vector<void *>::iterator i = v.begin();
15046 hex_iterator hex = begin_hex(), endh = end_hex();
15047 for (; hex != endh; ++hex, ++i)
15048 *i = get_user_pointer(hex);
15054template <
int dim,
int spacedim>
15057 const std::vector<void *> &v)
15059 Assert(v.size() == n_hexs(), ExcGridReadError());
15063 hex_iterator hex = begin_hex(), endh = end_hex();
15064 std::vector<void *>::const_iterator i = v.begin();
15065 for (; hex != endh; ++hex, ++i)
15066 set_user_pointer(hex, *i);
15076template <
int dim,
int spacedim>
15084 return begin_raw_line(level);
15086 return begin_raw_quad(level);
15088 return begin_raw_hex(level);
15091 return raw_cell_iterator();
15097template <
int dim,
int spacedim>
15105 return begin_line(level);
15107 return begin_quad(level);
15109 return begin_hex(level);
15112 return cell_iterator();
15118template <
int dim,
int spacedim>
15126 return begin_active_line(level);
15128 return begin_active_quad(level);
15130 return begin_active_hex(level);
15133 return active_cell_iterator();
15139template <
int dim,
int spacedim>
15144 const unsigned int level = levels.size() - 1;
15145 if (levels[level]->cells.n_objects() == 0)
15150 raw_cell_iterator ri(
const_cast<Triangulation<dim, spacedim> *
>(
this),
15152 levels[level]->cells.n_objects() - 1);
15155 if (ri->used() ==
true)
15158 if (ri->used() ==
true)
15165template <
int dim,
int spacedim>
15171 cell_iterator cell = last();
15176 if (cell->is_active() ==
true)
15179 if (cell->is_active() ==
true)
15187template <
int dim,
int spacedim>
15191 const CellId &cell_id)
const
15194 this->contains_cell(cell_id),
15196 "CellId is invalid for this triangulation.\n"
15197 "Either the provided CellId does not correspond to a cell in this "
15198 "triangulation object, or, in case you are using a parallel "
15199 "triangulation, may correspond to an artificial cell that is less "
15200 "refined on this processor. In the case of "
15201 "parallel::fullydistributed::Triangulation, the corresponding coarse "
15202 "cell might not be accessible by the current process."));
15204 cell_iterator cell(
15208 cell = cell->child(
static_cast<unsigned int>(child_index));
15215template <
int dim,
int spacedim>
15219 const auto coarse_cell_index =
15225 cell_iterator cell(
this, 0, coarse_cell_index);
15229 if (cell->has_children() ==
false)
15231 cell = cell->child(
static_cast<unsigned int>(child_index));
15239template <
int dim,
int spacedim>
15244 return cell_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
15251template <
int dim,
int spacedim>
15264 if (level >= levels.size())
15266 Assert(level < n_global_levels(),
15267 ExcInvalidLevel(level, n_global_levels()));
15273 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
15274 if (level < levels.size() - 1)
15275 return begin_raw(level + 1);
15281template <
int dim,
int spacedim>
15294 if (level >= levels.size())
15296 Assert(level < n_global_levels(),
15297 ExcInvalidLevel(level, n_global_levels()));
15303 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
15304 if (level < levels.size() - 1)
15305 return begin(level + 1);
15311template <
int dim,
int spacedim>
15324 if (level >= levels.size())
15326 Assert(level < n_global_levels(),
15327 ExcInvalidLevel(level, n_global_levels()));
15333 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
15334 return (level >= levels.size() - 1 ? active_cell_iterator(
end()) :
15335 begin_active(level + 1));
15340template <
int dim,
int spacedim>
15342IteratorRange<typename Triangulation<dim, spacedim>::
15346 return IteratorRange<typename Triangulation<dim, spacedim>::cell_iterator>(
15351template <
int dim,
int spacedim>
15353IteratorRange<typename Triangulation<dim, spacedim>::
15357 return IteratorRange<
15364template <
int dim,
int spacedim>
15366IteratorRange<typename Triangulation<dim, spacedim>::
15370 return IteratorRange<typename Triangulation<dim, spacedim>::cell_iterator>(
15376template <
int dim,
int spacedim>
15378IteratorRange<typename Triangulation<dim, spacedim>::
15382 return IteratorRange<
15384 begin_active(level), end_active(level));
15392template <
int dim,
int spacedim>
15401 return raw_face_iterator();
15403 return begin_line();
15405 return begin_quad();
15408 return face_iterator();
15414template <
int dim,
int spacedim>
15423 return raw_face_iterator();
15425 return begin_active_line();
15427 return begin_active_quad();
15430 return active_face_iterator();
15436template <
int dim,
int spacedim>
15445 return raw_face_iterator();
15452 return raw_face_iterator();
15458template <
int dim,
int spacedim>
15460IteratorRange<typename Triangulation<dim, spacedim>::
15464 return IteratorRange<
15466 begin_active_face(), end_face());
15472template <
int dim,
int spacedim>
15477 vertex_iterator i =
15478 raw_vertex_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this), 0, 0);
15482 while (i->used() ==
false)
15490template <
int dim,
int spacedim>
15496 return begin_vertex();
15501template <
int dim,
int spacedim>
15506 return raw_vertex_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
15518template <
int dim,
int spacedim>
15531 if (level >= levels.size())
15533 Assert(level < n_global_levels(),
15534 ExcInvalidLevel(level, n_global_levels()));
15543 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
15545 if (level >= levels.size() || levels[level]->cells.n_objects() == 0)
15548 return raw_line_iterator(
15549 const_cast<Triangulation<dim, spacedim> *
>(
this), level, 0);
15552 Assert(level == 0, ExcFacesHaveNoLevel());
15553 return raw_line_iterator(
15554 const_cast<Triangulation<dim, spacedim> *
>(
this), 0, 0);
15559template <
int dim,
int spacedim>
15565 raw_line_iterator ri = begin_raw_line(level);
15568 while (ri->used() ==
false)
15576template <
int dim,
int spacedim>
15580 const unsigned int level)
const
15583 line_iterator i = begin_line(level);
15586 while (i->has_children())
15594template <
int dim,
int spacedim>
15599 return raw_line_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
15610template <
int dim,
int spacedim>
15623 if (level >= levels.size())
15625 Assert(level < n_global_levels(),
15626 ExcInvalidLevel(level, n_global_levels()));
15634 return raw_hex_iterator();
15639 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
15641 if (level >= levels.size() || levels[level]->cells.n_objects() == 0)
15644 return raw_quad_iterator(
15645 const_cast<Triangulation<dim, spacedim> *
>(
this), level, 0);
15650 Assert(level == 0, ExcFacesHaveNoLevel());
15652 return raw_quad_iterator(
15653 const_cast<Triangulation<dim, spacedim> *
>(
this), 0, 0);
15659 return raw_hex_iterator();
15665template <
int dim,
int spacedim>
15671 raw_quad_iterator ri = begin_raw_quad(level);
15674 while (ri->used() ==
false)
15682template <
int dim,
int spacedim>
15686 const unsigned int level)
const
15689 quad_iterator i = begin_quad(level);
15692 while (i->has_children())
15700template <
int dim,
int spacedim>
15705 return raw_quad_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
15716template <
int dim,
int spacedim>
15729 if (level >= levels.size())
15731 Assert(level < n_global_levels(),
15732 ExcInvalidLevel(level, n_global_levels()));
15741 return raw_hex_iterator();
15746 Assert(level < levels.size(), ExcInvalidLevel(level, levels.size()));
15748 if (level >= levels.size() || levels[level]->cells.n_objects() == 0)
15751 return raw_hex_iterator(
15752 const_cast<Triangulation<dim, spacedim> *
>(
this), level, 0);
15757 return raw_hex_iterator();
15763template <
int dim,
int spacedim>
15769 raw_hex_iterator ri = begin_raw_hex(level);
15772 while (ri->used() ==
false)
15780template <
int dim,
int spacedim>
15786 hex_iterator i = begin_hex(level);
15789 while (i->has_children())
15797template <
int dim,
int spacedim>
15802 return raw_hex_iterator(
const_cast<Triangulation<dim, spacedim> *
>(
this),
15814 namespace TriangulationImplementation
15827 return c.n_active_lines;
15842 return c.n_active_quads;
15857 return c.n_active_hexes;
15864template <
int dim,
int spacedim>
15872template <
int dim,
int spacedim>
15879template <
int dim,
int spacedim>
15887template <
int dim,
int spacedim>
15895template <
int dim,
int spacedim>
15902 return n_used_vertices();
15914template <
int dim,
int spacedim>
15921 return n_vertices();
15923 return n_raw_lines();
15925 return n_raw_quads();
15933template <
int dim,
int spacedim>
15940 return n_used_vertices();
15942 return n_active_lines();
15944 return n_active_quads();
15952template <
int dim,
int spacedim>
15955 const unsigned int level)
const
15960 return n_raw_lines(level);
15962 return n_raw_quads(level);
15964 return n_raw_hexs(level);
15973template <
int dim,
int spacedim>
15976 const unsigned int level)
const
15981 return n_lines(level);
15983 return n_quads(level);
15985 return n_hexs(level);
15994template <
int dim,
int spacedim>
15997 const unsigned int level)
const
16002 return n_active_lines(level);
16004 return n_active_quads(level);
16006 return n_active_hexs(level);
16014template <
int dim,
int spacedim>
16018 if (anisotropic_refinement ==
false)
16020 for (
unsigned int lvl = 0; lvl < n_global_levels() - 1; ++lvl)
16026 for (
const auto &cell : active_cell_iterators())
16027 for (
const auto &i : cell->face_indices())
16028 if (cell->face(i)->has_children())
16035template <
int dim,
int spacedim>
16039 return number_cache.n_lines;
16044template <
int dim,
int spacedim>
16047 const unsigned int level)
const
16052 return levels[level]->cells.n_objects();
16055 Assert(
false, ExcFacesHaveNoLevel());
16060template <
int dim,
int spacedim>
16074template <
int dim,
int spacedim>
16077 const unsigned int level)
const
16080 Assert(dim == 1, ExcFacesHaveNoLevel());
16081 return number_cache.n_lines_level[level];
16085template <
int dim,
int spacedim>
16089 return number_cache.n_active_lines;
16093template <
int dim,
int spacedim>
16096 const unsigned int level)
const
16099 Assert(dim == 1, ExcFacesHaveNoLevel());
16101 return number_cache.n_active_lines_level[level];
16106template <
int dim,
int spacedim>
16110 if constexpr (dim == 1)
16113 return number_cache.n_quads;
16118template <
int dim,
int spacedim>
16121 const unsigned int level)
const
16124 if constexpr (dim == 1)
16126 else if constexpr (dim == 2)
16129 return number_cache.n_quads_level[level];
16131 else if constexpr (dim == 3)
16132 Assert(dim < 3, ExcFacesHaveNoLevel());
16141template <
int dim,
int spacedim>
16145 if constexpr (dim == 1)
16147 else if constexpr (dim == 2)
16149 else if constexpr (dim == 3)
16159template <
int dim,
int spacedim>
16162 const unsigned int level)
const
16165 if constexpr (dim == 1)
16167 else if constexpr (dim == 2)
16170 return levels[level]->cells.n_objects();
16172 else if constexpr (dim == 3)
16173 Assert(dim < 3, ExcFacesHaveNoLevel());
16182template <
int dim,
int spacedim>
16186 if constexpr (dim == 1)
16189 return number_cache.n_active_quads;
16194template <
int dim,
int spacedim>
16197 const unsigned int level)
const
16200 if constexpr (dim == 1)
16202 else if constexpr (dim == 2)
16205 return number_cache.n_active_quads_level[level];
16207 else if constexpr (dim == 3)
16208 Assert(dim < 3, ExcFacesHaveNoLevel());
16217template <
int dim,
int spacedim>
16221 if constexpr (dim == 3)
16222 return number_cache.n_hexes;
16229template <
int dim,
int spacedim>
16232 const unsigned int level)
const
16235 if constexpr (dim == 3)
16238 return number_cache.n_hexes_level[level];
16246template <
int dim,
int spacedim>
16249 const unsigned int level)
const
16252 if constexpr (dim == 3)
16255 return levels[level]->cells.n_objects();
16262template <
int dim,
int spacedim>
16266 if constexpr (dim == 3)
16267 return number_cache.n_active_hexes;
16274template <
int dim,
int spacedim>
16277 const unsigned int level)
const
16280 if constexpr (dim == 3)
16283 return number_cache.n_active_hexes_level[level];
16291template <
int dim,
int spacedim>
16295 return std::count(vertices_used.begin(), vertices_used.end(),
true);
16300template <
int dim,
int spacedim>
16304 return vertices_used;
16335template <
int dim,
int spacedim>
16339 cell_iterator cell =
begin(0),
16340 endc = (n_levels() > 1 ?
begin(1) : cell_iterator(
end()));
16343 unsigned int max_vertex_index = 0;
16344 for (; cell != endc; ++cell)
16346 if (cell->vertex_index(vertex) > max_vertex_index)
16347 max_vertex_index = cell->vertex_index(vertex);
16353 std::vector<unsigned short int> usage_count(max_vertex_index + 1, 0);
16357 for (cell =
begin(); cell != endc; ++cell)
16359 ++usage_count[cell->vertex_index(vertex)];
16362 static_cast<unsigned int>(
16363 *std::max_element(usage_count.begin(), usage_count.end())));
16368template <
int dim,
int spacedim>
16378template <
int dim,
int spacedim>
16387template <
int dim,
int spacedim>
16389const Triangulation<dim, spacedim>
16397template <
int dim,
int spacedim>
16400 const std::vector<GridTools::PeriodicFacePair<cell_iterator>>
16401 &periodicity_vector)
16403 periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(),
16404 periodicity_vector.begin(),
16405 periodicity_vector.end());
16408 update_periodic_face_map();
16413template <
int dim,
int spacedim>
16415const typename std::map<
16416 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
unsigned int>,
16417 std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
16422 return periodic_face_map;
16426template <
int dim,
int spacedim>
16433 if (
dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *
>(
16437 this->local_cell_relations.clear();
16440 for (
const auto &cell : this->active_cell_iterators())
16441 this->local_cell_relations.emplace_back(
16447template <
int dim,
int spacedim>
16451 if (
dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *
>(
16456 if (this->cell_attached_data.n_attached_data_sets > 0)
16458 this->local_cell_relations.clear();
16460 for (
const auto &cell : this->active_cell_iterators())
16463 (cell->level() > 0) && (cell->coarsen_flag_set());
16468 cell->refine_flag_set() ?
16472 local_cell_relations.emplace_back(cell, status);
16474 else if (cell->parent()->child_iterator_to_index(cell) == 0)
16476 local_cell_relations.emplace_back(
16480 for (
unsigned int i = 1; i < cell->parent()->n_children(); ++i)
16481 local_cell_relations.emplace_back(
16487 this->data_serializer.pack_data(
16488 this->local_cell_relations,
16489 this->cell_attached_data.pack_callbacks_fixed,
16490 this->cell_attached_data.pack_callbacks_variable,
16491 this->get_mpi_communicator());
16494 this->data_serializer.dest_data_fixed =
16495 this->data_serializer.src_data_fixed;
16496 this->data_serializer.dest_data_variable =
16497 this->data_serializer.src_data_variable;
16498 this->data_serializer.dest_sizes_variable =
16499 this->data_serializer.src_sizes_variable;
16505template <
int dim,
int spacedim>
16509 if (
dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *
>(
16514 if (this->cell_attached_data.n_attached_data_sets > 0)
16516 std::vector<typename internal::CellAttachedDataSerializer<dim, spacedim>::
16520 for (
const auto &cell : local_cell_relations)
16522 if (cell.first->has_children())
16527 temp.emplace_back(cell.first->child(0),
16531 temp.push_back(cell);
16534 this->local_cell_relations = std::move(temp);
16540template <
int dim,
int spacedim>
16554 if (smooth_grid & limit_level_difference_at_vertices)
16558 signals.pre_refinement();
16560 this->pack_data_serial();
16562 execute_coarsening();
16564 const DistortedCellList cells_with_distorted_children = execute_refinement();
16569 this->unpack_data_serial();
16571 reset_cell_vertex_indices_cache();
16575 line_to_adjacent_cells_map.reset();
16580 if (smooth_grid & limit_level_difference_at_vertices)
16585 this->policy->update_neighbors(*
this);
16586 reset_active_cell_indices();
16588 reset_global_cell_indices();
16591 signals.post_refinement();
16593 AssertThrow(cells_with_distorted_children.distorted_cells.empty(),
16594 cells_with_distorted_children);
16596 update_periodic_face_map();
16598 if (this->cell_attached_data.n_attached_data_sets == 0)
16599 this->update_cell_relations();
16625 for (
const auto &cell : cell_iterators())
16627 if (cell->has_children() && cell->reference_cell().is_hyper_cube())
16628 for (
const unsigned int f : cell->face_indices())
16629 if (cell->at_boundary(f) ==
false)
16631 for (
const auto &child : cell->child_iterators())
16634 child->at_boundary(f) ==
false,
16636 "We ended up with a triangulation whose child cells "
16637 "are not connected to their neighbors as expected. "
16638 "When you created the triangulation, did you forget "
16639 "to call GridTools::consistently_order_cells() "
16640 "before calling Triangulation::create_triangulation()?"));
16649template <
int dim,
int spacedim>
16653 unsigned int active_cell_index = 0;
16654 for (raw_cell_iterator cell = begin_raw(); cell !=
end(); ++cell)
16655 if ((cell->used() ==
false) || cell->has_children())
16659 cell->set_active_cell_index(active_cell_index);
16660 ++active_cell_index;
16668template <
int dim,
int spacedim>
16673 for (
unsigned int l = 0;
l < levels.size(); ++
l)
16676 for (
const auto &cell : cell_iterators_on_level(l))
16678 cell->set_global_level_cell_index(level_cell_index++);
16679 if (cell->is_active())
16680 cell->set_global_active_cell_index(active_cell_index++);
16688template <
int dim,
int spacedim>
16692 std::array<unsigned int, ReferenceCells::max_n_vertices<dim>()>
16694 for (
unsigned int l = 0;
l < levels.size(); ++
l)
16695 for (
const auto &cell : cell_iterators_on_level(l))
16700 const auto n_vertices = cell->n_vertices();
16703 cell, ArrayView<unsigned int>(cell_vertices.data(), n_vertices));
16704 for (
unsigned int vertex_no = 0; vertex_no < n_vertices; ++vertex_no)
16705 levels[l]->set_cached_vertex_index(cell->index(),
16707 cell_vertices[vertex_no]);
16713template <
int dim,
int spacedim>
16718 periodic_face_map.clear();
16720 typename std::vector<
16722 for (it = periodic_face_pairs_level_0.begin();
16723 it != periodic_face_pairs_level_0.end();
16726 update_periodic_face_map_recursively<dim, spacedim>(it->cell[0],
16731 periodic_face_map);
16733 const auto face_reference_cell =
16734 it->cell[0]->reference_cell().face_reference_cell(it->face_idx[0]);
16736 update_periodic_face_map_recursively<dim, spacedim>(
16741 face_reference_cell.get_inverse_combined_orientation(it->orientation),
16742 periodic_face_map);
16746 typename std::map<std::pair<cell_iterator, unsigned int>,
16747 std::pair<std::pair<cell_iterator, unsigned int>,
16750 for (it_test = periodic_face_map.begin(); it_test != periodic_face_map.end();
16754 it_test->first.first;
16756 it_test->second.first.first;
16761 Assert(periodic_face_map[it_test->second.first].first ==
16770template <
int dim,
int spacedim>
16772const std::vector<ReferenceCell<dim>>
16775 return this->reference_cells;
16780template <
int dim,
int spacedim>
16784 Assert(this->reference_cells.size() > 0,
16785 ExcMessage(
"You can't ask about the kinds of reference "
16786 "cells used by this triangulation if the "
16787 "triangulation doesn't yet have any cells in it."));
16788 return (this->reference_cells.size() == 1 &&
16789 this->reference_cells[0].is_hyper_cube());
16794template <
int dim,
int spacedim>
16798 Assert(this->reference_cells.size() > 0,
16799 ExcMessage(
"You can't ask about the kinds of reference "
16800 "cells used by this triangulation if the "
16801 "triangulation doesn't yet have any cells in it."));
16802 return (this->reference_cells.size() == 1 &&
16803 this->reference_cells[0].is_simplex());
16808template <
int dim,
int spacedim>
16812 Assert(this->reference_cells.size() > 0,
16813 ExcMessage(
"You can't ask about the kinds of reference "
16814 "cells used by this triangulation if the "
16815 "triangulation doesn't yet have any cells in it."));
16816 return reference_cells.size() > 1 ||
16817 ((reference_cells[0].is_hyper_cube() ==
false) &&
16818 (reference_cells[0].is_simplex() ==
false));
16823template <
int dim,
int spacedim>
16826 const std::function<std::vector<char>(
const cell_iterator &,
16827 const ::CellStatus)>
16829 const bool returns_variable_size_data)
16835 if (returns_variable_size_data)
16837 handle = 2 * this->cell_attached_data.pack_callbacks_variable.size();
16838 this->cell_attached_data.pack_callbacks_variable.push_back(pack_callback);
16842 handle = 2 * this->cell_attached_data.pack_callbacks_fixed.size() + 1;
16843 this->cell_attached_data.pack_callbacks_fixed.push_back(pack_callback);
16847 ++this->cell_attached_data.n_attached_data_sets;
16854template <
int dim,
int spacedim>
16857 const unsigned int handle,
16858 const std::function<
16859 void(
const cell_iterator &,
16860 const ::CellStatus,
16861 const boost::iterator_range<std::vector<char>::const_iterator> &)>
16865 this->data_serializer.unpack_data(this->local_cell_relations,
16870 --this->cell_attached_data.n_attached_data_sets;
16871 if (this->cell_attached_data.n_attached_deserialize > 0)
16872 --this->cell_attached_data.n_attached_deserialize;
16881 if (this->cell_attached_data.n_attached_data_sets == 0 &&
16882 this->cell_attached_data.n_attached_deserialize == 0)
16885 this->cell_attached_data.pack_callbacks_fixed.clear();
16886 this->cell_attached_data.pack_callbacks_variable.clear();
16887 this->data_serializer.clear();
16890 for (
auto &cell_rel : this->local_cell_relations)
16897template <
int dim,
int spacedim>
16900 const unsigned int global_first_cell,
16901 const unsigned int global_num_cells,
16902 const std::string &file_basename)
const
16905 auto tria =
const_cast<Triangulation<dim, spacedim> *
>(
this);
16908 for (
const auto &cell_rel : this->local_cell_relations)
16911 Assert((cell_rel.second ==
16916 if (this->cell_attached_data.n_attached_data_sets > 0)
16919 tria->data_serializer.pack_data(
16920 tria->local_cell_relations,
16921 tria->cell_attached_data.pack_callbacks_fixed,
16922 tria->cell_attached_data.pack_callbacks_variable,
16923 this->get_mpi_communicator());
16926 tria->data_serializer.save(global_first_cell,
16929 this->get_mpi_communicator());
16932 tria->data_serializer.clear();
16938 tria->cell_attached_data.n_attached_data_sets = 0;
16939 tria->cell_attached_data.pack_callbacks_fixed.clear();
16940 tria->cell_attached_data.pack_callbacks_variable.clear();
16945template <
int dim,
int spacedim>
16948 const unsigned int global_first_cell,
16949 const unsigned int global_num_cells,
16950 const unsigned int local_num_cells,
16951 const std::string &file_basename,
16952 const unsigned int n_attached_deserialize_fixed,
16953 const unsigned int n_attached_deserialize_variable)
16956 if (this->cell_attached_data.n_attached_deserialize > 0)
16958 this->data_serializer.load(global_first_cell,
16962 n_attached_deserialize_fixed,
16963 n_attached_deserialize_variable,
16964 this->get_mpi_communicator());
16966 this->data_serializer.unpack_cell_status(this->local_cell_relations);
16972 for (
const auto &cell_rel : this->local_cell_relations)
16974 Assert((cell_rel.second ==
16983template <
int dim,
int spacedim>
16991 vertices_used.clear();
16999 vertex_to_boundary_id_map_1d->clear();
17002 vertex_to_manifold_id_map_1d->clear();
17012 number_cache = internal::TriangulationImplementation::NumberCache<dim>();
17017template <
int dim,
int spacedim>
17019typename Triangulation<dim, spacedim>::DistortedCellList
17022 const DistortedCellList cells_with_distorted_children =
17023 this->policy->execute_refinement(*
this, check_for_distorted_cells);
17029 *
this, levels.size(), number_cache);
17033 for (
const auto &level : levels)
17040 for (
const auto &cell : this->cell_iterators())
17044 return cells_with_distorted_children;
17049template <
int dim,
int spacedim>
17055 const cell_iterator endc =
end();
17056 bool do_coarsen =
false;
17057 if (levels.size() >= 2)
17058 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
17059 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
17071 std::vector<unsigned int> line_cell_count(dim > 1 ? this->n_raw_lines() : 0);
17072 std::vector<unsigned int> quad_cell_count(dim > 2 ? this->n_raw_quads() : 0);
17074 for (
const auto &cell : this->cell_iterators())
17078 const auto line_indices = internal::TriaAccessorImplementation::
17079 Implementation::get_line_indices_of_cell(*cell);
17082 const unsigned int n_lines =
std::min(cell->n_lines(), 12u);
17083 for (
unsigned int l = 0;
l < n_lines; ++
l)
17084 ++line_cell_count[line_indices[l]];
17085 for (
const unsigned int q : cell->face_indices())
17086 ++quad_cell_count[cell->face_index(q)];
17089 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
17090 ++line_cell_count[cell->line(l)->index()];
17107 if (levels.size() >= 2)
17108 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
17109 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
17111 for (
unsigned int child = 0; child < cell->n_children(); ++child)
17113 Assert(cell->child(child)->coarsen_flag_set(),
17115 cell->child(child)->clear_coarsen_flag();
17118 signals.pre_coarsening_on_cell(cell);
17120 this->policy->delete_children(*
this,
17128 *
this, levels.size(), number_cache);
17133template <
int dim,
int spacedim>
17150 auto previous_coarsen_flags = internal::extract_raw_coarsen_flags(levels);
17152 bool continue_iterating =
true;
17156 if (smooth_grid & limit_level_difference_at_vertices)
17158 Assert(!anisotropic_refinement,
17159 ExcMessage(
"In case of anisotropic refinement the "
17160 "limit_level_difference_at_vertices flag for "
17161 "mesh smoothing must not be set!"));
17165 std::vector<int> vertex_level(vertices.size(), 0);
17166 for (
const auto &cell : this->active_cell_iterators())
17168 if (cell->refine_flag_set())
17169 for (
const unsigned int vertex : cell->vertex_indices())
17170 vertex_level[cell->vertex_index(vertex)] =
17171 std::max(vertex_level[cell->vertex_index(vertex)],
17172 cell->level() + 1);
17173 else if (!cell->coarsen_flag_set())
17174 for (
const unsigned int vertex : cell->vertex_indices())
17175 vertex_level[cell->vertex_index(vertex)] =
17176 std::max(vertex_level[cell->vertex_index(vertex)],
17187 for (
const unsigned int vertex : cell->vertex_indices())
17188 vertex_level[cell->vertex_index(vertex)] =
17189 std::max(vertex_level[cell->vertex_index(vertex)],
17190 cell->level() - 1);
17204 active_cell_iterator endc =
end();
17205 for (active_cell_iterator cell = last_active(); cell != endc; --cell)
17206 if (cell->refine_flag_set() ==
false)
17208 for (
const unsigned int vertex : cell->vertex_indices())
17209 if (vertex_level[cell->vertex_index(vertex)] >=
17213 cell->clear_coarsen_flag();
17218 if (vertex_level[cell->vertex_index(vertex)] >
17221 cell->set_refine_flag();
17223 for (
const unsigned int v : cell->vertex_indices())
17224 vertex_level[cell->vertex_index(v)] =
17225 std::max(vertex_level[cell->vertex_index(v)],
17226 cell->level() + 1);
17242 for (
const auto &acell : this->active_cell_iterators_on_level(0))
17243 acell->clear_coarsen_flag();
17245 const cell_iterator endc =
end();
17246 for (cell_iterator cell =
begin(n_levels() - 1); cell != endc; --cell)
17249 if (cell->is_active())
17252 const unsigned int n_children = cell->n_children();
17253 unsigned int flagged_children = 0;
17254 for (
unsigned int child = 0; child < n_children; ++child)
17256 const auto child_cell = cell->child(child);
17257 if (child_cell->is_active() && child_cell->coarsen_flag_set())
17259 ++flagged_children;
17261 child_cell->clear_coarsen_flag();
17267 if (flagged_children == n_children &&
17268 this->policy->coarsening_allowed(cell))
17269 for (
unsigned int c = 0; c < n_children; ++c)
17271 Assert(cell->child(c)->refine_flag_set() ==
false,
17274 cell->child(c)->set_coarsen_flag();
17280 auto current_coarsen_flags = internal::extract_raw_coarsen_flags(levels);
17282 continue_iterating = (current_coarsen_flags != previous_coarsen_flags);
17283 previous_coarsen_flags.swap(current_coarsen_flags);
17285 while (continue_iterating ==
true);
17297 const auto flags_before = internal::extract_raw_coarsen_flags(
levels);
17302 const auto flags_after = internal::extract_raw_coarsen_flags(
levels);
17304 return (flags_before != flags_after);
17315 const auto flags_before = internal::extract_raw_coarsen_flags(
levels);
17320 const auto flags_after = internal::extract_raw_coarsen_flags(
levels);
17322 return (flags_before != flags_after);
17333 const auto flags_before = internal::extract_raw_coarsen_flags(
levels);
17338 const auto flags_after = internal::extract_raw_coarsen_flags(
levels);
17340 return (flags_before != flags_after);
17351 template <
int dim,
int spacedim>
17353 possibly_do_not_produce_unrefined_islands(
17358 unsigned int n_neighbors = 0;
17361 unsigned int count = 0;
17369 if (face_will_be_refined_by_neighbor(cell, n))
17376 if (count == n_neighbors ||
17377 (count >= n_neighbors - 1 &&
17380 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
17386 (cell_will_be_coarsened(cell->
neighbor(face))))
17387 possibly_do_not_produce_unrefined_islands<dim, spacedim>(
17403 template <
int dim,
int spacedim>
17405 possibly_refine_unrefined_island(
17407 const bool allow_anisotropic_smoothing)
17416 const parallel::distributed::Triangulation<dim, spacedim> *
>(
17428 "The triangulation is trying to avoid unrefined islands "
17429 "during mesh refinement/coarsening, as you had requested "
17430 " by passing the appropriate 'smoothing flags' to the "
17431 "constructor of the triangulation. However, for objects "
17432 "of type parallel::distributed::Triangulation, control "
17433 "over which cells get refined rests with p4est, not the "
17434 "deal.II triangulation, and consequently it is not "
17435 "always possible to avoid unrefined islands in the mesh. "
17436 "Please remove the constructor argument to the triangulation "
17437 "object that requests mesh smoothing."));
17449 if (allow_anisotropic_smoothing ==
false)
17452 unsigned int refined_neighbors = 0, unrefined_neighbors = 0;
17456 if (face_will_be_refined_by_neighbor(cell, face))
17457 ++refined_neighbors;
17459 ++unrefined_neighbors;
17462 if (unrefined_neighbors < refined_neighbors)
17470 if (unrefined_neighbors > 0)
17473 (face_will_be_refined_by_neighbor(cell, face) ==
false) &&
17476 possibly_refine_unrefined_island<dim, spacedim>(
17477 cell->
neighbor(face), allow_anisotropic_smoothing);
17484 RefinementCase<dim> smoothing_cell_refinement_case =
17489 for (
unsigned int face_pair = 0;
17490 face_pair < GeometryInfo<dim>::faces_per_cell / 2;
17496 RefinementCase<dim> directional_cell_refinement_case =
17499 for (
unsigned int face_index = 0; face_index < 2; ++face_index)
17501 unsigned int face = 2 * face_pair + face_index;
17504 RefinementCase<dim - 1> expected_face_ref_case =
17505 RefinementCase<dim - 1>::no_refinement;
17508 face_will_be_refined_by_neighbor<dim, spacedim>(
17509 cell, face, expected_face_ref_case);
17523 directional_cell_refinement_case =
17524 (directional_cell_refinement_case &
17527 expected_face_ref_case,
17537 Assert(directional_cell_refinement_case <
17540 smoothing_cell_refinement_case =
17541 smoothing_cell_refinement_case | directional_cell_refinement_case;
17546 if (smoothing_cell_refinement_case)
17550 smoothing_cell_refinement_case);
17557template <
int dim,
int spacedim>
17563 const auto coarsen_flags_before = internal::extract_raw_coarsen_flags(levels);
17564 const auto refine_flags_before = internal::extract_raw_refine_flags(levels);
17582 auto coarsen_flags_before_loop = coarsen_flags_before;
17583 auto refine_flags_before_loop = refine_flags_before;
17660 if (((smooth_grid & coarsest_level_1) || (smooth_grid & patch_level_1)) &&
17663 for (
const auto &cell : active_cell_iterators_on_level(1))
17667 bool mesh_changed_in_this_loop =
false;
17676 if (smooth_grid & do_not_produce_unrefined_islands &&
17677 !(smooth_grid & patch_level_1))
17679 for (
const auto &cell : cell_iterators())
17683 if (!cell->
is_active() && cell_will_be_coarsened(cell))
17684 possibly_do_not_produce_unrefined_islands<dim, spacedim>(cell);
17708 if (smooth_grid & (eliminate_refined_inner_islands |
17709 eliminate_refined_boundary_islands) &&
17710 !(smooth_grid & patch_level_1))
17712 for (
const auto &cell : cell_iterators())
17722 bool all_children_active =
true;
17724 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
17729 all_children_active =
false;
17733 if (all_children_active)
17740 unsigned int unrefined_neighbors = 0, total_neighbors = 0;
17751 bool at_periodic_boundary =
false;
17755 const cell_iterator neighbor = cell->
neighbor(n);
17760 if (!face_will_be_refined_by_neighbor(cell, n))
17761 ++unrefined_neighbors;
17766 at_periodic_boundary =
true;
17781 if ((unrefined_neighbors == total_neighbors) &&
17783 (smooth_grid & eliminate_refined_inner_islands)) ||
17786 eliminate_refined_boundary_islands))) &&
17787 (total_neighbors != 0))
17790 for (
unsigned int c = 0; c < cell->
n_children(); ++c)
17811 if (smooth_grid & limit_level_difference_at_vertices)
17813 Assert(!anisotropic_refinement,
17814 ExcMessage(
"In case of anisotropic refinement the "
17815 "limit_level_difference_at_vertices flag for "
17816 "mesh smoothing must not be set!"));
17820 std::vector<int> vertex_level(vertices.size(), 0);
17821 for (
const auto &cell : active_cell_iterators())
17827 cell->
level() + 1);
17843 cell->
level() - 1);
17857 for (active_cell_iterator cell = last_active(); cell !=
end(); --cell)
17878 cell->
level() + 1);
17899 if (smooth_grid & eliminate_unrefined_islands)
17901 for (active_cell_iterator cell = last_active(); cell !=
end(); --cell)
17906 possibly_refine_unrefined_island<dim, spacedim>(
17907 cell, (smooth_grid & allow_anisotropic_smoothing) != 0);
17935 if (smooth_grid & patch_level_1)
17945 for (
const auto &cell : cell_iterators())
17957 RefinementCase<dim> combined_ref_case =
17959 for (
unsigned int i = 0; i < cell->
n_children(); ++i)
17960 combined_ref_case =
17963 for (
unsigned int i = 0; i < cell->
n_children(); ++i)
17965 cell_iterator child = cell->
child(i);
17967 child->clear_coarsen_flag();
17968 child->set_refine_flag(combined_ref_case);
17981 for (
const auto &cell : cell_iterators())
17994 const unsigned int n_children = cell->
n_children();
17995 bool has_active_grandchildren =
false;
17997 for (
unsigned int i = 0; i < n_children; ++i)
18000 has_active_grandchildren =
true;
18004 if (has_active_grandchildren ==
false)
18010 unsigned int n_grandchildren = 0;
18013 unsigned int n_coarsen_flags = 0;
18018 for (
unsigned int c = 0; c < n_children; ++c)
18023 cell_iterator child = cell->
child(c);
18025 const unsigned int nn_children = child->n_children();
18026 n_grandchildren += nn_children;
18031 if (child->child(0)->is_active())
18032 for (
unsigned int cc = 0; cc < nn_children; ++cc)
18033 if (child->child(cc)->coarsen_flag_set())
18048 if ((n_coarsen_flags != n_grandchildren) && (n_coarsen_flags > 0))
18049 for (
unsigned int c = 0; c < n_children; ++c)
18051 const cell_iterator child = cell->
child(c);
18052 if (child->child(0)->is_active())
18053 for (
unsigned int cc = 0; cc < child->n_children(); ++cc)
18054 child->child(cc)->clear_coarsen_flag();
18065 this->policy->prevent_distorted_boundary_cells(*
this);
18081 bool changed =
true;
18085 active_cell_iterator cell = last_active(), endc =
end();
18087 for (; cell != endc; --cell)
18088 if (cell->refine_flag_set())
18091 for (
const auto i : cell->face_indices())
18096 const bool has_periodic_neighbor =
18097 cell->has_periodic_neighbor(i);
18098 const bool has_neighbor_or_periodic_neighbor =
18099 !cell->at_boundary(i) || has_periodic_neighbor;
18100 if (has_neighbor_or_periodic_neighbor &&
18102 cell->refine_flag_set(), i) !=
18115 if (cell->neighbor_or_periodic_neighbor(i)->is_active())
18117 if ((!has_periodic_neighbor &&
18118 cell->neighbor_is_coarser(i)) ||
18119 (has_periodic_neighbor &&
18120 cell->periodic_neighbor_is_coarser(i)))
18122 if (cell->neighbor_or_periodic_neighbor(i)
18123 ->coarsen_flag_set())
18124 cell->neighbor_or_periodic_neighbor(i)
18125 ->clear_coarsen_flag();
18142 allow_anisotropic_smoothing)
18144 has_periodic_neighbor ?
18145 cell->periodic_neighbor(i)
18146 ->flag_for_face_refinement(
18148 ->periodic_neighbor_of_coarser_periodic_neighbor(
18151 RefinementCase<dim - 1>::cut_x) :
18153 ->flag_for_face_refinement(
18155 ->neighbor_of_coarser_neighbor(
18158 RefinementCase<dim - 1>::cut_x);
18162 ->neighbor_or_periodic_neighbor(
18164 ->refine_flag_set())
18166 cell->neighbor_or_periodic_neighbor(i)
18167 ->set_refine_flag();
18299 std::pair<unsigned int, unsigned int>
18301 has_periodic_neighbor ?
18303 ->periodic_neighbor_of_coarser_periodic_neighbor(
18305 cell->neighbor_of_coarser_neighbor(i);
18306 unsigned int refined_along_x = 0,
18307 refined_along_y = 0,
18308 to_be_refined_along_x = 0,
18309 to_be_refined_along_y = 0;
18311 const int this_face_index =
18312 cell->face_index(i);
18321 const auto parent_face = [&]() {
18322 if (has_periodic_neighbor)
18324 const auto neighbor =
18325 cell->periodic_neighbor(i);
18326 const auto parent_face_no =
18328 ->periodic_neighbor_of_periodic_neighbor(
18331 neighbor->periodic_neighbor(
18333 return parent->face(parent_face_no);
18336 return cell->neighbor(i)->face(
18340 if ((this_face_index ==
18341 parent_face->child_index(0)) ||
18342 (this_face_index ==
18343 parent_face->child_index(1)))
18351 RefinementCase<dim - 1> frc =
18352 parent_face->refinement_case();
18353 if (frc & RefinementCase<dim>::cut_x)
18355 if (frc & RefinementCase<dim>::cut_y)
18368 RefinementCase<dim - 1> flagged_frc =
18370 cell->refine_flag_set(),
18372 cell->face_orientation(i),
18373 cell->face_flip(i),
18374 cell->face_rotation(i));
18376 RefinementCase<dim>::cut_x)
18377 ++to_be_refined_along_x;
18379 RefinementCase<dim>::cut_y)
18380 ++to_be_refined_along_y;
18385 allow_anisotropic_smoothing) ||
18386 cell->neighbor_or_periodic_neighbor(i)
18387 ->refine_flag_set())
18389 if (refined_along_x +
18390 to_be_refined_along_x >
18394 ->neighbor_or_periodic_neighbor(i)
18395 ->flag_for_face_refinement(
18397 RefinementCase<dim -
18399 if (refined_along_y +
18400 to_be_refined_along_y >
18404 ->neighbor_or_periodic_neighbor(i)
18405 ->flag_for_face_refinement(
18407 RefinementCase<dim -
18413 ->neighbor_or_periodic_neighbor(i)
18414 ->refine_flag_set() !=
18418 cell->neighbor_or_periodic_neighbor(i)
18419 ->set_refine_flag();
18425 cell->neighbor_or_periodic_neighbor(i);
18426 RefinementCase<dim - 1> nb_frc =
18428 nb->refine_flag_set(),
18430 nb->face_orientation(nb_indices.first),
18431 nb->face_flip(nb_indices.first),
18432 nb->face_rotation(nb_indices.first));
18433 if ((nb_frc & RefinementCase<dim>::cut_x) &&
18434 !((refined_along_x != 0u) ||
18435 (to_be_refined_along_x != 0u)))
18436 changed |= cell->flag_for_face_refinement(
18439 if ((nb_frc & RefinementCase<dim>::cut_y) &&
18440 !((refined_along_y != 0u) ||
18441 (to_be_refined_along_y != 0u)))
18442 changed |= cell->flag_for_face_refinement(
18449 cell->neighbor_or_periodic_neighbor(i)
18450 ->clear_coarsen_flag();
18451 const unsigned int nb_nb =
18452 has_periodic_neighbor ?
18454 ->periodic_neighbor_of_periodic_neighbor(
18456 cell->neighbor_of_neighbor(i);
18457 const cell_iterator neighbor =
18458 cell->neighbor_or_periodic_neighbor(i);
18459 RefinementCase<dim - 1> face_ref_case =
18461 neighbor->refine_flag_set(),
18463 neighbor->face_orientation(nb_nb),
18464 neighbor->face_flip(nb_nb),
18465 neighbor->face_rotation(nb_nb));
18466 RefinementCase<dim - 1> needed_face_ref_case =
18468 cell->refine_flag_set(),
18470 cell->face_orientation(i),
18471 cell->face_flip(i),
18472 cell->face_rotation(i));
18477 if ((face_ref_case ==
18478 RefinementCase<dim>::cut_x &&
18479 needed_face_ref_case ==
18480 RefinementCase<dim>::cut_y) ||
18482 RefinementCase<dim>::cut_y &&
18483 needed_face_ref_case ==
18484 RefinementCase<dim>::cut_x))
18486 changed = cell->flag_for_face_refinement(
18488 neighbor->flag_for_face_refinement(
18489 nb_nb, needed_face_ref_case);
18495 RefinementCase<dim - 1>
18496 face_ref_case = cell->face(i)->refinement_case(),
18497 needed_face_ref_case =
18499 cell->refine_flag_set(),
18501 cell->face_orientation(i),
18502 cell->face_flip(i),
18503 cell->face_rotation(i));
18507 if ((face_ref_case == RefinementCase<dim>::cut_x &&
18508 needed_face_ref_case ==
18509 RefinementCase<dim>::cut_y) ||
18510 (face_ref_case == RefinementCase<dim>::cut_y &&
18511 needed_face_ref_case ==
18512 RefinementCase<dim>::cut_x))
18514 cell->flag_for_face_refinement(i,
18526 this->policy->prepare_refinement_dim_dependent(*
this);
18532 fix_coarsen_flags();
18535 auto coarsen_flags_after_loop =
18536 internal::extract_raw_coarsen_flags(levels);
18537 auto refine_flags_after_loop = internal::extract_raw_refine_flags(levels);
18540 mesh_changed_in_this_loop =
18541 ((coarsen_flags_before_loop != coarsen_flags_after_loop) ||
18542 (refine_flags_before_loop != refine_flags_after_loop));
18545 coarsen_flags_before_loop.swap(coarsen_flags_after_loop);
18546 refine_flags_before_loop.swap(refine_flags_after_loop);
18548 while (mesh_changed_in_this_loop);
18554 return ((coarsen_flags_before != coarsen_flags_before_loop) ||
18555 (refine_flags_before != refine_flags_before_loop));
18560template <
int dim,
int spacedim>
18563 const unsigned int magic_number1,
18564 const std::vector<bool> &v,
18565 const unsigned int magic_number2,
18568 const unsigned int N = v.size();
18569 unsigned char *flags =
new unsigned char[
N / 8 + 1];
18570 for (
unsigned int i = 0; i <
N / 8 + 1; ++i)
18573 for (
unsigned int position = 0; position <
N; ++position)
18574 flags[position / 8] |= (v[position] ? (1 << (position % 8)) : 0);
18583 out << magic_number1 <<
' ' <<
N << std::endl;
18584 for (
unsigned int i = 0; i <
N / 8 + 1; ++i)
18585 out <<
static_cast<unsigned int>(flags[i]) <<
' ';
18587 out << std::endl << magic_number2 << std::endl;
18595template <
int dim,
int spacedim>
18598 const unsigned int magic_number1,
18599 std::vector<bool> &v,
18600 const unsigned int magic_number2,
18605 unsigned int magic_number;
18606 in >> magic_number;
18607 AssertThrow(magic_number == magic_number1, ExcGridReadError());
18613 unsigned char *flags =
new unsigned char[
N / 8 + 1];
18614 unsigned short int tmp;
18615 for (
unsigned int i = 0; i <
N / 8 + 1; ++i)
18621 for (
unsigned int position = 0; position !=
N; ++position)
18622 v[position] = ((flags[position / 8] & (1 << (position % 8))) != 0);
18624 in >> magic_number;
18625 AssertThrow(magic_number == magic_number2, ExcGridReadError());
18634template <
int dim,
int spacedim>
18638 std::size_t mem = 0;
18639 mem +=
sizeof(MeshSmoothing);
18643 for (
const auto &level : levels)
18647 mem +=
sizeof(default_flat_manifold) +
sizeof(*default_flat_manifold);
18648 mem +=
sizeof(manifolds);
18649 mem +=
sizeof(smooth_grid);
18651 mem +=
sizeof(faces);
18660template <
int dim,
int spacedim>
18668#include "grid/tria.inst"