17 #include <deal.II/base/memory_consumption.h> 18 #include <deal.II/base/table.h> 19 #include <deal.II/base/geometry_info.h> 20 #include <deal.II/base/std_cxx11/bind.h> 22 #include <deal.II/grid/tria.h> 23 #include <deal.II/grid/tria_levels.h> 24 #include <deal.II/grid/tria_faces.h> 25 #include <deal.II/grid/manifold.h> 26 #include <deal.II/grid/tria_boundary.h> 27 #include <deal.II/grid/tria_accessor.h> 28 #include <deal.II/grid/tria_iterator.h> 29 #include <deal.II/grid/grid_tools.h> 30 #include <deal.II/grid/magic_numbers.h> 31 #include <deal.II/fe/mapping_q1.h> 32 #include <deal.II/lac/vector.h> 33 #include <deal.II/lac/full_matrix.h> 42 #include <deal.II/base/std_cxx11/array.h> 44 DEAL_II_NAMESPACE_OPEN
139 template <
int dim,
int spacedim>
142 Assert (cell->active() ==
false, ExcInternalError());
144 unsigned int n_active_children = 0;
145 for (
unsigned int i=0; i<cell->n_children(); ++i)
146 if (cell->child(i)->active())
149 return (n_active_children == 0) || (n_active_children == cell->n_children());
159 template <
int dim,
int spacedim>
165 if (cell->has_children())
167 unsigned int children_to_coarsen=0;
168 const unsigned int n_children=cell->n_children();
170 for (
unsigned int c=0; c<n_children; ++c)
171 if (cell->child(c)->active() &&
172 cell->child(c)->coarsen_flag_set())
173 ++children_to_coarsen;
174 if (children_to_coarsen==n_children)
177 for (
unsigned int c=0; c<n_children; ++c)
178 if (cell->child(c)->active())
179 cell->child(c)->clear_coarsen_flag();
189 Assert(cell->has_children(), ExcInternalError());
215 template <
int dim,
int spacedim>
218 const unsigned int face_no,
233 if (neighbor->has_children())
238 if (cell_will_be_coarsened(neighbor))
247 expected_face_ref_case=cell->face(face_no)->refinement_case();
259 const unsigned int neighbor_neighbor=cell->neighbor_face_no(face_no);
266 neighbor->face_orientation(neighbor_neighbor),
267 neighbor->face_flip(neighbor_neighbor),
268 neighbor->face_rotation(neighbor_neighbor));
272 const int this_face_index=cell->face_index(face_no);
278 if (neighbor_face->index()==this_face_index)
284 expected_face_ref_case = face_ref_case;
303 for (
unsigned int c=0; c<neighbor_face->n_children(); ++c)
304 if (neighbor_face->child_index(c)==this_face_index)
313 if ((neighbor_face->refinement_case() | face_ref_case)
314 == neighbor_face->refinement_case())
330 expected_face_ref_case = ~neighbor_face->refinement_case();
339 Assert(dim==3, ExcInternalError());
357 template <
int dim,
int spacedim>
360 const unsigned int face_no)
363 return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
370 template <
int dim,
int spacedim>
373 const unsigned int face_no,
376 return face_will_be_refined_by_neighbor_internal(cell, face_no,
377 expected_face_ref_case);
382 template <
int dim,
int spacedim>
386 std::vector<unsigned int> min_adjacent_cell_level (triangulation.
n_vertices(),
388 std::vector<unsigned int> max_adjacent_cell_level (triangulation.
n_vertices(),
393 cell != triangulation.
end(); ++cell)
394 for (
unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
396 min_adjacent_cell_level[cell->vertex_index(v)]
397 = std::min<unsigned int>
398 (min_adjacent_cell_level[cell->vertex_index(v)],
400 max_adjacent_cell_level[cell->vertex_index(v)]
401 = std::max<unsigned int> (min_adjacent_cell_level[cell->vertex_index(v)],
405 for (
unsigned int k=0; k<triangulation.
n_vertices(); ++k)
407 if (max_adjacent_cell_level[k] -
408 min_adjacent_cell_level[k] > 1)
421 template <
int dim,
int spacedim>
422 std::vector<unsigned int>
427 std::vector<unsigned int> line_cell_count(triangulation.
n_raw_lines(),0);
429 cell=triangulation.
begin(),
430 endc=triangulation.
end();
431 for (; cell!=endc; ++cell)
432 for (
unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
433 ++line_cell_count[cell->line_index(l)];
434 return line_cell_count;
437 return std::vector<unsigned int>();
449 template <
int dim,
int spacedim>
450 std::vector<unsigned int>
455 std::vector<unsigned int> quad_cell_count (triangulation.
n_raw_quads(),0);
457 cell=triangulation.
begin(),
458 endc=triangulation.
end();
459 for (; cell!=endc; ++cell)
460 for (
unsigned int q=0; q<GeometryInfo<dim>::faces_per_cell; ++q)
461 ++quad_cell_count[cell->quad_index(q)];
462 return quad_cell_count;
465 return std::vector<unsigned int>();
482 reorder_compatibility (
const std::vector<
CellData<1> > &,
491 reorder_compatibility (std::vector<
CellData<2> > &cells,
494 for (
unsigned int cell=0; cell<cells.size(); ++cell)
495 std::swap(cells[cell].vertices[2],cells[cell].vertices[3]);
500 reorder_compatibility (std::vector<
CellData<3> > &cells,
504 for (
unsigned int cell=0; cell<cells.size(); ++cell)
506 for (
unsigned int i=0; i<GeometryInfo<3>::vertices_per_cell; ++i)
507 tmp[i] = cells[cell].vertices[i];
508 for (
unsigned int i=0; i<GeometryInfo<3>::vertices_per_cell; ++i)
513 std::vector<CellData<2> >::iterator boundary_quad
515 std::vector<CellData<2> >::iterator end_quad
517 for (
unsigned int quad_no=0; boundary_quad!=end_quad; ++boundary_quad, ++quad_no)
518 std::swap(boundary_quad->vertices[2], boundary_quad->vertices[3]);
540 template <
int dim,
int spacedim>
542 middle_vertex_index(
const typename Triangulation<dim,spacedim>::line_iterator &line)
544 if (line->has_children())
545 return line->child(0)->vertex_index(1);
550 template <
int dim,
int spacedim>
552 middle_vertex_index(
const typename Triangulation<dim,spacedim>::quad_iterator &quad)
554 switch (static_cast<unsigned char> (quad->refinement_case()))
557 return middle_vertex_index<dim,spacedim>(quad->child(0)->line(1));
560 return middle_vertex_index<dim,spacedim>(quad->child(0)->line(3));
563 return quad->child(0)->vertex_index(3);
572 template <
int dim,
int spacedim>
574 middle_vertex_index(
const typename Triangulation<dim,spacedim>::hex_iterator &hex)
576 switch (static_cast<unsigned char> (hex->refinement_case()))
579 return middle_vertex_index<dim,spacedim>(hex->child(0)->quad(1));
582 return middle_vertex_index<dim,spacedim>(hex->child(0)->quad(3));
585 return middle_vertex_index<dim,spacedim>(hex->child(0)->quad(5));
588 return middle_vertex_index<dim,spacedim>(hex->child(0)->line(11));
591 return middle_vertex_index<dim,spacedim>(hex->child(0)->line(5));
594 return middle_vertex_index<dim,spacedim>(hex->child(0)->line(7));
597 return hex->child(0)->vertex_index(7);
618 template <
class TRIANGULATION>
620 typename TRIANGULATION::DistortedCellList
621 collect_distorted_coarse_cells (
const TRIANGULATION &)
623 return typename TRIANGULATION::DistortedCellList();
643 cell = triangulation.
begin(0); cell != triangulation.
end(0); ++cell)
646 for (
unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
647 vertices[i] = cell->vertex(i);
653 for (
unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
654 if (determinants[i] <= 1e-9 * std::pow (cell->diameter(),
657 distorted_cells.distorted_cells.push_back (cell);
662 return distorted_cells;
678 Assert (cell->has_children(), ExcInternalError());
680 for (
unsigned int c=0; c<cell->n_children(); ++c)
683 for (
unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
684 vertices[i] = cell->child(c)->vertex(i);
690 for (
unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
691 if (determinants[i] <= 1e-9 * std::pow (cell->child(c)->diameter(),
707 template <
int dim,
int spacedim>
722 template <
int spacedim>
729 template <
int dim,
int spacedim>
784 static const unsigned int left_right_offset[2][6][2] =
813 std::vector<typename Triangulation<dim,spacedim>::cell_iterator>
814 adjacent_cells(2*triangulation.
n_raw_faces(), dummy);
817 cell = triangulation.
begin(),
818 endc = triangulation.
end();
819 for (; cell != endc; ++cell)
820 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
826 offset = (cell->direction_flag()
828 left_right_offset[dim-2][f][cell->face_orientation(f)]
830 1-left_right_offset[dim-2][f][cell->face_orientation(f)]);
832 adjacent_cells[2*face->index() + offset] = cell;
842 if (cell->active() && face->has_children())
844 adjacent_cells[2*face->child(0)->index() + offset] = cell;
845 adjacent_cells[2*face->child(1)->index() + offset] = cell;
870 if (face->has_children() &&
875 for (
unsigned int c=0; c<face->n_children(); ++c)
876 adjacent_cells[2*face->child(c)->index() + offset] = cell;
877 if (face->child(0)->has_children())
879 adjacent_cells[2*face->child(0)->child(0)->index() + offset] = cell;
880 adjacent_cells[2*face->child(0)->child(1)->index() + offset] = cell;
882 if (face->child(1)->has_children())
884 adjacent_cells[2*face->child(1)->child(0)->index() + offset] = cell;
885 adjacent_cells[2*face->child(1)->child(1)->index() + offset] = cell;
896 for (cell=triangulation.
begin(); cell != endc; ++cell)
897 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
900 offset = (cell->direction_flag()
902 left_right_offset[dim-2][f][cell->face_orientation(f)]
904 1-left_right_offset[dim-2][f][cell->face_orientation(f)]);
905 cell->set_neighbor(f,
906 adjacent_cells[2*cell->face(f)->index() + 1 - offset]);
923 using ::Triangulation;
931 <<
"Something went wrong when making cell " << arg1
932 <<
". Read the docs and the source code " 933 <<
"for more information.");
940 <<
"Something went wrong upon construction of cell " 953 <<
"Cell " << arg1 <<
" has negative measure. This typically " 954 <<
"indicates some distortion in the cell, or a mistakenly " 955 <<
"swapped pair of vertices in the input to " 956 <<
"Triangulation::create_triangulation().");
966 <<
"Error while creating cell " << arg1
967 <<
": the vertex index " << arg2 <<
" must be between 0 and " 975 <<
"While trying to assign a boundary indicator to a line: " 976 <<
"the line with end vertices " << arg1 <<
" and " 977 << arg2 <<
" does not exist.");
984 <<
"While trying to assign a boundary indicator to a quad: " 985 <<
"the quad with bounding lines " << arg1 <<
", " << arg2
986 <<
", " << arg3 <<
", " << arg4 <<
" does not exist.");
994 <<
"The input data for creating a triangulation contained " 995 <<
"information about a line with indices " 996 << arg1 <<
" and " << arg2
997 <<
" that is described to have boundary indicator " 999 <<
". However, this is an internal line not located on the " 1000 <<
"boundary. You cannot assign a boundary indicator to it." 1003 <<
"If this happened at a place where you call " 1004 <<
"Triangulation::create_triangulation() yourself, you need " 1005 <<
"to check the SubCellData object you pass to this function." 1008 <<
"If this happened in a place where you are reading a mesh " 1009 <<
"from a file, then you need to investigate why such a line " 1010 <<
"ended up in the input file. A typical case is a geometry " 1011 <<
"that consisted of multiple parts and for which the mesh " 1012 <<
"generator program assumes that the interface between " 1013 <<
"two parts is a boundary when that isn't supposed to be " 1014 <<
"the case, or where the mesh generator simply assigns " 1015 <<
"'geometry indicators' to lines at the perimeter of " 1016 <<
"a part that are not supposed to be interpreted as " 1017 <<
"'boundary indicators'.");
1025 <<
"The input data for creating a triangulation contained " 1026 <<
"information about a quad with indices " 1027 << arg1 <<
", " << arg2 <<
", " << arg3 <<
", and " << arg4
1028 <<
" that is described to have boundary indicator " 1030 <<
". However, this is an internal quad not located on the " 1031 <<
"boundary. You cannot assign a boundary indicator to it." 1034 <<
"If this happened at a place where you call " 1035 <<
"Triangulation::create_triangulation() yourself, you need " 1036 <<
"to check the SubCellData object you pass to this function." 1039 <<
"If this happened in a place where you are reading a mesh " 1040 <<
"from a file, then you need to investigate why such a quad " 1041 <<
"ended up in the input file. A typical case is a geometry " 1042 <<
"that consisted of multiple parts and for which the mesh " 1043 <<
"generator program assumes that the interface between " 1044 <<
"two parts is a boundary when that isn't supposed to be " 1045 <<
"the case, or where the mesh generator simply assigns " 1046 <<
"'geometry indicators' to quads at the surface of " 1047 <<
"a part that are not supposed to be interpreted as " 1048 <<
"'boundary indicators'.");
1055 <<
"In SubCellData the line info of the line with vertex indices " 1056 << arg1 <<
" and " << arg2 <<
" appears more than once. " 1057 <<
"This is not allowed.");
1172 template <
int dim,
int spacedim>
1175 const unsigned int level_objects,
1179 typename Triangulation<dim,spacedim>::line_iterator line_iterator;
1181 typename Triangulation<dim,spacedim>::active_line_iterator active_line_iterator;
1184 if (level_objects > 0)
1188 for (
unsigned int level=0; level<level_objects; ++level)
1189 if (triangulation.
begin(level) !=
1190 triangulation.
end(level))
1211 for (
unsigned int level=0; level<number_cache.
n_levels; ++level)
1216 line_iterator line = triangulation.
begin_line (level),
1217 endc = (level == number_cache.
n_levels-1 ?
1218 line_iterator(triangulation.
end_line()) :
1220 for (; line!=endc; ++line)
1229 for (
unsigned int level=0; level<number_cache.
n_levels; ++level)
1236 for (; (line!=endc) && (line->level() ==
static_cast<signed int>(level)); ++line)
1248 line_iterator line = triangulation.
begin_line (),
1250 for (; line!=endc; ++line)
1257 for (; line!=endc; ++line)
1282 template <
int dim,
int spacedim>
1285 const unsigned int level_objects,
1289 compute_number_cache (triangulation,
1295 typename Triangulation<dim,spacedim>::quad_iterator quad_iterator;
1297 typename Triangulation<dim,spacedim>::active_quad_iterator active_quad_iterator;
1314 for (
unsigned int level=0; level<number_cache.
n_levels; ++level)
1319 quad_iterator quad = triangulation.
begin_quad (level),
1320 endc = (level == number_cache.
n_levels-1 ?
1321 quad_iterator(triangulation.
end_quad()) :
1323 for (; quad!=endc; ++quad)
1332 for (
unsigned int level=0; level<number_cache.
n_levels; ++level)
1339 for (; (quad!=endc) && (quad->level() ==
static_cast<signed int>(level)); ++quad)
1351 quad_iterator quad = triangulation.
begin_quad (),
1353 for (; quad!=endc; ++quad)
1360 for (; quad!=endc; ++quad)
1386 template <
int dim,
int spacedim>
1389 const unsigned int level_objects,
1393 compute_number_cache (triangulation,
1399 typename Triangulation<dim,spacedim>::hex_iterator hex_iterator;
1401 typename Triangulation<dim,spacedim>::active_hex_iterator active_hex_iterator;
1418 for (
unsigned int level=0; level<number_cache.
n_levels; ++level)
1423 hex_iterator hex = triangulation.
begin_hex (level),
1424 endc = (level == number_cache.
n_levels-1 ?
1425 hex_iterator(triangulation.
end_hex()) :
1427 for (; hex!=endc; ++hex)
1436 for (
unsigned int level=0; level<number_cache.
n_levels; ++level)
1442 endc = triangulation.
end_hex ();
1443 for (; (hex!=endc) && (hex->level() ==
static_cast<signed int>(level)); ++hex)
1455 hex_iterator hex = triangulation.
begin_hex (),
1456 endc = triangulation.
end_hex();
1457 for (; hex!=endc; ++hex)
1463 endc = triangulation.
end_hex();
1464 for (; hex!=endc; ++hex)
1478 template <
int spacedim>
1495 const unsigned int dim=1;
1499 triangulation.
vertices_used = std::vector<bool> (v.size(),
true);
1504 std::vector<std::vector<int> > lines_at_vertex (v.size());
1508 triangulation.
levels[0]->reserve_space (cells.size(), dim, spacedim);
1509 triangulation.
levels[0]->cells.reserve_space (0,cells.size());
1512 typename Triangulation<dim,spacedim>::raw_line_iterator
1514 for (
unsigned int cell=0; cell<cells.size(); ++cell)
1516 while (next_free_line->used())
1521 cells[cell].vertices[1]));
1522 next_free_line->set_used_flag ();
1523 next_free_line->set_material_id (cells[cell].material_id);
1524 next_free_line->set_manifold_id (cells[cell].manifold_id);
1526 next_free_line->set_subdomain_id (0);
1530 lines_at_vertex[cells[cell].vertices[0]].push_back (cell);
1531 lines_at_vertex[cells[cell].vertices[1]].push_back (cell);
1537 unsigned int boundary_nodes = 0;
1538 for (
unsigned int i=0; i<lines_at_vertex.size(); ++i)
1539 switch (lines_at_vertex[i].size())
1554 ExcMessage (
"You have a vertex in your triangulation " 1555 "at which more than two cells come together. " 1556 "(For one dimensional triangulation, cells are " 1559 "This is not currently supported because the " 1560 "Triangulation class makes the assumption that " 1561 "every cell has zero or one neighbors behind " 1562 "each face (here, behind each vertex), but in your " 1563 "situation there would be more than one." 1565 "Support for this is not currently implemented. " 1566 "If you need to work with triangulations where " 1567 "more than two cells come together at a vertex, " 1568 "duplicate the vertices once per cell (i.e., put " 1569 "multiple vertices at the same physical location, " 1570 "but using different vertex indices for each) " 1571 "and then ensure continuity of the solution by " 1572 "explicitly creating constraints that the degrees " 1573 "of freedom at these vertices have the same " 1574 "value, using the ConstraintMatrix class."));
1586 AssertThrow (((spacedim == 1) && (boundary_nodes == 2))
1589 ExcMessage(
"The Triangulation has too many end points"));
1595 typename Triangulation<dim,spacedim>::active_line_iterator
1598 for (; line!=triangulation.
end(); ++line)
1600 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
1606 if (lines_at_vertex[line->vertex_index(vertex)][0] == line->index())
1607 if (lines_at_vertex[line->vertex_index(vertex)].size() == 2)
1610 neighbor (&triangulation,
1612 lines_at_vertex[line->vertex_index(vertex)][1]);
1613 line->set_neighbor (vertex, neighbor);
1619 line->set_neighbor (vertex, triangulation.
end());
1626 neighbor (&triangulation,
1628 lines_at_vertex[line->vertex_index(vertex)][0]);
1629 line->set_neighbor (vertex, neighbor);
1640 cell != triangulation.
end(); ++cell)
1641 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
1647 if (cell->at_boundary(f))
1649 .vertex_to_boundary_id_map_1d)[cell->face(f)->vertex_index()]
1662 template <
int spacedim>
1673 const unsigned int dim=2;
1677 triangulation.
vertices_used = std::vector<bool> (v.size(),
true);
1692 std::map<std::pair<int,int>,
1693 typename Triangulation<dim,spacedim>::line_iterator> needed_lines;
1694 for (
unsigned int cell=0; cell<cells.size(); ++cell)
1696 for (
unsigned int vertex=0; vertex<4; ++vertex)
1698 ExcInvalidVertexIndex (cell, cells[cell].vertices[vertex],
1701 for (
unsigned int line=0; line<GeometryInfo<dim>::faces_per_cell; ++line)
1709 std::pair<int,int> line_vertices(
1739 AssertThrow (needed_lines.find(std::make_pair(line_vertices.second,
1740 line_vertices.first))
1743 ExcGridHasInvalidCell(cell));
1749 needed_lines[line_vertices] = triangulation.
end_line();
1757 std::vector<unsigned short int> vertex_touch_count (v.size(), 0);
1758 typename std::map<std::pair<int,int>,
1759 typename Triangulation<dim,spacedim>::line_iterator>::iterator i;
1760 for (i=needed_lines.begin(); i!=needed_lines.end(); i++)
1764 ++vertex_touch_count[i->first.first];
1765 ++vertex_touch_count[i->first.second];
1772 AssertThrow (* (std::min_element(vertex_touch_count.begin(),
1773 vertex_touch_count.end())) >= 2,
1774 ExcMessage(
"During creation of a triangulation, a part of the " 1775 "algorithm encountered a vertex that is part of only " 1776 "a single adjacent line. However, in 2d, every vertex " 1777 "needs to be at least part of two lines."));
1783 triangulation.
levels[0]->reserve_space (cells.size(), dim, spacedim);
1784 triangulation.
faces->lines.reserve_space (0,needed_lines.size());
1785 triangulation.
levels[0]->cells.reserve_space (0,cells.size());
1789 typename Triangulation<dim,spacedim>::raw_line_iterator
1791 typename std::map<std::pair<int,int>,
1792 typename Triangulation<dim,spacedim>::line_iterator>::iterator i;
1793 for (i = needed_lines.begin();
1794 line!=triangulation.
end_line(); ++line, ++i)
1798 line->set_used_flag ();
1799 line->clear_user_flag ();
1808 std::map<int,std::vector<typename Triangulation<dim,spacedim>::cell_iterator> >
1815 for (
unsigned int c=0; c<cells.size(); ++c, ++cell)
1817 typename Triangulation<dim,spacedim>::line_iterator
1819 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
1820 lines[line]=needed_lines[std::make_pair(
1827 lines[3]->index()));
1829 cell->set_used_flag ();
1830 cell->set_material_id (cells[c].material_id);
1831 cell->set_manifold_id (cells[c].manifold_id);
1832 cell->clear_user_data ();
1833 cell->set_subdomain_id (0);
1838 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
1839 adjacent_cells[lines[line]->index()].push_back (cell);
1844 for (
typename Triangulation<dim,spacedim>::line_iterator
1846 line!=triangulation.
end_line(); ++line)
1848 const unsigned int n_adj_cells = adjacent_cells[line->index()].size();
1861 ExcMessage (
"You have a line in your triangulation " 1862 "at which more than two cells come together. " 1864 "This is not currently supported because the " 1865 "Triangulation class makes the assumption that " 1866 "every cell has zero or one neighbors behind " 1867 "each face (here, behind each line), but in your " 1868 "situation there would be more than one." 1870 "Support for this is not currently implemented. " 1871 "If you need to work with triangulations where " 1872 "more than two cells come together at a line, " 1873 "duplicate the vertices once per cell (i.e., put " 1874 "multiple vertices at the same physical location, " 1875 "but using different vertex indices for each) " 1876 "and then ensure continuity of the solution by " 1877 "explicitly creating constraints that the degrees " 1878 "of freedom at these lines have the same " 1879 "value, using the ConstraintMatrix class."));
1885 if (n_adj_cells == 1)
1886 line->set_boundary_id (0);
1894 std::vector<CellData<1> >::const_iterator boundary_line
1896 std::vector<CellData<1> >::const_iterator end_boundary_line
1898 for (; boundary_line!=end_boundary_line; ++boundary_line)
1900 typename Triangulation<dim,spacedim>::line_iterator line;
1901 std::pair<int,int> line_vertices(std::make_pair(boundary_line->vertices[0],
1902 boundary_line->vertices[1]));
1903 if (needed_lines.find(line_vertices) != needed_lines.end())
1905 line = needed_lines[line_vertices];
1909 std::swap (line_vertices.first, line_vertices.second);
1910 if (needed_lines.find(line_vertices) != needed_lines.end())
1911 line = needed_lines[line_vertices];
1914 AssertThrow (
false, ExcLineInexistant(line_vertices.first,
1915 line_vertices.second));
1921 ExcMultiplySetLineInfoOfLine(line_vertices.first,
1922 line_vertices.second));
1932 ExcInteriorLineCantBeBoundary(line->vertex_index(0),
1933 line->vertex_index(1),
1934 boundary_line->boundary_id));
1935 line->set_boundary_id (boundary_line->boundary_id);
1938 line->set_manifold_id (boundary_line->manifold_id);
1944 cell=triangulation.
begin(); cell!=triangulation.
end(); ++cell)
1945 for (
unsigned int side=0; side<4; ++side)
1946 if (adjacent_cells[cell->line(side)->index()][0] == cell)
1950 if (adjacent_cells[cell->line(side)->index()].size() == 2)
1953 cell->set_neighbor (side,
1954 adjacent_cells[cell->line(side)->index()][1]);
1960 cell->set_neighbor (side,
1961 adjacent_cells[cell->line(side)->index()][0]);
2010 template <
int spacedim>
2021 const unsigned int dim=3;
2025 triangulation.
vertices_used = std::vector<bool> (v.size(),
true);
2035 for (
unsigned int cell_no = 0; cell_no<cells.size(); ++cell_no)
2037 cells[cell_no].vertices) >= 0,
2038 ExcGridHasInvalidCell(cell_no));
2060 typename std::map<std::pair<int,int>,
2061 typename Triangulation<dim,spacedim>::line_iterator> needed_lines;
2062 for (
unsigned int cell=0; cell<cells.size(); ++cell)
2066 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
2068 ExcInvalidVertexIndex (cell, cells[cell].vertices[vertex],
2071 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
2079 std::pair<int,int> line_vertices(
2086 if ( (needed_lines.find(std::make_pair(line_vertices.second,
2087 line_vertices.first))
2089 needed_lines.end()))
2095 needed_lines[line_vertices] = triangulation.
end_line();
2107 std::vector<unsigned short int> vertex_touch_count (v.size(), 0);
2108 typename std::map<std::pair<int,int>,
2109 typename Triangulation<dim,spacedim>::line_iterator>::iterator i;
2110 for (i=needed_lines.begin(); i!=needed_lines.end(); i++)
2114 ++vertex_touch_count[i->first.first];
2115 ++vertex_touch_count[i->first.second];
2122 AssertThrow (* (std::min_element(vertex_touch_count.begin(),
2123 vertex_touch_count.end())) >= 3,
2124 ExcMessage(
"During creation of a triangulation, a part of the " 2125 "algorithm encountered a vertex that is part of only " 2126 "one or two adjacent lines. However, in 3d, every vertex " 2127 "needs to be at least part of three lines."));
2137 triangulation.
levels[0]->reserve_space (cells.size(), dim, spacedim);
2138 triangulation.
faces->lines.reserve_space (0,needed_lines.size());
2142 typename Triangulation<dim,spacedim>::raw_line_iterator
2144 typename std::map<std::pair<int,int>,
2145 typename Triangulation<dim,spacedim>::line_iterator>::iterator i;
2146 for (i = needed_lines.begin(); line!=triangulation.
end_line(); ++line, ++i)
2150 line->set_used_flag ();
2151 line->clear_user_flag ();
2175 std::map<internal::Triangulation::TriaObject<2>,
2176 std::pair<typename Triangulation<dim,spacedim>::quad_iterator,
2177 std_cxx11::array<bool,GeometryInfo<dim>::lines_per_face> >,
2180 for (
unsigned int cell=0; cell<cells.size(); ++cell)
2208 std_cxx11::array<bool,GeometryInfo<dim>::lines_per_face> orientation;
2210 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
2212 line_list[line]=std::pair<int,int> (
2215 inverse_line_list[line]=std::pair<int,int> (
2220 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2230 for (
unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
2232 face_to_cell_lines(face,l)]) == needed_lines.end())
2236 orientation[l]=
true;
2242 orientation[l]=
false;
2247 quad(face_line_list[0],
2291 test_quad_1(quad.
face(2), quad.
face(3),
2293 test_quad_2(quad.
face(0), quad.
face(1),
2295 test_quad_3(quad.
face(3), quad.
face(2),
2297 test_quad_4(quad.
face(1), quad.
face(0),
2299 test_quad_5(quad.
face(2), quad.
face(3),
2301 test_quad_6(quad.
face(1), quad.
face(0),
2303 test_quad_7(quad.
face(3), quad.
face(2),
2305 if (needed_quads.find (test_quad_1) == needed_quads.end() &&
2306 needed_quads.find (test_quad_2) == needed_quads.end() &&
2307 needed_quads.find (test_quad_3) == needed_quads.end() &&
2308 needed_quads.find (test_quad_4) == needed_quads.end() &&
2309 needed_quads.find (test_quad_5) == needed_quads.end() &&
2310 needed_quads.find (test_quad_6) == needed_quads.end() &&
2311 needed_quads.find (test_quad_7) == needed_quads.end())
2312 needed_quads[quad] = std::make_pair(triangulation.
end_quad(),orientation);
2322 triangulation.
faces->quads.reserve_space (0,needed_quads.size());
2325 typename Triangulation<dim,spacedim>::raw_quad_iterator
2327 typename std::map<internal::Triangulation::TriaObject<2>,
2328 std::pair<typename Triangulation<dim,spacedim>::quad_iterator,
2329 std_cxx11::array<bool,GeometryInfo<dim>::lines_per_face> >,
2332 for (q = needed_quads.begin(); quad!=triangulation.
end_quad(); ++quad, ++q)
2334 quad->set (q->first);
2335 quad->set_used_flag ();
2336 quad->clear_user_flag ();
2339 quad->set_line_orientation(0,q->second.second[0]);
2340 quad->set_line_orientation(1,q->second.second[1]);
2341 quad->set_line_orientation(2,q->second.second[2]);
2342 quad->set_line_orientation(3,q->second.second[3]);
2347 q->second.first = quad;
2353 triangulation.
levels[0]->cells.reserve_space (cells.size());
2357 std::map<int,std::vector<typename Triangulation<dim,spacedim>::cell_iterator> >
2364 for (
unsigned int c=0; c<cells.size(); ++c, ++cell)
2379 unsigned int face_line_list[4];
2380 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
2382 line_list[line]=std::make_pair(
2385 inverse_line_list[line]=std::pair<int,int> (
2395 typename Triangulation<dim,spacedim>::quad_iterator
2400 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2402 for (
unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
2404 face_to_cell_lines(face,l)]) == needed_lines.end())
2406 face_to_cell_lines(face,l)]]->index();
2412 quad(face_line_list[0],
2417 if (needed_quads.find (quad) != needed_quads.end())
2428 face_iterator[face] = needed_quads[quad].first;
2429 face_orientation[face] =
true;
2430 face_flip[face]=
false;
2431 face_rotation[face]=
false;
2441 test_quad_1(quad.
face(2), quad.
face(3),
2443 test_quad_2(quad.
face(0), quad.
face(1),
2445 test_quad_3(quad.
face(3), quad.
face(2),
2447 test_quad_4(quad.
face(1), quad.
face(0),
2449 test_quad_5(quad.
face(2), quad.
face(3),
2451 test_quad_6(quad.
face(1), quad.
face(0),
2453 test_quad_7(quad.
face(3), quad.
face(2),
2455 if (needed_quads.find (test_quad_1) != needed_quads.end())
2457 face_iterator[face] = needed_quads[test_quad_1].first;
2458 face_orientation[face] =
false;
2459 face_flip[face]=
false;
2460 face_rotation[face]=
false;
2462 else if (needed_quads.find (test_quad_2) != needed_quads.end())
2464 face_iterator[face] = needed_quads[test_quad_2].first;
2465 face_orientation[face] =
false;
2466 face_flip[face]=
false;
2467 face_rotation[face]=
true;
2469 else if (needed_quads.find (test_quad_3) != needed_quads.end())
2471 face_iterator[face] = needed_quads[test_quad_3].first;
2472 face_orientation[face] =
false;
2473 face_flip[face]=
true;
2474 face_rotation[face]=
false;
2476 else if (needed_quads.find (test_quad_4) != needed_quads.end())
2478 face_iterator[face] = needed_quads[test_quad_4].first;
2479 face_orientation[face] =
false;
2480 face_flip[face]=
true;
2481 face_rotation[face]=
true;
2483 else if (needed_quads.find (test_quad_5) != needed_quads.end())
2485 face_iterator[face] = needed_quads[test_quad_5].first;
2486 face_orientation[face] =
true;
2487 face_flip[face]=
false;
2488 face_rotation[face]=
true;
2490 else if (needed_quads.find (test_quad_6) != needed_quads.end())
2492 face_iterator[face] = needed_quads[test_quad_6].first;
2493 face_orientation[face] =
true;
2494 face_flip[face]=
true;
2495 face_rotation[face]=
false;
2497 else if (needed_quads.find (test_quad_7) != needed_quads.end())
2499 face_iterator[face] = needed_quads[test_quad_7].first;
2500 face_orientation[face] =
true;
2501 face_flip[face]=
true;
2502 face_rotation[face]=
true;
2510 Assert(
false,ExcInternalError());
2519 face_iterator[1]->index(),
2520 face_iterator[2]->index(),
2521 face_iterator[3]->index(),
2522 face_iterator[4]->index(),
2523 face_iterator[5]->index()));
2525 cell->set_used_flag ();
2526 cell->set_material_id (cells[c].material_id);
2527 cell->set_manifold_id (cells[c].manifold_id);
2528 cell->clear_user_flag ();
2529 cell->clear_user_data ();
2530 cell->set_subdomain_id (0);
2534 for (
unsigned int quad=0; quad<GeometryInfo<dim>::faces_per_cell; ++quad)
2536 cell->set_face_orientation (quad, face_orientation[quad]);
2537 cell->set_face_flip (quad, face_flip[quad]);
2538 cell->set_face_rotation (quad, face_rotation[quad]);
2545 for (
unsigned int quad=0; quad<GeometryInfo<dim>::faces_per_cell; ++quad)
2546 adjacent_cells[face_iterator[quad]->index()].push_back (cell);
2563 std::multimap<unsigned int, std::pair<unsigned int, unsigned int> >
2565 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2566 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_face; ++line)
2567 cell_to_face_lines.insert(
2568 std::pair<
unsigned int, std::pair<unsigned int, unsigned int> > (
2570 std::pair<unsigned int, unsigned int> (face,line)));
2571 std::multimap<unsigned int, std::pair<unsigned int, unsigned int> >::const_iterator
2572 map_iter=cell_to_face_lines.begin();
2574 for (; map_iter!=cell_to_face_lines.end(); ++map_iter)
2576 const unsigned int cell_line=map_iter->first;
2577 const unsigned int face1=map_iter->second.first;
2578 const unsigned int line1=map_iter->second.second;
2580 Assert(map_iter!=cell_to_face_lines.end(), ExcInternalErrorOnCell(c));
2581 Assert(map_iter->first==cell_line, ExcInternalErrorOnCell(c));
2582 const unsigned int face2=map_iter->second.first;
2583 const unsigned int line2=map_iter->second.second;
2592 face_orientation[face1],
2594 face_rotation[face1])) ==
2597 face_orientation[face2],
2599 face_rotation[face2])),
2600 ExcInternalErrorOnCell(c));
2610 for (
typename Triangulation<dim,spacedim>::quad_iterator
2613 const unsigned int n_adj_cells = adjacent_cells[quad->index()].size();
2618 ExcInternalError());
2624 if (n_adj_cells == 1)
2625 quad->set_boundary_id (0);
2642 for (
typename Triangulation<dim,spacedim>::line_iterator
2667 for (
typename Triangulation<dim,spacedim>::quad_iterator
2669 if (quad->at_boundary())
2670 for (
unsigned int l=0; l<4; ++l)
2671 quad->line(l)->set_boundary_id (0);
2678 std::vector<CellData<1> >::const_iterator boundary_line
2680 std::vector<CellData<1> >::const_iterator end_boundary_line
2682 for (; boundary_line!=end_boundary_line; ++boundary_line)
2684 typename Triangulation<dim,spacedim>::line_iterator line;
2685 std::pair <int, int> line_vertices(std::make_pair(boundary_line->vertices[0],
2686 boundary_line->vertices[1]));
2687 if (needed_lines.find(line_vertices) != needed_lines.end())
2690 line = needed_lines[line_vertices];
2696 std::swap (line_vertices.first, line_vertices.second);
2697 if (needed_lines.find(line_vertices) != needed_lines.end())
2698 line = needed_lines[line_vertices];
2701 AssertThrow (
false, ExcLineInexistant(line_vertices.first,
2702 line_vertices.second));
2708 ExcInteriorLineCantBeBoundary(line->vertex_index(0),
2709 line->vertex_index(1),
2710 boundary_line->boundary_id));
2717 if (line->boundary_id() != 0)
2718 AssertThrow (line->boundary_id() == boundary_line->boundary_id,
2719 ExcMessage (
"Duplicate boundary lines are only allowed " 2720 "if they carry the same boundary indicator."));
2722 line->set_boundary_id (boundary_line->boundary_id);
2724 line->set_manifold_id(boundary_line->manifold_id);
2729 std::vector<CellData<2> >::const_iterator boundary_quad
2731 std::vector<CellData<2> >::const_iterator end_boundary_quad
2733 for (; boundary_quad!=end_boundary_quad; ++boundary_quad)
2735 typename Triangulation<dim,spacedim>::quad_iterator quad;
2736 typename Triangulation<dim,spacedim>::line_iterator line[4];
2745 for (
unsigned int i=0; i<4; ++i)
2747 std::pair<int, int> line_vertices(
2753 if (needed_lines.find(line_vertices) != needed_lines.end())
2754 line[i] = needed_lines[line_vertices];
2759 std::swap (line_vertices.first, line_vertices.second);
2760 if (needed_lines.find(line_vertices) != needed_lines.end())
2761 line[i] = needed_lines[line_vertices];
2765 AssertThrow (
false, ExcLineInexistant(line_vertices.first,
2766 line_vertices.second));
2784 quad_compare_1(line[0]->index(), line[1]->index(),
2785 line[2]->index(), line[3]->index());
2787 quad_compare_2(line[2]->index(), line[3]->index(),
2788 line[0]->index(), line[1]->index());
2801 static const unsigned int lex2cclock[4]= {3,1,0,2};
2807 typename Triangulation<dim,spacedim>::line_iterator
2808 line_counterclock[4];
2809 for (
unsigned int i=0; i<4; ++i)
2810 line_counterclock[lex2cclock[i]]=line[i];
2811 unsigned int n_rotations=0;
2812 bool not_found_quad_1;
2813 while ( (not_found_quad_1=(needed_quads.find(quad_compare_1) == needed_quads.end())) &&
2814 ( needed_quads.find(quad_compare_2) == needed_quads.end()) &&
2819 rotate(line_counterclock, line_counterclock+1, line_counterclock+4);
2823 for (
unsigned int i=0; i<4; ++i)
2825 quad_compare_1.
set_face(i, line_counterclock[lex2cclock[i]]->index());
2826 quad_compare_2.
set_face((i+2)%4, line_counterclock[lex2cclock[i]]->index());
2833 ExcQuadInexistant(line[0]->index(), line[1]->index(),
2834 line[2]->index(), line[3]->index()));
2836 if (not_found_quad_1)
2837 quad = needed_quads[quad_compare_2].first;
2839 quad = needed_quads[quad_compare_1].first;
2844 ExcInteriorQuadCantBeBoundary(quad->vertex_index(0),
2845 quad->vertex_index(1),
2846 quad->vertex_index(2),
2847 quad->vertex_index(3),
2848 boundary_quad->boundary_id));
2855 if (quad->boundary_id() != 0)
2856 AssertThrow (quad->boundary_id() == boundary_quad->boundary_id,
2857 ExcMessage (
"Duplicate boundary quads are only allowed " 2858 "if they carry the same boundary indicator."));
2860 quad->set_boundary_id (boundary_quad->boundary_id);
2861 quad->set_manifold_id (boundary_quad->manifold_id);
2868 cell=triangulation.
begin(); cell!=triangulation.
end(); ++cell)
2869 for (
unsigned int face=0; face<6; ++face)
2870 if (adjacent_cells[cell->quad(face)->index()][0] == cell)
2874 if (adjacent_cells[cell->quad(face)->index()].size() == 2)
2877 cell->set_neighbor (face,
2878 adjacent_cells[cell->quad(face)->index()][1]);
2884 cell->set_neighbor (face,
2885 adjacent_cells[cell->quad(face)->index()][0]);
2904 template <
int spacedim>
2909 std::vector<unsigned int> &,
2910 std::vector<unsigned int> &)
2912 const unsigned int dim = 1;
2924 Assert (!cell->child(0)->has_children() && !cell->child(1)->has_children(),
2925 ExcInternalError());
2930 if (cell->neighbor(0)->has_children())
2933 neighbor = cell->neighbor(0);
2934 Assert (neighbor->level() == cell->level(), ExcInternalError());
2937 neighbor = neighbor->child(1);
2940 Assert (neighbor->neighbor(1) == cell->child(0),
2941 ExcInternalError());
2942 neighbor->set_neighbor (1, cell);
2948 if (neighbor->has_children())
2949 neighbor = neighbor->child(1);
2958 if (cell->neighbor(1)->has_children())
2961 neighbor = cell->neighbor(1);
2962 Assert (neighbor->level() == cell->level(), ExcInternalError());
2965 neighbor = neighbor->child(0);
2968 Assert (neighbor->neighbor(0) == cell->child(1),
2969 ExcInternalError());
2970 neighbor->set_neighbor (0, cell);
2976 if (neighbor->has_children())
2977 neighbor = neighbor->child(0);
2987 triangulation.
vertices_used[cell->child(0)->vertex_index(1)] =
false;
2993 for (
unsigned int child=0; child<cell->n_children(); ++child)
2996 cell->child(child)->clear_user_flag();
2997 cell->child(child)->clear_used_flag();
3002 cell->clear_children ();
3003 cell->clear_user_flag();
3008 template <
int spacedim>
3013 std::vector<unsigned int> &line_cell_count,
3014 std::vector<unsigned int> &)
3016 const unsigned int dim=2;
3019 Assert(line_cell_count.size()==triangulation.
n_raw_lines(), ExcInternalError());
3023 std::vector<typename Triangulation<dim,spacedim>::line_iterator>
3026 lines_to_delete.reserve(4*2+4);
3031 for (
unsigned int c=0; c<cell->n_children(); ++c)
3034 child=cell->child(c);
3035 for (
unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
3036 --line_cell_count[child->line_index(l)];
3050 triangulation.
vertices_used[cell->child(0)->line(1)->vertex_index(1)] =
false;
3052 lines_to_delete.push_back(cell->child(0)->line(1));
3053 lines_to_delete.push_back(cell->child(0)->line(3));
3054 lines_to_delete.push_back(cell->child(3)->line(0));
3055 lines_to_delete.push_back(cell->child(3)->line(2));
3063 lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
3067 for (
unsigned int child=0; child<cell->n_children(); ++child)
3069 cell->child(child)->clear_user_data();
3070 cell->child(child)->clear_user_flag();
3071 cell->child(child)->clear_used_flag();
3076 cell->clear_children ();
3077 cell->clear_refinement_case();
3078 cell->clear_user_flag();
3084 for (
unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
3086 typename Triangulation<dim,spacedim>::line_iterator
3087 line=cell->line(line_no);
3089 if (line->has_children())
3094 Assert((line_cell_count[line->child_index(0)] == 0 &&
3095 line_cell_count[line->child_index(1)] == 0) ||
3096 (line_cell_count[line->child_index(0)] > 0 &&
3097 line_cell_count[line->child_index(1)] > 0),
3098 ExcInternalError());
3100 if (line_cell_count[line->child_index(0)]==0)
3102 for (
unsigned int c=0; c<2; ++c)
3103 Assert (!line->child(c)->has_children(),
3104 ExcInternalError());
3110 triangulation.
vertices_used[line->child(0)->vertex_index(1)] =
false;
3112 lines_to_delete.push_back(line->child(0));
3113 lines_to_delete.push_back(line->child(1));
3115 line->clear_children();
3127 typename std::vector<typename Triangulation<dim,spacedim>::line_iterator>::iterator
3128 line=lines_to_delete.
begin(),
3129 endline=lines_to_delete.end();
3130 for (; line!=endline; ++line)
3132 (*line)->clear_user_data();
3133 (*line)->clear_user_flag();
3134 (*line)->clear_used_flag();
3140 template <
int spacedim>
3145 std::vector<unsigned int> &line_cell_count,
3146 std::vector<unsigned int> &quad_cell_count)
3148 const unsigned int dim=3;
3150 Assert(line_cell_count.size()==triangulation.
n_raw_lines(), ExcInternalError());
3151 Assert(quad_cell_count.size()==triangulation.
n_raw_quads(), ExcInternalError());
3158 std::vector<typename Triangulation<dim,spacedim>::line_iterator>
3160 std::vector<typename Triangulation<dim,spacedim>::quad_iterator>
3163 lines_to_delete.reserve(12*2+6*4+6);
3164 quads_to_delete.reserve(6*4+12);
3168 for (
unsigned int c=0; c<cell->n_children(); ++c)
3171 child=cell->child(c);
3172 for (
unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
3173 --line_cell_count[child->line_index(l)];
3174 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
3175 --quad_cell_count[child->quad_index(f)];
3189 quads_to_delete.push_back(cell->child(0)->face(1));
3192 quads_to_delete.push_back(cell->child(0)->face(3));
3195 quads_to_delete.push_back(cell->child(0)->face(5));
3198 quads_to_delete.push_back(cell->child(0)->face(1));
3199 quads_to_delete.push_back(cell->child(0)->face(3));
3200 quads_to_delete.push_back(cell->child(3)->face(0));
3201 quads_to_delete.push_back(cell->child(3)->face(2));
3203 lines_to_delete.push_back(cell->child(0)->line(11));
3206 quads_to_delete.push_back(cell->child(0)->face(1));
3207 quads_to_delete.push_back(cell->child(0)->face(5));
3208 quads_to_delete.push_back(cell->child(3)->face(0));
3209 quads_to_delete.push_back(cell->child(3)->face(4));
3211 lines_to_delete.push_back(cell->child(0)->line(5));
3214 quads_to_delete.push_back(cell->child(0)->face(3));
3215 quads_to_delete.push_back(cell->child(0)->face(5));
3216 quads_to_delete.push_back(cell->child(3)->face(2));
3217 quads_to_delete.push_back(cell->child(3)->face(4));
3219 lines_to_delete.push_back(cell->child(0)->line(7));
3222 quads_to_delete.push_back(cell->child(0)->face(1));
3223 quads_to_delete.push_back(cell->child(2)->face(1));
3224 quads_to_delete.push_back(cell->child(4)->face(1));
3225 quads_to_delete.push_back(cell->child(6)->face(1));
3227 quads_to_delete.push_back(cell->child(0)->face(3));
3228 quads_to_delete.push_back(cell->child(1)->face(3));
3229 quads_to_delete.push_back(cell->child(4)->face(3));
3230 quads_to_delete.push_back(cell->child(5)->face(3));
3232 quads_to_delete.push_back(cell->child(0)->face(5));
3233 quads_to_delete.push_back(cell->child(1)->face(5));
3234 quads_to_delete.push_back(cell->child(2)->face(5));
3235 quads_to_delete.push_back(cell->child(3)->face(5));
3237 lines_to_delete.push_back(cell->child(0)->line(5));
3238 lines_to_delete.push_back(cell->child(0)->line(7));
3239 lines_to_delete.push_back(cell->child(0)->line(11));
3240 lines_to_delete.push_back(cell->child(7)->line(0));
3241 lines_to_delete.push_back(cell->child(7)->line(2));
3242 lines_to_delete.push_back(cell->child(7)->line(8));
3248 triangulation.
vertices_used[cell->child(0)->vertex_index(7)] =
false;
3253 Assert(
false, ExcInternalError());
3259 for (
unsigned int child=0; child<cell->n_children(); ++child)
3261 cell->child(child)->clear_user_data();
3262 cell->child(child)->clear_user_flag();
3264 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
3269 cell->child(child)->set_face_orientation (f,
true);
3270 cell->child(child)->set_face_flip(f,
false);
3271 cell->child(child)->set_face_rotation(f,
false);
3274 cell->child(child)->clear_used_flag();
3279 cell->clear_children ();
3280 cell->clear_refinement_case ();
3281 cell->clear_user_flag();
3291 for (
unsigned int quad_no=0; quad_no<GeometryInfo<dim>::faces_per_cell; ++quad_no)
3293 typename Triangulation<dim,spacedim>::quad_iterator
3294 quad=cell->face(quad_no);
3298 ExcInternalError());
3300 switch (quad->refinement_case())
3312 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3313 quad_cell_count[quad->child_index(1)] == 0) ||
3314 (quad_cell_count[quad->child_index(0)] > 0 &&
3315 quad_cell_count[quad->child_index(1)] > 0),
3316 ExcInternalError());
3321 unsigned int deleted_grandchildren=0;
3322 unsigned int number_of_child_refinements=0;
3324 for (
unsigned int c=0; c<2; ++c)
3325 if (quad->child(c)->has_children())
3327 ++number_of_child_refinements;
3331 Assert((quad_cell_count[quad->child(c)->child_index(0)] == 0 &&
3332 quad_cell_count[quad->child(c)->child_index(1)] == 0) ||
3333 (quad_cell_count[quad->child(c)->child_index(0)] > 0 &&
3334 quad_cell_count[quad->child(c)->child_index(1)] > 0),
3335 ExcInternalError());
3336 if (quad_cell_count[quad->child(c)->child_index(0)]==0)
3343 ExcInternalError());
3349 quads_to_delete.push_back(quad->child(c)->child(0));
3350 quads_to_delete.push_back(quad->child(c)->child(1));
3352 lines_to_delete.push_back(quad->child(c)->child(0)->line(1));
3354 lines_to_delete.push_back(quad->child(c)->child(0)->line(3));
3355 quad->child(c)->clear_children();
3356 quad->child(c)->clear_refinement_case();
3357 ++deleted_grandchildren;
3365 if (number_of_child_refinements>0 &&
3366 deleted_grandchildren==number_of_child_refinements)
3368 typename Triangulation<dim,spacedim>::line_iterator
3371 middle_line=quad->child(0)->line(1);
3373 middle_line=quad->child(0)->line(3);
3375 lines_to_delete.push_back(middle_line->child(0));
3376 lines_to_delete.push_back(middle_line->child(1));
3377 triangulation.
vertices_used[middle_vertex_index<dim,spacedim>(middle_line)]
3379 middle_line->clear_children();
3384 if (quad_cell_count[quad->child_index(0)]==0)
3390 quads_to_delete.push_back(quad->child(0));
3391 quads_to_delete.push_back(quad->child(1));
3393 lines_to_delete.push_back(quad->child(0)->line(1));
3395 lines_to_delete.push_back(quad->child(0)->line(3));
3415 if (quad->child(0)->has_children())
3448 const typename Triangulation<dim,spacedim>::quad_iterator
3449 switch_1=quad->child(0)->child(1),
3450 switch_2=quad->child(1)->child(0);
3452 Assert(!switch_1->has_children(), ExcInternalError());
3453 Assert(!switch_2->has_children(), ExcInternalError());
3455 const int switch_1_index=switch_1->index();
3456 const int switch_2_index=switch_2->index();
3457 for (
unsigned int l=0; l<triangulation.
levels.size(); ++l)
3458 for (
unsigned int h=0; h<triangulation.
levels[l]->cells.cells.size(); ++h)
3459 for (
unsigned int q=0; q<GeometryInfo<dim>::faces_per_cell; ++q)
3461 const int index=triangulation.
levels[l]->cells.cells[h].face(q);
3462 if (index==switch_1_index)
3463 triangulation.
levels[l]->cells.cells[h].set_face(q,switch_2_index);
3464 else if (index==switch_2_index)
3465 triangulation.
levels[l]->cells.cells[h].set_face(q,switch_1_index);
3470 const int switch_1_lines[4]=
3472 static_cast<signed int>(switch_1->line_index(0)),
3473 static_cast<signed int>(switch_1->line_index(1)),
3474 static_cast<signed int>(switch_1->line_index(2)),
3475 static_cast<signed int>(switch_1->line_index(3))
3477 const bool switch_1_line_orientations[4]=
3479 switch_1->line_orientation(0),
3480 switch_1->line_orientation(1),
3481 switch_1->line_orientation(2),
3482 switch_1->line_orientation(3)
3485 const unsigned int switch_1_user_index=switch_1->user_index();
3486 const bool switch_1_user_flag=switch_1->user_flag_set();
3489 switch_2->line_index(1),
3490 switch_2->line_index(2),
3491 switch_2->line_index(3)));
3492 switch_1->set_line_orientation(0, switch_2->line_orientation(0));
3493 switch_1->set_line_orientation(1, switch_2->line_orientation(1));
3494 switch_1->set_line_orientation(2, switch_2->line_orientation(2));
3495 switch_1->set_line_orientation(3, switch_2->line_orientation(3));
3496 switch_1->set_boundary_id(switch_2->boundary_id());
3497 switch_1->set_manifold_id(switch_2->manifold_id());
3498 switch_1->set_user_index(switch_2->user_index());
3499 if (switch_2->user_flag_set())
3500 switch_1->set_user_flag();
3502 switch_1->clear_user_flag();
3507 switch_1_lines[3]));
3508 switch_2->set_line_orientation(0, switch_1_line_orientations[0]);
3509 switch_2->set_line_orientation(1, switch_1_line_orientations[1]);
3510 switch_2->set_line_orientation(2, switch_1_line_orientations[2]);
3511 switch_2->set_line_orientation(3, switch_1_line_orientations[3]);
3512 switch_2->set_boundary_id(switch_1_boundary_id);
3513 switch_2->set_manifold_id(switch_1->manifold_id());
3514 switch_2->set_user_index(switch_1_user_index);
3515 if (switch_1_user_flag)
3516 switch_2->set_user_flag();
3518 switch_2->clear_user_flag();
3520 const unsigned int child_0=quad->child(0)->child_index(0);
3521 const unsigned int child_2=quad->child(1)->child_index(0);
3522 quad->clear_children();
3523 quad->clear_refinement_case();
3525 quad->set_children(0,child_0);
3526 quad->set_children(2,child_2);
3527 std::swap(quad_cell_count[child_0+1],quad_cell_count[child_2]);
3542 const unsigned int child_0=quad->child(0)->child_index(0);
3543 const unsigned int child_2=quad->child(1)->child_index(0);
3544 quad->clear_children();
3545 quad->clear_refinement_case();
3547 quad->set_children(0,child_0);
3548 quad->set_children(2,child_2);
3552 quad->clear_children();
3553 quad->clear_refinement_case();
3566 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3567 quad_cell_count[quad->child_index(1)] == 0 &&
3568 quad_cell_count[quad->child_index(2)] == 0 &&
3569 quad_cell_count[quad->child_index(3)] == 0) ||
3570 (quad_cell_count[quad->child_index(0)] > 0 &&
3571 quad_cell_count[quad->child_index(1)] > 0 &&
3572 quad_cell_count[quad->child_index(2)] > 0 &&
3573 quad_cell_count[quad->child_index(3)] > 0),
3574 ExcInternalError());
3576 if (quad_cell_count[quad->child_index(0)]==0)
3582 lines_to_delete.push_back(quad->child(0)->line(1));
3583 lines_to_delete.push_back(quad->child(3)->line(0));
3584 lines_to_delete.push_back(quad->child(0)->line(3));
3585 lines_to_delete.push_back(quad->child(3)->line(2));
3587 for (
unsigned int child=0; child<quad->n_children(); ++child)
3588 quads_to_delete.push_back(quad->child(child));
3590 triangulation.
vertices_used[quad->child(0)->vertex_index(3)] =
false;
3592 quad->clear_children();
3593 quad->clear_refinement_case();
3599 Assert(
false, ExcInternalError());
3613 for (
unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
3615 typename Triangulation<dim,spacedim>::line_iterator
3616 line=cell->line(line_no);
3620 ExcInternalError());
3622 if (line->has_children())
3627 Assert((line_cell_count[line->child_index(0)] == 0 &&
3628 line_cell_count[line->child_index(1)] == 0) ||
3629 (line_cell_count[line->child_index(0)] > 0 &&
3630 line_cell_count[line->child_index(1)] > 0),
3631 ExcInternalError());
3633 if (line_cell_count[line->child_index(0)]==0)
3635 for (
unsigned int c=0; c<2; ++c)
3636 Assert (!line->child(c)->has_children(),
3637 ExcInternalError());
3643 triangulation.
vertices_used[line->child(0)->vertex_index(1)] =
false;
3645 lines_to_delete.push_back(line->child(0));
3646 lines_to_delete.push_back(line->child(1));
3648 line->clear_children();
3660 typename std::vector<typename Triangulation<dim,spacedim>::line_iterator>::iterator
3661 line=lines_to_delete.begin(),
3662 endline=lines_to_delete.end();
3663 for (; line!=endline; ++line)
3665 (*line)->clear_user_data();
3666 (*line)->clear_user_flag();
3667 (*line)->clear_used_flag();
3670 typename std::vector<typename Triangulation<dim,spacedim>::quad_iterator>::iterator
3671 quad=quads_to_delete.begin(),
3672 endquad=quads_to_delete.end();
3673 for (; quad!=endquad; ++quad)
3675 (*quad)->clear_user_data();
3676 (*quad)->clear_children();
3677 (*quad)->clear_refinement_case();
3678 (*quad)->clear_user_flag();
3679 (*quad)->clear_used_flag();
3701 template <
int spacedim>
3705 unsigned int &next_unused_vertex,
3706 typename Triangulation<2,spacedim>::raw_line_iterator &next_unused_line,
3710 const unsigned int dim=2;
3713 cell->clear_refine_flag ();
3809 int new_vertices[9];
3810 for (
unsigned int vertex_no=0; vertex_no<4; ++vertex_no)
3811 new_vertices[vertex_no]=cell->vertex_index(vertex_no);
3812 for (
unsigned int line_no=0; line_no<4; ++line_no)
3813 if (cell->line(line_no)->has_children())
3814 new_vertices[4+line_no]=cell->line(line_no)->child(0)->vertex_index(1);
3824 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
3825 ++next_unused_vertex;
3827 ExcMessage(
"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."));
3830 new_vertices[8] = next_unused_vertex;
3842 if (dim == spacedim)
3845 triangulation.
vertices[next_unused_vertex] = cell->center(
true);
3855 if (cell->user_flag_set())
3858 cell->clear_user_flag();
3869 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
3870 if (cell->face(face)->at_boundary())
3891 std::vector<Point<spacedim> > ps(2);
3892 std::vector<double> ws(2, 0.5);
3893 ps[0] = cell->face(boundary_face)
3894 ->child(0)->vertex(1);
3896 ::opposite_face[boundary_face])
3897 ->child(0)->vertex(1);
3899 triangulation.
vertices[next_unused_vertex]
3900 = cell->get_manifold().get_new_point(qs);
3911 cell->clear_user_flag();
3916 Assert(cell->material_id()<= std::numeric_limits<types::material_id>::max(),
3917 ExcIndexRange(cell->material_id(),0,std::numeric_limits<types::material_id>::max()));
3921 triangulation.
vertices[next_unused_vertex] =
3928 typename Triangulation<dim,spacedim>::raw_line_iterator new_lines[12];
3929 unsigned int lmin=8;
3930 unsigned int lmax=12;
3937 for (
unsigned int l=lmin; l<lmax; ++l)
3939 while (next_unused_line->used() ==
true)
3941 new_lines[l] = next_unused_line;
3944 Assert (new_lines[l]->used() ==
false,
3945 ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3959 for (
unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
3960 for (
unsigned int c=0; c<2; ++c, ++l)
3961 new_lines[l]=cell->line(face_no)->child(c);
3962 Assert(l==8, ExcInternalError());
3980 new_lines[0]=cell->line(0);
3981 new_lines[1]=cell->line(1);
3982 new_lines[2]=cell->line(2)->child(0);
3983 new_lines[3]=cell->line(2)->child(1);
3984 new_lines[4]=cell->line(3)->child(0);
3985 new_lines[5]=cell->line(3)->child(1);
3997 new_lines[0]=cell->line(0)->child(0);
3998 new_lines[1]=cell->line(0)->child(1);
3999 new_lines[2]=cell->line(1)->child(0);
4000 new_lines[3]=cell->line(1)->child(1);
4001 new_lines[4]=cell->line(2);
4002 new_lines[5]=cell->line(3);
4007 for (
unsigned int l=lmin; l<lmax; ++l)
4009 new_lines[l]->set_used_flag();
4010 new_lines[l]->clear_user_flag();
4012 new_lines[l]->clear_children();
4015 new_lines[l]->set_manifold_id(cell->manifold_id());
4022 while (next_unused_cell->used() ==
true)
4025 const unsigned int n_children=
4027 for (
unsigned int i=0; i<n_children; ++i)
4029 Assert (next_unused_cell->used() ==
false,
4030 ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4031 subcells[i] = next_unused_cell;
4033 if (i%2==1 && i<n_children-1)
4034 while (next_unused_cell->used() ==
true)
4054 new_lines[8]->index(),
4055 new_lines[4]->index(),
4056 new_lines[10]->index()));
4059 new_lines[2]->index(),
4060 new_lines[5]->index(),
4061 new_lines[11]->index()));
4064 new_lines[9]->index(),
4065 new_lines[10]->index(),
4066 new_lines[6]->index()));
4069 new_lines[3]->index(),
4070 new_lines[11]->index(),
4071 new_lines[7]->index()));
4089 new_lines[6]->index(),
4090 new_lines[2]->index(),
4091 new_lines[4]->index()));
4094 new_lines[1]->index(),
4095 new_lines[3]->index(),
4096 new_lines[5]->index()));
4115 new_lines[2]->index(),
4116 new_lines[4]->index(),
4117 new_lines[6]->index()));
4120 new_lines[3]->index(),
4121 new_lines[6]->index(),
4122 new_lines[5]->index()));
4127 for (
unsigned int i=0; i<n_children; ++i)
4129 subcells[i]->set_used_flag();
4130 subcells[i]->clear_refine_flag();
4131 subcells[i]->clear_user_flag();
4132 subcells[i]->clear_user_data();
4133 subcells[i]->clear_children();
4136 subcells[i]->set_material_id (cell->material_id());
4137 subcells[i]->set_manifold_id (cell->manifold_id());
4138 subcells[i]->set_subdomain_id (subdomainid);
4141 subcells[i]->set_parent (cell->index ());
4149 for (
unsigned int i=0; i<n_children/2; ++i)
4150 cell->set_children (2*i, subcells[2*i]->index());
4152 cell->set_refinement_case(ref_case);
4160 for (
unsigned int c=0; c<n_children; ++c)
4161 cell->child(c)->set_direction_flag (cell->direction_flag());
4171 template <
int spacedim>
4177 const unsigned int dim = 1;
4185 endc = triangulation.
end();
4186 for (; cell != endc; ++cell)
4188 if (cell->refine_flag_set())
4201 unsigned int needed_vertices = 0;
4202 for (
int level=triangulation.
levels.size()-2; level>=0; --level)
4206 unsigned int flagged_cells = 0;
4210 for (; acell!=aendc; ++acell)
4211 if (acell->refine_flag_set())
4216 const unsigned int used_cells
4217 = std::count_if (triangulation.
levels[level+1]->cells.used.begin(),
4218 triangulation.
levels[level+1]->cells.used.end(),
4219 std::bind2nd (std::equal_to<bool>(),
true));
4224 triangulation.
levels[level+1]
4225 ->reserve_space(used_cells+
4232 triangulation.
levels[level+1]->cells
4237 needed_vertices += flagged_cells;
4242 needed_vertices += std::count_if (triangulation.
vertices_used.begin(),
4244 std::bind2nd (std::equal_to<bool>(),
4249 if (needed_vertices > triangulation.
vertices.size())
4251 triangulation.
vertices.resize (needed_vertices,
4253 triangulation.
vertices_used.resize (needed_vertices,
false);
4261 unsigned int next_unused_vertex = 0;
4263 for (
int level=triangulation.
levels.size()-2; level>=0; --level)
4270 next_unused_cell = triangulation.
begin_raw (level+1);
4272 for (; (cell!=endc) && (cell->level()==level); ++cell)
4273 if (cell->refine_flag_set())
4276 cell->clear_refine_flag ();
4280 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
4281 ++next_unused_vertex;
4283 ExcMessage(
"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."));
4288 triangulation.
vertices[next_unused_vertex] =
4298 while (next_unused_cell->used() ==
true)
4300 first_child = next_unused_cell;
4301 first_child->set_used_flag ();
4302 first_child->clear_user_data ();
4304 Assert (next_unused_cell->used() ==
false,
4305 ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4306 second_child = next_unused_cell;
4307 second_child->set_used_flag ();
4308 second_child->clear_user_data ();
4313 cell->set_children (0, first_child->index());
4314 first_child->clear_children ();
4317 next_unused_vertex));
4318 first_child->set_material_id (cell->material_id());
4319 first_child->set_manifold_id (cell->manifold_id());
4320 first_child->set_subdomain_id (subdomainid);
4321 first_child->set_direction_flag (cell->direction_flag());
4323 first_child->set_parent (cell->index ());
4327 first_child->face(1)->set_manifold_id(cell->manifold_id());
4332 first_child->set_neighbor (1, second_child);
4334 first_child->set_neighbor (0, cell->neighbor(0));
4335 else if (cell->neighbor(0)->active())
4341 Assert (cell->neighbor (0)->level () <= cell->level (),
4342 ExcInternalError ());
4343 first_child->set_neighbor (0, cell->neighbor(0));
4349 const unsigned int nbnb = cell->neighbor_of_neighbor (0);
4350 first_child->set_neighbor (0, cell->neighbor(0)->child(nbnb));
4355 left_neighbor = cell->neighbor(0);
4356 while (left_neighbor->has_children())
4358 left_neighbor = left_neighbor->child(nbnb);
4359 left_neighbor->set_neighbor (nbnb, first_child);
4364 second_child->clear_children ();
4367 cell->vertex_index(1)));
4368 second_child->set_neighbor (0, first_child);
4369 second_child->set_material_id (cell->material_id());
4370 second_child->set_manifold_id (cell->manifold_id());
4371 second_child->set_subdomain_id (subdomainid);
4372 second_child->set_direction_flag (cell->direction_flag());
4375 second_child->set_neighbor (1, cell->neighbor(1));
4376 else if (cell->neighbor(1)->active())
4378 Assert (cell->neighbor (1)->level () <= cell->level (),
4379 ExcInternalError ());
4380 second_child->set_neighbor (1, cell->neighbor(1));
4385 const unsigned int nbnb = cell->neighbor_of_neighbor (1);
4386 second_child->set_neighbor (1, cell->neighbor(1)->child(nbnb));
4389 right_neighbor = cell->neighbor(1);
4390 while (right_neighbor->has_children())
4392 right_neighbor = right_neighbor->child(nbnb);
4393 right_neighbor->set_neighbor (nbnb, second_child);
4397 triangulation.
signals.post_refinement_on_cell(cell);
4413 template <
int spacedim>
4417 const bool check_for_distorted_cells)
4419 const unsigned int dim = 2;
4428 endc = triangulation.
end();
4429 for (; cell != endc; ++cell)
4431 if (cell->refine_flag_set())
4441 for (
typename Triangulation<dim,spacedim>::line_iterator
4444 line->clear_user_flag();
4445 line->clear_user_data();
4450 unsigned int n_single_lines=0;
4454 unsigned int n_lines_in_pairs = 0;
4460 unsigned int needed_vertices = 0;
4461 for (
int level=triangulation.
levels.size()-2; level>=0; --level)
4465 unsigned int needed_cells = 0;
4470 for (; cell!=endc; ++cell)
4471 if (cell->refine_flag_set())
4482 n_single_lines += 4;
4494 n_single_lines += 1;
4503 for (
unsigned int line_no=0; line_no<GeometryInfo<dim>::faces_per_cell;
4509 typename Triangulation<dim,spacedim>::line_iterator
4510 line = cell->line(line_no);
4511 if (line->has_children() ==
false)
4513 line->set_user_flag ();
4521 if (line->at_boundary())
4524 line->set_user_index(line->boundary_id());
4529 line->set_user_index(cell->material_id());
4538 const unsigned int used_cells
4539 = std::count_if (triangulation.
levels[level+1]->cells.used.begin(),
4540 triangulation.
levels[level+1]->cells.used.end(),
4541 std::bind2nd (std::equal_to<bool>(),
true));
4547 triangulation.
levels[level+1]
4548 ->reserve_space (used_cells+needed_cells, 2, spacedim);
4552 triangulation.
levels[level+1]->cells.
4553 reserve_space (needed_cells,0);
4557 for (
typename Triangulation<dim,spacedim>::line_iterator
4559 if (line->user_flag_set())
4561 Assert (line->has_children() ==
false, ExcInternalError());
4562 n_lines_in_pairs += 2;
4563 needed_vertices += 1;
4572 triangulation.
faces->lines.
4573 reserve_space (n_lines_in_pairs, 0);
4577 std::bind2nd (std::equal_to<bool>(),
true));
4581 if (needed_vertices > triangulation.
vertices.size())
4584 triangulation.
vertices_used.resize (needed_vertices,
false);
4592 unsigned int next_unused_vertex = 0;
4599 typename Triangulation<dim,spacedim>::active_line_iterator
4602 typename Triangulation<dim,spacedim>::raw_line_iterator
4605 for (; line!=endl; ++line)
4606 if (line->user_flag_set())
4612 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
4613 ++next_unused_vertex;
4615 ExcMessage(
"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."));
4618 if (spacedim == dim)
4625 triangulation.
vertices[next_unused_vertex]
4626 = line->center(
true);
4635 triangulation.
vertices[next_unused_vertex]
4636 = triangulation.
get_manifold(line->user_index()).get_new_point_on_line (line);
4638 triangulation.
vertices[next_unused_vertex]
4639 = line->center(
true);
4644 bool pair_found=
false;
4646 for (; next_unused_line!=endl; ++next_unused_line)
4647 if (!next_unused_line->used() &&
4648 !(++next_unused_line)->used())
4656 Assert (pair_found, ExcInternalError());
4661 line->set_children (0, next_unused_line->index());
4664 const typename Triangulation<dim,spacedim>::raw_line_iterator
4665 children[2] = { next_unused_line,
4670 Assert (children[0]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4671 Assert (children[1]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4675 next_unused_vertex));
4678 line->vertex_index(1)));
4680 children[0]->set_used_flag();
4681 children[1]->set_used_flag();
4682 children[0]->clear_children();
4683 children[1]->clear_children();
4686 children[0]->clear_user_flag();
4687 children[1]->clear_user_flag();
4689 children[0]->set_boundary_id (line->boundary_id());
4690 children[1]->set_boundary_id (line->boundary_id());
4692 children[0]->set_manifold_id (line->manifold_id());
4693 children[1]->set_manifold_id (line->manifold_id());
4697 line->clear_user_flag ();
4706 triangulation.
faces->lines.
4707 reserve_space (0,n_single_lines);
4710 cells_with_distorted_children;
4714 typename Triangulation<dim,spacedim>::raw_line_iterator
4717 for (
int level=0; level<static_cast<int>(triangulation.
levels.size())-1; ++level)
4727 next_unused_cell = triangulation.
begin_raw (level+1);
4729 for (; cell!=endc; ++cell)
4730 if (cell->refine_flag_set())
4738 if (cell->at_boundary())
4739 cell->set_user_flag();
4743 create_children (triangulation,
4749 if ((check_for_distorted_cells ==
true)
4751 has_distorted_children (cell,
4756 triangulation.
signals.post_refinement_on_cell(cell);
4760 return cells_with_distorted_children;
4768 template <
int spacedim>
4772 const bool check_for_distorted_cells)
4774 const unsigned int dim = 3;
4780 Assert (spacedim == 3, ExcNotImplemented());
4789 endc = triangulation.
end();
4790 for (; cell != endc; ++cell)
4792 if (cell->refine_flag_set())
4802 triangulation.
faces->quads.clear_user_data();
4804 for (
typename Triangulation<dim,spacedim>::line_iterator
4806 line->clear_user_flag();
4807 for (
typename Triangulation<dim,spacedim>::quad_iterator
4809 quad->clear_user_flag();
4834 unsigned int needed_vertices = 0;
4835 unsigned int needed_lines_single = 0;
4836 unsigned int needed_quads_single = 0;
4837 unsigned int needed_lines_pair = 0;
4838 unsigned int needed_quads_pair = 0;
4839 for (
int level=triangulation.
levels.size()-2; level>=0; --level)
4843 unsigned int new_cells = 0;
4848 for (; acell!=aendc; ++acell)
4849 if (acell->refine_flag_set())
4859 ++needed_quads_single;
4867 ++needed_lines_single;
4868 needed_quads_single += 4;
4875 needed_lines_single += 6;
4876 needed_quads_single += 12;
4882 Assert(
false, ExcInternalError());
4890 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell;
4894 aface = acell->face(face);
4900 acell->face_orientation(face),
4901 acell->face_flip(face),
4902 acell->face_rotation(face));
4909 if (aface->number_of_children()<4)
4912 aface->set_user_flag();
4914 else if (aface->refinement_case()!=face_ref_case)
4927 ExcInternalError());
4928 aface->set_user_index(face_ref_case);
4934 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
4936 !acell->line(line)->has_children())
4937 acell->line(line)->set_user_flag();
4943 const unsigned int used_cells
4944 = std::count_if (triangulation.
levels[level+1]->cells.used.begin(),
4945 triangulation.
levels[level+1]->cells.used.end(),
4946 std::bind2nd (std::equal_to<bool>(),
true));
4952 triangulation.
levels[level+1]
4953 ->reserve_space (used_cells+new_cells, 3, spacedim);
4956 triangulation.
levels[level+1]->cells.reserve_space (new_cells);
4960 for (
typename Triangulation<dim,spacedim>::quad_iterator
4963 if (quad->user_flag_set())
4969 needed_quads_pair += 4;
4970 needed_lines_pair += 4;
4971 needed_vertices += 1;
4973 if (quad->user_index())
4977 needed_quads_pair += 2;
4978 needed_lines_single += 1;
4991 if (quad->has_children())
4993 Assert(quad->refinement_case()==
RefinementCase<dim-1>::isotropic_refinement, ExcInternalError());
4994 if ((face_refinement_cases[quad->user_index()]==
RefinementCase<dim-1>::cut_x
4995 && (quad->child(0)->line_index(1)+1!=quad->child(2)->line_index(1))) ||
4996 (face_refinement_cases[quad->user_index()]==
RefinementCase<dim-1>::cut_y
4997 && (quad->child(0)->line_index(3)+1!=quad->child(1)->line_index(3))))
4998 needed_lines_pair +=2;
5003 for (
typename Triangulation<dim,spacedim>::line_iterator
5005 if (line->user_flag_set())
5007 needed_lines_pair += 2;
5008 needed_vertices += 1;
5012 triangulation.
faces->lines.
5013 reserve_space (needed_lines_pair,needed_lines_single);
5015 triangulation.
faces->quads.
5016 reserve_space (needed_quads_pair,needed_quads_single);
5021 std::bind2nd (std::equal_to<bool>(),
true));
5025 if (needed_vertices > triangulation.
vertices.size())
5028 triangulation.
vertices_used.resize (needed_vertices,
false);
5045 if (!cell->refine_flag_set())
5046 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
5047 if (cell->line(line)->has_children())
5048 for (
unsigned int c=0; c<2; ++c)
5049 Assert (cell->line(line)->child(c)->user_flag_set() ==
false,
5050 ExcInternalError());
5061 unsigned int next_unused_vertex = 0;
5067 typename Triangulation<dim,spacedim>::active_line_iterator
5070 typename Triangulation<dim,spacedim>::raw_line_iterator
5073 for (; line!=endl; ++line)
5074 if (line->user_flag_set())
5080 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5081 ++next_unused_vertex;
5083 ExcMessage(
"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."));
5086 triangulation.
vertices[next_unused_vertex]
5087 = line->center(
true);
5092 next_unused_line=triangulation.
faces->lines.next_free_pair_object(triangulation);
5094 ExcInternalError());
5099 line->set_children (0, next_unused_line->index());
5102 const typename Triangulation<dim,spacedim>::raw_line_iterator
5103 children[2] = { next_unused_line,
5109 Assert (children[0]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5110 Assert (children[1]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5114 next_unused_vertex));
5117 line->vertex_index(1)));
5119 children[0]->set_used_flag();
5120 children[1]->set_used_flag();
5121 children[0]->clear_children();
5122 children[1]->clear_children();
5125 children[0]->clear_user_flag();
5126 children[1]->clear_user_flag();
5128 children[0]->set_boundary_id (line->boundary_id());
5129 children[1]->set_boundary_id (line->boundary_id());
5131 children[0]->set_manifold_id (line->manifold_id());
5132 children[1]->set_manifold_id (line->manifold_id());
5137 line->clear_user_flag ();
5172 typename Triangulation<dim,spacedim>::quad_iterator
5175 typename Triangulation<dim,spacedim>::raw_line_iterator
5177 typename Triangulation<dim,spacedim>::raw_quad_iterator
5180 for (; quad!=endq; ++quad)
5182 if (quad->user_index())
5184 RefinementCase<dim-1> aniso_quad_ref_case=face_refinement_cases[quad->user_index()];
5193 if (aniso_quad_ref_case == quad->refinement_case())
5198 ExcInternalError());
5203 ExcInternalError());
5206 typename Triangulation<dim,spacedim>::raw_line_iterator new_line;
5208 new_line=triangulation.
faces->lines.next_free_single_object(triangulation);
5209 Assert (new_line->used() ==
false,
5210 ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5225 unsigned int vertex_indices[2];
5228 vertex_indices[0]=quad->line(2)->child(0)->vertex_index(1);
5229 vertex_indices[1]=quad->line(3)->child(0)->vertex_index(1);
5233 vertex_indices[0]=quad->line(0)->child(0)->vertex_index(1);
5234 vertex_indices[1]=quad->line(1)->child(0)->vertex_index(1);
5239 new_line->set_used_flag();
5240 new_line->clear_user_flag();
5242 new_line->clear_children();
5243 new_line->set_boundary_id(quad->boundary_id());
5244 new_line->set_manifold_id(quad->manifold_id());
5252 const unsigned int index[2][2]=
5260 typename Triangulation<dim,spacedim>::raw_quad_iterator new_quads[2];
5262 next_unused_quad=triangulation.
faces->quads.next_free_pair_object(triangulation);
5263 new_quads[0] = next_unused_quad;
5264 Assert (new_quads[0]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5267 new_quads[1] = next_unused_quad;
5268 Assert (new_quads[1]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5276 quad->line(2)->child(index[0][quad->line_orientation(2)])->index(),
5277 quad->line(3)->child(index[0][quad->line_orientation(3)])->index()));
5280 quad->line_index(1),
5281 quad->line(2)->child(index[1][quad->line_orientation(2)])->index(),
5282 quad->line(3)->child(index[1][quad->line_orientation(3)])->index()));
5287 ::
TriaObject<2>(quad->line(0)->child(index[0][quad->line_orientation(0)])->index(),
5288 quad->line(1)->child(index[0][quad->line_orientation(1)])->index(),
5289 quad->line_index(2),
5290 new_line->index()));
5292 ::
TriaObject<2>(quad->line(0)->child(index[1][quad->line_orientation(0)])->index(),
5293 quad->line(1)->child(index[1][quad->line_orientation(1)])->index(),
5295 quad->line_index(3)));
5298 for (
unsigned int i=0; i<2; ++i)
5300 new_quads[i]->set_used_flag();
5301 new_quads[i]->clear_user_flag();
5303 new_quads[i]->clear_children();
5304 new_quads[i]->set_boundary_id (quad->boundary_id());
5305 new_quads[i]->set_manifold_id (quad->manifold_id());
5309 for (
unsigned int j=0; j<GeometryInfo<dim>::lines_per_face; ++j)
5310 new_quads[i]->set_line_orientation(j,
true);
5316 new_quads[0]->set_line_orientation(0,quad->line_orientation(0));
5317 new_quads[0]->set_line_orientation(2,quad->line_orientation(2));
5318 new_quads[1]->set_line_orientation(1,quad->line_orientation(1));
5319 new_quads[1]->set_line_orientation(3,quad->line_orientation(3));
5322 new_quads[0]->set_line_orientation(3,quad->line_orientation(3));
5323 new_quads[1]->set_line_orientation(2,quad->line_orientation(2));
5327 new_quads[0]->set_line_orientation(1,quad->line_orientation(1));
5328 new_quads[1]->set_line_orientation(0,quad->line_orientation(0));
5354 typename Triangulation<dim,spacedim>::line_iterator old_child[2];
5357 old_child[0]=quad->child(0)->line(1);
5358 old_child[1]=quad->child(2)->line(1);
5364 old_child[0]=quad->child(0)->line(3);
5365 old_child[1]=quad->child(1)->line(3);
5368 if (old_child[0]->index()+1 != old_child[1]->index())
5373 typename Triangulation<dim,spacedim>::raw_line_iterator new_child[2];
5375 new_child[0]=new_child[1]=triangulation.
faces->lines.next_free_pair_object(triangulation);
5378 new_child[0]->set_used_flag();
5379 new_child[1]->set_used_flag();
5381 const int old_index_0=old_child[0]->index(),
5382 old_index_1=old_child[1]->index(),
5383 new_index_0=new_child[0]->index(),
5384 new_index_1=new_child[1]->index();
5388 for (
unsigned int q=0; q<triangulation.
faces->quads.cells.size(); ++q)
5389 for (
unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
5391 const int this_index=triangulation.
faces->quads.cells[q].face(l);
5392 if (this_index==old_index_0)
5393 triangulation.
faces->quads.cells[q].set_face(l,new_index_0);
5394 else if (this_index==old_index_1)
5395 triangulation.
faces->quads.cells[q].set_face(l,new_index_1);
5399 for (
unsigned int i=0; i<2; ++i)
5401 Assert(!old_child[i]->has_children(), ExcInternalError());
5404 old_child[i]->vertex_index(1)));
5405 new_child[i]->set_boundary_id(old_child[i]->boundary_id());
5406 new_child[i]->set_manifold_id(old_child[i]->manifold_id());
5407 new_child[i]->set_user_index(old_child[i]->user_index());
5408 if (old_child[i]->user_flag_set())
5409 new_child[i]->set_user_flag();
5411 new_child[i]->clear_user_flag();
5413 new_child[i]->clear_children();
5415 old_child[i]->clear_user_flag();
5416 old_child[i]->clear_user_index();
5417 old_child[i]->clear_used_flag();
5425 new_line->set_children(0, quad->child(0)->line_index(1));
5426 Assert(new_line->child(1)==quad->child(2)->line(1),
5427 ExcInternalError());
5462 const typename Triangulation<dim,spacedim>::quad_iterator
5463 switch_1=quad->child(1),
5464 switch_2=quad->child(2);
5465 const int switch_1_index=switch_1->index();
5466 const int switch_2_index=switch_2->index();
5467 for (
unsigned int l=0; l<triangulation.
levels.size(); ++l)
5468 for (
unsigned int h=0; h<triangulation.
levels[l]->cells.cells.size(); ++h)
5469 for (
unsigned int q=0; q<GeometryInfo<dim>::faces_per_cell; ++q)
5471 const int face_index=triangulation.
levels[l]->cells.cells[h].face(q);
5472 if (face_index==switch_1_index)
5473 triangulation.
levels[l]->cells.cells[h].set_face(q,switch_2_index);
5474 else if (face_index==switch_2_index)
5475 triangulation.
levels[l]->cells.cells[h].set_face(q,switch_1_index);
5479 const unsigned int switch_1_lines[4]=
5481 switch_1->line_index(0),
5482 switch_1->line_index(1),
5483 switch_1->line_index(2),
5484 switch_1->line_index(3)
5486 const bool switch_1_line_orientations[4]=
5488 switch_1->line_orientation(0),
5489 switch_1->line_orientation(1),
5490 switch_1->line_orientation(2),
5491 switch_1->line_orientation(3)
5494 const unsigned int switch_1_user_index=switch_1->user_index();
5495 const bool switch_1_user_flag=switch_1->user_flag_set();
5496 const RefinementCase<dim-1> switch_1_refinement_case=switch_1->refinement_case();
5497 const int switch_1_first_child_pair=(switch_1_refinement_case ? switch_1->child_index(0) : -1);
5498 const int switch_1_second_child_pair=(switch_1_refinement_case==
RefinementCase<dim-1>::cut_xy ? switch_1->child_index(2) : -1);
5501 switch_2->line_index(1),
5502 switch_2->line_index(2),
5503 switch_2->line_index(3)));
5504 switch_1->set_line_orientation(0, switch_2->line_orientation(0));
5505 switch_1->set_line_orientation(1, switch_2->line_orientation(1));
5506 switch_1->set_line_orientation(2, switch_2->line_orientation(2));
5507 switch_1->set_line_orientation(3, switch_2->line_orientation(3));
5508 switch_1->set_boundary_id(switch_2->boundary_id());
5509 switch_1->set_manifold_id(switch_2->manifold_id());
5510 switch_1->set_user_index(switch_2->user_index());
5511 if (switch_2->user_flag_set())
5512 switch_1->set_user_flag();
5514 switch_1->clear_user_flag();
5515 switch_1->clear_refinement_case();
5516 switch_1->set_refinement_case(switch_2->refinement_case());
5517 switch_1->clear_children();
5518 if (switch_2->refinement_case())
5519 switch_1->set_children(0, switch_2->child_index(0));
5521 switch_1->set_children(2, switch_2->child_index(2));
5526 switch_1_lines[3]));
5527 switch_2->set_line_orientation(0, switch_1_line_orientations[0]);
5528 switch_2->set_line_orientation(1, switch_1_line_orientations[1]);
5529 switch_2->set_line_orientation(2, switch_1_line_orientations[2]);
5530 switch_2->set_line_orientation(3, switch_1_line_orientations[3]);
5531 switch_2->set_boundary_id(switch_1_boundary_id);
5532 switch_2->set_manifold_id(switch_1->manifold_id());
5533 switch_2->set_user_index(switch_1_user_index);
5534 if (switch_1_user_flag)
5535 switch_2->set_user_flag();
5537 switch_2->clear_user_flag();
5538 switch_2->clear_refinement_case();
5539 switch_2->set_refinement_case(switch_1_refinement_case);
5540 switch_2->clear_children();
5541 switch_2->set_children(0, switch_1_first_child_pair);
5542 switch_2->set_children(2, switch_1_second_child_pair);
5545 new_quads[0]->set_children(0, quad->child_index(0));
5547 new_quads[1]->set_children(0, quad->child_index(2));
5552 new_quads[0]->set_children(0, quad->child_index(0));
5554 new_quads[1]->set_children(0, quad->child_index(2));
5555 new_line->set_children(0, quad->child(0)->line_index(3));
5556 Assert(new_line->child(1)==quad->child(1)->line(3),
5557 ExcInternalError());
5559 quad->clear_children();
5563 quad->set_children (0, new_quads[0]->index());
5565 quad->set_refinement_case(aniso_quad_ref_case);
5572 if (quad->user_flag_set())
5582 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
5583 ++next_unused_vertex;
5585 ExcMessage(
"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."));
5592 const RefinementCase<dim-1> quad_ref_case=quad->refinement_case();
5601 typename Triangulation<dim,spacedim>::line_iterator middle_line;
5603 middle_line=quad->child(0)->line(1);
5605 middle_line=quad->child(0)->line(3);
5612 if (!middle_line->has_children())
5619 triangulation.
vertices[next_unused_vertex]
5620 = middle_line->center(
true);
5625 next_unused_line=triangulation.
faces->lines.next_free_pair_object(triangulation);
5629 middle_line->set_children (0, next_unused_line->index());
5632 const typename Triangulation<dim,spacedim>::raw_line_iterator
5633 children[2] = { next_unused_line,
5640 Assert (children[0]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5641 Assert (children[1]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5645 next_unused_vertex));
5648 middle_line->vertex_index(1)));
5650 children[0]->set_used_flag();
5651 children[1]->set_used_flag();
5652 children[0]->clear_children();
5653 children[1]->clear_children();
5656 children[0]->clear_user_flag();
5657 children[1]->clear_user_flag();
5659 children[0]->set_boundary_id (middle_line->boundary_id());
5660 children[1]->set_boundary_id (middle_line->boundary_id());
5662 children[0]->set_manifold_id (middle_line->manifold_id());
5663 children[1]->set_manifold_id (middle_line->manifold_id());
5669 quad->clear_user_flag();
5680 if (quad->at_boundary() ||
5682 triangulation.
vertices[next_unused_vertex]
5683 = quad->center(
true);
5707 triangulation.
vertices[next_unused_vertex] =
5708 quad->center(
true,
true);
5714 typename Triangulation<dim,spacedim>::raw_line_iterator new_lines[4];
5716 for (
unsigned int i=0; i<4; ++i)
5725 next_unused_line=triangulation.
faces->lines.next_free_pair_object(triangulation);
5727 new_lines[i] = next_unused_line;
5730 Assert (new_lines[i]->used() ==
false,
5731 ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5751 const unsigned int vertex_indices[5]
5752 = { quad->line(0)->child(0)->vertex_index(1),
5753 quad->line(1)->child(0)->vertex_index(1),
5754 quad->line(2)->child(0)->vertex_index(1),
5755 quad->line(3)->child(0)->vertex_index(1),
5768 for (
unsigned int i=0; i<4; ++i)
5770 new_lines[i]->set_used_flag();
5771 new_lines[i]->clear_user_flag();
5773 new_lines[i]->clear_children();
5774 new_lines[i]->set_boundary_id(quad->boundary_id());
5775 new_lines[i]->set_manifold_id(quad->manifold_id());
5794 const unsigned int index[2][2]=
5800 const int line_indices[12]
5801 = { quad->line(0)->child(index[0][quad->line_orientation(0)])->index(),
5802 quad->line(0)->child(index[1][quad->line_orientation(0)])->index(),
5803 quad->line(1)->child(index[0][quad->line_orientation(1)])->index(),
5804 quad->line(1)->child(index[1][quad->line_orientation(1)])->index(),
5805 quad->line(2)->child(index[0][quad->line_orientation(2)])->index(),
5806 quad->line(2)->child(index[1][quad->line_orientation(2)])->index(),
5807 quad->line(3)->child(index[0][quad->line_orientation(3)])->index(),
5808 quad->line(3)->child(index[1][quad->line_orientation(3)])->index(),
5809 new_lines[0]->index(),
5810 new_lines[1]->index(),
5811 new_lines[2]->index(),
5812 new_lines[3]->index()
5818 typename Triangulation<dim,spacedim>::raw_quad_iterator new_quads[4];
5820 next_unused_quad=triangulation.
faces->quads.next_free_pair_object(triangulation);
5822 new_quads[0] = next_unused_quad;
5823 Assert (new_quads[0]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5826 new_quads[1] = next_unused_quad;
5827 Assert (new_quads[1]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5829 next_unused_quad=triangulation.
faces->quads.next_free_pair_object(triangulation);
5830 new_quads[2] = next_unused_quad;
5831 Assert (new_quads[2]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5834 new_quads[3] = next_unused_quad;
5835 Assert (new_quads[3]->used() ==
false,
ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5838 quad->set_children (0, new_quads[0]->index());
5839 quad->set_children (2, new_quads[2]->index());
5868 for (
unsigned int i=0; i<4; ++i)
5870 new_quads[i]->set_used_flag();
5871 new_quads[i]->clear_user_flag();
5873 new_quads[i]->clear_children();
5874 new_quads[i]->set_boundary_id (quad->boundary_id());
5875 new_quads[i]->set_manifold_id (quad->manifold_id());
5879 for (
unsigned int j=0; j<GeometryInfo<dim>::lines_per_face; ++j)
5880 new_quads[i]->set_line_orientation(j,
true);
5886 new_quads[0]->set_line_orientation(0,quad->line_orientation(0));
5887 new_quads[0]->set_line_orientation(2,quad->line_orientation(2));
5888 new_quads[1]->set_line_orientation(1,quad->line_orientation(1));
5889 new_quads[1]->set_line_orientation(2,quad->line_orientation(2));
5890 new_quads[2]->set_line_orientation(0,quad->line_orientation(0));
5891 new_quads[2]->set_line_orientation(3,quad->line_orientation(3));
5892 new_quads[3]->set_line_orientation(1,quad->line_orientation(1));
5893 new_quads[3]->set_line_orientation(3,quad->line_orientation(3));
5897 quad->clear_user_flag ();
5908 cells_with_distorted_children;
5910 for (
unsigned int level=0; level!=triangulation.
levels.size()-1; ++level)
5915 typename Triangulation<dim,spacedim>::active_hex_iterator
5918 typename Triangulation<dim,spacedim>::raw_hex_iterator
5921 for (; hex!=endh; ++hex)
5922 if (hex->refine_flag_set())
5930 hex->clear_refine_flag ();
5931 hex->set_refinement_case(ref_case);
5942 unsigned int n_new_lines=0;
5943 unsigned int n_new_quads=0;
5944 unsigned int n_new_hexes=0;
5967 Assert(
false, ExcInternalError());
5973 std::vector<typename Triangulation<dim,spacedim>::raw_line_iterator>
5974 new_lines(n_new_lines);
5975 for (
unsigned int i=0; i<n_new_lines; ++i)
5977 new_lines[i] = triangulation.
faces->lines.next_free_single_object(triangulation);
5979 Assert (new_lines[i]->used() ==
false,
5980 ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5981 new_lines[i]->set_used_flag();
5982 new_lines[i]->clear_user_flag();
5983 new_lines[i]->clear_user_data();
5984 new_lines[i]->clear_children();
5988 new_lines[i]->set_manifold_id(hex->manifold_id());
5993 std::vector<typename Triangulation<dim,spacedim>::raw_quad_iterator>
5994 new_quads(n_new_quads);
5995 for (
unsigned int i=0; i<n_new_quads; ++i)
5997 new_quads[i] = triangulation.
faces->quads.next_free_single_object(triangulation);
5999 Assert (new_quads[i]->used() ==
false,
6000 ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6001 new_quads[i]->set_used_flag();
6002 new_quads[i]->clear_user_flag();
6003 new_quads[i]->clear_user_data();
6004 new_quads[i]->clear_children();
6008 new_quads[i]->set_manifold_id (hex->manifold_id());
6011 for (
unsigned int j=0; j<GeometryInfo<dim>::lines_per_face; ++j)
6012 new_quads[i]->set_line_orientation(j,
true);
6019 std::vector<typename Triangulation<dim,spacedim>::raw_hex_iterator>
6020 new_hexes(n_new_hexes);
6021 for (
unsigned int i=0; i<n_new_hexes; ++i)
6024 next_unused_hex=triangulation.
levels[level+1]->cells.next_free_hex(triangulation,level+1);
6028 new_hexes[i]=next_unused_hex;
6030 Assert (new_hexes[i]->used() ==
false,
6031 ExcMessage(
"Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6032 new_hexes[i]->set_used_flag();
6033 new_hexes[i]->clear_user_flag();
6034 new_hexes[i]->clear_user_data();
6035 new_hexes[i]->clear_children();
6038 new_hexes[i]->set_material_id (hex->material_id());
6039 new_hexes[i]->set_manifold_id (hex->manifold_id());
6040 new_hexes[i]->set_subdomain_id (subdomainid);
6043 new_hexes[i]->set_parent (hex->index ());
6055 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
6057 new_hexes[i]->set_face_orientation(f,
true);
6058 new_hexes[i]->set_face_flip(f,
false);
6059 new_hexes[i]->set_face_rotation(f,
false);
6063 for (
unsigned int i=0; i<n_new_hexes/2; ++i)
6064 hex->set_children (2*i, new_hexes[2*i]->index());
6072 = { hex->face_orientation (0),
6073 hex->face_orientation (1),
6074 hex->face_orientation (2),
6075 hex->face_orientation (3),
6076 hex->face_orientation (4),
6077 hex->face_orientation (5)
6082 = { hex->face_flip (0),
6092 = { hex->face_rotation (0),
6093 hex->face_rotation (1),
6094 hex->face_rotation (2),
6095 hex->face_rotation (3),
6096 hex->face_rotation (4),
6097 hex->face_rotation (5)
6102 const unsigned int *vertex_indices=0;
6103 const typename Triangulation<dim,spacedim>::raw_line_iterator
6105 const unsigned int *line_indices=0;
6106 const bool *line_orientation=0;
6107 const int *quad_indices=0;
6117 const unsigned int child_at_origin[2][2][2]=
6203 const typename Triangulation<dim,spacedim>::raw_line_iterator
6207 hex->face(2)->child(0)
6209 hex->face(3)->child(0)
6211 hex->face(4)->child(0)
6213 hex->face(5)->child(0)
6219 unsigned int line_indices_x[4];
6221 for (
unsigned int i=0; i<4; ++i)
6222 line_indices_x[i]=lines[i]->index();
6223 line_indices=&line_indices_x[0];
6230 bool line_orientation_x[4];
6236 const unsigned int middle_vertices[2]=
6238 hex->line(2)->child(0)->vertex_index(1),
6239 hex->line(7)->child(0)->vertex_index(1)
6242 for (
unsigned int i=0; i<4; ++i)
6243 if (lines[i]->vertex_index(i%2)==middle_vertices[i%2])
6244 line_orientation_x[i]=
true;
6249 Assert(lines[i]->vertex_index((i+1)%2)==middle_vertices[i%2],
6250 ExcInternalError());
6251 line_orientation_x[i]=
false;
6254 line_orientation=&line_orientation_x[0];
6264 new_quads[0]->set_line_orientation(0,line_orientation[0]);
6265 new_quads[0]->set_line_orientation(1,line_orientation[1]);
6266 new_quads[0]->set_line_orientation(2,line_orientation[2]);
6267 new_quads[0]->set_line_orientation(3,line_orientation[3]);
6300 const int quad_indices_x[11]
6303 new_quads[0]->index(),
6305 hex->face(0)->index(),
6307 hex->face(1)->index(),
6309 hex->face(2)->child_index( child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]),
6310 hex->face(2)->child_index(1-child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]),
6312 hex->face(3)->child_index( child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]),
6313 hex->face(3)->child_index(1-child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]),
6315 hex->face(4)->child_index( child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]),
6316 hex->face(4)->child_index(1-child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]),
6318 hex->face(5)->child_index( child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]),
6319 hex->face(5)->child_index(1-child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]])
6322 quad_indices=&quad_indices_x[0];
6404 const typename Triangulation<dim,spacedim>::raw_line_iterator
6408 hex->face(0)->child(0)
6410 hex->face(1)->child(0)
6412 hex->face(4)->child(0)
6414 hex->face(5)->child(0)
6420 unsigned int line_indices_y[4];
6422 for (
unsigned int i=0; i<4; ++i)
6423 line_indices_y[i]=lines[i]->index();
6424 line_indices=&line_indices_y[0];
6431 bool line_orientation_y[4];
6437 const unsigned int middle_vertices[2]=
6439 hex->line(0)->child(0)->vertex_index(1),
6440 hex->line(5)->child(0)->vertex_index(1)
6443 for (
unsigned int i=0; i<4; ++i)
6444 if (lines[i]->vertex_index(i%2)==middle_vertices[i%2])
6445 line_orientation_y[i]=
true;
6449 Assert(lines[i]->vertex_index((i+1)%2)==middle_vertices[i%2],
6450 ExcInternalError());
6451 line_orientation_y[i]=
false;
6454 line_orientation=&line_orientation_y[0];
6464 new_quads[0]->set_line_orientation(0,line_orientation[2]);
6465 new_quads[0]->set_line_orientation(1,line_orientation[3]);
6466 new_quads[0]->set_line_orientation(2,line_orientation[0]);
6467 new_quads[0]->set_line_orientation(3,line_orientation[1]);
6500 const int quad_indices_y[11]
6503 new_quads[0]->index(),
6505 hex->face(0)->child_index( child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]),
6506 hex->face(0)->child_index(1-child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]),
6508 hex->face(1)->child_index( child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]),
6509 hex->face(1)->child_index(1-child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]),
6511 hex->face(2)->index(),
6513 hex->face(3)->index(),
6515 hex->face(4)->child_index( child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]),
6516 hex->face(4)->child_index(1-child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]),
6518 hex->face(5)->child_index( child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]),
6519 hex->face(5)->child_index(1-child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]])
6522 quad_indices=&quad_indices_y[0];
6606 const typename Triangulation<dim,spacedim>::raw_line_iterator
6610 hex->face(0)->child(0)
6612 hex->face(1)->child(0)
6614 hex->face(2)->child(0)
6616 hex->face(3)->child(0)
6622 unsigned int line_indices_z[4];
6624 for (
unsigned int i=0; i<4; ++i)
6625 line_indices_z[i]=lines[i]->index();
6626 line_indices=&line_indices_z[0];
6633 bool line_orientation_z[4];
6639 const unsigned int middle_vertices[2]=
6641 middle_vertex_index<dim,spacedim>(hex->line(8)),
6642 middle_vertex_index<dim,spacedim>(hex->line(11))
6645 for (
unsigned int i=0; i<4; ++i)
6646 if (lines[i]->vertex_index(i%2)==middle_vertices[i%2])
6647 line_orientation_z[i]=
true;
6651 Assert(lines[i]->vertex_index((i+1)%2)==middle_vertices[i%2],
6652 ExcInternalError());
6653 line_orientation_z[i]=
false;
6656 line_orientation=&line_orientation_z[0];
6666 new_quads[0]->set_line_orientation(0,line_orientation[0]);
6667 new_quads[0]->set_line_orientation(1,line_orientation[1]);
6668 new_quads[0]->set_line_orientation(2,line_orientation[2]);
6669 new_quads[0]->set_line_orientation(3,line_orientation[3]);
6703 const int quad_indices_z[11]
6706 new_quads[0]->index(),
6708 hex->face(0)->child_index( child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]),
6709 hex->face(0)->child_index(1-child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]),
6711 hex->face(1)->child_index( child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]),
6712 hex->face(1)->child_index(1-child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]),
6714 hex->face(2)->child_index( child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]),
6715 hex->face(2)->child_index(1-child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]),
6717 hex->face(3)->child_index( child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]),
6718 hex->face(3)->child_index(1-child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]),
6720 hex->face(4)->index(),
6722 hex->face(5)->index()
6724 quad_indices=&quad_indices_z[0];
6765 TriaObject<1>(middle_vertex_index<dim,spacedim>(hex->face(4)),
6766 middle_vertex_index<dim,spacedim>(hex->face(5))));
6832 const typename Triangulation<dim,spacedim>::raw_line_iterator
6836 hex->face(0)->child(0)
6838 hex->face(1)->child(0)
6840 hex->face(2)->child(0)
6842 hex->face(3)->child(0)
6868 unsigned int line_indices_xy[13];
6870 for (
unsigned int i=0; i<13; ++i)
6871 line_indices_xy[i]=lines[i]->index();
6872 line_indices=&line_indices_xy[0];
6879 bool line_orientation_xy[13];
6883 const unsigned int middle_vertices[4]=
6885 hex->line(0)->child(0)->vertex_index(1),
6886 hex->line(1)->child(0)->vertex_index(1),
6887 hex->line(2)->child(0)->vertex_index(1),
6888 hex->line(3)->child(0)->vertex_index(1),
6894 for (
unsigned int i=0; i<4; ++i)
6895 if (lines[i]->vertex_index(0)==middle_vertices[i])
6896 line_orientation_xy[i]=
true;
6900 Assert(lines[i]->vertex_index(1)==middle_vertices[i],
6901 ExcInternalError());
6902 line_orientation_xy[i]=
false;
6911 for (
unsigned int i=4; i<12; ++i)
6912 if (lines[i]->vertex_index((i+1)%2) ==
6913 middle_vertex_index<dim,spacedim>(hex->face(3+i/4)))
6914 line_orientation_xy[i]=
true;
6919 Assert(lines[i]->vertex_index(i%2) ==
6920 (middle_vertex_index<dim,spacedim>(hex->face(3+i/4))),
6921 ExcInternalError());
6922 line_orientation_xy[i]=
false;
6928 line_orientation_xy[12]=
true;
6929 line_orientation=&line_orientation_xy[0];
6975 new_quads[0]->set_line_orientation(0,line_orientation[2]);
6976 new_quads[0]->set_line_orientation(2,line_orientation[4]);
6977 new_quads[0]->set_line_orientation(3,line_orientation[8]);
6979 new_quads[1]->set_line_orientation(1,line_orientation[3]);
6980 new_quads[1]->set_line_orientation(2,line_orientation[5]);
6981 new_quads[1]->set_line_orientation(3,line_orientation[9]);
6983 new_quads[2]->set_line_orientation(0,line_orientation[6]);
6984 new_quads[2]->set_line_orientation(1,line_orientation[10]);
6985 new_quads[2]->set_line_orientation(2,line_orientation[0]);
6987 new_quads[3]->set_line_orientation(0,line_orientation[7]);
6988 new_quads[3]->set_line_orientation(1,line_orientation[11]);
6989 new_quads[3]->set_line_orientation(3,line_orientation[1]);
7022 const int quad_indices_xy[20]
7025 new_quads[0]->index(),
7026 new_quads[1]->index(),
7027 new_quads[2]->index(),
7028 new_quads[3]->index(),
7030 hex->face(0)->child_index( child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]),
7031 hex->face(0)->child_index(1-child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]),
7033 hex->face(1)->child_index( child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]),
7034 hex->face(1)->child_index(1-child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]),
7036 hex->face(2)->child_index( child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]),
7037 hex->face(2)->child_index(1-child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]),
7039 hex->face(3)->child_index( child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]),
7040 hex->face(3)->child_index(1-child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]),
7052 quad_indices=&quad_indices_xy[0];
7107 TriaObject<1>(middle_vertex_index<dim,spacedim>(hex->face(2)),
7108 middle_vertex_index<dim,spacedim>(hex->face(3))));
7174 const typename Triangulation<dim,spacedim>::raw_line_iterator
7178 hex->face(0)->child(0)
7180 hex->face(1)->child(0)
7182 hex->face(4)->child(0)
7184 hex->face(5)->child(0)
7210 unsigned int line_indices_xz[13];
7212 for (
unsigned int i=0; i<13; ++i)
7213 line_indices_xz[i]=lines[i]->index();
7214 line_indices=&line_indices_xz[0];
7221 bool line_orientation_xz[13];
7225 const unsigned int middle_vertices[4]=
7227 hex->line(8)->child(0)->vertex_index(1),
7228 hex->line(9)->child(0)->vertex_index(1),
7229 hex->line(2)->child(0)->vertex_index(1),
7230 hex->line(6)->child(0)->vertex_index(1),
7235 for (
unsigned int i=0; i<4; ++i)
7236 if (lines[i]->vertex_index(0)==middle_vertices[i])
7237 line_orientation_xz[i]=
true;
7241 Assert(lines[i]->vertex_index(1)==middle_vertices[i],
7242 ExcInternalError());
7243 line_orientation_xz[i]=
false;
7252 for (
unsigned int i=4; i<12; ++i)
7253 if (lines[i]->vertex_index((i+1)%2) ==
7254 middle_vertex_index<dim,spacedim>(hex->face(1+i/4)))
7255 line_orientation_xz[i]=
true;
7260 Assert(lines[i]->vertex_index(i%2) ==
7261 (middle_vertex_index<dim,spacedim>(hex->face(1+i/4))),
7262 ExcInternalError());
7263 line_orientation_xz[i]=
false;
7269 line_orientation_xz[12]=
true;
7270 line_orientation=&line_orientation_xz[0];
7317 new_quads[0]->set_line_orientation(0,line_orientation[0]);
7318 new_quads[0]->set_line_orientation(2,line_orientation[6]);
7319 new_quads[0]->set_line_orientation(3,line_orientation[10]);
7321 new_quads[1]->set_line_orientation(1,line_orientation[1]);
7322 new_quads[1]->set_line_orientation(2,line_orientation[7]);
7323 new_quads[1]->set_line_orientation(3,line_orientation[11]);
7325 new_quads[2]->set_line_orientation(0,line_orientation[4]);
7326 new_quads[2]->set_line_orientation(1,line_orientation[8]);
7327 new_quads[2]->set_line_orientation(2,line_orientation[2]);
7329 new_quads[3]->set_line_orientation(0,line_orientation[5]);
7330 new_quads[3]->set_line_orientation(1,line_orientation[9]);
7331 new_quads[3]->set_line_orientation(3,line_orientation[3]);
7365 const int quad_indices_xz[20]
7368 new_quads[0]->index(),
7369 new_quads[1]->index(),
7370 new_quads[2]->index(),
7371 new_quads[3]->index(),
7373 hex->face(0)->child_index( child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]),
7374 hex->face(0)->child_index(1-child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]),
7376 hex->face(1)->child_index( child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]),
7377 hex->face(1)->child_index(1-child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]),
7389 hex->face(4)->child_index( child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]),
7390 hex->face(4)->child_index(1-child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]),
7392 hex->face(5)->child_index( child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]),
7393 hex->face(5)->child_index(1-child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]])
7395 quad_indices=&quad_indices_xz[0];
7460 TriaObject<1>(middle_vertex_index<dim,spacedim>(hex->face(0)),
7461 middle_vertex_index<dim,spacedim>(hex->face(1))));
7528 const typename Triangulation<dim,spacedim>::raw_line_iterator
7532 hex->face(2)->child(0)
7534 hex->face(3)->child(0)
7536 hex->face(4)->child(0)
7538 hex->face(5)->child(0)
7564 unsigned int line_indices_yz[13];
7566 for (
unsigned int i=0; i<13; ++i)
7567 line_indices_yz[i]=lines[i]->index();
7568 line_indices=&line_indices_yz[0];
7575 bool line_orientation_yz[13];
7579 const unsigned int middle_vertices[4]=
7581 hex->line(8)->child(0)->vertex_index(1),
7582 hex->line(10)->child(0)->vertex_index(1),
7583 hex->line(0)->child(0)->vertex_index(1),
7584 hex->line(4)->child(0)->vertex_index(1),
7589 for (
unsigned int i=0; i<4; ++i)
7590 if (lines[i]->vertex_index(0)==middle_vertices[i])
7591 line_orientation_yz[i]=
true;
7595 Assert(lines[i]->vertex_index(1)==middle_vertices[i],
7596 ExcInternalError());
7597 line_orientation_yz[i]=
false;
7606 for (
unsigned int i=4; i<12; ++i)
7607 if (lines[i]->vertex_index((i+1)%2) ==
7608 middle_vertex_index<dim,spacedim>(hex->face(i/4-1)))
7609 line_orientation_yz[i]=
true;
7614 Assert(lines[i]->vertex_index(i%2) ==
7615 (middle_vertex_index<dim,spacedim>(hex->face(i/4-1))),
7616 ExcInternalError());
7617 line_orientation_yz[i]=
false;
7623 line_orientation_yz[12]=
true;
7624 line_orientation=&line_orientation_yz[0];
7665 new_quads[0]->set_line_orientation(0,line_orientation[6]);
7666 new_quads[0]->set_line_orientation(1,line_orientation[10]);
7667 new_quads[0]->set_line_orientation(2,line_orientation[0]);
7669 new_quads[1]->set_line_orientation(0,line_orientation[7]);
7670 new_quads[1]->set_line_orientation(1,line_orientation[11]);
7671 new_quads[1]->set_line_orientation(3,line_orientation[1]);
7673 new_quads[2]->set_line_orientation(0,line_orientation[2]);
7674 new_quads[2]->set_line_orientation(2,line_orientation[4]);
7675 new_quads[2]->set_line_orientation(3,line_orientation[8]);
7677 new_quads[3]->set_line_orientation(1,line_orientation[3]);
7678 new_quads[3]->set_line_orientation(2,line_orientation[5]);
7679 new_quads[3]->set_line_orientation(3,line_orientation[9]);
7713 const int quad_indices_yz[20]
7716 new_quads[0]->index(),
7717 new_quads[1]->index(),
7718 new_quads[2]->index(),
7719 new_quads[3]->index(),
7731 hex->face(2)->child_index( child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]),
7732 hex->face(2)->child_index(1-child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]),
7734 hex->face(3)->child_index( child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]),
7735 hex->face(3)->child_index(1-child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]),
7737 hex->face(4)->child_index( child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]),
7738 hex->face(4)->child_index(1-child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]),
7740 hex->face(5)->child_index( child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]),
7741 hex->face(5)->child_index(1-child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]])
7743 quad_indices=&quad_indices_yz[0];
7800 while (triangulation.
vertices_used[next_unused_vertex] ==
true)
7801 ++next_unused_vertex;
7803 ExcMessage(
"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."));
7812 triangulation.
vertices[next_unused_vertex] =
7813 hex->center(
true,
true);
7834 const unsigned int vertex_indices_xyz[7]
7835 = { middle_vertex_index<dim,spacedim>(hex->face(0)),
7836 middle_vertex_index<dim,spacedim>(hex->face(1)),
7837 middle_vertex_index<dim,spacedim>(hex->face(2)),
7838 middle_vertex_index<dim,spacedim>(hex->face(3)),
7839 middle_vertex_index<dim,spacedim>(hex->face(4)),
7840 middle_vertex_index<dim,spacedim>(hex->face(5)),
7843 vertex_indices=&vertex_indices_xyz[0];
7923 const typename Triangulation<dim,spacedim>::raw_line_iterator
7989 lines=&lines_xyz[0];
7991 unsigned int line_indices_xyz[30];
7992 for (
unsigned int i=0; i<30; ++i)
7993 line_indices_xyz[i]=lines[i]->index();
7994 line_indices=&line_indices_xyz[0];
8001 bool line_orientation_xyz[30];
8009 for (
unsigned int i=0; i<24; ++i)
8010 if (lines[i]->vertex_index((i+1)%2)==vertex_indices[i/4])
8011 line_orientation_xyz[i]=
true;
8016 Assert(lines[i]->vertex_index(i%2)==vertex_indices[i/4],
8017 ExcInternalError());
8018 line_orientation_xyz[i]=
false;
8023 for (
unsigned int i=24; i<30; ++i)
8024 line_orientation_xyz[i]=
true;
8025 line_orientation=&line_orientation_xyz[0];
8121 new_quads[0]->set_line_orientation(0,line_orientation[10]);
8122 new_quads[0]->set_line_orientation(2,line_orientation[16]);
8124 new_quads[1]->set_line_orientation(1,line_orientation[14]);
8125 new_quads[1]->set_line_orientation(2,line_orientation[17]);
8127 new_quads[2]->set_line_orientation(0,line_orientation[11]);
8128 new_quads[2]->set_line_orientation(3,line_orientation[20]);
8130 new_quads[3]->set_line_orientation(1,line_orientation[15]);
8131 new_quads[3]->set_line_orientation(3,line_orientation[21]);
8133 new_quads[4]->set_line_orientation(0,line_orientation[18]);
8134 new_quads[4]->set_line_orientation(2,line_orientation[0]);
8136 new_quads[5]->set_line_orientation(1,line_orientation[22]);
8137 new_quads[5]->set_line_orientation(2,line_orientation[1]);
8139 new_quads[6]->set_line_orientation(0,line_orientation[19]);
8140 new_quads[6]->set_line_orientation(3,line_orientation[4]);
8142 new_quads[7]->set_line_orientation(1,line_orientation[23]);
8143 new_quads[7]->set_line_orientation(3,line_orientation[5]);
8145 new_quads[8]->set_line_orientation(0,line_orientation[2]);
8146 new_quads[8]->set_line_orientation(2,line_orientation[8]);
8148 new_quads[9]->set_line_orientation(1,line_orientation[6]);
8149 new_quads[9]->set_line_orientation(2,line_orientation[9]);
8151 new_quads[10]->set_line_orientation(0,line_orientation[3]);
8152 new_quads[10]->set_line_orientation(3,line_orientation[12]);
8154 new_quads[11]->set_line_orientation(1,line_orientation[7]);
8155 new_quads[11]->set_line_orientation(3,line_orientation[13]);
8196 const int quad_indices_xyz[36]
8199 new_quads[0]->index(),
8200 new_quads[1]->index(),
8201 new_quads[2]->index(),
8202 new_quads[3]->index(),
8203 new_quads[4]->index(),
8204 new_quads[5]->index(),
8205 new_quads[6]->index(),
8206 new_quads[7]->index(),
8207 new_quads[8]->index(),
8208 new_quads[9]->index(),
8209 new_quads[10]->index(),
8210 new_quads[11]->index(),
8242 quad_indices=&quad_indices_xyz[0];
8312 Assert(
false, ExcInternalError());
8335 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
8336 for (
unsigned int s=0;
8341 const unsigned int current_child
8353 new_hexes[current_child]->set_face_orientation (f, f_or[f]);
8354 new_hexes[current_child]->set_face_flip (f, f_fl[f]);
8355 new_hexes[current_child]->set_face_rotation (f, f_ro[f]);
8360 if ((check_for_distorted_cells ==
true)
8362 has_distorted_children (hex,
8371 triangulation.
signals.post_refinement_on_cell(hex);
8380 triangulation.
faces->quads.clear_user_data();
8383 return cells_with_distorted_children;
8399 template <
int spacedim>
8405 template <
int dim,
int spacedim>
8412 if (spacedim>dim)
return;
8415 cell=triangulation.
begin(); cell!=triangulation.
end(); ++cell)
8416 if (cell->at_boundary() &&
8417 cell->refine_flag_set() &&
8424 for (
unsigned int face_no=0;
8425 face_no<GeometryInfo<dim>::faces_per_cell;
8427 if (cell->face(face_no)->at_boundary())
8437 face = cell->face(face_no);
8441 = face->center(
true);
8446 transform_real_to_unit_cell(cell,
8457 const double dist = std::fabs(new_unit[face_no/2] - face_no%2);
8462 const double allowed=0.25;
8465 cell->flag_for_face_refinement(face_no);
8484 template <
int dim,
int spacedim>
8490 ExcMessage (
"Wrong function called -- there should " 8491 "be a specialization."));
8495 template <
int spacedim>
8500 const unsigned int dim = 3;
8512 bool mesh_changed =
false;
8516 mesh_changed =
false;
8527 cell=triangulation.
begin(); cell!=triangulation.
end(); ++cell)
8528 if (cell->refine_flag_set())
8530 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
8535 cell->line(line)->set_user_flag();
8537 else if (cell->has_children() && !cell->child(0)->coarsen_flag_set())
8539 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
8544 cell->line(line)->set_user_flag();
8546 else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
8547 cell->set_user_flag();
8555 cell=triangulation.
last_active(); cell!=triangulation.
end(); --cell)
8556 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
8558 if (cell->line(line)->has_children())
8567 bool offending_line_found =
false;
8569 for (
unsigned int c=0; c<2; ++c)
8571 Assert (cell->line(line)->child(c)->has_children() ==
false,
8572 ExcInternalError());
8574 if (cell->line(line)->child(c)->user_flag_set () &&
8580 cell->clear_coarsen_flag ();
8587 cell->flag_for_line_refinement(line);
8589 cell->set_refine_flag();
8591 for (
unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
8595 cell->line(l)->set_user_flag();
8598 offending_line_found =
true;
8604 for (
unsigned int l=0;
8605 l<GeometryInfo<dim>::lines_per_cell; ++l)
8606 if (!cell->line(l)->has_children() &&
8610 cell->line(l)->set_user_flag();
8616 if (offending_line_found)
8618 mesh_changed =
true;
8636 cell=triangulation.
last(); cell!=triangulation.
end(); --cell)
8638 if (cell->user_flag_set())
8639 for (
unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
8640 if (cell->line(line)->has_children() &&
8641 (cell->line(line)->child(0)->user_flag_set() ||
8642 cell->line(line)->child(1)->user_flag_set()))
8644 for (
unsigned int c=0; c<cell->n_children(); ++c)
8645 cell->child(c)->clear_coarsen_flag ();
8646 cell->clear_user_flag();
8647 for (
unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
8652 cell->line(l)->set_user_flag();
8653 mesh_changed =
true;
8658 while (mesh_changed ==
true);
8669 template <
int dim,
int spacedim>
8680 for (
unsigned int n=0; n<GeometryInfo<dim>::faces_per_cell; ++n)
8689 const unsigned int n_subfaces
8692 if (n_subfaces == 0 || cell->at_boundary(n))
8694 for (
unsigned int c=0; c<n_subfaces; ++c)
8698 child_cell_on_face(ref_case,
8702 child_neighbor = child->neighbor(n);
8703 if (!child->neighbor_is_coarser(n))
8716 if ((child_neighbor->has_children() &&
8717 !child_neighbor->user_flag_set())||
8726 child_neighbor->refine_flag_set())
8737 template <
int dim,
int spacedim>
8743 template <
int dim,
int spacedim>
8749 template <
int dim,
int spacedim>
8752 const bool check_for_distorted_cells)
8754 smooth_grid(smooth_grid),
8756 anisotropic_refinement(false),
8757 check_for_distorted_cells(check_for_distorted_cells),
8758 vertex_to_boundary_id_map_1d (0),
8759 vertex_to_manifold_id_map_1d (0)
8764 =
new std::map<unsigned int, types::boundary_id>();
8766 =
new std::map<unsigned int, types::manifold_id>();
8776 template <
int dim,
int spacedim>
8789 "because copying Triangulation objects is not " 8790 "allowed. Use Triangulation::copy_from() instead."));
8795 template <
int dim,
int spacedim>
8798 for (
unsigned int i=0; i<
levels.size(); ++i)
8808 ExcInternalError());
8815 ExcInternalError());
8821 template <
int dim,
int spacedim>
8830 template <
int dim,
int spacedim>
8841 template <
int dim,
int spacedim>
8849 template <
int dim,
int spacedim>
8857 manifold[m_number] = &manifold_object;
8861 template <
int dim,
int spacedim>
8869 template <
int dim,
int spacedim>
8880 template <
int dim,
int spacedim>
8887 for (; cell != endc; ++cell)
8888 cell->set_all_manifold_ids(m_number);
8892 template <
int dim,
int spacedim>
8899 for (; cell != endc; ++cell)
8900 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
8901 if (cell->face(f)->at_boundary())
8902 cell->face(f)->set_all_manifold_ids(m_number);
8906 template <
int dim,
int spacedim>
8911 bool boundary_found =
false;
8915 for (; cell != endc; ++cell)
8918 for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
8919 if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id()==b_id)
8921 boundary_found =
true;
8922 cell->face(f)->set_manifold_id(m_number);
8927 for (
unsigned int e=0; e<GeometryInfo<dim>::lines_per_cell; ++e)
8928 if (cell->line(e)->at_boundary() && cell->line(e)->boundary_id()==b_id)
8930 boundary_found =
true;
8931 cell->line(e)->set_manifold_id(m_number);
8935 (void)boundary_found;
8936 Assert(boundary_found, ExcBoundaryIdNotFound(b_id));
8940 template <
int dim,
int spacedim>
8947 ExcMessage(
"You tried to get a Boundary, but I only have a Manifold."));
8953 template <
int dim,
int spacedim>
8965 return *(it->second);
8978 template <
int dim,
int spacedim>
8979 std::vector<types::boundary_id>
8986 std::vector<types::boundary_id> boundary_ids;
8987 for (std::map<unsigned int, types::boundary_id>::const_iterator
8991 boundary_ids.push_back (p->second);
8993 return boundary_ids;
8997 std::set<types::boundary_id> b_ids;
8999 for (; cell!=
end(); ++cell)
9000 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
9001 if (cell->at_boundary(face))
9002 b_ids.insert(cell->face(face)->boundary_id());
9003 std::vector<types::boundary_id> boundary_ids(b_ids.begin(), b_ids.end());
9004 return boundary_ids;
9010 template <
int dim,
int spacedim>
9011 std::vector<types::boundary_id>
9019 template <
int dim,
int spacedim>
9020 std::vector<types::manifold_id>
9023 std::set<types::manifold_id> m_ids;
9025 for (; cell!=
end(); ++cell)
9027 m_ids.insert(cell->manifold_id());
9029 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
9030 if (cell->at_boundary(face))
9031 m_ids.insert(cell->face(face)->manifold_id());
9033 std::vector<types::manifold_id> manifold_indicators(m_ids.begin(), m_ids.end());
9034 return manifold_indicators;
9040 template <
int dim,
int spacedim>
9051 (dim == 1 || old_tria.
faces != NULL),
9052 ExcMessage(
"When calling Triangulation::copy_triangulation(), " 9053 "the target triangulation must be empty but the source " 9054 "triangulation (the argument to this function) must contain " 9055 "something. Here, it seems like the source does not " 9056 "contain anything at all."));
9069 bdry_iterator = old_tria.
manifold.begin();
9070 for (; bdry_iterator != old_tria.
manifold.end() ; bdry_iterator++)
9071 manifold[bdry_iterator->first] = bdry_iterator->second;
9075 for (
unsigned int level=0; level<old_tria.
levels.size(); ++level)
9078 TriaLevel<dim>(*old_tria.
levels[level]));
9086 = (
new std::map<unsigned int, types::boundary_id>
9091 = (
new std::map<unsigned int, types::manifold_id>
9107 template <
int dim,
int spacedim>
9114 std::vector<CellData<dim> > reordered_cells (cells);
9118 reorder_compatibility (reordered_cells, reordered_subcelldata);
9127 template <
int dim,
int spacedim>
9166 DistortedCellList distorted_cells = collect_distorted_coarse_cells (*
this);
9170 AssertThrow (distorted_cells.distorted_cells.size() == 0,
9211 bool values [][2] = {{
false,
true},
9214 for (
unsigned int i=0; i< GeometryInfo< dim >::faces_per_cell; ++i)
9215 for (
unsigned int j=0; j< GeometryInfo< dim >::faces_per_cell; ++j)
9216 correct(i,j) = ( values[i][j]);
9221 bool values [][4]= {{
false,true ,true ,
false},
9222 {true ,
false,
false,
true },
9223 {true ,
false,
false,
true },
9224 {
false,true ,true ,
false}
9226 for (
unsigned int i=0; i< GeometryInfo< dim >::faces_per_cell; ++i)
9227 for (
unsigned int j=0; j< GeometryInfo< dim >::faces_per_cell; ++j)
9228 correct(i,j) = ( values[i][j]);
9232 Assert (
false, ExcNotImplemented());
9236 std::list<active_cell_iterator> this_round, next_round;
9243 while (this_round.size() > 0)
9245 for (
typename std::list<active_cell_iterator>::iterator cell = this_round.begin();
9246 cell != this_round.end(); ++cell)
9248 for (
unsigned int i = 0; i < GeometryInfo< dim >::faces_per_cell; ++i)
9250 if ( !((*cell)->face(i)->at_boundary()) )
9252 neighbor = (*cell)->neighbor(i);
9254 unsigned int cf = (*cell)->face_index(i);
9256 while (neighbor->face_index(j) != cf)
9261 if ( (correct(i,j) && !(*cell)->direction_flag())
9263 (!correct(i,j) && (*cell)->direction_flag()) )
9265 if (neighbor->user_flag_set() ==
false)
9267 neighbor->set_direction_flag (
false);
9268 neighbor->set_user_flag ();
9269 next_round.push_back (neighbor);
9272 Assert (neighbor->direction_flag() ==
false,
9273 ExcNonOrientableTriangulation());
9284 if (next_round.size() == 0)
9286 cell !=
end(); ++cell)
9287 if (cell->user_flag_set() ==
false)
9289 next_round.push_back (cell);
9290 cell->set_direction_flag (
true);
9291 cell->set_user_flag ();
9295 this_round = next_round;
9307 template <
int dim,
int spacedim>
9314 cell !=
end(); ++cell)
9315 cell->set_direction_flag (!cell->direction_flag());
9320 template <
int dim,
int spacedim>
9327 for (; cell != endc; ++cell)
9329 cell->clear_coarsen_flag();
9330 cell->set_refine_flag ();
9336 template <
int dim,
int spacedim>
9339 for (
unsigned int i=0; i<times; ++i)
9352 template <
int dim,
int spacedim>
9356 std::vector<bool>::iterator i = v.begin();
9359 for (; cell!=endc; ++cell)
9360 for (
unsigned int j=0; j<dim; ++j,++i)
9361 if (cell->refine_flag_set() & (1<<j) )
9364 Assert (i == v.end(), ExcInternalError());
9369 template <
int dim,
int spacedim>
9372 std::vector<bool> v;
9380 template <
int dim,
int spacedim>
9383 std::vector<bool> v;
9391 template <
int dim,
int spacedim>
9398 std::vector<bool>::const_iterator i = v.begin();
9399 for (; cell!=endc; ++cell)
9401 unsigned int ref_case=0;
9403 for (
unsigned int j=0; j<dim; ++j, ++i)
9407 ExcGridReadError());
9411 cell->clear_refine_flag();
9414 Assert (i == v.end(), ExcInternalError());
9419 template <
int dim,
int spacedim>
9423 std::vector<bool>::iterator i = v.begin();
9426 for (; cell!=endc; ++cell, ++i)
9427 *i = cell->coarsen_flag_set();
9429 Assert (i == v.end(), ExcInternalError());
9434 template <
int dim,
int spacedim>
9437 std::vector<bool> v;
9445 template <
int dim,
int spacedim>
9448 std::vector<bool> v;
9449 read_bool_vector (mn_tria_coarsen_flags_begin, v, mn_tria_coarsen_flags_end,
9456 template <
int dim,
int spacedim>
9463 std::vector<bool>::const_iterator i = v.begin();
9464 for (; cell!=endc; ++cell, ++i)
9466 cell->set_coarsen_flag();
9468 cell->clear_coarsen_flag();
9470 Assert (i == v.end(), ExcInternalError());
9474 template <
int dim,
int spacedim>
9491 for (
unsigned int level=0; level<
levels.size(); ++level)
9492 levels[level]->cells.clear_user_data();
9505 faces->
lines.clear_user_data();
9511 faces->
lines.clear_user_data();
9517 template <
int dim,
int spacedim>
9532 for (
unsigned int level=0; level<
levels.size(); ++level)
9533 levels[level]->cells.clear_user_flags();
9540 faces->lines.clear_user_flags();
9545 template <
int dim,
int spacedim>
9564 for (
unsigned int level=0; level<
levels.size(); ++level)
9565 levels[level]->cells.clear_user_flags();
9572 faces->quads.clear_user_flags();
9577 template <
int dim,
int spacedim>
9603 for (
unsigned int level=0; level<
levels.size(); ++level)
9604 levels[level]->cells.clear_user_flags();
9609 template <
int dim,
int spacedim>
9617 template <
int dim,
int spacedim>
9627 template <
int dim,
int spacedim>
9639 Assert (
false, ExcNotImplemented());
9644 template <
int dim,
int spacedim>
9651 std::vector<bool> tmp;
9654 v.insert (v.end(), tmp.begin(), tmp.end());
9659 v.insert (v.end(), tmp.begin(), tmp.end());
9665 v.insert (v.end(), tmp.begin(), tmp.end());
9669 Assert (
false, ExcNotImplemented());
9674 template <
int dim,
int spacedim>
9686 Assert (
false, ExcNotImplemented());
9691 template <
int dim,
int spacedim>
9695 std::vector<bool> tmp;
9699 tmp.insert (tmp.end(),
9700 v.begin(), v.begin()+
n_lines());
9707 tmp.insert (tmp.end(),
9715 tmp.insert (tmp.end(),
9721 Assert (
false, ExcNotImplemented());
9726 template <
int dim,
int spacedim>
9730 std::vector<bool>::iterator i = v.begin();
9733 for (; line!=endl; ++line, ++i)
9734 *i = line->user_flag_set();
9736 Assert (i == v.end(), ExcInternalError());
9741 template <
int dim,
int spacedim>
9744 std::vector<bool> v;
9746 write_bool_vector (mn_tria_line_user_flags_begin, v, mn_tria_line_user_flags_end,
9752 template <
int dim,
int spacedim>
9755 std::vector<bool> v;
9756 read_bool_vector (mn_tria_line_user_flags_begin, v, mn_tria_line_user_flags_end,
9763 template <
int dim,
int spacedim>
9770 std::vector<bool>::const_iterator i = v.begin();
9771 for (; line!=endl; ++line, ++i)
9773 line->set_user_flag();
9775 line->clear_user_flag();
9777 Assert (i == v.end(), ExcInternalError());
9783 template <
typename Iterator>
9784 bool get_user_flag (
const Iterator &i)
9786 return i->user_flag_set();
9791 template <
int structdim,
int dim,
int spacedim>
9794 Assert (
false, ExcInternalError());
9800 template <
typename Iterator>
9801 void set_user_flag (
const Iterator &i)
9808 template <
int structdim,
int dim,
int spacedim>
9811 Assert (
false, ExcInternalError());
9816 template <
typename Iterator>
9817 void clear_user_flag (
const Iterator &i)
9819 i->clear_user_flag();
9824 template <
int structdim,
int dim,
int spacedim>
9827 Assert (
false, ExcInternalError());
9832 template <
int dim,
int spacedim>
9839 std::vector<bool>::iterator i = v.begin();
9842 for (; quad!=endq; ++quad, ++i)
9843 *i = get_user_flag (quad);
9845 Assert (i == v.end(), ExcInternalError());
9851 template <
int dim,
int spacedim>
9854 std::vector<bool> v;
9856 write_bool_vector (mn_tria_quad_user_flags_begin, v, mn_tria_quad_user_flags_end,
9862 template <
int dim,
int spacedim>
9865 std::vector<bool> v;
9866 read_bool_vector (mn_tria_quad_user_flags_begin, v, mn_tria_quad_user_flags_end,
9873 template <
int dim,
int spacedim>
9882 std::vector<bool>::const_iterator i = v.begin();
9883 for (; quad!=endq; ++quad, ++i)
9885 set_user_flag(quad);
9887 clear_user_flag(quad);
9889 Assert (i == v.end(), ExcInternalError());
9895 template <
int dim,
int spacedim>
9898 v.resize (
n_hexs(),
false);
9902 std::vector<bool>::iterator i = v.begin();
9905 for (; hex!=endh; ++hex, ++i)
9906 *i = get_user_flag (hex);
9908 Assert (i == v.end(), ExcInternalError());
9914 template <
int dim,
int spacedim>
9917 std::vector<bool> v;
9925 template <
int dim,
int spacedim>
9928 std::vector<bool> v;
9929 read_bool_vector (mn_tria_hex_user_flags_begin, v, mn_tria_hex_user_flags_end,
9936 template <
int dim,
int spacedim>
9945 std::vector<bool>::const_iterator i = v.begin();
9946 for (; hex!=endh; ++hex, ++i)
9950 clear_user_flag(hex);
9952 Assert (i == v.end(), ExcInternalError());
9958 template <
int dim,
int spacedim>
9965 std::vector<unsigned int> tmp;
9968 v.insert (v.end(), tmp.begin(), tmp.end());
9973 v.insert (v.end(), tmp.begin(), tmp.end());
9979 v.insert (v.end(), tmp.begin(), tmp.end());
9983 Assert (
false, ExcNotImplemented());
9988 template <
int dim,
int spacedim>
9992 std::vector<unsigned int> tmp;
9996 tmp.insert (tmp.end(),
9997 v.begin(), v.begin()+
n_lines());
10004 tmp.insert (tmp.end(),
10012 tmp.insert (tmp.end(),
10018 Assert (
false, ExcNotImplemented());
10025 template <
typename Iterator>
10026 unsigned int get_user_index (
const Iterator &i)
10028 return i->user_index();
10033 template <
int structdim,
int dim,
int spacedim>
10036 Assert (
false, ExcInternalError());
10042 template <
typename Iterator>
10043 void set_user_index (
const Iterator &i,
10044 const unsigned int x)
10046 i->set_user_index(x);
10051 template <
int structdim,
int dim,
int spacedim>
10053 const unsigned int)
10055 Assert (
false, ExcInternalError());
10060 template <
int dim,
int spacedim>
10064 std::vector<unsigned int>::iterator i = v.begin();
10067 for (; line!=endl; ++line, ++i)
10068 *i = line->user_index();
10073 template <
int dim,
int spacedim>
10080 std::vector<unsigned int>::const_iterator i = v.begin();
10081 for (; line!=endl; ++line, ++i)
10082 line->set_user_index(*i);
10086 template <
int dim,
int spacedim>
10093 std::vector<unsigned int>::iterator i = v.begin();
10096 for (; quad!=endq; ++quad, ++i)
10097 *i = get_user_index(quad);
10103 template <
int dim,
int spacedim>
10112 std::vector<unsigned int>::const_iterator i = v.begin();
10113 for (; quad!=endq; ++quad, ++i)
10114 set_user_index(quad, *i);
10119 template <
int dim,
int spacedim>
10126 std::vector<unsigned int>::iterator i = v.begin();
10129 for (; hex!=endh; ++hex, ++i)
10130 *i = get_user_index(hex);
10136 template <
int dim,
int spacedim>
10145 std::vector<unsigned int>::const_iterator i = v.begin();
10146 for (; hex!=endh; ++hex, ++i)
10147 set_user_index(hex, *i);
10158 template <
typename Iterator>
10159 void *get_user_pointer (
const Iterator &i)
10161 return i->user_pointer();
10166 template <
int structdim,
int dim,
int spacedim>
10169 Assert (
false, ExcInternalError());
10175 template <
typename Iterator>
10176 void set_user_pointer (
const Iterator &i,
10179 i->set_user_pointer(x);
10184 template <
int structdim,
int dim,
int spacedim>
10188 Assert (
false, ExcInternalError());
10193 template <
int dim,
int spacedim>
10200 std::vector<void *> tmp;
10203 v.insert (v.end(), tmp.begin(), tmp.end());
10208 v.insert (v.end(), tmp.begin(), tmp.end());
10214 v.insert (v.end(), tmp.begin(), tmp.end());
10218 Assert (
false, ExcNotImplemented());
10223 template <
int dim,
int spacedim>
10227 std::vector<void *> tmp;
10231 tmp.insert (tmp.end(),
10232 v.begin(), v.begin()+
n_lines());
10239 tmp.insert (tmp.end(),
10247 tmp.insert (tmp.end(),
10253 Assert (
false, ExcNotImplemented());
10258 template <
int dim,
int spacedim>
10262 std::vector<void *>::iterator i = v.begin();
10265 for (; line!=endl; ++line, ++i)
10266 *i = line->user_pointer();
10271 template <
int dim,
int spacedim>
10278 std::vector<void *>::const_iterator i = v.begin();
10279 for (; line!=endl; ++line, ++i)
10280 line->set_user_pointer(*i);
10285 template <
int dim,
int spacedim>
10292 std::vector<void *>::iterator i = v.begin();
10295 for (; quad!=endq; ++quad, ++i)
10296 *i = get_user_pointer(quad);
10302 template <
int dim,
int spacedim>
10311 std::vector<void *>::const_iterator i = v.begin();
10312 for (; quad!=endq; ++quad, ++i)
10313 set_user_pointer(quad, *i);
10318 template <
int dim,
int spacedim>
10325 std::vector<void *>::iterator i = v.begin();
10328 for (; hex!=endh; ++hex, ++i)
10329 *i = get_user_pointer(hex);
10335 template <
int dim,
int spacedim>
10344 std::vector<void *>::const_iterator i = v.begin();
10345 for (; hex!=endh; ++hex, ++i)
10346 set_user_pointer(hex, *i);
10355 template <
int dim,
int spacedim>
10368 Assert (
false, ExcNotImplemented());
10375 template <
int dim,
int spacedim>
10388 Assert (
false, ExcImpossibleInDim(dim));
10395 template <
int dim,
int spacedim>
10408 Assert (
false, ExcNotImplemented());
10415 template <
int dim,
int spacedim>
10419 const unsigned int level =
levels.size()-1;
10422 if (
levels[level]->cells.cells.size() ==0)
10429 levels[level]->cells.cells.size()-1);
10432 if (ri->used()==
true)
10435 if (ri->used()==
true)
10442 template <
int dim,
int spacedim>
10452 if (cell->active()==
true)
10455 if (cell->active()==
true)
10463 template <
int dim,
int spacedim>
10474 template <
int dim,
int spacedim>
10479 if (level <
levels.size()-1)
10486 template <
int dim,
int spacedim>
10490 if (level <
levels.size()-1)
10491 return begin (level+1);
10497 template <
int dim,
int spacedim>
10502 return (level >=
levels.size()-1 ?
10509 template <
int dim,
int spacedim>
10519 template <
int dim,
int spacedim>
10530 template <
int dim,
int spacedim>
10541 template <
int dim,
int spacedim>
10554 template <
int dim,
int spacedim>
10561 Assert (
false, ExcImpossibleInDim(1));
10568 Assert (
false, ExcNotImplemented());
10575 template <
int dim,
int spacedim>
10582 Assert (
false, ExcImpossibleInDim(1));
10589 Assert (
false, ExcNotImplemented());
10596 template <
int dim,
int spacedim>
10603 Assert (
false, ExcImpossibleInDim(1));
10610 Assert (
false, ExcNotImplemented());
10619 template <
int dim,
int spacedim>
10620 typename Triangulation<dim,spacedim>::vertex_iterator
10627 Assert(
false, ExcNotImplemented());
10628 return raw_vertex_iterator();
10638 while (i->used() ==
false)
10647 template <
int dim,
int spacedim>
10648 typename Triangulation<dim,spacedim>::active_vertex_iterator
10656 template <
int dim,
int spacedim>
10657 typename Triangulation<dim,spacedim>::vertex_iterator
10662 Assert(
false, ExcNotImplemented());
10663 return raw_vertex_iterator();
10678 template <
int dim,
int spacedim>
10679 typename Triangulation<dim, spacedim>::raw_line_iterator
10687 if (level >=
levels.size() ||
levels[level]->cells.cells.size() == 0)
10695 Assert (level == 0, ExcFacesHaveNoLevel());
10703 template <
int dim,
int spacedim>
10704 typename Triangulation<dim, spacedim>::line_iterator
10711 while (ri->used() ==
false)
10719 template <
int dim,
int spacedim>
10720 typename Triangulation<dim, spacedim>::active_line_iterator
10727 while (i->has_children())
10735 template <
int dim,
int spacedim>
10736 typename Triangulation<dim, spacedim>::line_iterator
10749 template <
int dim,
int spacedim>
10750 typename Triangulation<dim,spacedim>::raw_quad_iterator
10756 Assert (
false, ExcImpossibleInDim(1));
10757 return raw_hex_iterator();
10762 if (level >=
levels.size() ||
levels[level]->cells.cells.size() == 0)
10772 Assert (level == 0, ExcFacesHaveNoLevel());
10781 Assert (
false, ExcNotImplemented());
10782 return raw_hex_iterator();
10788 template <
int dim,
int spacedim>
10789 typename Triangulation<dim,spacedim>::quad_iterator
10796 while (ri->used() ==
false)
10804 template <
int dim,
int spacedim>
10805 typename Triangulation<dim,spacedim>::active_quad_iterator
10812 while (i->has_children())
10820 template <
int dim,
int spacedim>
10821 typename Triangulation<dim,spacedim>::quad_iterator
10833 template <
int dim,
int spacedim>
10834 typename Triangulation<dim,spacedim>::raw_hex_iterator
10841 Assert (
false, ExcImpossibleInDim(1));
10842 return raw_hex_iterator();
10847 if (level >=
levels.size() ||
levels[level]->cells.cells.size() == 0)
10856 Assert (
false, ExcNotImplemented());
10857 return raw_hex_iterator();
10863 template <
int dim,
int spacedim>
10864 typename Triangulation<dim,spacedim>::hex_iterator
10871 while (ri->used() ==
false)
10879 template <
int dim,
int spacedim>
10880 typename Triangulation<dim, spacedim>::active_hex_iterator
10887 while (i->has_children())
10895 template <
int dim,
int spacedim>
10896 typename Triangulation<dim, spacedim>::hex_iterator
10965 template <
int dim,
int spacedim>
10968 return internal::Triangulation::n_cells (
number_cache);
10972 template <
int dim,
int spacedim>
10975 return internal::Triangulation::n_active_cells (
number_cache);
10978 template <
int dim,
int spacedim>
10986 template <
int dim,
int spacedim>
10998 Assert (
false, ExcNotImplemented());
11004 template <
int dim,
int spacedim>
11014 Assert (
false, ExcNotImplemented());
11020 template <
int dim,
int spacedim>
11032 Assert (
false, ExcNotImplemented());
11038 template <
int dim,
int spacedim>
11050 Assert (
false, ExcNotImplemented());
11057 template <
int dim,
int spacedim>
11069 Assert (
false, ExcNotImplemented());
11076 template <
int dim,
int spacedim>
11088 Assert (
false, ExcNotImplemented());
11094 template <
int dim,
int spacedim>
11105 template <
int dim,
int spacedim>
11117 return levels[level]->cells.cells.size();
11124 Assert(
false, ExcNotImplemented());
11134 return levels[level]->cells.cells.size();
11141 Assert(
false, ExcNotImplemented());
11150 return levels[level]->cells.cells.size();
11156 Assert(
false, ExcNotImplemented());
11162 template <
int dim,
int spacedim>
11165 Assert(
false, ExcFacesHaveNoLevel());
11170 template <
int dim,
int spacedim>
11173 return faces->lines.cells.size();
11177 template <
int dim,
int spacedim>
11181 ExcIndexRange (level, 0,
number_cache.n_lines_level.size()));
11182 Assert (dim == 1, ExcFacesHaveNoLevel());
11187 template <
int dim,
int spacedim>
11194 template <
int dim,
int spacedim>
11198 ExcIndexRange (level, 0,
number_cache.n_lines_level.size()));
11199 Assert (dim == 1, ExcFacesHaveNoLevel());
11334 template <
int dim,
int spacedim>
11341 template <
int dim,
int spacedim>
11344 Assert (dim == 2, ExcFacesHaveNoLevel());
11346 ExcIndexRange (level, 0,
number_cache.n_quads_level.size()));
11356 return levels[level]->cells.cells.size();
11365 return levels[level]->cells.cells.size();
11372 Assert(
false, ExcFacesHaveNoLevel());
11380 template <
int dim,
int spacedim>
11383 Assert (
false, ExcNotImplemented());
11392 return faces->quads.cells.size();
11397 template <
int dim,
int spacedim>
11404 template <
int dim,
int spacedim>
11408 ExcIndexRange (level, 0,
number_cache.n_quads_level.size()));
11409 Assert (dim == 2, ExcFacesHaveNoLevel());
11415 template <
int dim,
int spacedim>
11423 template <
int dim,
int spacedim>
11431 template <
int dim,
int spacedim>
11438 template <
int dim,
int spacedim>
11446 template <
int dim,
int spacedim>
11465 ExcIndexRange (level, 0,
number_cache.n_hexes_level.size()));
11476 return levels[level]->cells.cells.size();
11492 ExcIndexRange (level, 0,
number_cache.n_hexes_level.size()));
11499 template <
int dim,
int spacedim>
11504 std::bind2nd (std::equal_to<bool>(),
true));
11509 template <
int dim,
int spacedim>
11510 const std::vector<bool> &
11541 template <
int dim,
int spacedim>
11548 unsigned int max_vertex_index = 0;
11549 for (; cell!=endc; ++cell)
11550 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
11551 if (cell->vertex_index(vertex) > max_vertex_index)
11552 max_vertex_index = cell->vertex_index(vertex);
11558 std::vector<unsigned short int> usage_count (max_vertex_index+1, 0);
11562 for (cell=
begin(); cell!=endc; ++cell)
11563 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
11564 ++usage_count[cell->vertex_index(vertex)];
11567 static_cast<unsigned int>(*std::max_element (usage_count.begin(),
11568 usage_count.end())));
11573 template <
int dim,
int spacedim>
11582 template <
int dim,
int spacedim>
11591 template <
int dim,
int spacedim>
11600 template <
int dim,
int spacedim>
11610 Assert (satisfies_level1_at_vertex_rule (*
this) ==
true,
11611 ExcInternalError());
11618 const DistortedCellList
11624 if (
smooth_grid & limit_level_difference_at_vertices)
11625 Assert (satisfies_level1_at_vertex_rule (*
this) ==
true,
11626 ExcInternalError());
11630 update_neighbors(*
this);
11636 AssertThrow (cells_with_distorted_children.distorted_cells.size() == 0,
11637 cells_with_distorted_children);
11642 template <
int dim,
int spacedim>
11646 unsigned int active_cell_index = 0;
11648 if ((cell->used() ==
false) || cell->has_children())
11652 cell->set_active_cell_index (active_cell_index);
11653 ++active_cell_index;
11661 template<
int dim,
int spacedim>
11668 for (
unsigned int i=0; i<
levels.size(); ++i)
11684 template <
int dim,
int spacedim>
11688 const DistortedCellList
11689 cells_with_distorted_children
11701 for (
unsigned int level=0; level<
levels.size(); ++level)
11702 levels[level]->cells.monitor_memory (dim);
11715 while (cell != endc)
11716 Assert (!(cell++)->refine_flag_set(), ExcInternalError ());
11719 return cells_with_distorted_children;
11724 template <
int dim,
int spacedim>
11732 std::vector<unsigned int> line_cell_count = count_cells_bounded_by_line (*
this);
11733 std::vector<unsigned int> quad_cell_count = count_cells_bounded_by_quad (*
this);
11755 for (; cell!=endc; ++cell)
11756 if (!cell->active())
11757 if (cell->child(0)->coarsen_flag_set())
11759 cell->set_user_flag();
11760 for (
unsigned int child=0; child<cell->n_children(); ++child)
11762 Assert (cell->child(child)->coarsen_flag_set(),
11763 ExcInternalError());
11764 cell->child(child)->clear_coarsen_flag();
11786 for (cell =
last(); cell!=endc; --cell)
11787 if (cell->level()<=
static_cast<int>(
levels.size()-2) && cell->user_flag_set())
11807 for (cell=
begin(); cell!=endc; ++cell)
11808 Assert (cell->user_flag_set() ==
false, ExcInternalError());
11814 template <
int dim,
int spacedim>
11833 std::vector<int> vertex_level (
vertices.size(), 0);
11835 bool continue_iterating =
true;
11842 ExcMessage(
"In case of anisotropic refinement the " 11843 "limit_level_difference_at_vertices flag for " 11844 "mesh smoothing must not be set!"));
11848 std::fill (vertex_level.begin(), vertex_level.end(), 0);
11851 for (; cell!=endc; ++cell)
11853 if (cell->refine_flag_set())
11854 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
11856 vertex_level[cell->vertex_index(vertex)]
11857 = std::max (vertex_level[cell->vertex_index(vertex)],
11859 else if (!cell->coarsen_flag_set())
11860 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
11862 vertex_level[cell->vertex_index(vertex)]
11863 = std::max (vertex_level[cell->vertex_index(vertex)],
11873 Assert (cell->coarsen_flag_set(), ExcInternalError());
11874 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
11876 vertex_level[cell->vertex_index(vertex)]
11877 = std::max (vertex_level[cell->vertex_index(vertex)],
11893 if (cell->refine_flag_set() ==
false)
11895 for (
unsigned int vertex=0;
11896 vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
11897 if (vertex_level[cell->vertex_index(vertex)] >=
11901 cell->clear_coarsen_flag();
11906 if (vertex_level[cell->vertex_index(vertex)] >
11909 cell->set_refine_flag();
11911 for (
unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell;
11913 vertex_level[cell->vertex_index(v)]
11914 = std::max (vertex_level[cell->vertex_index(v)],
11942 for (; acell!=end_ac; ++acell)
11943 acell->clear_coarsen_flag();
11947 for (; cell!=endc; ++cell)
11950 if (cell->active())
11953 const unsigned int n_children=cell->n_children();
11954 unsigned int flagged_children=0;
11955 for (
unsigned int child=0; child<n_children; ++child)
11956 if (cell->child(child)->active() &&
11957 cell->child(child)->coarsen_flag_set())
11959 ++flagged_children;
11961 cell->child(child)->clear_coarsen_flag();
11966 if (flagged_children == n_children)
11967 cell->set_user_flag();
11973 for (cell=
begin(); cell!=endc; ++cell)
11974 Assert (cell->coarsen_flag_set() ==
false, ExcInternalError());
11994 for (cell=
last(); cell!=endc; --cell)
11995 if (cell->user_flag_set())
11998 if (internal::Triangulation::Implementation::template coarsening_allowed<dim,spacedim>(cell))
11999 for (
unsigned int c=0; c<cell->n_children(); ++c)
12001 Assert (cell->child(c)->refine_flag_set()==
false,
12002 ExcInternalError());
12004 cell->child(c)->set_coarsen_flag();
12017 continue_iterating = (current_coarsen_flags != previous_coarsen_flags);
12018 previous_coarsen_flags = current_coarsen_flags;
12020 while (continue_iterating ==
true);
12030 std::vector<bool> flags_before;
12036 std::vector<bool> flags_after;
12039 return (flags_before != flags_after);
12048 std::vector<bool> flags_before;
12054 std::vector<bool> flags_after;
12057 return (flags_before != flags_after);
12066 std::vector<bool> flags_before;
12072 std::vector<bool> flags_after;
12075 return (flags_before != flags_after);
12088 template <
int dim,
int spacedim>
12090 possibly_do_not_produce_unrefined_islands(
12093 Assert (cell->has_children(), ExcInternalError());
12095 unsigned int n_neighbors=0;
12098 unsigned int count=0;
12099 for (
unsigned int n=0; n<GeometryInfo<dim>::faces_per_cell; ++n)
12105 if (face_will_be_refined_by_neighbor(cell,n))
12112 if (count==n_neighbors ||
12113 (count>=n_neighbors-1 &&
12116 for (
unsigned int c=0; c<cell->n_children(); ++c)
12117 cell->child(c)->clear_coarsen_flag();
12119 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
12120 if (!cell->at_boundary(face)
12122 ( !cell->neighbor(face)->active() )
12123 && (cell_will_be_coarsened(cell->neighbor(face))) )
12124 possibly_do_not_produce_unrefined_islands<dim,spacedim>( cell->neighbor(face) );
12139 template <
int dim,
int spacedim>
12141 possibly_refine_unrefined_island
12145 Assert (cell->has_children() ==
false, ExcInternalError());
12146 Assert (cell->refine_flag_set() ==
false, ExcInternalError());
12158 if (allow_anisotropic_smoothing ==
false)
12161 unsigned int refined_neighbors = 0,
12162 unrefined_neighbors = 0;
12163 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
12164 if (!cell->at_boundary(face))
12166 if (face_will_be_refined_by_neighbor(cell,face))
12167 ++refined_neighbors;
12169 ++unrefined_neighbors;
12172 if (unrefined_neighbors < refined_neighbors)
12174 cell->clear_coarsen_flag();
12175 cell->set_refine_flag ();
12180 if (unrefined_neighbors > 0)
12181 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
12182 if (!cell->at_boundary(face)
12184 (face_will_be_refined_by_neighbor(cell,face) ==
false)
12186 (cell->neighbor(face)->has_children() ==
false)
12188 (cell->neighbor(face)->refine_flag_set() ==
false))
12189 possibly_refine_unrefined_island<dim,spacedim>
12190 (cell->neighbor(face),
12203 for (
unsigned int face_pair=0;
12204 face_pair<GeometryInfo<dim>::faces_per_cell/2; ++face_pair)
12212 for (
unsigned int face_index=0; face_index<2; ++face_index)
12214 unsigned int face=2*face_pair+face_index;
12221 face_will_be_refined_by_neighbor<dim,spacedim>(cell,face,expected_face_ref_case);
12235 directional_cell_refinement_case
12236 = (directional_cell_refinement_case &
12238 expected_face_ref_case,
12240 cell->face_orientation(face),
12241 cell->face_flip(face),
12242 cell->face_rotation(face)));
12248 Assert(directional_cell_refinement_case <
12250 ExcInternalError());
12251 smoothing_cell_refinement_case = smoothing_cell_refinement_case |
12252 directional_cell_refinement_case;
12257 if (smoothing_cell_refinement_case)
12259 cell->clear_coarsen_flag();
12260 cell->set_refine_flag(cell->refine_flag_set() |
12261 smoothing_cell_refinement_case);
12268 template <
int dim,
int spacedim>
12273 std::vector<bool> flags_before[2];
12293 std::vector<bool> flags_before_loop[2] = {flags_before[0],
12379 for (; cell!=endc; ++cell)
12380 cell->clear_coarsen_flag();
12383 bool mesh_changed_in_this_loop =
false;
12398 for (cell=
begin(); cell!=endc; ++cell)
12402 if (!cell->active() && cell_will_be_coarsened(cell))
12403 possibly_do_not_produce_unrefined_islands<dim,spacedim>(cell);
12434 for (cell=
begin(); cell!=endc; ++cell)
12435 if (!cell->active() ||
12437 cell->refine_flag_set() &&
12438 cell->is_locally_owned()))
12445 bool all_children_active =
true;
12446 if (!cell->active())
12447 for (
unsigned int c=0; c<cell->n_children(); ++c)
12448 if (!cell->child(c)->active() ||
12449 cell->child(c)->is_ghost() ||
12450 cell->child(c)->is_artificial())
12452 all_children_active =
false;
12456 if (all_children_active)
12463 unsigned int unrefined_neighbors = 0,
12464 total_neighbors = 0;
12466 for (
unsigned int n=0; n<GeometryInfo<dim>::faces_per_cell; ++n)
12473 if (!face_will_be_refined_by_neighbor(cell,n))
12474 ++unrefined_neighbors;
12490 if ((unrefined_neighbors == total_neighbors)
12497 (total_neighbors != 0))
12499 if (!cell->active())
12500 for (
unsigned int c=0; c<cell->n_children(); ++c)
12502 cell->child(c)->clear_refine_flag ();
12503 cell->child(c)->set_coarsen_flag ();
12506 cell->clear_refine_flag();
12524 ExcMessage(
"In case of anisotropic refinement the " 12525 "limit_level_difference_at_vertices flag for " 12526 "mesh smoothing must not be set!"));
12530 std::vector<int> vertex_level (
vertices.size(), 0);
12533 for (; cell!=endc; ++cell)
12535 if (cell->refine_flag_set())
12536 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
12538 vertex_level[cell->vertex_index(vertex)]
12539 = std::max (vertex_level[cell->vertex_index(vertex)],
12541 else if (!cell->coarsen_flag_set())
12542 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
12544 vertex_level[cell->vertex_index(vertex)]
12545 = std::max (vertex_level[cell->vertex_index(vertex)],
12553 Assert (cell->coarsen_flag_set(), ExcInternalError());
12554 for (
unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
12556 vertex_level[cell->vertex_index(vertex)]
12557 = std::max (vertex_level[cell->vertex_index(vertex)],
12573 if (cell->refine_flag_set() ==
false)
12575 for (
unsigned int vertex=0;
12576 vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
12577 if (vertex_level[cell->vertex_index(vertex)] >=
12581 cell->clear_coarsen_flag();
12586 if (vertex_level[cell->vertex_index(vertex)] >
12589 cell->set_refine_flag();
12591 for (
unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell;
12593 vertex_level[cell->vertex_index(v)]
12594 = std::max (vertex_level[cell->vertex_index(v)],
12621 for (; cell != endc; --cell)
12625 possibly_refine_unrefined_island<dim,spacedim>
12667 if (!cell->active())
12672 Assert (cell_is_patch_level_1(cell),
12673 ExcInternalError());
12674 if (cell->child(0)->has_children() ==
true)
12680 for (
unsigned int i=0; i<cell->n_children(); ++i)
12681 combined_ref_case = combined_ref_case |
12682 cell->child(i)->refine_flag_set();
12684 for (
unsigned int i=0; i<cell->n_children(); ++i)
12688 child->clear_coarsen_flag();
12689 child->set_refine_flag(combined_ref_case);
12710 cell->child(0)->active())
12717 const unsigned int n_children=cell->n_children();
12718 bool has_active_grandchildren =
false;
12720 for (
unsigned int i=0; i<n_children; ++i)
12721 if (cell->child(i)->child(0)->active())
12723 has_active_grandchildren =
true;
12727 if (has_active_grandchildren ==
false)
12733 unsigned int n_grandchildren=0;
12736 unsigned int n_coarsen_flags=0;
12741 for (
unsigned int c=0; c<n_children; ++c)
12748 const unsigned int nn_children=child->n_children();
12749 n_grandchildren += nn_children;
12754 if (child->child(0)->active())
12755 for (
unsigned int cc=0; cc<nn_children; ++cc)
12756 if (child->child(cc)->coarsen_flag_set())
12771 if ((n_coarsen_flags != n_grandchildren)
12773 (n_coarsen_flags > 0))
12774 for (
unsigned int c=0; c<n_children; ++c)
12777 if (child->child(0)->active())
12778 for (
unsigned int cc=0; cc<child->n_children(); ++cc)
12779 child->child(cc)->clear_coarsen_flag();
12806 bool changed =
true;
12813 for (; cell != endc; --cell)
12814 if (cell->refine_flag_set())
12817 for (
unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
12837 if (cell->neighbor(i)->active())
12839 if (cell->neighbor_is_coarser(i))
12841 if (cell->neighbor(i)->coarsen_flag_set())
12842 cell->neighbor(i)->clear_coarsen_flag();
12859 changed=cell->neighbor(i)->flag_for_face_refinement(cell->neighbor_of_coarser_neighbor(i).first,
12863 if (!cell->neighbor(i)->refine_flag_set())
12865 cell->neighbor(i)->set_refine_flag();
12988 std::pair<unsigned int, unsigned int> nb_indices
12989 =cell->neighbor_of_coarser_neighbor(i);
12990 unsigned int refined_along_x=0,
12992 to_be_refined_along_x=0,
12993 to_be_refined_along_y=0;
12995 const int this_face_index=cell->face_index(i);
12998 if ((this_face_index
12999 == cell->neighbor(i)->face(nb_indices.first)->child_index(0)) ||
13001 == cell->neighbor(i)->face(nb_indices.first)->child_index(1)))
13009 RefinementCase<dim-1> frc=cell->neighbor(i)->face(nb_indices.first)->refinement_case();
13027 cell->face_orientation(i),
13028 cell->face_flip(i),
13029 cell->face_rotation(i));
13031 ++to_be_refined_along_x;
13033 ++to_be_refined_along_y;
13037 cell->neighbor(i)->refine_flag_set())
13039 if (refined_along_x + to_be_refined_along_x > 1)
13040 changed |= cell->neighbor(i)->flag_for_face_refinement(nb_indices.first,
13042 if (refined_along_y + to_be_refined_along_y > 1)
13043 changed |= cell->neighbor(i)->flag_for_face_refinement(nb_indices.first,
13050 cell->neighbor(i)->set_refine_flag();
13058 nb->face_orientation(nb_indices.first),
13059 nb->face_flip(nb_indices.first),
13060 nb->face_rotation(nb_indices.first));
13061 if ((nb_frc & RefinementCase<dim>::cut_x) &&
13062 !(refined_along_x || to_be_refined_along_x))
13064 if ((nb_frc & RefinementCase<dim>::cut_y) &&
13065 !(refined_along_y || to_be_refined_along_y))
13071 cell->neighbor(i)->clear_coarsen_flag();
13072 const unsigned int nb_nb=cell->neighbor_of_neighbor(i);
13077 neighbor->face_orientation(nb_nb),
13078 neighbor->face_flip(nb_nb),
13079 neighbor->face_rotation(nb_nb));
13083 cell->face_orientation(i),
13084 cell->face_flip(i),
13085 cell->face_rotation(i));
13093 changed=cell->flag_for_face_refinement(i, face_ref_case);
13094 neighbor->flag_for_face_refinement(nb_nb, needed_face_ref_case);
13100 RefinementCase<dim-1> face_ref_case = cell->face(i)->refinement_case(),
13103 cell->face_orientation(i),
13104 cell->face_flip(i),
13105 cell->face_rotation(i));
13111 changed=cell->flag_for_face_refinement(i, face_ref_case);
13133 std::vector<bool> flags_after_loop[2];
13139 mesh_changed_in_this_loop
13140 = ((flags_before_loop[0] != flags_after_loop[0]) ||
13141 (flags_before_loop[1] != flags_after_loop[1]));
13145 flags_before_loop[0].swap(flags_after_loop[0]);
13146 flags_before_loop[1].swap(flags_after_loop[1]);
13148 while (mesh_changed_in_this_loop);
13154 return ((flags_before[0] != flags_before_loop[0]) ||
13155 (flags_before[1] != flags_before_loop[1]));
13161 template <
int dim,
int spacedim>
13163 const std::vector<bool> &v,
13164 const unsigned int magic_number2,
13167 const unsigned int N = v.size();
13168 unsigned char *flags =
new unsigned char[N/8+1];
13169 for (
unsigned int i=0; i<N/8+1; ++i) flags[i]=0;
13171 for (
unsigned int position=0; position<N; ++position)
13172 flags[position/8] |= (v[position] ? (1<<(position%8)) : 0);
13181 out << magic_number1 <<
' ' << N << std::endl;
13182 for (
unsigned int i=0; i<N/8+1; ++i)
13183 out << static_cast<unsigned int>(flags[i]) <<
' ';
13185 out << std::endl << magic_number2 << std::endl;
13193 template <
int dim,
int spacedim>
13195 std::vector<bool> &v,
13196 const unsigned int magic_number2,
13201 unsigned int magic_number;
13202 in >> magic_number;
13203 AssertThrow (magic_number==magic_number1, ExcGridReadError());
13209 unsigned char *flags =
new unsigned char[N/8+1];
13210 unsigned short int tmp;
13211 for (
unsigned int i=0; i<N/8+1; ++i)
13217 for (
unsigned int position=0; position!=N; ++position)
13218 v[position] = (flags[position/8] & (1<<(position%8)));
13220 in >> magic_number;
13221 AssertThrow (magic_number==magic_number2, ExcGridReadError());
13230 template <
int dim,
int spacedim>
13234 std::size_t mem = 0;
13236 for (
unsigned int i=0; i<
levels.size(); ++i)
13243 mem +=
sizeof (
faces);
13252 template<
int dim,
int spacedim>
13264 Assert(
false, ExcImpossibleInDim(1));
13275 Assert(
false, ExcImpossibleInDim(1));
13286 Assert(
false, ExcImpossibleInDim(2));
13295 #include "tria.inst" 13297 DEAL_II_NAMESPACE_CLOSE
std::vector< CellData< 1 > > boundary_lines
std::vector<::internal::Triangulation::TriaLevel< dim > * > levels
static void compute_number_cache(const Triangulation< dim, spacedim > &triangulation, const unsigned int level_objects, internal::Triangulation::NumberCache< 1 > &number_cache)
TriaActiveIterator< CellAccessor< dim, spacedim > > active_cell_iterator
boost::signals2::signal< void()> any_change
unsigned int n_active_cells() const
void set_boundary(const types::manifold_id number, const Boundary< dim, spacedim > &boundary_object)
unsigned int n_vertices() const
std::map< unsigned int, types::boundary_id > * vertex_to_boundary_id_map_1d
unsigned int n_used_vertices() const
const types::manifold_id flat_manifold_id
static const unsigned int invalid_unsigned_int
void load_user_flags_line(std::istream &in)
virtual ~DistortedCellList()
hex_iterator begin_hex(const unsigned int level=0) const
const types::subdomain_id invalid_subdomain_id
active_face_iterator begin_active_face() const
::internal::Triangulation::NumberCache< dim > number_cache
void save_user_flags_quad(std::ostream &out) const
virtual bool has_hanging_nodes() const
TriaRawIterator< CellAccessor< dim, spacedim > > raw_cell_iterator
unsigned int n_raw_cells(const unsigned int level) const
cell_iterator last() const
void save_user_pointers_hex(std::vector< void *> &v) const
vertex_iterator begin_vertex() const
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
unsigned int n_cells() const
static void create_triangulation(const std::vector< Point< spacedim > > &v, const std::vector< CellData< 2 > > &cells, const SubCellData &subcelldata, Triangulation< 2, spacedim > &triangulation)
int face(const unsigned int i) const
unsigned int n_raw_hexs() const
static void delete_children(Triangulation< 1, spacedim > &triangulation, typename Triangulation< 1, spacedim >::cell_iterator &cell, std::vector< unsigned int > &, std::vector< unsigned int > &)
std::list< typename Triangulation< dim, spacedim >::cell_iterator > distorted_cells
std::vector< Point< spacedim > > vertices
boost::signals2::signal< void()> clear
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
static void compute_number_cache(const Triangulation< dim, spacedim > &triangulation, const unsigned int level_objects, internal::Triangulation::NumberCache< 3 > &number_cache)
IteratorRange< active_cell_iterator > active_cell_iterators() const
std::vector< types::boundary_id > get_boundary_indicators() const DEAL_II_DEPRECATED
void set_all_refine_flags()
void load_user_flags(std::istream &in)
::ExceptionBase & ExcMessage(std::string arg1)
static void compute_number_cache(const Triangulation< dim, spacedim > &triangulation, const unsigned int level_objects, internal::Triangulation::NumberCache< 2 > &number_cache)
bool anisotropic_refinement
std::vector< unsigned int > n_active_lines_level
std::vector< unsigned int > n_hexes_level
DeclException1(ExcGridHasInvalidCell, int,<< "Something went wrong when making cell "<< arg1<< ". Read the docs and the source code "<< "for more information.")
boost::signals2::signal< void(const typename Triangulation< dim, spacedim >::cell_iterator &cell)> pre_coarsening_on_cell
static void read_bool_vector(const unsigned int magic_number1, std::vector< bool > &v, const unsigned int magic_number2, std::istream &in)
line_iterator begin_line(const unsigned int level=0) const
line_iterator end_line() const
virtual void copy_triangulation(const Triangulation< dim, spacedim > &old_tria)
void clear_user_flags_quad()
unsigned int n_faces() const
std::vector< unsigned int > n_quads_level
void flip_all_direction_flags()
static RefinementCase< dim-1 > face_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static Triangulation< 3, spacedim >::DistortedCellList execute_refinement(Triangulation< 3, spacedim > &triangulation, const bool check_for_distorted_cells)
active_cell_iterator begin_active(const unsigned int level=0) const
#define AssertThrow(cond, exc)
static Triangulation< 1, spacedim >::DistortedCellList execute_refinement(Triangulation< 1, spacedim > &triangulation, const bool)
static void create_children(Triangulation< 2, spacedim > &triangulation, unsigned int &next_unused_vertex, typename Triangulation< 2, spacedim >::raw_line_iterator &next_unused_line, typename Triangulation< 2, spacedim >::raw_cell_iterator &next_unused_cell, typename Triangulation< 2, spacedim >::cell_iterator &cell)
static void write_bool_vector(const unsigned int magic_number1, const std::vector< bool > &v, const unsigned int magic_number2, std::ostream &out)
unsigned int n_active_hexs() const
void save_user_pointers_quad(std::vector< void *> &v) const
std::vector< unsigned int > n_active_quads_level
void execute_coarsening()
const Boundary< dim, spacedim > & get_boundary(const types::manifold_id number) const
cell_iterator begin(const unsigned int level=0) const
void load_user_flags_quad(std::istream &in)
unsigned int n_levels() const
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
raw_hex_iterator begin_raw_hex(const unsigned int level=0) const
void load_user_pointers_line(const std::vector< void *> &v)
unsigned int n_active_faces() const
raw_quad_iterator begin_raw_quad(const unsigned int level=0) const
cell_iterator end() const
static void create_triangulation(const std::vector< Point< spacedim > > &v, const std::vector< CellData< 3 > > &cells, const SubCellData &subcelldata, Triangulation< 3, spacedim > &triangulation)
void load_coarsen_flags(std::istream &out)
void load_user_flags_hex(std::istream &in)
DeclException5(ExcInteriorQuadCantBeBoundary, int, int, int, int, types::boundary_id,<< "The input data for creating a triangulation contained "<< "information about a quad with indices "<< arg1<< ", "<< arg2<< ", "<< arg3<< ", and "<< arg4<< " that is described to have boundary indicator "<<(int) arg5<< ". However, this is an internal quad not located on the "<< "boundary. You cannot assign a boundary indicator to it."<< std::endl<< std::endl<< "If this happened at a place where you call "<< "Triangulation::create_triangulation() yourself, you need "<< "to check the SubCellData object you pass to this function."<< std::endl<< std::endl<< "If this happened in a place where you are reading a mesh "<< "from a file, then you need to investigate why such a quad "<< "ended up in the input file. A typical case is a geometry "<< "that consisted of multiple parts and for which the mesh "<< "generator program assumes that the interface between "<< "two parts is a boundary when that isn't supposed to be "<< "the case, or where the mesh generator simply assigns "<< "'geometry indicators' to quads at the surface of "<< "a part that are not supposed to be interpreted as "<< "'boundary indicators'.")
active_quad_iterator begin_active_quad(const unsigned int level=0) const
virtual void execute_coarsening_and_refinement()
unsigned int n_active_lines() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
virtual bool prepare_coarsening_and_refinement()
bool check_consistency(const unsigned int dim) const
void reset_active_cell_indices()
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
DistortedCellList execute_refinement()
const bool check_for_distorted_cells
boost::signals2::signal< void()> pre_refinement
void save_user_flags_line(std::ostream &out) const
void save_user_pointers_line(std::vector< void *> &v) const
static bool coarsening_allowed(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
Triangulation(const MeshSmoothing smooth_grid=none, const bool check_for_distorted_cells=false)
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
void clear_user_flags_hex()
void load_user_indices_line(const std::vector< unsigned int > &v)
unsigned int global_dof_index
static unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
vertex_iterator end_vertex() const
unsigned int n_raw_quads() const
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim-1 > &face_refinement_case=RefinementCase< dim-1 >::isotropic_refinement)
#define Assert(cond, exc)
active_hex_iterator begin_active_hex(const unsigned int level=0) const
void save_user_indices_line(std::vector< unsigned int > &v) const
std::map< types::manifold_id, SmartPointer< const Manifold< dim, spacedim >, Triangulation< dim, spacedim > > > manifold
void set_all_manifold_ids(const types::manifold_id number)
boost::signals2::signal< void()> create
static const StraightBoundary< dim, spacedim > straight_boundary
unsigned int max_adjacent_cells() const
bool get_anisotropic_refinement_flag() const
unsigned int n_quads() const
void clear_user_data(const unsigned int i)
std::vector< bool > vertices_used
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std_cxx11::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std_cxx11::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std_cxx11::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
DeclException4(ExcQuadInexistant, int, int, int, int,<< "While trying to assign a boundary indicator to a quad: "<< "the quad with bounding lines "<< arg1<< ", "<< arg2<< ", "<< arg3<< ", "<< arg4<< " does not exist.")
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
void save_coarsen_flags(std::ostream &out) const
void save_user_indices(std::vector< unsigned int > &v) const
TriaObjects< TriaObject< 1 > > lines
void load_user_indices(const std::vector< unsigned int > &v)
void save_refine_flags(std::ostream &out) const
unsigned int n_lines() const
void clear_despite_subscriptions()
raw_cell_iterator begin_raw(const unsigned int level=0) const
unsigned int n_raw_lines() const
virtual std::size_t memory_consumption() const
std::map< unsigned int, types::manifold_id > * vertex_to_manifold_id_map_1d
virtual void set_mesh_smoothing(const MeshSmoothing mesh_smoothing)
void load_user_pointers_hex(const std::vector< void *> &v)
void save_user_flags_hex(std::ostream &out) const
face_iterator begin_face() const
unsigned int subdomain_id
boost::signals2::signal< void()> post_refinement
virtual types::global_dof_index n_global_active_cells() const
unsigned int n_hexs() const
DeclException2(ExcLineInexistant, int, int,<< "While trying to assign a boundary indicator to a line: "<< "the line with end vertices "<< arg1<< " and "<< arg2<< " does not exist.")
void set_face(const unsigned int i, const int index)
unsigned int n_active_quads() const
TriaObjects< TriaObject< 1 > > lines
active_line_iterator begin_active_line(const unsigned int level=0) const
raw_cell_iterator end_raw(const unsigned int level) const
unsigned int n_raw_faces() const
void load_user_indices_hex(const std::vector< unsigned int > &v)
active_vertex_iterator begin_active_vertex() const
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim-dim, spacedim >(&forms)[vertices_per_cell])
DeclException3(ExcInvalidVertexIndex, int, int, int,<< "Error while creating cell "<< arg1<< ": the vertex index "<< arg2<< " must be between 0 and "<< arg3<< ".")
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::vector< types::boundary_id > get_boundary_ids() const
const types::manifold_id invalid_manifold_id
std::vector< unsigned int > n_active_hexes_level
void load_user_pointers_quad(const std::vector< void *> &v)
active_cell_iterator end_active(const unsigned int level) const
static void prevent_distorted_boundary_cells(const Triangulation< 1, spacedim > &)
::internal::Triangulation::TriaFaces< dim > * faces
void save_user_indices_quad(std::vector< unsigned int > &v) const
unsigned int n_active_quads
const std::vector< bool > & get_used_vertices() const
virtual unsigned int n_global_levels() const
hex_iterator end_hex() const
std::vector< CellData< 2 > > boundary_quads
MeshSmoothing smooth_grid
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
void load_user_pointers(const std::vector< void *> &v)
void refine_global(const unsigned int times=1)
void save_user_pointers(std::vector< void *> &v) const
void save_user_flags(std::ostream &out) const
Iterator points to a valid object.
quad_iterator end_quad() const
IteratorRange< cell_iterator > cell_iterators() const
unsigned char boundary_id
face_iterator end_face() const
void set_all_manifold_ids_on_boundary(const types::manifold_id number)
active_cell_iterator last_active() const
const types::boundary_id internal_face_boundary_id
void clear_user_flags_line()
IteratorState::IteratorStates state() const
static Triangulation< 2, spacedim >::DistortedCellList execute_refinement(Triangulation< 2, spacedim > &triangulation, const bool check_for_distorted_cells)
quad_iterator begin_quad(const unsigned int level=0) const
static void prepare_refinement_dim_dependent(const Triangulation< dim, spacedim > &)
std::vector< types::manifold_id > get_manifold_ids() const
bool vertex_used(const unsigned int index) const
unsigned int n_active_hexes
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
virtual void create_triangulation_compatibility(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
static RefinementCase< dim > min_cell_refinement_case_for_face_refinement(const RefinementCase< dim-1 > &face_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
void load_user_indices_quad(const std::vector< unsigned int > &v)
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_lines
raw_line_iterator begin_raw_line(const unsigned int level=0) const
boost::signals2::signal< void(const Triangulation< dim, spacedim > &destination_tria)> copy
void load_refine_flags(std::istream &in)
std::vector< unsigned int > n_lines_level
Triangulation< dim, spacedim > & get_triangulation()
static void create_triangulation(const std::vector< Point< spacedim > > &v, const std::vector< CellData< 1 > > &cells, const SubCellData &, Triangulation< 1, spacedim > &triangulation)
void save_user_indices_hex(std::vector< unsigned int > &v) const