Reference documentation for deal.II version Git bed997f895 2020-09-22 11:49:20 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
19 
20 #include <deal.II/fe/mapping_q1.h>
21 
25 #include <deal.II/grid/manifold.h>
26 #include <deal.II/grid/tria.h>
31 
33 #include <deal.II/lac/vector.h>
34 
35 #include <algorithm>
36 #include <array>
37 #include <cmath>
38 #include <functional>
39 #include <list>
40 #include <map>
41 #include <memory>
42 #include <numeric>
43 
44 
46 
47 
48 namespace internal
49 {
50  namespace TriangulationImplementation
51  {
53  : n_levels(0)
54  , n_lines(0)
55  , n_active_lines(0)
56  // all other fields are
57  // default constructed
58  {}
59 
60 
61 
62  std::size_t
64  {
65  return (MemoryConsumption::memory_consumption(n_levels) +
69  MemoryConsumption::memory_consumption(n_active_lines_level));
70  }
71 
72 
74  : n_quads(0)
75  , n_active_quads(0)
76  // all other fields are
77  // default constructed
78  {}
79 
80 
81 
82  std::size_t
84  {
89  MemoryConsumption::memory_consumption(n_active_quads_level));
90  }
91 
92 
93 
95  : n_hexes(0)
96  , n_active_hexes(0)
97  // all other fields are
98  // default constructed
99  {}
100 
101 
102 
103  std::size_t
105  {
109  MemoryConsumption::memory_consumption(n_active_hexes) +
110  MemoryConsumption::memory_consumption(n_active_hexes_level));
111  }
112  } // namespace TriangulationImplementation
113 } // namespace internal
114 
115 // anonymous namespace for internal helper functions
116 namespace
117 {
118  // return whether the given cell is
119  // patch_level_1, i.e. determine
120  // whether either all or none of
121  // its children are further
122  // refined. this function can only
123  // be called for non-active cells.
124  template <int dim, int spacedim>
125  bool
126  cell_is_patch_level_1(
128  {
129  Assert(cell->is_active() == false, ExcInternalError());
130 
131  unsigned int n_active_children = 0;
132  for (unsigned int i = 0; i < cell->n_children(); ++i)
133  if (cell->child(i)->is_active())
134  ++n_active_children;
135 
136  return (n_active_children == 0) ||
137  (n_active_children == cell->n_children());
138  }
139 
140 
141 
142  // return, whether a given @p cell will be
143  // coarsened, which is the case if all
144  // children are active and have their coarsen
145  // flag set. In case only part of the coarsen
146  // flags are set, remove them.
147  template <int dim, int spacedim>
148  bool
149  cell_will_be_coarsened(
151  {
152  // only cells with children should be
153  // considered for coarsening
154 
155  if (cell->has_children())
156  {
157  unsigned int children_to_coarsen = 0;
158  const unsigned int n_children = cell->n_children();
159 
160  for (unsigned int c = 0; c < n_children; ++c)
161  if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
162  ++children_to_coarsen;
163  if (children_to_coarsen == n_children)
164  return true;
165  else
166  for (unsigned int c = 0; c < n_children; ++c)
167  if (cell->child(c)->is_active())
168  cell->child(c)->clear_coarsen_flag();
169  }
170  // no children, so no coarsening
171  // possible. however, no children also
172  // means that this cell will be in the same
173  // state as if it had children and was
174  // coarsened. So, what should we return -
175  // false or true?
176  // make sure we do not have to do this at
177  // all...
178  Assert(cell->has_children(), ExcInternalError());
179  // ... and then simply return false
180  return false;
181  }
182 
183 
184  // return, whether the face @p face_no of the
185  // given @p cell will be refined after the
186  // current refinement step, considering
187  // refine and coarsen flags and considering
188  // only those refinemnts that will be caused
189  // by the neighboring cell.
190 
191  // this function is used on both active cells
192  // and cells with children. on cells with
193  // children it also of interest to know 'how'
194  // the face will be refined. thus there is an
195  // additional third argument @p
196  // expected_face_ref_case returning just
197  // that. be aware, that this variable will
198  // only contain useful information if this
199  // function is called for an active cell.
200  //
201  // thus, this is an internal function, users
202  // should call one of the two alternatives
203  // following below.
204  template <int dim, int spacedim>
205  bool
206  face_will_be_refined_by_neighbor_internal(
208  const unsigned int face_no,
209  RefinementCase<dim - 1> &expected_face_ref_case)
210  {
211  // first of all: set the default value for
212  // expected_face_ref_case, which is no
213  // refinement at all
214  expected_face_ref_case = RefinementCase<dim - 1>::no_refinement;
215 
216  const typename Triangulation<dim, spacedim>::cell_iterator neighbor =
217  cell->neighbor(face_no);
218 
219  // If we are at the boundary, there is no
220  // neighbor which could refine the face
221  if (neighbor.state() != IteratorState::valid)
222  return false;
223 
224  if (neighbor->has_children())
225  {
226  // if the neighbor is refined, it may be
227  // coarsened. if so, then it won't refine
228  // the face, no matter what else happens
229  if (cell_will_be_coarsened(neighbor))
230  return false;
231  else
232  // if the neighbor is refined, then it
233  // is also refined at our current
234  // face. He will stay so without
235  // coarsening, so return true in that
236  // case.
237  {
238  expected_face_ref_case = cell->face(face_no)->refinement_case();
239  return true;
240  }
241  }
242 
243  // now, the neighbor is not refined, but
244  // perhaps it will be
245  const RefinementCase<dim> nb_ref_flag = neighbor->refine_flag_set();
246  if (nb_ref_flag != RefinementCase<dim>::no_refinement)
247  {
248  // now we need to know, which of the
249  // neighbors faces points towards us
250  const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
251  // check, whether the cell will be
252  // refined in a way that refines our
253  // face
254  const RefinementCase<dim - 1> face_ref_case =
256  nb_ref_flag,
257  neighbor_neighbor,
258  neighbor->face_orientation(neighbor_neighbor),
259  neighbor->face_flip(neighbor_neighbor),
260  neighbor->face_rotation(neighbor_neighbor));
261  if (face_ref_case != RefinementCase<dim - 1>::no_refinement)
262  {
264  neighbor_face = neighbor->face(neighbor_neighbor);
265  const int this_face_index = cell->face_index(face_no);
266 
267  // there are still two basic
268  // possibilities here: the neighbor
269  // might be coarser or as coarse
270  // as we are
271  if (neighbor_face->index() == this_face_index)
272  // the neighbor is as coarse as
273  // we are and will be refined at
274  // the face of consideration, so
275  // return true
276  {
277  expected_face_ref_case = face_ref_case;
278  return true;
279  }
280  else
281  {
282  // the neighbor is coarser.
283  // this is the most complicated
284  // case. It might be, that the
285  // neighbor's face will be
286  // refined, but that we will
287  // not see this, as we are
288  // refined in a similar way.
289 
290  // so, the neighbor's face must
291  // have children. check, if our
292  // cell's face is one of these
293  // (it could also be a
294  // grand_child)
295  for (unsigned int c = 0; c < neighbor_face->n_children(); ++c)
296  if (neighbor_face->child_index(c) == this_face_index)
297  {
298  // if the flagged refine
299  // case of the face is a
300  // subset or the same as
301  // the current refine case,
302  // then the face, as seen
303  // from our cell, won't be
304  // refined by the neighbor
305  if ((neighbor_face->refinement_case() | face_ref_case) ==
306  neighbor_face->refinement_case())
307  return false;
308  else
309  {
310  // if we are active, we
311  // must be an
312  // anisotropic child
313  // and the coming
314  // face_ref_case is
315  // isotropic. Thus,
316  // from our cell we
317  // will see exactly the
318  // opposite refine case
319  // that the face has
320  // now...
321  Assert(
322  face_ref_case ==
324  ExcInternalError());
325  expected_face_ref_case =
326  ~neighbor_face->refinement_case();
327  return true;
328  }
329  }
330 
331  // so, obviously we were not
332  // one of the children, but a
333  // grandchild. This is only
334  // possible in 3d.
335  Assert(dim == 3, ExcInternalError());
336  // In that case, however, no
337  // matter what the neighbor
338  // does, it won't be finer
339  // after the next refinement
340  // step.
341  return false;
342  }
343  } // if face will be refined
344  } // if neighbor is flagged for refinement
345 
346  // no cases left, so the neighbor will not
347  // refine the face
348  return false;
349  }
350 
351  // version of above function for both active
352  // and non-active cells
353  template <int dim, int spacedim>
354  bool
355  face_will_be_refined_by_neighbor(
357  const unsigned int face_no)
358  {
359  RefinementCase<dim - 1> dummy = RefinementCase<dim - 1>::no_refinement;
360  return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
361  }
362 
363  // version of above function for active cells
364  // only. Additionally returning the refine
365  // case (to come) of the face under
366  // consideration
367  template <int dim, int spacedim>
368  bool
369  face_will_be_refined_by_neighbor(
371  const unsigned int face_no,
372  RefinementCase<dim - 1> &expected_face_ref_case)
373  {
374  return face_will_be_refined_by_neighbor_internal(cell,
375  face_no,
376  expected_face_ref_case);
377  }
378 
379 
380 
381  template <int dim, int spacedim>
382  bool
383  satisfies_level1_at_vertex_rule(
385  {
386  std::vector<unsigned int> min_adjacent_cell_level(
387  triangulation.n_vertices(), triangulation.n_levels());
388  std::vector<unsigned int> max_adjacent_cell_level(
389  triangulation.n_vertices(), 0);
390 
391  for (const auto &cell : triangulation.active_cell_iterators())
392  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
393  {
394  min_adjacent_cell_level[cell->vertex_index(v)] =
395  std::min<unsigned int>(
396  min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
397  max_adjacent_cell_level[cell->vertex_index(v)] =
398  std::max<unsigned int>(
399  min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
400  }
401 
402  for (unsigned int k = 0; k < triangulation.n_vertices(); ++k)
403  if (triangulation.vertex_used(k))
404  if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
405  return false;
406  return true;
407  }
408 
409 
410 
417  template <int dim, int spacedim>
418  std::vector<unsigned int>
419  count_cells_bounded_by_line(const Triangulation<dim, spacedim> &triangulation)
420  {
421  if (dim >= 2)
422  {
423  std::vector<unsigned int> line_cell_count(triangulation.n_raw_lines(),
424  0);
425  for (const auto &cell : triangulation.cell_iterators())
426  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
427  ++line_cell_count[cell->line_index(l)];
428  return line_cell_count;
429  }
430  else
431  return std::vector<unsigned int>();
432  }
433 
434 
435 
442  template <int dim, int spacedim>
443  std::vector<unsigned int>
444  count_cells_bounded_by_quad(const Triangulation<dim, spacedim> &triangulation)
445  {
446  if (dim >= 3)
447  {
448  std::vector<unsigned int> quad_cell_count(triangulation.n_raw_quads(),
449  0);
450  for (const auto &cell : triangulation.cell_iterators())
451  for (unsigned int q : GeometryInfo<dim>::face_indices())
452  ++quad_cell_count[cell->quad_index(q)];
453  return quad_cell_count;
454  }
455  else
456  return std::vector<unsigned int>();
457  }
458 
459 
460 
472  void
473  reorder_compatibility(const std::vector<CellData<1>> &, const SubCellData &)
474  {
475  // nothing to do here: the format
476  // hasn't changed for 1d
477  }
478 
479 
480  void reorder_compatibility(std::vector<CellData<2>> &cells,
481  const SubCellData &)
482  {
483  for (auto &cell : cells)
484  std::swap(cell.vertices[2], cell.vertices[3]);
485  }
486 
487 
488  void reorder_compatibility(std::vector<CellData<3>> &cells,
489  SubCellData & subcelldata)
490  {
491  unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
492  for (auto &cell : cells)
493  {
494  for (const unsigned int i : GeometryInfo<3>::vertex_indices())
495  tmp[i] = cell.vertices[i];
496  for (const unsigned int i : GeometryInfo<3>::vertex_indices())
497  cell.vertices[GeometryInfo<3>::ucd_to_deal[i]] = tmp[i];
498  }
499 
500  // now points in boundary quads
501  std::vector<CellData<2>>::iterator boundary_quad =
502  subcelldata.boundary_quads.begin();
503  std::vector<CellData<2>>::iterator end_quad =
504  subcelldata.boundary_quads.end();
505  for (unsigned int quad_no = 0; boundary_quad != end_quad;
506  ++boundary_quad, ++quad_no)
507  std::swap(boundary_quad->vertices[2], boundary_quad->vertices[3]);
508  }
509 
510 
511 
529  template <int dim, int spacedim>
530  unsigned int
531  middle_vertex_index(
532  const typename Triangulation<dim, spacedim>::line_iterator &line)
533  {
534  if (line->has_children())
535  return line->child(0)->vertex_index(1);
537  }
538 
539 
540  template <int dim, int spacedim>
541  unsigned int
542  middle_vertex_index(
543  const typename Triangulation<dim, spacedim>::quad_iterator &quad)
544  {
545  switch (static_cast<unsigned char>(quad->refinement_case()))
546  {
548  return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
549  break;
551  return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
552  break;
554  return quad->child(0)->vertex_index(3);
555  break;
556  default:
557  break;
558  }
560  }
561 
562 
563  template <int dim, int spacedim>
564  unsigned int
565  middle_vertex_index(
566  const typename Triangulation<dim, spacedim>::hex_iterator &hex)
567  {
568  switch (static_cast<unsigned char>(hex->refinement_case()))
569  {
571  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
572  break;
574  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
575  break;
577  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
578  break;
580  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
581  break;
583  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
584  break;
586  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
587  break;
589  return hex->child(0)->vertex_index(7);
590  break;
591  default:
592  break;
593  }
595  }
596 
597 
610  template <class TRIANGULATION>
611  inline typename TRIANGULATION::DistortedCellList
612  collect_distorted_coarse_cells(const TRIANGULATION &)
613  {
614  return typename TRIANGULATION::DistortedCellList();
615  }
616 
617 
618 
627  template <int dim>
629  collect_distorted_coarse_cells(const Triangulation<dim, dim> &triangulation)
630  {
631  typename Triangulation<dim, dim>::DistortedCellList distorted_cells;
632  for (const auto &cell : triangulation.cell_iterators_on_level(0))
633  {
635  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
636  vertices[i] = cell->vertex(i);
637 
639  GeometryInfo<dim>::alternating_form_at_vertices(vertices, determinants);
640 
641  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
642  if (determinants[i] <= 1e-9 * std::pow(cell->diameter(), 1. * dim))
643  {
644  distorted_cells.distorted_cells.push_back(cell);
645  break;
646  }
647  }
648 
649  return distorted_cells;
650  }
651 
652 
659  template <int dim>
660  bool
661  has_distorted_children(
662  const typename Triangulation<dim, dim>::cell_iterator &cell,
663  std::integral_constant<int, dim>,
664  std::integral_constant<int, dim>)
665  {
666  Assert(cell->has_children(), ExcInternalError());
667 
668  for (unsigned int c = 0; c < cell->n_children(); ++c)
669  {
671  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
672  vertices[i] = cell->child(c)->vertex(i);
673 
675  GeometryInfo<dim>::alternating_form_at_vertices(vertices, determinants);
676 
677  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
678  if (determinants[i] <=
679  1e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
680  return true;
681  }
682 
683  return false;
684  }
685 
686 
694  template <int dim, int spacedim>
695  bool
696  has_distorted_children(
698  std::integral_constant<int, dim>,
699  std::integral_constant<int, spacedim>)
700  {
701  return false;
702  }
703 
704 
705 
710  template <int spacedim>
711  void update_neighbors(Triangulation<1, spacedim> &)
712  {}
713 
714 
715  template <int dim, int spacedim>
716  void
717  update_neighbors(Triangulation<dim, spacedim> &triangulation)
718  {
719  // each face can be neighbored on two sides
720  // by cells. according to the face's
721  // intrinsic normal we define the left
722  // neighbor as the one for which the face
723  // normal points outward, and store that
724  // one first; the second one is then
725  // the right neighbor for which the
726  // face normal points inward. This
727  // information depends on the type of cell
728  // and local number of face for the
729  // 'standard ordering and orientation' of
730  // faces and then on the face_orientation
731  // information for the real mesh. Set up a
732  // table to have fast access to those
733  // offsets (0 for left and 1 for
734  // right). Some of the values are invalid
735  // as they reference too large face
736  // numbers, but we just leave them at a
737  // zero value.
738  //
739  // Note, that in 2d for lines as faces the
740  // normal direction given in the
741  // GeometryInfo class is not consistent. We
742  // thus define here that the normal for a
743  // line points to the right if the line
744  // points upwards.
745  //
746  // There is one more point to
747  // consider, however: if we have
748  // dim<spacedim, then we may have
749  // cases where cells are
750  // inverted. In effect, both
751  // cells think they are the left
752  // neighbor of an edge, for
753  // example, which leads us to
754  // forget neighborship
755  // information (a case that shows
756  // this is
757  // codim_one/hanging_nodes_02). We
758  // store whether a cell is
759  // inverted using the
760  // direction_flag, so if a cell
761  // has a false direction_flag,
762  // then we need to invert our
763  // selection whether we are a
764  // left or right neighbor in all
765  // following computations.
766  //
767  // first index: dimension (minus 2)
768  // second index: local face index
769  // third index: face_orientation (false and true)
770  static const unsigned int left_right_offset[2][6][2] = {
771  // quadrilateral
772  {{0, 1}, // face 0, face_orientation = false and true
773  {1, 0}, // face 1, face_orientation = false and true
774  {1, 0}, // face 2, face_orientation = false and true
775  {0, 1}, // face 3, face_orientation = false and true
776  {0, 0}, // face 4, invalid face
777  {0, 0}}, // face 5, invalid face
778  // hexahedron
779  {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
780 
781  // now create a vector of the two active
782  // neighbors (left and right) for each face
783  // and fill it by looping over all cells. For
784  // cases with anisotropic refinement and more
785  // then one cell neighboring at a given side
786  // of the face we will automatically get the
787  // active one on the highest level as we loop
788  // over cells from lower levels first.
789  const typename Triangulation<dim, spacedim>::cell_iterator dummy;
790  std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
791  adjacent_cells(2 * triangulation.n_raw_faces(), dummy);
792 
793  for (const auto &cell : triangulation.cell_iterators())
794  for (auto f : GeometryInfo<dim>::face_indices())
795  {
796  const typename Triangulation<dim, spacedim>::face_iterator face =
797  cell->face(f);
798 
799  const unsigned int offset =
800  (cell->direction_flag() ?
801  left_right_offset[dim - 2][f][cell->face_orientation(f)] :
802  1 - left_right_offset[dim - 2][f][cell->face_orientation(f)]);
803 
804  adjacent_cells[2 * face->index() + offset] = cell;
805 
806  // if this cell is not refined, but the
807  // face is, then we'll have to set our
808  // cell as neighbor for the child faces
809  // as well. Fortunately the normal
810  // orientation of children will be just
811  // the same.
812  if (dim == 2)
813  {
814  if (cell->is_active() && face->has_children())
815  {
816  adjacent_cells[2 * face->child(0)->index() + offset] = cell;
817  adjacent_cells[2 * face->child(1)->index() + offset] = cell;
818  }
819  }
820  else // -> dim == 3
821  {
822  // We need the same as in 2d
823  // here. Furthermore, if the face is
824  // refined with cut_x or cut_y then
825  // those children again in the other
826  // direction, and if this cell is
827  // refined isotropically (along the
828  // face) then the neighbor will
829  // (probably) be refined as cut_x or
830  // cut_y along the face. For those
831  // neighboring children cells, their
832  // neighbor will be the current,
833  // inactive cell, as our children are
834  // too fine to be neighbors. Catch that
835  // case by also acting on inactive
836  // cells with isotropic refinement
837  // along the face. If the situation
838  // described is not present, the data
839  // will be overwritten later on when we
840  // visit cells on finer levels, so no
841  // harm will be done.
842  if (face->has_children() &&
843  (cell->is_active() ||
845  cell->refinement_case(), f) ==
847  {
848  for (unsigned int c = 0; c < face->n_children(); ++c)
849  adjacent_cells[2 * face->child(c)->index() + offset] = cell;
850  if (face->child(0)->has_children())
851  {
852  adjacent_cells[2 * face->child(0)->child(0)->index() +
853  offset] = cell;
854  adjacent_cells[2 * face->child(0)->child(1)->index() +
855  offset] = cell;
856  }
857  if (face->child(1)->has_children())
858  {
859  adjacent_cells[2 * face->child(1)->child(0)->index() +
860  offset] = cell;
861  adjacent_cells[2 * face->child(1)->child(1)->index() +
862  offset] = cell;
863  }
864  } // if cell active and face refined
865  } // else -> dim==3
866  } // for all faces of all cells
867 
868  // now loop again over all cells and set the
869  // corresponding neighbor cell. Note, that we
870  // have to use the opposite of the
871  // left_right_offset in this case as we want
872  // the offset of the neighbor, not our own.
873  for (const auto &cell : triangulation.cell_iterators())
874  for (auto f : GeometryInfo<dim>::face_indices())
875  {
876  const unsigned int offset =
877  (cell->direction_flag() ?
878  left_right_offset[dim - 2][f][cell->face_orientation(f)] :
879  1 - left_right_offset[dim - 2][f][cell->face_orientation(f)]);
880  cell->set_neighbor(
881  f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
882  }
883  }
884 
885 
886  template <int dim, int spacedim>
887  void
888  update_periodic_face_map_recursively(
889  const typename Triangulation<dim, spacedim>::cell_iterator &cell_1,
890  const typename Triangulation<dim, spacedim>::cell_iterator &cell_2,
891  unsigned int n_face_1,
892  unsigned int n_face_2,
893  const std::bitset<3> & orientation,
894  typename std::map<
896  unsigned int>,
897  std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
898  unsigned int>,
899  std::bitset<3>>> &periodic_face_map)
900  {
901  using FaceIterator = typename Triangulation<dim, spacedim>::face_iterator;
902  const FaceIterator face_1 = cell_1->face(n_face_1);
903  const FaceIterator face_2 = cell_2->face(n_face_2);
904 
905  const bool face_orientation = orientation[0];
906  const bool face_flip = orientation[1];
907  const bool face_rotation = orientation[2];
908 
909  Assert((dim != 1) || (face_orientation == true && face_flip == false &&
910  face_rotation == false),
911  ExcMessage("The supplied orientation "
912  "(face_orientation, face_flip, face_rotation) "
913  "is invalid for 1D"));
914 
915  Assert((dim != 2) || (face_orientation == true && face_rotation == false),
916  ExcMessage("The supplied orientation "
917  "(face_orientation, face_flip, face_rotation) "
918  "is invalid for 2D"));
919 
920  Assert(face_1 != face_2, ExcMessage("face_1 and face_2 are equal!"));
921 
922  Assert(face_1->at_boundary() && face_2->at_boundary(),
923  ExcMessage("Periodic faces must be on the boundary"));
924 
925  // Check if the requirement that each edge can only have at most one hanging
926  // node, and as a consequence neighboring cells can differ by at most
927  // one refinement level is enforced. In 1d, there are no hanging nodes and
928  // so neighboring cells can differ by more than one refinement level.
929  Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
930  ExcInternalError());
931 
932  // insert periodic face pair for both cells
933  using CellFace =
934  std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
935  unsigned int>;
936  const CellFace cell_face_1(cell_1, n_face_1);
937  const CellFace cell_face_2(cell_2, n_face_2);
938  const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
939  cell_face_2, orientation);
940 
941  const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
942  periodic_faces(cell_face_1, cell_face_orientation_2);
943 
944  // Only one periodic neighbor is allowed
945  Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
946  periodic_face_map.insert(periodic_faces);
947 
948  if (dim == 1)
949  {
950  if (cell_1->has_children())
951  {
952  if (cell_2->has_children())
953  {
954  update_periodic_face_map_recursively<dim, spacedim>(
955  cell_1->child(n_face_1),
956  cell_2->child(n_face_2),
957  n_face_1,
958  n_face_2,
959  orientation,
960  periodic_face_map);
961  }
962  else // only face_1 has children
963  {
964  update_periodic_face_map_recursively<dim, spacedim>(
965  cell_1->child(n_face_1),
966  cell_2,
967  n_face_1,
968  n_face_2,
969  orientation,
970  periodic_face_map);
971  }
972  }
973  }
974  else // dim == 2 || dim == 3
975  {
976  // A lookup table on how to go through the child cells depending on the
977  // orientation:
978  // see Documentation of GeometryInfo for details
979 
980  static const int lookup_table_2d[2][2] =
981  // flip:
982  {
983  {0, 1}, // false
984  {1, 0} // true
985  };
986 
987  static const int lookup_table_3d[2][2][2][4] =
988  // orientation flip rotation
989  {{{
990  {0, 2, 1, 3}, // false false false
991  {2, 3, 0, 1} // false false true
992  },
993  {
994  {3, 1, 2, 0}, // false true false
995  {1, 0, 3, 2} // false true true
996  }},
997  {{
998  {0, 1, 2, 3}, // true false false
999  {1, 3, 0, 2} // true false true
1000  },
1001  {
1002  {3, 2, 1, 0}, // true true false
1003  {2, 0, 3, 1} // true true true
1004  }}};
1005 
1006  if (cell_1->has_children())
1007  {
1008  if (cell_2->has_children())
1009  {
1010  // In the case that both faces have children, we loop over all
1011  // children and apply update_periodic_face_map_recursively
1012  // recursively:
1013 
1014  Assert(face_1->n_children() ==
1016  face_2->n_children() ==
1018  ExcNotImplemented());
1019 
1020  for (unsigned int i = 0;
1021  i < GeometryInfo<dim>::max_children_per_face;
1022  ++i)
1023  {
1024  // Lookup the index for the second face
1025  unsigned int j = 0;
1026  switch (dim)
1027  {
1028  case 2:
1029  j = lookup_table_2d[face_flip][i];
1030  break;
1031  case 3:
1032  j = lookup_table_3d[face_orientation][face_flip]
1033  [face_rotation][i];
1034  break;
1035  default:
1036  AssertThrow(false, ExcNotImplemented());
1037  }
1038 
1039  // find subcell ids that belong to the subface indices
1040  unsigned int child_cell_1 =
1042  cell_1->refinement_case(),
1043  n_face_1,
1044  i,
1045  cell_1->face_orientation(n_face_1),
1046  cell_1->face_flip(n_face_1),
1047  cell_1->face_rotation(n_face_1),
1048  face_1->refinement_case());
1049  unsigned int child_cell_2 =
1051  cell_2->refinement_case(),
1052  n_face_2,
1053  j,
1054  cell_2->face_orientation(n_face_2),
1055  cell_2->face_flip(n_face_2),
1056  cell_2->face_rotation(n_face_2),
1057  face_2->refinement_case());
1058 
1059  Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
1060  face_1->child(i),
1061  ExcInternalError());
1062  Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
1063  face_2->child(j),
1064  ExcInternalError());
1065 
1066  // precondition: subcell has the same orientation as cell
1067  // (so that the face numbers coincide) recursive call
1068  update_periodic_face_map_recursively<dim, spacedim>(
1069  cell_1->child(child_cell_1),
1070  cell_2->child(child_cell_2),
1071  n_face_1,
1072  n_face_2,
1073  orientation,
1074  periodic_face_map);
1075  }
1076  }
1077  else // only face_1 has children
1078  {
1079  for (unsigned int i = 0;
1080  i < GeometryInfo<dim>::max_children_per_face;
1081  ++i)
1082  {
1083  // find subcell ids that belong to the subface indices
1084  unsigned int child_cell_1 =
1086  cell_1->refinement_case(),
1087  n_face_1,
1088  i,
1089  cell_1->face_orientation(n_face_1),
1090  cell_1->face_flip(n_face_1),
1091  cell_1->face_rotation(n_face_1),
1092  face_1->refinement_case());
1093 
1094  // recursive call
1095  update_periodic_face_map_recursively<dim, spacedim>(
1096  cell_1->child(child_cell_1),
1097  cell_2,
1098  n_face_1,
1099  n_face_2,
1100  orientation,
1101  periodic_face_map);
1102  }
1103  }
1104  }
1105  }
1106  }
1107 
1108 
1109 } // end of anonymous namespace
1110 
1111 
1112 namespace internal
1113 {
1114  namespace TriangulationImplementation
1115  {
1116  // make sure that if in the following we
1117  // write Triangulation<dim,spacedim>
1118  // we mean the *class*
1119  // ::Triangulation, not the
1120  // enclosing namespace
1121  // internal::TriangulationImplementation
1123 
1129  int,
1130  << "Something went wrong when making cell " << arg1
1131  << ". Read the docs and the source code "
1132  << "for more information.");
1138  int,
1139  << "Something went wrong upon construction of cell "
1140  << arg1);
1151  int,
1152  << "Cell " << arg1
1153  << " has negative measure. This typically "
1154  << "indicates some distortion in the cell, or a mistakenly "
1155  << "swapped pair of vertices in the input to "
1156  << "Triangulation::create_triangulation().");
1165  int,
1166  int,
1167  int,
1168  << "Error while creating cell " << arg1
1169  << ": the vertex index " << arg2 << " must be between 0 and "
1170  << arg3 << ".");
1176  int,
1177  int,
1178  << "While trying to assign a boundary indicator to a line: "
1179  << "the line with end vertices " << arg1 << " and " << arg2
1180  << " does not exist.");
1186  int,
1187  int,
1188  int,
1189  int,
1190  << "While trying to assign a boundary indicator to a quad: "
1191  << "the quad with bounding lines " << arg1 << ", " << arg2
1192  << ", " << arg3 << ", " << arg4 << " does not exist.");
1199  int,
1200  int,
1202  << "The input data for creating a triangulation contained "
1203  << "information about a line with indices " << arg1 << " and " << arg2
1204  << " that is described to have boundary indicator "
1205  << static_cast<int>(arg3)
1206  << ". However, this is an internal line not located on the "
1207  << "boundary. You cannot assign a boundary indicator to it." << std::endl
1208  << std::endl
1209  << "If this happened at a place where you call "
1210  << "Triangulation::create_triangulation() yourself, you need "
1211  << "to check the SubCellData object you pass to this function."
1212  << std::endl
1213  << std::endl
1214  << "If this happened in a place where you are reading a mesh "
1215  << "from a file, then you need to investigate why such a line "
1216  << "ended up in the input file. A typical case is a geometry "
1217  << "that consisted of multiple parts and for which the mesh "
1218  << "generator program assumes that the interface between "
1219  << "two parts is a boundary when that isn't supposed to be "
1220  << "the case, or where the mesh generator simply assigns "
1221  << "'geometry indicators' to lines at the perimeter of "
1222  << "a part that are not supposed to be interpreted as "
1223  << "'boundary indicators'.");
1230  int,
1231  int,
1232  int,
1233  int,
1235  << "The input data for creating a triangulation contained "
1236  << "information about a quad with indices " << arg1 << ", " << arg2
1237  << ", " << arg3 << ", and " << arg4
1238  << " that is described to have boundary indicator "
1239  << static_cast<int>(arg5)
1240  << ". However, this is an internal quad not located on the "
1241  << "boundary. You cannot assign a boundary indicator to it." << std::endl
1242  << std::endl
1243  << "If this happened at a place where you call "
1244  << "Triangulation::create_triangulation() yourself, you need "
1245  << "to check the SubCellData object you pass to this function."
1246  << std::endl
1247  << std::endl
1248  << "If this happened in a place where you are reading a mesh "
1249  << "from a file, then you need to investigate why such a quad "
1250  << "ended up in the input file. A typical case is a geometry "
1251  << "that consisted of multiple parts and for which the mesh "
1252  << "generator program assumes that the interface between "
1253  << "two parts is a boundary when that isn't supposed to be "
1254  << "the case, or where the mesh generator simply assigns "
1255  << "'geometry indicators' to quads at the surface of "
1256  << "a part that are not supposed to be interpreted as "
1257  << "'boundary indicators'.");
1264  int,
1265  int,
1266  << "In SubCellData the line info of the line with vertex indices " << arg1
1267  << " and " << arg2 << " appears more than once. "
1268  << "This is not allowed.");
1275  int,
1276  int,
1277  std::string,
1278  << "In SubCellData the line info of the line with vertex indices " << arg1
1279  << " and " << arg2 << " appears multiple times with different (valid) "
1280  << arg3 << ". This is not allowed.");
1287  int,
1288  int,
1289  int,
1290  int,
1291  std::string,
1292  << "In SubCellData the quad info of the quad with line indices " << arg1
1293  << ", " << arg2 << ", " << arg3 << " and " << arg4
1294  << " appears multiple times with different (valid) " << arg5
1295  << ". This is not allowed.");
1296 
1297  /*
1298  * Reserve space for TriaFaces. Details:
1299  *
1300  * Reserve space for line_orientations.
1301  *
1302  * @note Used only for dim=3.
1303  */
1304  void
1305  reserve_space(TriaFaces & tria_faces,
1306  const unsigned int new_quads_in_pairs,
1307  const unsigned int new_quads_single)
1308  {
1309  AssertDimension(tria_faces.dim, 3);
1310 
1311  Assert(new_quads_in_pairs % 2 == 0, ExcInternalError());
1312 
1313  unsigned int next_free_single = 0;
1314  unsigned int next_free_pair = 0;
1315 
1316  // count the number of objects, of unused single objects and of
1317  // unused pairs of objects
1318  unsigned int n_quads = 0;
1319  unsigned int n_unused_pairs = 0;
1320  unsigned int n_unused_singles = 0;
1321  for (unsigned int i = 0; i < tria_faces.quads.used.size(); ++i)
1322  {
1323  if (tria_faces.quads.used[i])
1324  ++n_quads;
1325  else if (i + 1 < tria_faces.quads.used.size())
1326  {
1327  if (tria_faces.quads.used[i + 1])
1328  {
1329  ++n_unused_singles;
1330  if (next_free_single == 0)
1331  next_free_single = i;
1332  }
1333  else
1334  {
1335  ++n_unused_pairs;
1336  if (next_free_pair == 0)
1337  next_free_pair = i;
1338  ++i;
1339  }
1340  }
1341  else
1342  ++n_unused_singles;
1343  }
1344  Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
1345  tria_faces.quads.used.size(),
1346  ExcInternalError());
1347 
1348  // how many single quads are needed in addition to n_unused_quads?
1349  const int additional_single_quads = new_quads_single - n_unused_singles;
1350 
1351  unsigned int new_size =
1352  tria_faces.quads.used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
1353  if (additional_single_quads > 0)
1354  new_size += additional_single_quads;
1355 
1356  // see above...
1357  if (new_size > tria_faces.quads.n_objects())
1358  {
1359  // reserve the field of the derived class
1360  tria_faces.quads_line_orientations.reserve(
1361  new_size * GeometryInfo<2>::lines_per_cell);
1362  tria_faces.quads_line_orientations.insert(
1363  tria_faces.quads_line_orientations.end(),
1364  new_size * GeometryInfo<2>::lines_per_cell -
1365  tria_faces.quads_line_orientations.size(),
1366  true);
1367 
1368  tria_faces.quad_reference_cell_type.reserve(new_size);
1369  tria_faces.quad_reference_cell_type.insert(
1370  tria_faces.quad_reference_cell_type.end(),
1371  new_size - tria_faces.quad_reference_cell_type.size(),
1373  }
1374  }
1375 
1376 
1377 
1391  void
1392  reserve_space(TriaLevel & tria_level,
1393  const unsigned int total_cells,
1394  const unsigned int dimension,
1395  const unsigned int space_dimension)
1396  {
1397  // we need space for total_cells cells. Maybe we have more already
1398  // with those cells which are unused, so only allocate new space if
1399  // needed.
1400  //
1401  // note that all arrays should have equal sizes (checked by
1402  // @p{monitor_memory}
1403  if (total_cells > tria_level.refine_flags.size())
1404  {
1405  tria_level.refine_flags.reserve(total_cells);
1406  tria_level.refine_flags.insert(tria_level.refine_flags.end(),
1407  total_cells -
1408  tria_level.refine_flags.size(),
1409  /*RefinementCase::no_refinement=*/0);
1410 
1411  tria_level.coarsen_flags.reserve(total_cells);
1412  tria_level.coarsen_flags.insert(tria_level.coarsen_flags.end(),
1413  total_cells -
1414  tria_level.coarsen_flags.size(),
1415  false);
1416 
1417  tria_level.active_cell_indices.reserve(total_cells);
1418  tria_level.active_cell_indices.insert(
1419  tria_level.active_cell_indices.end(),
1420  total_cells - tria_level.active_cell_indices.size(),
1422 
1423  tria_level.subdomain_ids.reserve(total_cells);
1424  tria_level.subdomain_ids.insert(tria_level.subdomain_ids.end(),
1425  total_cells -
1426  tria_level.subdomain_ids.size(),
1427  0);
1428 
1429  tria_level.level_subdomain_ids.reserve(total_cells);
1430  tria_level.level_subdomain_ids.insert(
1431  tria_level.level_subdomain_ids.end(),
1432  total_cells - tria_level.level_subdomain_ids.size(),
1433  0);
1434 
1435  if (dimension < space_dimension)
1436  {
1437  tria_level.direction_flags.reserve(total_cells);
1438  tria_level.direction_flags.insert(
1439  tria_level.direction_flags.end(),
1440  total_cells - tria_level.direction_flags.size(),
1441  true);
1442  }
1443  else
1444  tria_level.direction_flags.clear();
1445 
1446  tria_level.parents.reserve((total_cells + 1) / 2);
1447  tria_level.parents.insert(tria_level.parents.end(),
1448  (total_cells + 1) / 2 -
1449  tria_level.parents.size(),
1450  -1);
1451 
1452  tria_level.neighbors.reserve(total_cells * (2 * dimension));
1453  tria_level.neighbors.insert(tria_level.neighbors.end(),
1454  total_cells * (2 * dimension) -
1455  tria_level.neighbors.size(),
1456  std::make_pair(-1, -1));
1457 
1458 
1459  if (tria_level.dim == 3)
1460  {
1461  tria_level.face_orientations.reserve(
1462  total_cells * GeometryInfo<3>::faces_per_cell);
1463  tria_level.face_orientations.insert(
1464  tria_level.face_orientations.end(),
1465  total_cells * GeometryInfo<3>::faces_per_cell -
1466  tria_level.face_orientations.size(),
1467  true);
1468  }
1469 
1470  if (tria_level.dim == 2 || tria_level.dim == 3)
1471  {
1472  tria_level.reference_cell_type.reserve(total_cells);
1473  tria_level.reference_cell_type.insert(
1474  tria_level.reference_cell_type.end(),
1475  total_cells - tria_level.reference_cell_type.size(),
1476  tria_level.dim == 2 ? ReferenceCell::Type::Quad :
1478  }
1479  }
1480  }
1481 
1482 
1483 
1488  int,
1489  int,
1490  << "The containers have sizes " << arg1 << " and " << arg2
1491  << ", which is not as expected.");
1492 
1498  void
1499  monitor_memory(const TriaLevel & tria_level,
1500  const unsigned int true_dimension)
1501  {
1502  (void)tria_level;
1503  (void)true_dimension;
1504  Assert(2 * true_dimension * tria_level.refine_flags.size() ==
1505  tria_level.neighbors.size(),
1506  ExcMemoryInexact(tria_level.refine_flags.size(),
1507  tria_level.neighbors.size()));
1508  Assert(2 * true_dimension * tria_level.coarsen_flags.size() ==
1509  tria_level.neighbors.size(),
1510  ExcMemoryInexact(tria_level.coarsen_flags.size(),
1511  tria_level.neighbors.size()));
1512  }
1513 
1514 
1515 
1528  void
1529  reserve_space(TriaObjects & tria_objects,
1530  const unsigned int new_objects_in_pairs,
1531  const unsigned int new_objects_single = 0)
1532  {
1533  if (tria_objects.structdim <= 2)
1534  {
1535  Assert(new_objects_in_pairs % 2 == 0, ExcInternalError());
1536 
1537  tria_objects.next_free_single = 0;
1538  tria_objects.next_free_pair = 0;
1539  tria_objects.reverse_order_next_free_single = false;
1540 
1541  // count the number of objects, of unused single objects and of
1542  // unused pairs of objects
1543  unsigned int n_objects = 0;
1544  unsigned int n_unused_pairs = 0;
1545  unsigned int n_unused_singles = 0;
1546  for (unsigned int i = 0; i < tria_objects.used.size(); ++i)
1547  {
1548  if (tria_objects.used[i])
1549  ++n_objects;
1550  else if (i + 1 < tria_objects.used.size())
1551  {
1552  if (tria_objects.used[i + 1])
1553  {
1554  ++n_unused_singles;
1555  if (tria_objects.next_free_single == 0)
1556  tria_objects.next_free_single = i;
1557  }
1558  else
1559  {
1560  ++n_unused_pairs;
1561  if (tria_objects.next_free_pair == 0)
1562  tria_objects.next_free_pair = i;
1563  ++i;
1564  }
1565  }
1566  else
1567  ++n_unused_singles;
1568  }
1569  Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
1570  tria_objects.used.size(),
1571  ExcInternalError());
1572 
1573  // how many single objects are needed in addition to
1574  // n_unused_objects?
1575  const int additional_single_objects =
1576  new_objects_single - n_unused_singles;
1577 
1578  unsigned int new_size = tria_objects.used.size() +
1579  new_objects_in_pairs - 2 * n_unused_pairs;
1580  if (additional_single_objects > 0)
1581  new_size += additional_single_objects;
1582 
1583  // only allocate space if necessary
1584  if (new_size > tria_objects.n_objects())
1585  {
1586  unsigned int faces_per_cell = 1;
1587  unsigned int max_children_per_cell = 1;
1588 
1589  if (tria_objects.structdim == 1)
1590  faces_per_cell = GeometryInfo<1>::faces_per_cell;
1591  else if (tria_objects.structdim == 2)
1592  faces_per_cell = GeometryInfo<2>::faces_per_cell;
1593  else if (tria_objects.structdim == 3)
1594  faces_per_cell = GeometryInfo<3>::faces_per_cell;
1595  else
1596  AssertThrow(false, ExcNotImplemented());
1597 
1598  if (tria_objects.structdim == 1)
1599  max_children_per_cell = GeometryInfo<1>::max_children_per_cell;
1600  else if (tria_objects.structdim == 2)
1601  max_children_per_cell = GeometryInfo<2>::max_children_per_cell;
1602  else if (tria_objects.structdim == 3)
1603  max_children_per_cell = GeometryInfo<3>::max_children_per_cell;
1604  else
1605  AssertThrow(false, ExcNotImplemented());
1606 
1607  tria_objects.cells.reserve(new_size * faces_per_cell);
1608  tria_objects.cells.insert(tria_objects.cells.end(),
1609  (new_size - tria_objects.n_objects()) *
1610  faces_per_cell,
1611  -1);
1612 
1613  tria_objects.used.reserve(new_size);
1614  tria_objects.used.insert(tria_objects.used.end(),
1615  new_size - tria_objects.used.size(),
1616  false);
1617 
1618  tria_objects.user_flags.reserve(new_size);
1619  tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1620  new_size -
1621  tria_objects.user_flags.size(),
1622  false);
1623 
1624  const unsigned int factor = max_children_per_cell / 2;
1625  tria_objects.children.reserve(factor * new_size);
1626  tria_objects.children.insert(tria_objects.children.end(),
1627  factor * new_size -
1628  tria_objects.children.size(),
1629  -1);
1630 
1631  if (tria_objects.structdim > 1)
1632  {
1633  tria_objects.refinement_cases.reserve(new_size);
1634  tria_objects.refinement_cases.insert(
1635  tria_objects.refinement_cases.end(),
1636  new_size - tria_objects.refinement_cases.size(),
1637  /*RefinementCase::no_refinement=*/0);
1638  }
1639 
1640  // first reserve, then resize. Otherwise the std library can
1641  // decide to allocate more entries.
1642  tria_objects.boundary_or_material_id.reserve(new_size);
1643  tria_objects.boundary_or_material_id.resize(new_size);
1644 
1645  tria_objects.user_data.reserve(new_size);
1646  tria_objects.user_data.resize(new_size);
1647 
1648  tria_objects.manifold_id.reserve(new_size);
1649  tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1650  new_size -
1651  tria_objects.manifold_id.size(),
1653  }
1654 
1655  if (n_unused_singles == 0)
1656  {
1657  tria_objects.next_free_single = new_size - 1;
1658  tria_objects.reverse_order_next_free_single = true;
1659  }
1660  }
1661  else
1662  {
1663  const unsigned int new_hexes = new_objects_in_pairs;
1664 
1665  const unsigned int new_size =
1666  new_hexes + std::count(tria_objects.used.begin(),
1667  tria_objects.used.end(),
1668  true);
1669 
1670  // see above...
1671  if (new_size > tria_objects.n_objects())
1672  {
1673  unsigned int faces_per_cell = 1;
1674 
1675  if (tria_objects.structdim == 1)
1676  faces_per_cell = GeometryInfo<1>::faces_per_cell;
1677  else if (tria_objects.structdim == 2)
1678  faces_per_cell = GeometryInfo<2>::faces_per_cell;
1679  else if (tria_objects.structdim == 3)
1680  faces_per_cell = GeometryInfo<3>::faces_per_cell;
1681  else
1682  AssertThrow(false, ExcNotImplemented());
1683 
1684  tria_objects.cells.reserve(new_size * faces_per_cell);
1685  tria_objects.cells.insert(tria_objects.cells.end(),
1686  (new_size - tria_objects.n_objects()) *
1687  faces_per_cell,
1688  -1);
1689 
1690  tria_objects.used.reserve(new_size);
1691  tria_objects.used.insert(tria_objects.used.end(),
1692  new_size - tria_objects.used.size(),
1693  false);
1694 
1695  tria_objects.user_flags.reserve(new_size);
1696  tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1697  new_size -
1698  tria_objects.user_flags.size(),
1699  false);
1700 
1701  tria_objects.children.reserve(4 * new_size);
1702  tria_objects.children.insert(tria_objects.children.end(),
1703  4 * new_size -
1704  tria_objects.children.size(),
1705  -1);
1706 
1707  // for the following fields, we know exactly how many elements
1708  // we need, so first reserve then resize (resize itself, at least
1709  // with some compiler libraries, appears to round up the size it
1710  // actually reserves)
1711  tria_objects.boundary_or_material_id.reserve(new_size);
1712  tria_objects.boundary_or_material_id.resize(new_size);
1713 
1714  tria_objects.manifold_id.reserve(new_size);
1715  tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1716  new_size -
1717  tria_objects.manifold_id.size(),
1719 
1720  tria_objects.user_data.reserve(new_size);
1721  tria_objects.user_data.resize(new_size);
1722 
1723  tria_objects.refinement_cases.reserve(new_size);
1724  tria_objects.refinement_cases.insert(
1725  tria_objects.refinement_cases.end(),
1726  new_size - tria_objects.refinement_cases.size(),
1727  /*RefinementCase::no_refinement=*/0);
1728  }
1729  tria_objects.next_free_single = tria_objects.next_free_pair = 0;
1730  }
1731  }
1732 
1733 
1734 
1740  void
1741  monitor_memory(const TriaObjects &tria_object, const unsigned int)
1742  {
1743  Assert(tria_object.n_objects() == tria_object.used.size(),
1744  ExcMemoryInexact(tria_object.n_objects(),
1745  tria_object.used.size()));
1746  Assert(tria_object.n_objects() == tria_object.user_flags.size(),
1747  ExcMemoryInexact(tria_object.n_objects(),
1748  tria_object.user_flags.size()));
1749  Assert(tria_object.n_objects() ==
1750  tria_object.boundary_or_material_id.size(),
1751  ExcMemoryInexact(tria_object.n_objects(),
1752  tria_object.boundary_or_material_id.size()));
1753  Assert(tria_object.n_objects() == tria_object.manifold_id.size(),
1754  ExcMemoryInexact(tria_object.n_objects(),
1755  tria_object.manifold_id.size()));
1756  Assert(tria_object.n_objects() == tria_object.user_data.size(),
1757  ExcMemoryInexact(tria_object.n_objects(),
1758  tria_object.user_data.size()));
1759 
1760  if (tria_object.structdim == 1)
1761  {
1762  Assert(1 * tria_object.n_objects() == tria_object.children.size(),
1763  ExcMemoryInexact(tria_object.n_objects(),
1764  tria_object.children.size()));
1765  }
1766  else if (tria_object.structdim == 2)
1767  {
1768  Assert(2 * tria_object.n_objects() == tria_object.children.size(),
1769  ExcMemoryInexact(tria_object.n_objects(),
1770  tria_object.children.size()));
1771  }
1772  else if (tria_object.structdim == 3)
1773  {
1774  Assert(4 * tria_object.n_objects() == tria_object.children.size(),
1775  ExcMemoryInexact(tria_object.n_objects(),
1776  tria_object.children.size()));
1777  }
1778  }
1779 
1876  {
1888  template <int dim, int spacedim>
1889  static void
1892  const unsigned int level_objects,
1894  {
1895  using line_iterator =
1897 
1898  number_cache.n_levels = 0;
1899  if (level_objects > 0)
1900  // find the last level on which there are used cells
1901  for (unsigned int level = 0; level < level_objects; ++level)
1902  if (triangulation.begin(level) != triangulation.end(level))
1903  number_cache.n_levels = level + 1;
1904 
1905  // no cells at all?
1906  Assert(number_cache.n_levels > 0, ExcInternalError());
1907 
1909  // update the number of lines on the different levels in the
1910  // cache
1911  number_cache.n_lines = 0;
1912  number_cache.n_active_lines = 0;
1913 
1914  // for 1d, lines have levels so take count the objects per
1915  // level and globally
1916  if (dim == 1)
1917  {
1918  number_cache.n_lines_level.resize(number_cache.n_levels);
1919  number_cache.n_active_lines_level.resize(number_cache.n_levels);
1920 
1921  for (unsigned int level = 0; level < number_cache.n_levels; ++level)
1922  {
1923  // count lines on this level
1924  number_cache.n_lines_level[level] = 0;
1925  number_cache.n_active_lines_level[level] = 0;
1926 
1927  line_iterator line = triangulation.begin_line(level),
1928  endc =
1929  (level == number_cache.n_levels - 1 ?
1930  line_iterator(triangulation.end_line()) :
1931  triangulation.begin_line(level + 1));
1932  for (; line != endc; ++line)
1933  {
1934  ++number_cache.n_lines_level[level];
1935  if (line->has_children() == false)
1936  ++number_cache.n_active_lines_level[level];
1937  }
1938 
1939  // update total number of lines
1940  number_cache.n_lines += number_cache.n_lines_level[level];
1941  number_cache.n_active_lines +=
1942  number_cache.n_active_lines_level[level];
1943  }
1944  }
1945  else
1946  {
1947  // for dim>1, there are no levels for lines
1948  number_cache.n_lines_level.clear();
1949  number_cache.n_active_lines_level.clear();
1950 
1951  line_iterator line = triangulation.begin_line(),
1952  endc = triangulation.end_line();
1953  for (; line != endc; ++line)
1954  {
1955  ++number_cache.n_lines;
1956  if (line->has_children() == false)
1957  ++number_cache.n_active_lines;
1958  }
1959  }
1960  }
1961 
1976  template <int dim, int spacedim>
1977  static void
1980  const unsigned int level_objects,
1982  {
1983  // update lines and n_levels in number_cache. since we don't
1984  // access any of these numbers, we can do this in the
1985  // background
1986  Threads::Task<void> update_lines = Threads::new_task(
1987  static_cast<
1988  void (*)(const Triangulation<dim, spacedim> &,
1989  const unsigned int,
1991  &compute_number_cache<dim, spacedim>),
1992  triangulation,
1993  level_objects,
1995  number_cache));
1996 
1997  using quad_iterator =
1999 
2001  // update the number of quads on the different levels in the
2002  // cache
2003  number_cache.n_quads = 0;
2004  number_cache.n_active_quads = 0;
2005 
2006  // for 2d, quads have levels so take count the objects per
2007  // level and globally
2008  if (dim == 2)
2009  {
2010  // count the number of levels; the function we called above
2011  // on a separate Task for lines also does this and puts it into
2012  // number_cache.n_levels, but this datum may not yet be
2013  // available as we call the function on a separate task
2014  unsigned int n_levels = 0;
2015  if (level_objects > 0)
2016  // find the last level on which there are used cells
2017  for (unsigned int level = 0; level < level_objects; ++level)
2018  if (triangulation.begin(level) != triangulation.end(level))
2019  n_levels = level + 1;
2020 
2021  number_cache.n_quads_level.resize(n_levels);
2022  number_cache.n_active_quads_level.resize(n_levels);
2023 
2024  for (unsigned int level = 0; level < n_levels; ++level)
2025  {
2026  // count quads on this level
2027  number_cache.n_quads_level[level] = 0;
2028  number_cache.n_active_quads_level[level] = 0;
2029 
2030  quad_iterator quad = triangulation.begin_quad(level),
2031  endc =
2032  (level == n_levels - 1 ?
2033  quad_iterator(triangulation.end_quad()) :
2034  triangulation.begin_quad(level + 1));
2035  for (; quad != endc; ++quad)
2036  {
2037  ++number_cache.n_quads_level[level];
2038  if (quad->has_children() == false)
2039  ++number_cache.n_active_quads_level[level];
2040  }
2041 
2042  // update total number of quads
2043  number_cache.n_quads += number_cache.n_quads_level[level];
2044  number_cache.n_active_quads +=
2045  number_cache.n_active_quads_level[level];
2046  }
2047  }
2048  else
2049  {
2050  // for dim>2, there are no levels for quads
2051  number_cache.n_quads_level.clear();
2052  number_cache.n_active_quads_level.clear();
2053 
2054  quad_iterator quad = triangulation.begin_quad(),
2055  endc = triangulation.end_quad();
2056  for (; quad != endc; ++quad)
2057  {
2058  ++number_cache.n_quads;
2059  if (quad->has_children() == false)
2060  ++number_cache.n_active_quads;
2061  }
2062  }
2063 
2064  // wait for the background computation for lines
2065  update_lines.join();
2066  }
2067 
2083  template <int dim, int spacedim>
2084  static void
2087  const unsigned int level_objects,
2089  {
2090  // update quads, lines and n_levels in number_cache. since we
2091  // don't access any of these numbers, we can do this in the
2092  // background
2093  Threads::Task<void> update_quads_and_lines = Threads::new_task(
2094  static_cast<
2095  void (*)(const Triangulation<dim, spacedim> &,
2096  const unsigned int,
2098  &compute_number_cache<dim, spacedim>),
2099  triangulation,
2100  level_objects,
2102  number_cache));
2103 
2104  using hex_iterator =
2106 
2108  // update the number of hexes on the different levels in the
2109  // cache
2110  number_cache.n_hexes = 0;
2111  number_cache.n_active_hexes = 0;
2112 
2113  // for 3d, hexes have levels so take count the objects per
2114  // level and globally
2115  if (dim == 3)
2116  {
2117  // count the number of levels; the function we called
2118  // above on a separate Task for quads (recursively, via
2119  // the lines function) also does this and puts it into
2120  // number_cache.n_levels, but this datum may not yet be
2121  // available as we call the function on a separate task
2122  unsigned int n_levels = 0;
2123  if (level_objects > 0)
2124  // find the last level on which there are used cells
2125  for (unsigned int level = 0; level < level_objects; ++level)
2126  if (triangulation.begin(level) != triangulation.end(level))
2127  n_levels = level + 1;
2128 
2129  number_cache.n_hexes_level.resize(n_levels);
2130  number_cache.n_active_hexes_level.resize(n_levels);
2131 
2132  for (unsigned int level = 0; level < n_levels; ++level)
2133  {
2134  // count hexes on this level
2135  number_cache.n_hexes_level[level] = 0;
2136  number_cache.n_active_hexes_level[level] = 0;
2137 
2138  hex_iterator hex = triangulation.begin_hex(level),
2139  endc = (level == n_levels - 1 ?
2140  hex_iterator(triangulation.end_hex()) :
2141  triangulation.begin_hex(level + 1));
2142  for (; hex != endc; ++hex)
2143  {
2144  ++number_cache.n_hexes_level[level];
2145  if (hex->has_children() == false)
2146  ++number_cache.n_active_hexes_level[level];
2147  }
2148 
2149  // update total number of hexes
2150  number_cache.n_hexes += number_cache.n_hexes_level[level];
2151  number_cache.n_active_hexes +=
2152  number_cache.n_active_hexes_level[level];
2153  }
2154  }
2155  else
2156  {
2157  // for dim>3, there are no levels for hexes
2158  number_cache.n_hexes_level.clear();
2159  number_cache.n_active_hexes_level.clear();
2160 
2161  hex_iterator hex = triangulation.begin_hex(),
2162  endc = triangulation.end_hex();
2163  for (; hex != endc; ++hex)
2164  {
2165  ++number_cache.n_hexes;
2166  if (hex->has_children() == false)
2167  ++number_cache.n_active_hexes;
2168  }
2169  }
2170 
2171  // wait for the background computation for quads
2172  update_quads_and_lines.join();
2173  }
2174 
2175 
2179  template <int dim, int spacedim>
2180  static void
2182  const std::vector<CellData<dim>> & cells,
2183  const SubCellData & subcelldata,
2185  {
2186  AssertThrow(vertices.size() > 0, ExcMessage("No vertices given"));
2187  AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
2188 
2189  // Check that all cells have positive volume.
2190 #ifndef _MSC_VER
2191  // TODO: The following code does not compile with MSVC. Find a way
2192  // around it
2193  if (dim == spacedim)
2194  for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2195  {
2196  // If we should check for distorted cells, then we permit them
2197  // to exist. If a cell has negative measure, then it must be
2198  // distorted (the converse is not necessarily true); hence
2199  // throw an exception if no such cells should exist.
2200  if (tria.check_for_distorted_cells)
2201  {
2202  const double cell_measure =
2203  GridTools::cell_measure<dim>(vertices,
2205  cells[cell_no].vertices));
2206  AssertThrow(cell_measure > 0, ExcGridHasInvalidCell(cell_no));
2207  }
2208  }
2209 #endif
2210 
2211  // clear old content
2212  tria.levels.clear();
2213  tria.levels.push_back(
2214  std::make_unique<
2216 
2217  if (dim > 1)
2218  tria.faces = std::make_unique<
2220 
2221  // copy vertices
2222  tria.vertices = vertices;
2223  tria.vertices_used.assign(vertices.size(), true);
2224 
2225  // compute connectivity
2226  const auto connectivity = build_connectivity<unsigned int>(cells);
2227  const unsigned int n_cell = cells.size();
2228 
2229  // TriaObjects: lines
2230  if (dim >= 2)
2231  {
2232  auto &lines_0 = tria.faces->lines; // data structure to be filled
2233 
2234  // get connectivity between quads and lines
2235  const auto & crs = connectivity.entity_to_entities(1, 0);
2236  const unsigned int n_lines = crs.ptr.size() - 1;
2237 
2238  // allocate memory
2239  reserve_space_(lines_0, n_lines);
2240 
2241  // loop over lines
2242  for (unsigned int line = 0; line < n_lines; ++line)
2243  for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
2244  ++i, ++j)
2245  lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
2246  crs.col[i]; // set vertex indices
2247  }
2248 
2249  // TriaObjects: quads
2250  if (dim == 3)
2251  {
2252  auto &quads_0 = tria.faces->quads; // data structures to be filled
2253  auto &faces = *tria.faces;
2254 
2255  // get connectivity between quads and lines
2256  const auto & crs = connectivity.entity_to_entities(2, 1);
2257  const unsigned int n_quads = crs.ptr.size() - 1;
2258 
2259  // allocate memory
2260  reserve_space_(quads_0, n_quads);
2261  reserve_space_(faces, 2 /*structdim*/, n_quads);
2262 
2263  // loop over all quads -> entity type, line indices/orientations
2264  for (unsigned int q = 0, k = 0; q < n_quads; ++q)
2265  {
2266  // set entity type of quads
2267  faces.quad_reference_cell_type[q] =
2268  connectivity.entity_types(2)[q];
2269 
2270  // loop over all its lines
2271  for (unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
2272  ++i, ++j, ++k)
2273  {
2274  // set line index
2275  quads_0.cells[q * GeometryInfo<2>::faces_per_cell + j] =
2276  crs.col[i];
2277 
2278  // set line orientations
2279  faces.quads_line_orientations
2281  connectivity.entity_orientations(1)[k];
2282  }
2283  }
2284  }
2285 
2286  // TriaObjects/TriaLevel: cell
2287  {
2288  auto &cells_0 = tria.levels[0]->cells; // data structure to be filled
2289  auto &level = *tria.levels[0];
2290 
2291  // get connectivity between cells/faces and cells/cells
2292  const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2293  const auto &nei = connectivity.entity_to_entities(dim, dim);
2294 
2295  // in 2D optional: since in in pure QUAD meshes same line
2296  // orientations can be guaranteed
2297  const bool orientation_needed =
2298  dim == 3 ||
2299  (dim == 2 &&
2300  std::any_of(connectivity.entity_orientations(1).begin(),
2301  connectivity.entity_orientations(1).end(),
2302  [](const auto &i) { return i == 0; }));
2303 
2304  // allocate memory
2305  reserve_space_(cells_0, n_cell);
2306  reserve_space_(level, spacedim, n_cell, orientation_needed);
2307 
2308  // loop over all cells
2309  for (unsigned int cell = 0; cell < n_cell; ++cell)
2310  {
2311  // set material ids
2312  cells_0.boundary_or_material_id[cell].material_id =
2313  cells[cell].material_id;
2314 
2315  // set manifold ids
2316  cells_0.manifold_id[cell] = cells[cell].manifold_id;
2317 
2318  // set entity types
2319  level.reference_cell_type[cell] =
2320  connectivity.entity_types(dim)[cell];
2321 
2322  // loop over faces
2323  for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2324  ++i, ++j)
2325  {
2326  // set neighbor if not at boundary
2327  if (nei.col[i] != static_cast<unsigned int>(-1))
2328  level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
2329  j] = {0, nei.col[i]};
2330 
2331  // set face indices
2332  cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
2333  crs.col[i];
2334 
2335  // set face orientation if needed
2336  if (orientation_needed)
2337  level.face_orientations
2338  [cell * GeometryInfo<dim>::faces_per_cell + j] =
2339  connectivity.entity_orientations(dim - 1)[i];
2340  }
2341  }
2342  }
2343 
2344  // TriaFaces: boundary id of boundary faces
2345  if (dim > 1)
2346  {
2347  auto &bids_face = dim == 3 ?
2348  tria.faces->quads.boundary_or_material_id :
2349  tria.faces->lines.boundary_or_material_id;
2350 
2351  // count number of cells a face is belonging to
2352  std::vector<unsigned int> count(bids_face.size(), 0);
2353 
2354  // get connectivity between cells/faces
2355  const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2356 
2357  // count how many cells are adjacent to the same face
2358  for (unsigned int cell = 0; cell < cells.size(); ++cell)
2359  for (unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
2360  count[crs.col[i]]++;
2361 
2362  // loop over all faces
2363  for (unsigned int face = 0; face < count.size(); ++face)
2364  {
2365  if (count[face] != 1) // inner face
2366  continue;
2367 
2368  // boundary faces ...
2369  bids_face[face].boundary_id = 0;
2370 
2371  if (dim != 3)
2372  continue;
2373 
2374  // ... and the lines of quads in 3D
2375  const auto &crs = connectivity.entity_to_entities(2, 1);
2376  for (unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
2377  tria.faces->lines.boundary_or_material_id[crs.col[i]]
2378  .boundary_id = 0;
2379  }
2380  }
2381  else // 1D
2382  {
2383  static const unsigned int t_tba = static_cast<unsigned int>(-1);
2384  static const unsigned int t_inner = static_cast<unsigned int>(-2);
2385 
2386  std::vector<unsigned int> type(vertices.size(), t_tba);
2387 
2388  const auto &crs = connectivity.entity_to_entities(1, 0);
2389 
2390  for (unsigned int cell = 0; cell < cells.size(); ++cell)
2391  for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2392  ++i, ++j)
2393  if (type[crs.col[i]] != t_inner)
2394  type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
2395 
2396  for (unsigned int face = 0; face < type.size(); ++face)
2397  {
2398  // note: we also treat manifolds here!?
2399  (*tria.vertex_to_manifold_id_map_1d)[face] =
2401  if (type[face] != t_inner && type[face] != t_tba)
2402  (*tria.vertex_to_boundary_id_map_1d)[face] = type[face];
2403  }
2404  }
2405 
2406  // SubCellData: line
2407  if (dim >= 2)
2408  process_subcelldata(connectivity.entity_to_entities(1, 0),
2409  tria.faces->lines,
2410  subcelldata.boundary_lines);
2411 
2412  // SubCellData: quad
2413  if (dim == 3)
2414  process_subcelldata(connectivity.entity_to_entities(2, 0),
2415  tria.faces->quads,
2416  subcelldata.boundary_quads);
2417  }
2418 
2419 
2420  template <int structdim, typename T>
2421  static void
2423  const CRS<T> & crs,
2424  TriaObjects & obj,
2425  const std::vector<CellData<structdim>> &boundary_objects_in)
2426  {
2427  AssertDimension(obj.structdim, structdim);
2428 
2429  if (boundary_objects_in.size() == 0)
2430  return; // empty subcelldata -> nothing to do
2431 
2432  // pre-sort subcelldata
2433  auto boundary_objects = boundary_objects_in;
2434 
2435  // ... sort vertices
2436  for (auto &boundary_object : boundary_objects)
2437  std::sort(boundary_object.vertices.begin(),
2438  boundary_object.vertices.end());
2439 
2440  // ... sort cells
2441  std::sort(boundary_objects.begin(),
2442  boundary_objects.end(),
2443  [](const auto &a, const auto &b) {
2444  return a.vertices < b.vertices;
2445  });
2446 
2447  unsigned int counter = 0;
2448 
2449  std::vector<unsigned int> key;
2451 
2452  for (unsigned int o = 0; o < obj.n_objects(); ++o)
2453  {
2454  auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
2455  auto &manifold_id = obj.manifold_id[o];
2456 
2457  // assert that object has not been visited yet and its value
2458  // has not been modified yet
2459  AssertThrow(boundary_id == 0 ||
2461  ExcNotImplemented());
2463  ExcNotImplemented());
2464 
2465  // create key
2466  key.assign(crs.col.data() + crs.ptr[o],
2467  crs.col.data() + crs.ptr[o + 1]);
2468  std::sort(key.begin(), key.end());
2469 
2470  // is subcelldata provided? -> binary search
2471  const auto subcell_object =
2472  std::lower_bound(boundary_objects.begin(),
2473  boundary_objects.end(),
2474  key,
2475  [&](const auto &cell, const auto &key) {
2476  return cell.vertices < key;
2477  });
2478 
2479  // no subcelldata provided for this object
2480  if (subcell_object == boundary_objects.end() ||
2481  subcell_object->vertices != key)
2482  continue;
2483 
2484  counter++;
2485 
2486  // set manifold id
2487  manifold_id = subcell_object->manifold_id;
2488 
2489  // set boundary id
2490  if (subcell_object->boundary_id !=
2492  {
2494  ExcNotImplemented());
2495  boundary_id = subcell_object->boundary_id;
2496  }
2497  }
2498 
2499  // make sure that all subcelldata entries have been processed
2500  // TODO: this is not guaranteed, why?
2501  // AssertDimension(counter, boundary_objects_in.size());
2502  }
2503 
2504 
2505 
2506  static void
2508  const unsigned structdim,
2509  const unsigned int size)
2510  {
2511  const unsigned int dim = faces.dim;
2512 
2513  const unsigned int faces_per_cell =
2514  structdim == 1 ? GeometryInfo<1>::faces_per_cell :
2515  (structdim == 2 ? GeometryInfo<2>::faces_per_cell :
2517 
2518  if (dim == 3 && structdim == 2)
2519  {
2520  // quad entity types
2521  faces.quad_reference_cell_type.assign(size,
2523 
2524  // quad line orientations
2525  faces.quads_line_orientations.assign(size * faces_per_cell, -1);
2526  }
2527  }
2528 
2529 
2530 
2531  static void
2533  const unsigned int spacedim,
2534  const unsigned int size,
2535  const bool orientation_needed)
2536  {
2537  const unsigned int dim = level.dim;
2538 
2539  const unsigned int faces_per_cell =
2540  dim == 1 ? GeometryInfo<1>::faces_per_cell :
2541  (dim == 2 ? GeometryInfo<2>::faces_per_cell :
2543 
2544  level.active_cell_indices.assign(size, -1);
2545  level.subdomain_ids.assign(size, 0);
2546  level.level_subdomain_ids.assign(size, 0);
2547 
2548  level.refine_flags.assign(size, false);
2549  level.coarsen_flags.assign(size, false);
2550 
2551  level.parents.assign((size + 1) / 2, -1);
2552 
2553  if (dim < spacedim)
2554  level.direction_flags.assign(size, true);
2555 
2556  level.neighbors.assign(size * faces_per_cell, {-1, -1});
2557 
2559 
2560  if (orientation_needed)
2561  level.face_orientations.assign(size * faces_per_cell, -1);
2562  }
2563 
2564 
2565 
2566  static void
2567  reserve_space_(TriaObjects &obj, const unsigned int size)
2568  {
2569  const unsigned int structdim = obj.structdim;
2570 
2571  const unsigned int max_children_per_cell =
2572  structdim == 1 ?
2574  (structdim == 2 ? GeometryInfo<2>::max_children_per_cell :
2576  const unsigned int faces_per_cell =
2577  structdim == 1 ? GeometryInfo<1>::faces_per_cell :
2578  (structdim == 2 ? GeometryInfo<2>::faces_per_cell :
2580 
2581  obj.used.assign(size, true);
2582  obj.boundary_or_material_id.assign(
2583  size,
2585  BoundaryOrMaterialId());
2586  obj.manifold_id.assign(size, -1);
2587  obj.user_flags.assign(size, false);
2588  obj.user_data.resize(size);
2589 
2590  if (structdim > 1) // TODO: why?
2591  obj.refinement_cases.assign(size, 0);
2592 
2593  obj.children.assign(max_children_per_cell / 2 * size, -1);
2594 
2595  obj.cells.assign(faces_per_cell * size, -1);
2596 
2597  if (structdim <= 2)
2598  {
2599  obj.next_free_single = size - 1;
2600  obj.next_free_pair = 0;
2601  obj.reverse_order_next_free_single = true;
2602  }
2603  else
2604  {
2605  obj.next_free_single = obj.next_free_pair = 0;
2606  }
2607  }
2608 
2609 
2625  template <int spacedim>
2626  static void delete_children(
2629  std::vector<unsigned int> &,
2630  std::vector<unsigned int> &)
2631  {
2632  const unsigned int dim = 1;
2633 
2634  // first we need to reset the
2635  // neighbor pointers of the
2636  // neighbors of this cell's
2637  // children to this cell. This is
2638  // different for one dimension,
2639  // since there neighbors can have a
2640  // refinement level differing from
2641  // that of this cell's children by
2642  // more than one level.
2643 
2644  Assert(!cell->child(0)->has_children() &&
2645  !cell->child(1)->has_children(),
2646  ExcInternalError());
2647 
2648  // first do it for the cells to the
2649  // left
2650  if (cell->neighbor(0).state() == IteratorState::valid)
2651  if (cell->neighbor(0)->has_children())
2652  {
2654  cell->neighbor(0);
2655  Assert(neighbor->level() == cell->level(), ExcInternalError());
2656 
2657  // right child
2658  neighbor = neighbor->child(1);
2659  while (true)
2660  {
2661  Assert(neighbor->neighbor(1) == cell->child(0),
2662  ExcInternalError());
2663  neighbor->set_neighbor(1, cell);
2664 
2665  // move on to further
2666  // children on the
2667  // boundary between this
2668  // cell and its neighbor
2669  if (neighbor->has_children())
2670  neighbor = neighbor->child(1);
2671  else
2672  break;
2673  }
2674  }
2675 
2676  // now do it for the cells to the
2677  // left
2678  if (cell->neighbor(1).state() == IteratorState::valid)
2679  if (cell->neighbor(1)->has_children())
2680  {
2682  cell->neighbor(1);
2683  Assert(neighbor->level() == cell->level(), ExcInternalError());
2684 
2685  // left child
2686  neighbor = neighbor->child(0);
2687  while (true)
2688  {
2689  Assert(neighbor->neighbor(0) == cell->child(1),
2690  ExcInternalError());
2691  neighbor->set_neighbor(0, cell);
2692 
2693  // move on to further
2694  // children on the
2695  // boundary between this
2696  // cell and its neighbor
2697  if (neighbor->has_children())
2698  neighbor = neighbor->child(0);
2699  else
2700  break;
2701  }
2702  }
2703 
2704 
2705  // delete the vertex which will not
2706  // be needed anymore. This vertex
2707  // is the second of the first child
2708  triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
2709 
2710  // invalidate children. clear user
2711  // pointers, to avoid that they may
2712  // appear at unwanted places later
2713  // on...
2714  for (unsigned int child = 0; child < cell->n_children(); ++child)
2715  {
2716  cell->child(child)->clear_user_data();
2717  cell->child(child)->clear_user_flag();
2718  cell->child(child)->clear_used_flag();
2719  }
2720 
2721 
2722  // delete pointer to children
2723  cell->clear_children();
2724  cell->clear_user_flag();
2725  }
2726 
2727 
2728 
2729  template <int spacedim>
2730  static void delete_children(
2733  std::vector<unsigned int> & line_cell_count,
2734  std::vector<unsigned int> &)
2735  {
2736  const unsigned int dim = 2;
2737  const RefinementCase<dim> ref_case = cell->refinement_case();
2738 
2739  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2740  ExcInternalError());
2741 
2742  // vectors to hold all lines which
2743  // may be deleted
2744  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2745  lines_to_delete(0);
2746 
2747  lines_to_delete.reserve(4 * 2 + 4);
2748 
2749  // now we decrease the counters for
2750  // lines contained in the child
2751  // cells
2752  for (unsigned int c = 0; c < cell->n_children(); ++c)
2753  {
2755  cell->child(c);
2756  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2757  --line_cell_count[child->line_index(l)];
2758  }
2759 
2760 
2761  // delete the vertex which will not
2762  // be needed anymore. This vertex
2763  // is the second of the second line
2764  // of the first child, if the cell
2765  // is refined with cut_xy, else there
2766  // is no inner vertex.
2767  // additionally delete unneeded inner
2768  // lines
2769  if (ref_case == RefinementCase<dim>::cut_xy)
2770  {
2771  triangulation
2772  .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
2773 
2774  lines_to_delete.push_back(cell->child(0)->line(1));
2775  lines_to_delete.push_back(cell->child(0)->line(3));
2776  lines_to_delete.push_back(cell->child(3)->line(0));
2777  lines_to_delete.push_back(cell->child(3)->line(2));
2778  }
2779  else
2780  {
2781  unsigned int inner_face_no =
2782  ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
2783 
2784  // the inner line will not be
2785  // used any more
2786  lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2787  }
2788 
2789  // invalidate children
2790  for (unsigned int child = 0; child < cell->n_children(); ++child)
2791  {
2792  cell->child(child)->clear_user_data();
2793  cell->child(child)->clear_user_flag();
2794  cell->child(child)->clear_used_flag();
2795  }
2796 
2797 
2798  // delete pointer to children
2799  cell->clear_children();
2800  cell->clear_refinement_case();
2801  cell->clear_user_flag();
2802 
2803  // look at the refinement of outer
2804  // lines. if nobody needs those
2805  // anymore we can add them to the
2806  // list of lines to be deleted.
2807  for (unsigned int line_no = 0;
2808  line_no < GeometryInfo<dim>::lines_per_cell;
2809  ++line_no)
2810  {
2812  cell->line(line_no);
2813 
2814  if (line->has_children())
2815  {
2816  // if one of the cell counters is
2817  // zero, the other has to be as well
2818 
2819  Assert((line_cell_count[line->child_index(0)] == 0 &&
2820  line_cell_count[line->child_index(1)] == 0) ||
2821  (line_cell_count[line->child_index(0)] > 0 &&
2822  line_cell_count[line->child_index(1)] > 0),
2823  ExcInternalError());
2824 
2825  if (line_cell_count[line->child_index(0)] == 0)
2826  {
2827  for (unsigned int c = 0; c < 2; ++c)
2828  Assert(!line->child(c)->has_children(),
2829  ExcInternalError());
2830 
2831  // we may delete the line's
2832  // children and the middle vertex
2833  // as no cell references them
2834  // anymore
2835  triangulation
2836  .vertices_used[line->child(0)->vertex_index(1)] = false;
2837 
2838  lines_to_delete.push_back(line->child(0));
2839  lines_to_delete.push_back(line->child(1));
2840 
2841  line->clear_children();
2842  }
2843  }
2844  }
2845 
2846  // finally, delete unneeded lines
2847 
2848  // clear user pointers, to avoid that
2849  // they may appear at unwanted places
2850  // later on...
2851  // same for user flags, then finally
2852  // delete the lines
2853  typename std::vector<
2855  line = lines_to_delete.begin(),
2856  endline = lines_to_delete.end();
2857  for (; line != endline; ++line)
2858  {
2859  (*line)->clear_user_data();
2860  (*line)->clear_user_flag();
2861  (*line)->clear_used_flag();
2862  }
2863  }
2864 
2865 
2866 
2867  template <int spacedim>
2868  static void delete_children(
2871  std::vector<unsigned int> & line_cell_count,
2872  std::vector<unsigned int> & quad_cell_count)
2873  {
2874  const unsigned int dim = 3;
2875 
2876  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2877  ExcInternalError());
2878  Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
2879  ExcInternalError());
2880 
2881  // first of all, we store the RefineCase of
2882  // this cell
2883  const RefinementCase<dim> ref_case = cell->refinement_case();
2884  // vectors to hold all lines and quads which
2885  // may be deleted
2886  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2887  lines_to_delete(0);
2888  std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
2889  quads_to_delete(0);
2890 
2891  lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
2892  quads_to_delete.reserve(6 * 4 + 12);
2893 
2894  // now we decrease the counters for lines and
2895  // quads contained in the child cells
2896  for (unsigned int c = 0; c < cell->n_children(); ++c)
2897  {
2899  cell->child(c);
2900  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2901  --line_cell_count[child->line_index(l)];
2902  for (auto f : GeometryInfo<dim>::face_indices())
2903  --quad_cell_count[child->quad_index(f)];
2904  }
2905 
2907  // delete interior quads and lines and the
2908  // interior vertex, depending on the
2909  // refinement case of the cell
2910  //
2911  // for append quads and lines: only append
2912  // them to the list of objects to be deleted
2913 
2914  switch (ref_case)
2915  {
2917  quads_to_delete.push_back(cell->child(0)->face(1));
2918  break;
2920  quads_to_delete.push_back(cell->child(0)->face(3));
2921  break;
2923  quads_to_delete.push_back(cell->child(0)->face(5));
2924  break;
2926  quads_to_delete.push_back(cell->child(0)->face(1));
2927  quads_to_delete.push_back(cell->child(0)->face(3));
2928  quads_to_delete.push_back(cell->child(3)->face(0));
2929  quads_to_delete.push_back(cell->child(3)->face(2));
2930 
2931  lines_to_delete.push_back(cell->child(0)->line(11));
2932  break;
2934  quads_to_delete.push_back(cell->child(0)->face(1));
2935  quads_to_delete.push_back(cell->child(0)->face(5));
2936  quads_to_delete.push_back(cell->child(3)->face(0));
2937  quads_to_delete.push_back(cell->child(3)->face(4));
2938 
2939  lines_to_delete.push_back(cell->child(0)->line(5));
2940  break;
2942  quads_to_delete.push_back(cell->child(0)->face(3));
2943  quads_to_delete.push_back(cell->child(0)->face(5));
2944  quads_to_delete.push_back(cell->child(3)->face(2));
2945  quads_to_delete.push_back(cell->child(3)->face(4));
2946 
2947  lines_to_delete.push_back(cell->child(0)->line(7));
2948  break;
2950  quads_to_delete.push_back(cell->child(0)->face(1));
2951  quads_to_delete.push_back(cell->child(2)->face(1));
2952  quads_to_delete.push_back(cell->child(4)->face(1));
2953  quads_to_delete.push_back(cell->child(6)->face(1));
2954 
2955  quads_to_delete.push_back(cell->child(0)->face(3));
2956  quads_to_delete.push_back(cell->child(1)->face(3));
2957  quads_to_delete.push_back(cell->child(4)->face(3));
2958  quads_to_delete.push_back(cell->child(5)->face(3));
2959 
2960  quads_to_delete.push_back(cell->child(0)->face(5));
2961  quads_to_delete.push_back(cell->child(1)->face(5));
2962  quads_to_delete.push_back(cell->child(2)->face(5));
2963  quads_to_delete.push_back(cell->child(3)->face(5));
2964 
2965  lines_to_delete.push_back(cell->child(0)->line(5));
2966  lines_to_delete.push_back(cell->child(0)->line(7));
2967  lines_to_delete.push_back(cell->child(0)->line(11));
2968  lines_to_delete.push_back(cell->child(7)->line(0));
2969  lines_to_delete.push_back(cell->child(7)->line(2));
2970  lines_to_delete.push_back(cell->child(7)->line(8));
2971  // delete the vertex which will not
2972  // be needed anymore. This vertex
2973  // is the vertex at the heart of
2974  // this cell, which is the sixth of
2975  // the first child
2976  triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
2977  false;
2978  break;
2979  default:
2980  // only remaining case is
2981  // no_refinement, thus an error
2982  Assert(false, ExcInternalError());
2983  break;
2984  }
2985 
2986 
2987  // invalidate children
2988  for (unsigned int child = 0; child < cell->n_children(); ++child)
2989  {
2990  cell->child(child)->clear_user_data();
2991  cell->child(child)->clear_user_flag();
2992 
2993  for (auto f : GeometryInfo<dim>::face_indices())
2994  {
2995  // set flags denoting deviations from
2996  // standard orientation of faces back
2997  // to initialization values
2998  cell->child(child)->set_face_orientation(f, true);
2999  cell->child(child)->set_face_flip(f, false);
3000  cell->child(child)->set_face_rotation(f, false);
3001  }
3002 
3003  cell->child(child)->clear_used_flag();
3004  }
3005 
3006 
3007  // delete pointer to children
3008  cell->clear_children();
3009  cell->clear_refinement_case();
3010  cell->clear_user_flag();
3011 
3012  // so far we only looked at inner quads,
3013  // lines and vertices. Now we have to
3014  // consider outer ones as well. here, we have
3015  // to check, whether there are other cells
3016  // still needing these objects. otherwise we
3017  // can delete them. first for quads (and
3018  // their inner lines).
3019 
3020  for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
3021  {
3023  cell->face(quad_no);
3024 
3025  Assert(
3026  (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
3027  quad->has_children()) ||
3028  GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
3029  RefinementCase<dim - 1>::no_refinement,
3030  ExcInternalError());
3031 
3032  switch (quad->refinement_case())
3033  {
3034  case RefinementCase<dim - 1>::no_refinement:
3035  // nothing to do as the quad
3036  // is not refined
3037  break;
3038  case RefinementCase<dim - 1>::cut_x:
3039  case RefinementCase<dim - 1>::cut_y:
3040  {
3041  // if one of the cell counters is
3042  // zero, the other has to be as
3043  // well
3044  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3045  quad_cell_count[quad->child_index(1)] == 0) ||
3046  (quad_cell_count[quad->child_index(0)] > 0 &&
3047  quad_cell_count[quad->child_index(1)] > 0),
3048  ExcInternalError());
3049  // it might be, that the quad is
3050  // refined twice anisotropically,
3051  // first check, whether we may
3052  // delete possible grand_children
3053  unsigned int deleted_grandchildren = 0;
3054  unsigned int number_of_child_refinements = 0;
3055 
3056  for (unsigned int c = 0; c < 2; ++c)
3057  if (quad->child(c)->has_children())
3058  {
3059  ++number_of_child_refinements;
3060  // if one of the cell counters is
3061  // zero, the other has to be as
3062  // well
3063  Assert(
3064  (quad_cell_count[quad->child(c)->child_index(0)] ==
3065  0 &&
3066  quad_cell_count[quad->child(c)->child_index(1)] ==
3067  0) ||
3068  (quad_cell_count[quad->child(c)->child_index(0)] >
3069  0 &&
3070  quad_cell_count[quad->child(c)->child_index(1)] >
3071  0),
3072  ExcInternalError());
3073  if (quad_cell_count[quad->child(c)->child_index(0)] ==
3074  0)
3075  {
3076  // Assert, that the two
3077  // anisotropic
3078  // refinements add up to
3079  // isotropic refinement
3080  Assert(quad->refinement_case() +
3081  quad->child(c)->refinement_case() ==
3083  ExcInternalError());
3084  // we may delete the
3085  // quad's children and
3086  // the inner line as no
3087  // cell references them
3088  // anymore
3089  quads_to_delete.push_back(
3090  quad->child(c)->child(0));
3091  quads_to_delete.push_back(
3092  quad->child(c)->child(1));
3093  if (quad->child(c)->refinement_case() ==
3095  lines_to_delete.push_back(
3096  quad->child(c)->child(0)->line(1));
3097  else
3098  lines_to_delete.push_back(
3099  quad->child(c)->child(0)->line(3));
3100  quad->child(c)->clear_children();
3101  quad->child(c)->clear_refinement_case();
3102  ++deleted_grandchildren;
3103  }
3104  }
3105  // if no grandchildren are left, we
3106  // may as well delete the
3107  // refinement of the inner line
3108  // between our children and the
3109  // corresponding vertex
3110  if (number_of_child_refinements > 0 &&
3111  deleted_grandchildren == number_of_child_refinements)
3112  {
3114  middle_line;
3115  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3116  middle_line = quad->child(0)->line(1);
3117  else
3118  middle_line = quad->child(0)->line(3);
3119 
3120  lines_to_delete.push_back(middle_line->child(0));
3121  lines_to_delete.push_back(middle_line->child(1));
3122  triangulation
3123  .vertices_used[middle_vertex_index<dim, spacedim>(
3124  middle_line)] = false;
3125  middle_line->clear_children();
3126  }
3127 
3128  // now consider the direct children
3129  // of the given quad
3130  if (quad_cell_count[quad->child_index(0)] == 0)
3131  {
3132  // we may delete the quad's
3133  // children and the inner line
3134  // as no cell references them
3135  // anymore
3136  quads_to_delete.push_back(quad->child(0));
3137  quads_to_delete.push_back(quad->child(1));
3138  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3139  lines_to_delete.push_back(quad->child(0)->line(1));
3140  else
3141  lines_to_delete.push_back(quad->child(0)->line(3));
3142 
3143  // if the counters just dropped
3144  // to zero, otherwise the
3145  // children would have been
3146  // deleted earlier, then this
3147  // cell's children must have
3148  // contained the anisotropic
3149  // quad children. thus, if
3150  // those have again anisotropic
3151  // children, which are in
3152  // effect isotropic children of
3153  // the original quad, those are
3154  // still needed by a
3155  // neighboring cell and we
3156  // cannot delete them. instead,
3157  // we have to reset this quad's
3158  // refine case to isotropic and
3159  // set the children
3160  // accordingly.
3161  if (quad->child(0)->has_children())
3162  if (quad->refinement_case() ==
3164  {
3165  // now evereything is
3166  // quite complicated. we
3167  // have the children
3168  // numbered according to
3169  //
3170  // *---*---*
3171  // |n+1|m+1|
3172  // *---*---*
3173  // | n | m |
3174  // *---*---*
3175  //
3176  // from the original
3177  // anisotropic
3178  // refinement. we have to
3179  // reorder them as
3180  //
3181  // *---*---*
3182  // | m |m+1|
3183  // *---*---*
3184  // | n |n+1|
3185  // *---*---*
3186  //
3187  // for isotropic refinement.
3188  //
3189  // this is a bit ugly, of
3190  // course: loop over all
3191  // cells on all levels
3192  // and look for faces n+1
3193  // (switch_1) and m
3194  // (switch_2).
3195  const typename Triangulation<dim, spacedim>::
3196  quad_iterator switch_1 =
3197  quad->child(0)->child(1),
3198  switch_2 =
3199  quad->child(1)->child(0);
3200 
3201  Assert(!switch_1->has_children(),
3202  ExcInternalError());
3203  Assert(!switch_2->has_children(),
3204  ExcInternalError());
3205 
3206  const int switch_1_index = switch_1->index();
3207  const int switch_2_index = switch_2->index();
3208  for (unsigned int l = 0;
3209  l < triangulation.levels.size();
3210  ++l)
3211  for (unsigned int h = 0;
3212  h <
3213  triangulation.levels[l]->cells.n_objects();
3214  ++h)
3215  for (const unsigned int q :
3217  {
3218  const int index =
3219  triangulation.levels[l]
3220  ->cells.get_bounding_object_indices(
3221  h)[q];
3222  if (index == switch_1_index)
3223  triangulation.levels[l]
3224  ->cells.get_bounding_object_indices(
3225  h)[q] = switch_2_index;
3226  else if (index == switch_2_index)
3227  triangulation.levels[l]
3228  ->cells.get_bounding_object_indices(
3229  h)[q] = switch_1_index;
3230  }
3231  // now we have to copy
3232  // all information of the
3233  // two quads
3234  const int switch_1_lines[4] = {
3235  static_cast<signed int>(
3236  switch_1->line_index(0)),
3237  static_cast<signed int>(
3238  switch_1->line_index(1)),
3239  static_cast<signed int>(
3240  switch_1->line_index(2)),
3241  static_cast<signed int>(
3242  switch_1->line_index(3))};
3243  const bool switch_1_line_orientations[4] = {
3244  switch_1->line_orientation(0),
3245  switch_1->line_orientation(1),
3246  switch_1->line_orientation(2),
3247  switch_1->line_orientation(3)};
3248  const types::boundary_id switch_1_boundary_id =
3249  switch_1->boundary_id();
3250  const unsigned int switch_1_user_index =
3251  switch_1->user_index();
3252  const bool switch_1_user_flag =
3253  switch_1->user_flag_set();
3254 
3255  switch_1->set_bounding_object_indices(
3256  {switch_2->line_index(0),
3257  switch_2->line_index(1),
3258  switch_2->line_index(2),
3259  switch_2->line_index(3)});
3260  switch_1->set_line_orientation(
3261  0, switch_2->line_orientation(0));
3262  switch_1->set_line_orientation(
3263  1, switch_2->line_orientation(1));
3264  switch_1->set_line_orientation(
3265  2, switch_2->line_orientation(2));
3266  switch_1->set_line_orientation(
3267  3, switch_2->line_orientation(3));
3268  switch_1->set_boundary_id_internal(
3269  switch_2->boundary_id());
3270  switch_1->set_manifold_id(
3271  switch_2->manifold_id());
3272  switch_1->set_user_index(switch_2->user_index());
3273  if (switch_2->user_flag_set())
3274  switch_1->set_user_flag();
3275  else
3276  switch_1->clear_user_flag();
3277 
3278  switch_2->set_bounding_object_indices(
3279  {switch_1_lines[0],
3280  switch_1_lines[1],
3281  switch_1_lines[2],
3282  switch_1_lines[3]});
3283  switch_2->set_line_orientation(
3284  0, switch_1_line_orientations[0]);
3285  switch_2->set_line_orientation(
3286  1, switch_1_line_orientations[1]);
3287  switch_2->set_line_orientation(
3288  2, switch_1_line_orientations[2]);
3289  switch_2->set_line_orientation(
3290  3, switch_1_line_orientations[3]);
3291  switch_2->set_boundary_id_internal(
3292  switch_1_boundary_id);
3293  switch_2->set_manifold_id(
3294  switch_1->manifold_id());
3295  switch_2->set_user_index(switch_1_user_index);
3296  if (switch_1_user_flag)
3297  switch_2->set_user_flag();
3298  else
3299  switch_2->clear_user_flag();
3300 
3301  const unsigned int child_0 =
3302  quad->child(0)->child_index(0);
3303  const unsigned int child_2 =
3304  quad->child(1)->child_index(0);
3305  quad->clear_children();
3306  quad->clear_refinement_case();
3307  quad->set_refinement_case(
3309  quad->set_children(0, child_0);
3310  quad->set_children(2, child_2);
3311  std::swap(quad_cell_count[child_0 + 1],
3312  quad_cell_count[child_2]);
3313  }
3314  else
3315  {
3316  // the face was refined
3317  // with cut_y, thus the
3318  // children are already
3319  // in correct order. we
3320  // only have to set them
3321  // correctly, deleting
3322  // the indirection of two
3323  // anisotropic refinement
3324  // and going directly
3325  // from the quad to
3326  // isotropic children
3327  const unsigned int child_0 =
3328  quad->child(0)->child_index(0);
3329  const unsigned int child_2 =
3330  quad->child(1)->child_index(0);
3331  quad->clear_children();
3332  quad->clear_refinement_case();
3333  quad->set_refinement_case(
3335  quad->set_children(0, child_0);
3336  quad->set_children(2, child_2);
3337  }
3338  else
3339  {
3340  quad->clear_children();
3341  quad->clear_refinement_case();
3342  }
3343  }
3344  break;
3345  }
3346  case RefinementCase<dim - 1>::cut_xy:
3347  {
3348  // if one of the cell counters is
3349  // zero, the others have to be as
3350  // well
3351 
3352  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3353  quad_cell_count[quad->child_index(1)] == 0 &&
3354  quad_cell_count[quad->child_index(2)] == 0 &&
3355  quad_cell_count[quad->child_index(3)] == 0) ||
3356  (quad_cell_count[quad->child_index(0)] > 0 &&
3357  quad_cell_count[quad->child_index(1)] > 0 &&
3358  quad_cell_count[quad->child_index(2)] > 0 &&
3359  quad_cell_count[quad->child_index(3)] > 0),
3360  ExcInternalError());
3361 
3362  if (quad_cell_count[quad->child_index(0)] == 0)
3363  {
3364  // we may delete the quad's
3365  // children, the inner lines
3366  // and the middle vertex as no
3367  // cell references them anymore
3368  lines_to_delete.push_back(quad->child(0)->line(1));
3369  lines_to_delete.push_back(quad->child(3)->line(0));
3370  lines_to_delete.push_back(quad->child(0)->line(3));
3371  lines_to_delete.push_back(quad->child(3)->line(2));
3372 
3373  for (unsigned int child = 0; child < quad->n_children();
3374  ++child)
3375  quads_to_delete.push_back(quad->child(child));
3376 
3377  triangulation
3378  .vertices_used[quad->child(0)->vertex_index(3)] =
3379  false;
3380 
3381  quad->clear_children();
3382  quad->clear_refinement_case();
3383  }
3384  }
3385  break;
3386 
3387  default:
3388  Assert(false, ExcInternalError());
3389  break;
3390  }
3391  }
3392 
3393  // now we repeat a similar procedure
3394  // for the outer lines of this cell.
3395 
3396  // if in debug mode: check that each
3397  // of the lines for which we consider
3398  // deleting the children in fact has
3399  // children (the bits/coarsening_3d
3400  // test tripped over this initially)
3401  for (unsigned int line_no = 0;
3402  line_no < GeometryInfo<dim>::lines_per_cell;
3403  ++line_no)
3404  {
3406  cell->line(line_no);
3407 
3408  Assert(
3409  (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
3410  line->has_children()) ||
3411  GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
3413  ExcInternalError());
3414 
3415  if (line->has_children())
3416  {
3417  // if one of the cell counters is
3418  // zero, the other has to be as well
3419 
3420  Assert((line_cell_count[line->child_index(0)] == 0 &&
3421  line_cell_count[line->child_index(1)] == 0) ||
3422  (line_cell_count[line->child_index(0)] > 0 &&
3423  line_cell_count[line->child_index(1)] > 0),
3424  ExcInternalError());
3425 
3426  if (line_cell_count[line->child_index(0)] == 0)
3427  {
3428  for (unsigned int c = 0; c < 2; ++c)
3429  Assert(!line->child(c)->has_children(),
3430  ExcInternalError());
3431 
3432  // we may delete the line's
3433  // children and the middle vertex
3434  // as no cell references them
3435  // anymore
3436  triangulation
3437  .vertices_used[line->child(0)->vertex_index(1)] = false;
3438 
3439  lines_to_delete.push_back(line->child(0));
3440  lines_to_delete.push_back(line->child(1));
3441 
3442  line->clear_children();
3443  }
3444  }
3445  }
3446 
3447  // finally, delete unneeded quads and lines
3448 
3449  // clear user pointers, to avoid that
3450  // they may appear at unwanted places
3451  // later on...
3452  // same for user flags, then finally
3453  // delete the quads and lines
3454  typename std::vector<
3456  line = lines_to_delete.begin(),
3457  endline = lines_to_delete.end();
3458  for (; line != endline; ++line)
3459  {
3460  (*line)->clear_user_data();
3461  (*line)->clear_user_flag();
3462  (*line)->clear_used_flag();
3463  }
3464 
3465  typename std::vector<
3467  quad = quads_to_delete.begin(),
3468  endquad = quads_to_delete.end();
3469  for (; quad != endquad; ++quad)
3470  {
3471  (*quad)->clear_user_data();
3472  (*quad)->clear_children();
3473  (*quad)->clear_refinement_case();
3474  (*quad)->clear_user_flag();
3475  (*quad)->clear_used_flag();
3476  }
3477  }
3478 
3479 
3497  template <int spacedim>
3498  static void create_children(
3500  unsigned int & next_unused_vertex,
3502  &next_unused_line,
3504  &next_unused_cell,
3505  const typename Triangulation<2, spacedim>::cell_iterator &cell)
3506  {
3507  const unsigned int dim = 2;
3508  // clear refinement flag
3509  const RefinementCase<dim> ref_case = cell->refine_flag_set();
3510  cell->clear_refine_flag();
3511 
3512  /* For the refinement process: since we go the levels up from the
3513  lowest, there are (unlike above) only two possibilities: a neighbor
3514  cell is on the same level or one level up (in both cases, it may or
3515  may not be refined later on, but we don't care here).
3516 
3517  First:
3518  Set up an array of the 3x3 vertices, which are distributed on the
3519  cell (the array consists of indices into the @p{vertices} std::vector
3520 
3521  2--7--3
3522  | | |
3523  4--8--5
3524  | | |
3525  0--6--1
3526 
3527  note: in case of cut_x or cut_y not all these vertices are needed for
3528  the new cells
3529 
3530  Second:
3531  Set up an array of the new lines (the array consists of iterator
3532  pointers into the lines arrays)
3533 
3534  .-6-.-7-. The directions are: .->-.->-.
3535  1 9 3 ^ ^ ^
3536  .-10.11-. .->-.->-.
3537  0 8 2 ^ ^ ^
3538  .-4-.-5-. .->-.->-.
3539 
3540  cut_x:
3541  .-4-.-5-.
3542  | | |
3543  0 6 1
3544  | | |
3545  .-2-.-3-.
3546 
3547  cut_y:
3548  .---5---.
3549  1 3
3550  .---6---.
3551  0 2
3552  .---4---.
3553 
3554 
3555  Third:
3556  Set up an array of neighbors:
3557 
3558  6 7
3559  .--.--.
3560  1| | |3
3561  .--.--.
3562  0| | |2
3563  .--.--.
3564  4 5
3565 
3566  We need this array for two reasons: first to get the lines which will
3567  bound the four subcells (if the neighboring cell is refined, these
3568  lines already exist), and second to update neighborship information.
3569  Since if a neighbor is not refined, its neighborship record only
3570  points to the present, unrefined, cell rather than the children we
3571  are presently creating, we only need the neighborship information
3572  if the neighbor cells are refined. In all other cases, we store
3573  the unrefined neighbor address
3574 
3575  We also need for every neighbor (if refined) which number among its
3576  neighbors the present (unrefined) cell has, since that number is to
3577  be replaced and because that also is the number of the subline which
3578  will be the interface between that neighbor and the to be created
3579  cell. We will store this number (between 0 and 3) in the field
3580  @p{neighbors_neighbor}.
3581 
3582  It would be sufficient to use the children of the common line to the
3583  neighbor, if we only wanted to get the new sublines and the new
3584  vertex, but because we need to update the neighborship information of
3585  the two refined subcells of the neighbor, we need to search these
3586  anyway.
3587 
3588  Convention:
3589  The created children are numbered like this:
3590 
3591  .--.--.
3592  |2 . 3|
3593  .--.--.
3594  |0 | 1|
3595  .--.--.
3596  */
3597  // collect the
3598  // indices of the
3599  // eight
3600  // surrounding
3601  // vertices
3602  // 2--7--3
3603  // | | |
3604  // 4--9--5
3605  // | | |
3606  // 0--6--1
3607  int new_vertices[9];
3608  for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3609  new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3610  for (unsigned int line_no = 0; line_no < 4; ++line_no)
3611  if (cell->line(line_no)->has_children())
3612  new_vertices[4 + line_no] =
3613  cell->line(line_no)->child(0)->vertex_index(1);
3614 
3615  if (ref_case == RefinementCase<dim>::cut_xy)
3616  {
3617  // find the next
3618  // unused vertex and
3619  // allocate it for
3620  // the new vertex we
3621  // need here
3622  while (triangulation.vertices_used[next_unused_vertex] == true)
3623  ++next_unused_vertex;
3624  Assert(
3625  next_unused_vertex < triangulation.vertices.size(),
3626  ExcMessage(
3627  "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."));
3628  triangulation.vertices_used[next_unused_vertex] = true;
3629 
3630  new_vertices[8] = next_unused_vertex;
3631 
3632  // if this quad lives
3633  // in 2d, then we can
3634  // compute the new
3635  // central vertex
3636  // location just from
3637  // the surrounding
3638  // ones. If this is
3639  // not the case, then
3640  // we need to ask a
3641  // boundary object
3642  if (dim == spacedim)
3643  {
3644  // triangulation.vertices[next_unused_vertex] = new_point;
3645  triangulation.vertices[next_unused_vertex] = cell->center(true);
3646 
3647  // if the user_flag is set, i.e. if the cell is at the
3648  // boundary, use a different calculation of the middle vertex
3649  // here. this is of advantage if the boundary is strongly
3650  // curved (whereas the cell is not) and the cell has a high
3651  // aspect ratio.
3652  if (cell->user_flag_set())
3653  {
3654  // first reset the user_flag and then refine
3655  cell->clear_user_flag();
3656  triangulation.vertices[next_unused_vertex] =
3657  cell->center(true, true);
3658  }
3659  }
3660  else
3661  {
3662  // if this quad lives in a higher dimensional space
3663  // then we don't need to worry if it is at the
3664  // boundary of the manifold -- we always have to use
3665  // the boundary object anyway; so ignore whether the
3666  // user flag is set or not
3667  cell->clear_user_flag();
3668 
3669  // determine middle vertex by transfinite interpolation to be
3670  // consistent with what happens to quads in a Triangulation<3,
3671  // 3> when they are refined
3672  triangulation.vertices[next_unused_vertex] =
3673  cell->center(true, true);
3674  }
3675  }
3676 
3677 
3678  // Now the lines:
3679  typename Triangulation<dim, spacedim>::raw_line_iterator new_lines[12];
3680  unsigned int lmin = 8;
3681  unsigned int lmax = 12;
3682  if (ref_case != RefinementCase<dim>::cut_xy)
3683  {
3684  lmin = 6;
3685  lmax = 7;
3686  }
3687 
3688  for (unsigned int l = lmin; l < lmax; ++l)
3689  {
3690  while (next_unused_line->used() == true)
3691  ++next_unused_line;
3692  new_lines[l] = next_unused_line;
3693  ++next_unused_line;
3694 
3695  Assert(
3696  new_lines[l]->used() == false,
3697  ExcMessage(
3698  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3699  }
3700 
3701  if (ref_case == RefinementCase<dim>::cut_xy)
3702  {
3703  // .-6-.-7-.
3704  // 1 9 3
3705  // .-10.11-.
3706  // 0 8 2
3707  // .-4-.-5-.
3708 
3709  // lines 0-7 already exist, create only the four interior
3710  // lines 8-11
3711  unsigned int l = 0;
3712  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
3713  for (unsigned int c = 0; c < 2; ++c, ++l)
3714  new_lines[l] = cell->line(face_no)->child(c);
3715  Assert(l == 8, ExcInternalError());
3716 
3717  new_lines[8]->set_bounding_object_indices(
3718  {new_vertices[6], new_vertices[8]});
3719  new_lines[9]->set_bounding_object_indices(
3720  {new_vertices[8], new_vertices[7]});
3721  new_lines[10]->set_bounding_object_indices(
3722  {new_vertices[4], new_vertices[8]});
3723  new_lines[11]->set_bounding_object_indices(
3724  {new_vertices[8], new_vertices[5]});
3725  }
3726  else if (ref_case == RefinementCase<dim>::cut_x)
3727  {
3728  // .-4-.-5-.
3729  // | | |
3730  // 0 6 1
3731  // | | |
3732  // .-2-.-3-.
3733  new_lines[0] = cell->line(0);
3734  new_lines[1] = cell->line(1);
3735  new_lines[2] = cell->line(2)->child(0);
3736  new_lines[3] = cell->line(2)->child(1);
3737  new_lines[4] = cell->line(3)->child(0);
3738  new_lines[5] = cell->line(3)->child(1);
3739  new_lines[6]->set_bounding_object_indices(
3740  {new_vertices[6], new_vertices[7]});
3741  }
3742  else
3743  {
3745  // .---5---.
3746  // 1 3
3747  // .---6---.
3748  // 0 2
3749  // .---4---.
3750  new_lines[0] = cell->line(0)->child(0);
3751  new_lines[1] = cell->line(0)->child(1);
3752  new_lines[2] = cell->line(1)->child(0);
3753  new_lines[3] = cell->line(1)->child(1);
3754  new_lines[4] = cell->line(2);
3755  new_lines[5] = cell->line(3);
3756  new_lines[6]->set_bounding_object_indices(
3757  {new_vertices[4], new_vertices[5]});
3758  }
3759 
3760  for (unsigned int l = lmin; l < lmax; ++l)
3761  {
3762  new_lines[l]->set_used_flag();
3763  new_lines[l]->clear_user_flag();
3764  new_lines[l]->clear_user_data();
3765  new_lines[l]->clear_children();
3766  // interior line
3767  new_lines[l]->set_boundary_id_internal(
3769  new_lines[l]->set_manifold_id(cell->manifold_id());
3770  }
3771 
3772  // Now add the four (two)
3773  // new cells!
3776  while (next_unused_cell->used() == true)
3777  ++next_unused_cell;
3778 
3779  const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
3780  for (unsigned int i = 0; i < n_children; ++i)
3781  {
3782  Assert(
3783  next_unused_cell->used() == false,
3784  ExcMessage(
3785  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3786  subcells[i] = next_unused_cell;
3787  ++next_unused_cell;
3788  if (i % 2 == 1 && i < n_children - 1)
3789  while (next_unused_cell->used() == true)
3790  ++next_unused_cell;
3791  }
3792 
3793  if (ref_case == RefinementCase<dim>::cut_xy)
3794  {
3795  // children:
3796  // .--.--.
3797  // |2 . 3|
3798  // .--.--.
3799  // |0 | 1|
3800  // .--.--.
3801  // lines:
3802  // .-6-.-7-.
3803  // 1 9 3
3804  // .-10.11-.
3805  // 0 8 2
3806  // .-4-.-5-.
3807  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3808  new_lines[8]->index(),
3809  new_lines[4]->index(),
3810  new_lines[10]->index()});
3811  subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3812  new_lines[2]->index(),
3813  new_lines[5]->index(),
3814  new_lines[11]->index()});
3815  subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3816  new_lines[9]->index(),
3817  new_lines[10]->index(),
3818  new_lines[6]->index()});
3819  subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3820  new_lines[3]->index(),
3821  new_lines[11]->index(),
3822  new_lines[7]->index()});
3823  }
3824  else if (ref_case == RefinementCase<dim>::cut_x)
3825  {
3826  // children:
3827  // .--.--.
3828  // | . |
3829  // .0 . 1.
3830  // | | |
3831  // .--.--.
3832  // lines:
3833  // .-4-.-5-.
3834  // | | |
3835  // 0 6 1
3836  // | | |
3837  // .-2-.-3-.
3838  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3839  new_lines[6]->index(),
3840  new_lines[2]->index(),
3841  new_lines[4]->index()});
3842  subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3843  new_lines[1]->index(),
3844  new_lines[3]->index(),
3845  new_lines[5]->index()});
3846  }
3847  else
3848  {
3850  // children:
3851  // .-----.
3852  // | 1 |
3853  // .-----.
3854  // | 0 |
3855  // .-----.
3856  // lines:
3857  // .---5---.
3858  // 1 3
3859  // .---6---.
3860  // 0 2
3861  // .---4---.
3862  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3863  new_lines[2]->index(),
3864  new_lines[4]->index(),
3865  new_lines[6]->index()});
3866  subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
3867  new_lines[3]->index(),
3868  new_lines[6]->index(),
3869  new_lines[5]->index()});
3870  }
3871 
3872  types::subdomain_id subdomainid = cell->subdomain_id();
3873 
3874  for (unsigned int i = 0; i < n_children; ++i)
3875  {
3876  subcells[i]->set_used_flag();
3877  subcells[i]->clear_refine_flag();
3878  subcells[i]->clear_user_flag();
3879  subcells[i]->clear_user_data();
3880  subcells[i]->clear_children();
3881  // inherit material
3882  // properties
3883  subcells[i]->set_material_id(cell->material_id());
3884  subcells[i]->set_manifold_id(cell->manifold_id());
3885  subcells[i]->set_subdomain_id(subdomainid);
3886 
3887  if (i % 2 == 0)
3888  subcells[i]->set_parent(cell->index());
3889  }
3890 
3891 
3892 
3893  // set child index for
3894  // even children children
3895  // i=0,2 (0)
3896  for (unsigned int i = 0; i < n_children / 2; ++i)
3897  cell->set_children(2 * i, subcells[2 * i]->index());
3898  // set the refine case
3899  cell->set_refinement_case(ref_case);
3900 
3901  // note that the
3902  // refinement flag was
3903  // already cleared at the
3904  // beginning of this function
3905 
3906  if (dim < spacedim)
3907  for (unsigned int c = 0; c < n_children; ++c)
3908  cell->child(c)->set_direction_flag(cell->direction_flag());
3909  }
3910 
3911 
3912 
3917  template <int spacedim>
3920  const bool /*check_for_distorted_cells*/)
3921  {
3922  const unsigned int dim = 1;
3923 
3924  // check whether a new level is needed we have to check for
3925  // this on the highest level only (on this, all used cells are
3926  // also active, so we only have to check for this)
3927  {
3929  cell = triangulation.begin_active(triangulation.levels.size() - 1),
3930  endc = triangulation.end();
3931  for (; cell != endc; ++cell)
3932  if (cell->used())
3933  if (cell->refine_flag_set())
3934  {
3935  triangulation.levels.push_back(
3936  std::make_unique<
3938  break;
3939  }
3940  }
3941 
3942 
3943  // check how much space is needed on every level we need not
3944  // check the highest level since either - on the highest level
3945  // no cells are flagged for refinement - there are, but
3946  // prepare_refinement added another empty level
3947  unsigned int needed_vertices = 0;
3948  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
3949  {
3950  // count number of flagged
3951  // cells on this level
3952  unsigned int flagged_cells = 0;
3953 
3954  for (const auto &acell :
3955  triangulation.active_cell_iterators_on_level(level))
3956  if (acell->refine_flag_set())
3957  ++flagged_cells;
3958 
3959  // count number of used cells
3960  // on the next higher level
3961  const unsigned int used_cells =
3962  std::count(triangulation.levels[level + 1]->cells.used.begin(),
3963  triangulation.levels[level + 1]->cells.used.end(),
3964  true);
3965 
3966  // reserve space for the used_cells cells already existing
3967  // on the next higher level as well as for the
3968  // 2*flagged_cells that will be created on that level
3969  reserve_space(*triangulation.levels[level + 1],
3971  flagged_cells,
3972  1,
3973  spacedim);
3974  // reserve space for 2*flagged_cells new lines on the next
3975  // higher level
3976  reserve_space(triangulation.levels[level + 1]->cells,
3978  flagged_cells,
3979  0);
3980 
3981  needed_vertices += flagged_cells;
3982  }
3983 
3984  // add to needed vertices how many
3985  // vertices are already in use
3986  needed_vertices += std::count(triangulation.vertices_used.begin(),
3987  triangulation.vertices_used.end(),
3988  true);
3989  // if we need more vertices: create them, if not: leave the
3990  // array as is, since shrinking is not really possible because
3991  // some of the vertices at the end may be in use
3992  if (needed_vertices > triangulation.vertices.size())
3993  {
3994  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
3995  triangulation.vertices_used.resize(needed_vertices, false);
3996  }
3997 
3998 
3999  // Do REFINEMENT on every level; exclude highest level as
4000  // above
4001 
4002  // index of next unused vertex
4003  unsigned int next_unused_vertex = 0;
4004 
4005  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4006  {
4008  next_unused_cell = triangulation.begin_raw(level + 1);
4009 
4010  for (const auto &cell :
4011  triangulation.active_cell_iterators_on_level(level))
4012  if (cell->refine_flag_set())
4013  {
4014  // clear refinement flag
4015  cell->clear_refine_flag();
4016 
4017  // search for next unused
4018  // vertex
4019  while (triangulation.vertices_used[next_unused_vertex] ==
4020  true)
4021  ++next_unused_vertex;
4022  Assert(
4023  next_unused_vertex < triangulation.vertices.size(),
4024  ExcMessage(
4025  "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."));
4026 
4027  // Now we always ask the cell itself where to put
4028  // the new point. The cell in turn will query the
4029  // manifold object internally.
4030  triangulation.vertices[next_unused_vertex] =
4031  cell->center(true);
4032 
4033  triangulation.vertices_used[next_unused_vertex] = true;
4034 
4035  // search for next two unused cell (++ takes care of
4036  // the end of the vector)
4038  first_child,
4039  second_child;
4040  while (next_unused_cell->used() == true)
4041  ++next_unused_cell;
4042  first_child = next_unused_cell;
4043  first_child->set_used_flag();
4044  first_child->clear_user_data();
4045  ++next_unused_cell;
4046  Assert(
4047  next_unused_cell->used() == false,
4048  ExcMessage(
4049  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4050  second_child = next_unused_cell;
4051  second_child->set_used_flag();
4052  second_child->clear_user_data();
4053 
4054  types::subdomain_id subdomainid = cell->subdomain_id();
4055 
4056  // insert first child
4057  cell->set_children(0, first_child->index());
4058  first_child->clear_children();
4059  first_child->set_bounding_object_indices(
4060  {cell->vertex_index(0), next_unused_vertex});
4061  first_child->set_material_id(cell->material_id());
4062  first_child->set_manifold_id(cell->manifold_id());
4063  first_child->set_subdomain_id(subdomainid);
4064  first_child->set_direction_flag(cell->direction_flag());
4065 
4066  first_child->set_parent(cell->index());
4067 
4068  // Set manifold id of the right face. Only do this
4069  // on the first child.
4070  first_child->face(1)->set_manifold_id(cell->manifold_id());
4071 
4072  // reset neighborship info (refer to
4073  // internal::TriangulationImplementation::TriaLevel<0> for
4074  // details)
4075  first_child->set_neighbor(1, second_child);
4076  if (cell->neighbor(0).state() != IteratorState::valid)
4077  first_child->set_neighbor(0, cell->neighbor(0));
4078  else if (cell->neighbor(0)->is_active())
4079  {
4080  // since the neighbors level is always <=level,
4081  // if the cell is active, then there are no
4082  // cells to the left which may want to know
4083  // about this new child cell.
4084  Assert(cell->neighbor(0)->level() <= cell->level(),
4085  ExcInternalError());
4086  first_child->set_neighbor(0, cell->neighbor(0));
4087  }
4088  else
4089  // left neighbor is refined
4090  {
4091  // set neighbor to cell on same level
4092  const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4093  first_child->set_neighbor(0,
4094  cell->neighbor(0)->child(nbnb));
4095 
4096  // reset neighbor info of all right descendant
4097  // of the left neighbor of cell
4099  left_neighbor = cell->neighbor(0);
4100  while (left_neighbor->has_children())
4101  {
4102  left_neighbor = left_neighbor->child(nbnb);
4103  left_neighbor->set_neighbor(nbnb, first_child);
4104  }
4105  }
4106 
4107  // insert second child
4108  second_child->clear_children();
4109  second_child->set_bounding_object_indices(
4110  {next_unused_vertex, cell->vertex_index(1)});
4111  second_child->set_neighbor(0, first_child);
4112  second_child->set_material_id(cell->material_id());
4113  second_child->set_manifold_id(cell->manifold_id());
4114  second_child->set_subdomain_id(subdomainid);
4115  second_child->set_direction_flag(cell->direction_flag());
4116 
4117  if (cell->neighbor(1).state() != IteratorState::valid)
4118  second_child->set_neighbor(1, cell->neighbor(1));
4119  else if (cell->neighbor(1)->is_active())
4120  {
4121  Assert(cell->neighbor(1)->level() <= cell->level(),
4122  ExcInternalError());
4123  second_child->set_neighbor(1, cell->neighbor(1));
4124  }
4125  else
4126  // right neighbor is refined same as above
4127  {
4128  const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4129  second_child->set_neighbor(
4130  1, cell->neighbor(1)->child(nbnb));
4131 
4133  right_neighbor = cell->neighbor(1);
4134  while (right_neighbor->has_children())
4135  {
4136  right_neighbor = right_neighbor->child(nbnb);
4137  right_neighbor->set_neighbor(nbnb, second_child);
4138  }
4139  }
4140  // inform all listeners that cell refinement is done
4141  triangulation.signals.post_refinement_on_cell(cell);
4142  }
4143  }
4144 
4145  // in 1d, we can not have distorted children unless the parent
4146  // was already distorted (that is because we don't use
4147  // boundary information for 1d triangulations). so return an
4148  // empty list
4150  }
4151 
4152 
4157  template <int spacedim>
4160  const bool check_for_distorted_cells)
4161  {
4162  const unsigned int dim = 2;
4163 
4164  // check whether a new level is needed we have to check for
4165  // this on the highest level only (on this, all used cells are
4166  // also active, so we only have to check for this)
4167  {
4169  cell = triangulation.begin_active(triangulation.levels.size() - 1),
4170  endc = triangulation.end();
4171  for (; cell != endc; ++cell)
4172  if (cell->used())
4173  if (cell->refine_flag_set())
4174  {
4175  triangulation.levels.push_back(
4176  std::make_unique<
4178  break;
4179  }
4180  }
4181 
4182  // TODO[WB]: we clear user flags and pointers of lines; we're going
4183  // to use them to flag which lines need refinement
4184  for (typename Triangulation<dim, spacedim>::line_iterator line =
4185  triangulation.begin_line();
4186  line != triangulation.end_line();
4187  ++line)
4188  {
4189  line->clear_user_flag();
4190  line->clear_user_data();
4191  }
4192  // running over all cells and lines count the number
4193  // n_single_lines of lines which can be stored as single
4194  // lines, e.g. inner lines
4195  unsigned int n_single_lines = 0;
4196 
4197  // New lines to be created: number lines which are stored in
4198  // pairs (the children of lines must be stored in pairs)
4199  unsigned int n_lines_in_pairs = 0;
4200 
4201  // check how much space is needed on every level we need not
4202  // check the highest level since either - on the highest level
4203  // no cells are flagged for refinement - there are, but
4204  // prepare_refinement added another empty level
4205  unsigned int needed_vertices = 0;
4206  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4207  {
4208  // count number of flagged cells on this level and compute
4209  // how many new vertices and new lines will be needed
4210  unsigned int needed_cells = 0;
4211 
4212  for (const auto &cell :
4213  triangulation.active_cell_iterators_on_level(level))
4214  if (cell->refine_flag_set())
4215  {
4216  if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
4217  {
4218  needed_cells += 4;
4219 
4220  // new vertex at center of cell is needed in any
4221  // case
4222  ++needed_vertices;
4223 
4224  // the four inner lines can be stored as singles
4225  n_single_lines += 4;
4226  }
4227  else // cut_x || cut_y
4228  {
4229  // set the flag showing that anisotropic
4230  // refinement is used for at least one cell
4231  triangulation.anisotropic_refinement = true;
4232 
4233  needed_cells += 2;
4234  // no vertex at center
4235 
4236  // the inner line can be stored as single
4237  n_single_lines += 1;
4238  }
4239 
4240  // mark all faces (lines) for refinement; checking
4241  // locally whether the neighbor would also like to
4242  // refine them is rather difficult for lines so we
4243  // only flag them and after visiting all cells, we
4244  // decide which lines need refinement;
4245  for (const unsigned int line_no :
4247  {
4249  cell->refine_flag_set(), line_no) ==
4251  {
4253  line = cell->line(line_no);
4254  if (line->has_children() == false)
4255  line->set_user_flag();
4256  }
4257  }
4258  }
4259 
4260 
4261  // count number of used cells on the next higher level
4262  const unsigned int used_cells =
4263  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4264  triangulation.levels[level + 1]->cells.used.end(),
4265  true);
4266 
4267 
4268  // reserve space for the used_cells cells already existing
4269  // on the next higher level as well as for the
4270  // needed_cells that will be created on that level
4271  reserve_space(*triangulation.levels[level + 1],
4272  used_cells + needed_cells,
4273  2,
4274  spacedim);
4275 
4276  // reserve space for needed_cells new quads on the next
4277  // higher level
4278  reserve_space(triangulation.levels[level + 1]->cells,
4279  needed_cells,
4280  0);
4281  }
4282 
4283  // now count the lines which were flagged for refinement
4284  for (typename Triangulation<dim, spacedim>::line_iterator line =
4285  triangulation.begin_line();
4286  line != triangulation.end_line();
4287  ++line)
4288  if (line->user_flag_set())
4289  {
4290  Assert(line->has_children() == false, ExcInternalError());
4291  n_lines_in_pairs += 2;
4292  needed_vertices += 1;
4293  }
4294  // reserve space for n_lines_in_pairs new lines. note, that
4295  // we can't reserve space for the single lines here as well,
4296  // as all the space reserved for lines in pairs would be
4297  // counted as unused and we would end up with too little space
4298  // to store all lines. memory reservation for n_single_lines
4299  // can only be done AFTER we refined the lines of the current
4300  // cells
4301  reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4302 
4303  // add to needed vertices how many vertices are already in use
4304  needed_vertices += std::count(triangulation.vertices_used.begin(),
4305  triangulation.vertices_used.end(),
4306  true);
4307  // if we need more vertices: create them, if not: leave the
4308  // array as is, since shrinking is not really possible because
4309  // some of the vertices at the end may be in use
4310  if (needed_vertices > triangulation.vertices.size())
4311  {
4312  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4313  triangulation.vertices_used.resize(needed_vertices, false);
4314  }
4315 
4316 
4317  // Do REFINEMENT on every level; exclude highest level as
4318  // above
4319 
4320  // index of next unused vertex
4321  unsigned int next_unused_vertex = 0;
4322 
4323  // first the refinement of lines. children are stored
4324  // pairwise
4325  {
4326  // only active objects can be refined further
4328  line = triangulation.begin_active_line(),
4329  endl = triangulation.end_line();
4331  next_unused_line = triangulation.begin_raw_line();
4332 
4333  for (; line != endl; ++line)
4334  if (line->user_flag_set())
4335  {
4336  // this line needs to be refined
4337 
4338  // find the next unused vertex and set it
4339  // appropriately
4340  while (triangulation.vertices_used[next_unused_vertex] == true)
4341  ++next_unused_vertex;
4342  Assert(
4343  next_unused_vertex < triangulation.vertices.size(),
4344  ExcMessage(
4345  "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."));
4346  triangulation.vertices_used[next_unused_vertex] = true;
4347 
4348  triangulation.vertices[next_unused_vertex] = line->center(true);
4349 
4350  // now that we created the right point, make up the
4351  // two child lines. To this end, find a pair of
4352  // unused lines
4353  bool pair_found = false;
4354  (void)pair_found;
4355  for (; next_unused_line != endl; ++next_unused_line)
4356  if (!next_unused_line->used() &&
4357  !(++next_unused_line)->used())
4358  {
4359  // go back to the first of the two unused
4360  // lines
4361  --next_unused_line;
4362  pair_found = true;
4363  break;
4364  }
4365  Assert(pair_found, ExcInternalError());
4366 
4367  // there are now two consecutive unused lines, such
4368  // that the children of a line will be consecutive.
4369  // then set the child pointer of the present line
4370  line->set_children(0, next_unused_line->index());
4371 
4372  // set the two new lines
4374  children[2] = {next_unused_line, ++next_unused_line};
4375  // some tests; if any of the iterators should be
4376  // invalid, then already dereferencing will fail
4377  Assert(
4378  children[0]->used() == false,
4379  ExcMessage(
4380  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4381  Assert(
4382  children[1]->used() == false,
4383  ExcMessage(
4384  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4385 
4386  children[0]->set_bounding_object_indices(
4387  {line->vertex_index(0), next_unused_vertex});
4388  children[1]->set_bounding_object_indices(
4389  {next_unused_vertex, line->vertex_index(1)});
4390 
4391  children[0]->set_used_flag();
4392  children[1]->set_used_flag();
4393  children[0]->clear_children();
4394  children[1]->clear_children();
4395  children[0]->clear_user_data();
4396  children[1]->clear_user_data();
4397  children[0]->clear_user_flag();
4398  children[1]->clear_user_flag();
4399 
4400 
4401  children[0]->set_boundary_id_internal(line->boundary_id());
4402  children[1]->set_boundary_id_internal(line->boundary_id());
4403 
4404  children[0]->set_manifold_id(line->manifold_id());
4405  children[1]->set_manifold_id(line->manifold_id());
4406 
4407  // finally clear flag indicating the need for
4408  // refinement
4409  line->clear_user_flag();
4410  }
4411  }
4412 
4413 
4414  // Now set up the new cells
4415 
4416  // reserve space for inner lines (can be stored as single
4417  // lines)
4418  reserve_space(triangulation.faces->lines, 0, n_single_lines);
4419 
4421  cells_with_distorted_children;
4422 
4423  // reset next_unused_line, as now also single empty places in
4424  // the vector can be used
4426  next_unused_line = triangulation.begin_raw_line();
4427 
4428  for (int level = 0;
4429  level < static_cast<int>(triangulation.levels.size()) - 1;
4430  ++level)
4431  {
4433  next_unused_cell = triangulation.begin_raw(level + 1);
4434 
4435  for (const auto &cell :
4436  triangulation.active_cell_iterators_on_level(level))
4437  if (cell->refine_flag_set())
4438  {
4439  // set the user flag to indicate, that at least one
4440  // line is at the boundary
4441 
4442  // TODO[Tobias Leicht] find a better place to set
4443  // this flag, so that we do not need so much time to
4444  // check each cell here
4445  if (cell->at_boundary())
4446  cell->set_user_flag();
4447 
4448  // actually set up the children and update neighbor
4449  // information
4450  create_children(triangulation,
4451  next_unused_vertex,
4452  next_unused_line,
4453  next_unused_cell,
4454  cell);
4455 
4456  if (check_for_distorted_cells &&
4457  has_distorted_children(
4458  cell,
4459  std::integral_constant<int, dim>(),
4460  std::integral_constant<int, spacedim>()))
4461  cells_with_distorted_children.distorted_cells.push_back(
4462  cell);
4463  // inform all listeners that cell refinement is done
4464  triangulation.signals.post_refinement_on_cell(cell);
4465  }
4466  }
4467 
4468  return cells_with_distorted_children;
4469  }
4470 
4471 
4476  template <int spacedim>
4479  const bool check_for_distorted_cells)
4480  {
4481  const unsigned int dim = 3;
4482 
4483  // this function probably also works for spacedim>3 but it
4484  // isn't tested. it will probably be necessary to pull new
4485  // vertices onto the manifold just as we do for the other
4486  // functions above.
4487  Assert(spacedim == 3, ExcNotImplemented());
4488 
4489  // check whether a new level is needed we have to check for
4490  // this on the highest level only (on this, all used cells are
4491  // also active, so we only have to check for this)
4492  {
4494  cell = triangulation.begin_active(triangulation.levels.size() - 1),
4495  endc = triangulation.end();
4496  for (; cell != endc; ++cell)
4497  if (cell->used())
4498  if (cell->refine_flag_set())
4499  {
4500  triangulation.levels.push_back(
4501  std::make_unique<
4503  break;
4504  }
4505  }
4506 
4507 
4508  // first clear user flags for quads and lines; we're going to
4509  // use them to flag which lines and quads need refinement
4510  triangulation.faces->quads.clear_user_data();
4511 
4512  for (typename Triangulation<dim, spacedim>::line_iterator line =
4513  triangulation.begin_line();
4514  line != triangulation.end_line();
4515  ++line)
4516  line->clear_user_flag();
4517  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
4518  triangulation.begin_quad();
4519  quad != triangulation.end_quad();
4520  ++quad)
4521  quad->clear_user_flag();
4522 
4523  // create an array of face refine cases. User indices of faces
4524  // will be set to values corresponding with indices in this
4525  // array.
4526  const RefinementCase<dim - 1> face_refinement_cases[4] = {
4527  RefinementCase<dim - 1>::no_refinement,
4528  RefinementCase<dim - 1>::cut_x,
4529  RefinementCase<dim - 1>::cut_y,
4530  RefinementCase<dim - 1>::cut_xy};
4531 
4532  // check how much space is needed on every level we need not
4533  // check the highest level since either
4534  // - on the highest level no cells are flagged for refinement
4535  // - there are, but prepare_refinement added another empty
4536  // level which then is the highest level
4537 
4538  // variables to hold the number of newly to be created
4539  // vertices, lines and quads. as these are stored globally,
4540  // declare them outside the loop over al levels. we need lines
4541  // and quads in pairs for refinement of old ones and lines and
4542  // quads, that can be stored as single ones, as they are newly
4543  // created in the inside of an existing cell
4544  unsigned int needed_vertices = 0;
4545  unsigned int needed_lines_single = 0;
4546  unsigned int needed_quads_single = 0;
4547  unsigned int needed_lines_pair = 0;
4548  unsigned int needed_quads_pair = 0;
4549  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4550  {
4551  // count number of flagged cells on this level and compute
4552  // how many new vertices and new lines will be needed
4553  unsigned int new_cells = 0;
4554 
4555  for (const auto &acell :
4556  triangulation.active_cell_iterators_on_level(level))
4557  if (acell->refine_flag_set())
4558  {
4559  RefinementCase<dim> ref_case = acell->refine_flag_set();
4560 
4561  // now for interior vertices, lines and quads, which
4562  // are needed in any case
4563  if (ref_case == RefinementCase<dim>::cut_x ||
4564  ref_case == RefinementCase<dim>::cut_y ||
4565  ref_case == RefinementCase<dim>::cut_z)
4566  {
4567  ++needed_quads_single;
4568  new_cells += 2;
4569  triangulation.anisotropic_refinement = true;
4570  }
4571  else if (ref_case == RefinementCase<dim>::cut_xy ||
4572  ref_case == RefinementCase<dim>::cut_xz ||
4573  ref_case == RefinementCase<dim>::cut_yz)
4574  {
4575  ++needed_lines_single;
4576  needed_quads_single += 4;
4577  new_cells += 4;
4578  triangulation.anisotropic_refinement = true;
4579  }
4580  else if (ref_case == RefinementCase<dim>::cut_xyz)
4581  {
4582  ++needed_vertices;
4583  needed_lines_single += 6;
4584  needed_quads_single += 12;
4585  new_cells += 8;
4586  }
4587  else
4588  {
4589  // we should never get here
4590  Assert(false, ExcInternalError());
4591  }
4592 
4593  // mark all faces for refinement; checking locally
4594  // if and how the neighbor would like to refine
4595  // these is difficult so we only flag them and after
4596  // visiting all cells, we decide which faces need
4597  // which refinement;
4598  for (const unsigned int face :
4600  {
4602  aface = acell->face(face);
4603  // get the RefineCase this faces has for the
4604  // given RefineCase of the cell
4605  RefinementCase<dim - 1> face_ref_case =
4607  ref_case,
4608  face,
4609  acell->face_orientation(face),
4610  acell->face_flip(face),
4611  acell->face_rotation(face));
4612  // only do something, if this face has to be
4613  // refined
4614  if (face_ref_case)
4615  {
4616  if (face_ref_case ==
4618  {
4619  if (aface->number_of_children() < 4)
4620  // we use user_flags to denote needed
4621  // isotropic refinement
4622  aface->set_user_flag();
4623  }
4624  else if (aface->refinement_case() != face_ref_case)
4625  // we use user_indices to denote needed
4626  // anisotropic refinement. note, that we
4627  // can have at most one anisotropic
4628  // refinement case for this face, as
4629  // otherwise prepare_refinement() would
4630  // have changed one of the cells to yield
4631  // isotropic refinement at this
4632  // face. therefore we set the user_index
4633  // uniquely
4634  {
4635  Assert(aface->refinement_case() ==
4637  dim - 1>::isotropic_refinement ||
4638  aface->refinement_case() ==
4639  RefinementCase<dim - 1>::no_refinement,
4640  ExcInternalError());
4641  aface->set_user_index(face_ref_case);
4642  }
4643  }
4644  } // for all faces
4645 
4646  // flag all lines, that have to be refined
4647  for (unsigned int line = 0;
4648  line < GeometryInfo<dim>::lines_per_cell;
4649  ++line)
4651  line) &&
4652  !acell->line(line)->has_children())
4653  acell->line(line)->set_user_flag();
4654 
4655  } // if refine_flag set and for all cells on this level
4656 
4657 
4658  // count number of used cells on the next higher level
4659  const unsigned int used_cells =
4660  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4661  triangulation.levels[level + 1]->cells.used.end(),
4662  true);
4663 
4664 
4665  // reserve space for the used_cells cells already existing
4666  // on the next higher level as well as for the
4667  // 8*flagged_cells that will be created on that level
4668  reserve_space(*triangulation.levels[level + 1],
4669  used_cells + new_cells,
4670  3,
4671  spacedim);
4672  // reserve space for 8*flagged_cells new hexes on the next
4673  // higher level
4674  reserve_space(triangulation.levels[level + 1]->cells, new_cells);
4675  } // for all levels
4676  // now count the quads and lines which were flagged for
4677  // refinement
4678  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
4679  triangulation.begin_quad();
4680  quad != triangulation.end_quad();
4681  ++quad)
4682  {
4683  if (quad->user_flag_set())
4684  {
4685  // isotropic refinement: 1 interior vertex, 4 quads
4686  // and 4 interior lines. we store the interior lines
4687  // in pairs in case the face is already or will be
4688  // refined anisotropically
4689  needed_quads_pair += 4;
4690  needed_lines_pair += 4;
4691  needed_vertices += 1;
4692  }
4693  if (quad->user_index())
4694  {
4695  // anisotropic refinement: 1 interior
4696  // line and two quads
4697  needed_quads_pair += 2;
4698  needed_lines_single += 1;
4699  // there is a kind of complicated situation here which
4700  // requires our attention. if the quad is refined
4701  // isotropcally, two of the interior lines will get a
4702  // new mother line - the interior line of our
4703  // anisotropically refined quad. if those two lines
4704  // are not consecutive, we cannot do so and have to
4705  // replace them by two lines that are consecutive. we
4706  // try to avoid that situation, but it may happen
4707  // nevertheless through repeated refinement and
4708  // coarsening. thus we have to check here, as we will
4709  // need some additional space to store those new lines
4710  // in case we need them...
4711  if (quad->has_children())
4712  {
4713  Assert(quad->refinement_case() ==
4714  RefinementCase<dim - 1>::isotropic_refinement,
4715  ExcInternalError());
4716  if ((face_refinement_cases[quad->user_index()] ==
4717  RefinementCase<dim - 1>::cut_x &&
4718  (quad->child(0)->line_index(1) + 1 !=
4719  quad->child(2)->line_index(1))) ||
4720  (face_refinement_cases[quad->user_index()] ==
4721  RefinementCase<dim - 1>::cut_y &&
4722  (quad->child(0)->line_index(3) + 1 !=
4723  quad->child(1)->line_index(3))))
4724  needed_lines_pair += 2;
4725  }
4726  }
4727  }
4728 
4729  for (typename Triangulation<dim, spacedim>::line_iterator line =
4730  triangulation.begin_line();
4731  line != triangulation.end_line();
4732  ++line)
4733  if (line->user_flag_set())
4734  {
4735  needed_lines_pair += 2;
4736  needed_vertices += 1;
4737  }
4738 
4739  // reserve space for needed_lines new lines stored in pairs
4740  reserve_space(triangulation.faces->lines,
4741  needed_lines_pair,
4742  needed_lines_single);
4743  // reserve space for needed_quads new quads stored in pairs
4744  reserve_space(*triangulation.faces,
4745  needed_quads_pair,
4746  needed_quads_single);
4747  reserve_space(triangulation.faces->quads,
4748  needed_quads_pair,
4749  needed_quads_single);
4750 
4751 
4752  // add to needed vertices how many vertices are already in use
4753  needed_vertices += std::count(triangulation.vertices_used.begin(),
4754  triangulation.vertices_used.end(),
4755  true);
4756  // if we need more vertices: create them, if not: leave the
4757  // array as is, since shrinking is not really possible because
4758  // some of the vertices at the end may be in use
4759  if (needed_vertices > triangulation.vertices.size())
4760  {
4761  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4762  triangulation.vertices_used.resize(needed_vertices, false);
4763  }
4764 
4765 
4767  // Before we start with the actual refinement, we do some
4768  // sanity checks if in debug mode. especially, we try to catch
4769  // the notorious problem with lines being twice refined,
4770  // i.e. there are cells adjacent at one line ("around the
4771  // edge", but not at a face), with two cells differing by more
4772  // than one refinement level
4773  //
4774  // this check is very simple to implement here, since we have
4775  // all lines flagged if they shall be refined
4776 #ifdef DEBUG
4777  for (const auto &cell : triangulation.active_cell_iterators())
4778  if (!cell->refine_flag_set())
4779  for (unsigned int line = 0;
4780  line < GeometryInfo<dim>::lines_per_cell;
4781  ++line)
4782  if (cell->line(line)->has_children())
4783  for (unsigned int c = 0; c < 2; ++c)
4784  Assert(cell->line(line)->child(c)->user_flag_set() == false,
4785  ExcInternalError());
4786 #endif
4787 
4789  // Do refinement on every level
4790  //
4791  // To make life a bit easier, we first refine those lines and
4792  // quads that were flagged for refinement and then compose the
4793  // newly to be created cells.
4794  //
4795  // index of next unused vertex
4796  unsigned int next_unused_vertex = 0;
4797 
4798  // first for lines
4799  {
4800  // only active objects can be refined further
4802  line = triangulation.begin_active_line(),
4803  endl = triangulation.end_line();
4805  next_unused_line = triangulation.begin_raw_line();
4806 
4807  for (; line != endl; ++line)
4808  if (line->user_flag_set())
4809  {
4810  // this line needs to be refined
4811 
4812  // find the next unused vertex and set it
4813  // appropriately
4814  while (triangulation.vertices_used[next_unused_vertex] == true)
4815  ++next_unused_vertex;
4816  Assert(
4817  next_unused_vertex < triangulation.vertices.size(),
4818  ExcMessage(
4819  "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."));
4820  triangulation.vertices_used[next_unused_vertex] = true;
4821 
4822  triangulation.vertices[next_unused_vertex] = line->center(true);
4823 
4824  // now that we created the right point, make up the
4825  // two child lines (++ takes care of the end of the
4826  // vector)
4827  next_unused_line =
4828  triangulation.faces->lines.template next_free_pair_object<1>(
4829  triangulation);
4830  Assert(next_unused_line.state() == IteratorState::valid,
4831  ExcInternalError());
4832 
4833  // now we found two consecutive unused lines, such
4834  // that the children of a line will be consecutive.
4835  // then set the child pointer of the present line
4836  line->set_children(0, next_unused_line->index());
4837 
4838  // set the two new lines
4840  children[2] = {next_unused_line, ++next_unused_line};
4841 
4842  // some tests; if any of the iterators should be
4843  // invalid, then already dereferencing will fail
4844  Assert(
4845  children[0]->used() == false,
4846  ExcMessage(
4847  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4848  Assert(
4849  children[1]->used() == false,
4850  ExcMessage(
4851  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4852 
4853  children[0]->set_bounding_object_indices(
4854  {line->vertex_index(0), next_unused_vertex});
4855  children[1]->set_bounding_object_indices(
4856  {next_unused_vertex, line->vertex_index(1)});
4857 
4858  children[0]->set_used_flag();
4859  children[1]->set_used_flag();
4860  children[0]->clear_children();
4861  children[1]->clear_children();
4862  children[0]->clear_user_data();
4863  children[1]->clear_user_data();
4864  children[0]->clear_user_flag();
4865  children[1]->clear_user_flag();
4866 
4867  children[0]->set_boundary_id_internal(line->boundary_id());
4868  children[1]->set_boundary_id_internal(line->boundary_id());
4869 
4870  children[0]->set_manifold_id(line->manifold_id());
4871  children[1]->set_manifold_id(line->manifold_id());
4872 
4873  // finally clear flag
4874  // indicating the need
4875  // for refinement
4876  line->clear_user_flag();
4877  }
4878  }
4879 
4880 
4882  // now refine marked quads
4884 
4885  // here we encounter several cases:
4886 
4887  // a) the quad is unrefined and shall be refined isotropically
4888 
4889  // b) the quad is unrefined and shall be refined
4890  // anisotropically
4891 
4892  // c) the quad is unrefined and shall be refined both
4893  // anisotropically and isotropically (this is reduced to case
4894  // b) and then case b) for the children again)
4895 
4896  // d) the quad is refined anisotropically and shall be refined
4897  // isotropically (this is reduced to case b) for the
4898  // anisotropic children)
4899 
4900  // e) the quad is refined isotropically and shall be refined
4901  // anisotropically (this is transformed to case c), however we
4902  // might have to renumber/rename children...)
4903 
4904  // we need a loop in cases c) and d), as the anisotropic
4905  // children migt have a lower index than the mother quad
4906  for (unsigned int loop = 0; loop < 2; ++loop)
4907  {
4908  // usually, only active objects can be refined
4909  // further. however, in cases d) and e) that is not true,
4910  // so we have to use 'normal' iterators here
4912  quad = triangulation.begin_quad(),
4913  endq = triangulation.end_quad();
4915  next_unused_line = triangulation.begin_raw_line();
4917  next_unused_quad = triangulation.begin_raw_quad();
4918 
4919  for (; quad != endq; ++quad)
4920  {
4921  if (quad->user_index())
4922  {
4923  RefinementCase<dim - 1> aniso_quad_ref_case =
4924  face_refinement_cases[quad->user_index()];
4925  // there is one unlikely event here, where we
4926  // already have refind the face: if the face was
4927  // refined anisotropically and we want to refine
4928  // it isotropically, both children are flagged for
4929  // anisotropic refinement. however, if those
4930  // children were already flagged for anisotropic
4931  // refinement, they might already be processed and
4932  // refined.
4933  if (aniso_quad_ref_case == quad->refinement_case())
4934  continue;
4935 
4936  Assert(quad->refinement_case() ==
4937  RefinementCase<dim - 1>::cut_xy ||
4938  quad->refinement_case() ==
4939  RefinementCase<dim - 1>::no_refinement,
4940  ExcInternalError());
4941 
4942  // this quad needs to be refined anisotropically
4943  Assert(quad->user_index() ==
4944  RefinementCase<dim - 1>::cut_x ||
4945  quad->user_index() ==
4946  RefinementCase<dim - 1>::cut_y,
4947  ExcInternalError());
4948 
4949  // make the new line interior to the quad
4951  new_line;
4952 
4953  new_line =
4954  triangulation.faces->lines
4955  .template next_free_single_object<1>(triangulation);
4956  Assert(
4957  new_line->used() == false,
4958  ExcMessage(
4959  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4960 
4961  // first collect the
4962  // indices of the vertices:
4963  // *--1--*
4964  // | | |
4965  // | | | cut_x
4966  // | | |
4967  // *--0--*
4968  //
4969  // *-----*
4970  // | |
4971  // 0-----1 cut_y
4972  // | |
4973  // *-----*
4974  unsigned int vertex_indices[2];
4975  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
4976  {
4977  vertex_indices[0] =
4978  quad->line(2)->child(0)->vertex_index(1);
4979  vertex_indices[1] =
4980  quad->line(3)->child(0)->vertex_index(1);
4981  }
4982  else
4983  {
4984  vertex_indices[0] =
4985  quad->line(0)->child(0)->vertex_index(1);
4986  vertex_indices[1] =
4987  quad->line(1)->child(0)->vertex_index(1);
4988  }
4989 
4990  new_line->set_bounding_object_indices(
4991  {vertex_indices[0], vertex_indices[1]});
4992  new_line->set_used_flag();
4993  new_line->clear_user_flag();
4994  new_line->clear_user_data();
4995  new_line->clear_children();
4996  new_line->set_boundary_id_internal(quad->boundary_id());
4997  new_line->set_manifold_id(quad->manifold_id());
4998 
4999  // child 0 and 1 of a line are switched if the
5000  // line orientation is false. set up a miniature
5001  // table, indicating which child to take for line
5002  // orientations false and true. first index: child
5003  // index in standard orientation, second index:
5004  // line orientation
5005  const unsigned int index[2][2] = {
5006  {1, 0}, // child 0, line_orientation=false and true
5007  {0, 1}}; // child 1, line_orientation=false and true
5008 
5009  // find some space (consecutive) for the two newly
5010  // to be created quads.
5012  new_quads[2];
5013 
5014  next_unused_quad =
5015  triangulation.faces->quads
5016  .template next_free_pair_object<2>(triangulation);
5017  new_quads[0] = next_unused_quad;
5018  Assert(
5019  new_quads[0]->used() == false,
5020  ExcMessage(
5021  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5022 
5023  ++next_unused_quad;
5024  new_quads[1] = next_unused_quad;
5025  Assert(
5026  new_quads[1]->used() == false,
5027  ExcMessage(
5028  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5029 
5030 
5031  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5032  {
5033  new_quads[0]->set_bounding_object_indices(
5034  {static_cast<int>(quad->line_index(0)),
5035  new_line->index(),
5036  quad->line(2)
5037  ->child(index[0][quad->line_orientation(2)])
5038  ->index(),
5039  quad->line(3)
5040  ->child(index[0][quad->line_orientation(3)])
5041  ->index()});
5042  new_quads[1]->set_bounding_object_indices(
5043  {new_line->index(),
5044  static_cast<int>(quad->line_index(1)),
5045  quad->line(2)
5046  ->child(index[1][quad->line_orientation(2)])
5047  ->index(),
5048  quad->line(3)
5049  ->child(index[1][quad->line_orientation(3)])
5050  ->index()});
5051  }
5052  else
5053  {
5054  new_quads[0]->set_bounding_object_indices(
5055  {quad->line(0)
5056  ->child(index[0][quad->line_orientation(0)])
5057  ->index(),
5058  quad->line(1)
5059  ->child(index[0][quad->line_orientation(1)])
5060  ->index(),
5061  static_cast<int>(quad->line_index(2)),
5062  new_line->index()});
5063  new_quads[1]->set_bounding_object_indices(
5064  {quad->line(0)
5065  ->child(index[1][quad->line_orientation(0)])
5066  ->index(),
5067  quad->line(1)
5068  ->child(index[1][quad->line_orientation(1)])
5069  ->index(),
5070  new_line->index(),
5071  static_cast<int>(quad->line_index(3))});
5072  }
5073 
5074  for (const auto &new_quad : new_quads)
5075  {
5076  new_quad->set_used_flag();
5077  new_quad->clear_user_flag();
5078  new_quad->clear_user_data();
5079  new_quad->clear_children();
5080  new_quad->set_boundary_id_internal(quad->boundary_id());
5081  new_quad->set_manifold_id(quad->manifold_id());
5082  // set all line orientations to true, change
5083  // this after the loop, as we have to consider
5084  // different lines for each child
5085  for (unsigned int j = 0;
5086  j < GeometryInfo<dim>::lines_per_face;
5087  ++j)
5088  new_quad->set_line_orientation(j, true);
5089  }
5090  // now set the line orientation of children of
5091  // outer lines correctly, the lines in the
5092  // interior of the refined quad are automatically
5093  // oriented conforming to the standard
5094  new_quads[0]->set_line_orientation(
5095  0, quad->line_orientation(0));
5096  new_quads[0]->set_line_orientation(
5097  2, quad->line_orientation(2));
5098  new_quads[1]->set_line_orientation(
5099  1, quad->line_orientation(1));
5100  new_quads[1]->set_line_orientation(
5101  3, quad->line_orientation(3));
5102  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5103  {
5104  new_quads[0]->set_line_orientation(
5105  3, quad->line_orientation(3));
5106  new_quads[1]->set_line_orientation(
5107  2, quad->line_orientation(2));
5108  }
5109  else
5110  {
5111  new_quads[0]->set_line_orientation(
5112  1, quad->line_orientation(1));
5113  new_quads[1]->set_line_orientation(
5114  0, quad->line_orientation(0));
5115  }
5116 
5117  // test, whether this face is refined
5118  // isotropically already. if so, set the correct
5119  // children pointers.
5120  if (quad->refinement_case() ==
5121  RefinementCase<dim - 1>::cut_xy)
5122  {
5123  // we will put a new refinemnt level of
5124  // anisotropic refinement between the
5125  // unrefined and isotropically refined quad
5126  // ending up with the same fine quads but
5127  // introducing anisotropically refined ones as
5128  // children of the unrefined quad and mother
5129  // cells of the original fine ones.
5130 
5131  // this process includes the creation of a new
5132  // middle line which we will assign as the
5133  // mother line of two of the existing inner
5134  // lines. If those inner lines are not
5135  // consecutive in memory, we won't find them
5136  // later on, so we have to create new ones
5137  // instead and replace all occurrences of the
5138  // old ones with those new ones. As this is
5139  // kind of ugly, we hope we don't have to do
5140  // it often...
5142  old_child[2];
5143  if (aniso_quad_ref_case ==
5145  {
5146  old_child[0] = quad->child(0)->line(1);
5147  old_child[1] = quad->child(2)->line(1);
5148  }
5149  else
5150  {
5151  Assert(aniso_quad_ref_case ==
5153  ExcInternalError());
5154 
5155  old_child[0] = quad->child(0)->line(3);
5156  old_child[1] = quad->child(1)->line(3);
5157  }
5158 
5159  if (old_child[0]->index() + 1 != old_child[1]->index())
5160  {
5161  // this is exactly the ugly case we taked
5162  // about. so, no coimplaining, lets get
5163  // two new lines and copy all info
5164  typename Triangulation<dim,
5165  spacedim>::raw_line_iterator
5166  new_child[2];
5167 
5168  new_child[0] = new_child[1] =
5169  triangulation.faces->lines
5170  .template next_free_pair_object<1>(
5171  triangulation);
5172  ++new_child[1];
5173 
5174  new_child[0]->set_used_flag();
5175  new_child[1]->set_used_flag();
5176 
5177  const int old_index_0 = old_child[0]->index(),
5178  old_index_1 = old_child[1]->index(),
5179  new_index_0 = new_child[0]->index(),
5180  new_index_1 = new_child[1]->index();
5181 
5182  // loop over all quads and replace the old
5183  // lines
5184  for (unsigned int q = 0;
5185  q < triangulation.faces->quads.n_objects();
5186  ++q)
5187  for (unsigned int l = 0;
5188  l < GeometryInfo<dim>::lines_per_face;
5189  ++l)
5190  {
5191  const int this_index =
5192  triangulation.faces->quads
5193  .get_bounding_object_indices(q)[l];
5194  if (this_index == old_index_0)
5195  triangulation.faces->quads
5196  .get_bounding_object_indices(q)[l] =
5197  new_index_0;
5198  else if (this_index == old_index_1)
5199  triangulation.faces->quads
5200  .get_bounding_object_indices(q)[l] =
5201  new_index_1;
5202  }
5203  // now we have to copy all information of
5204  // the two lines
5205  for (unsigned int i = 0; i < 2; ++i)
5206  {
5207  Assert(!old_child[i]->has_children(),
5208  ExcInternalError());
5209 
5210  new_child[i]->set_bounding_object_indices(
5211  {old_child[i]->vertex_index(0),
5212  old_child[i]->vertex_index(1)});
5213  new_child[i]->set_boundary_id_internal(
5214  old_child[i]->boundary_id());
5215  new_child[i]->set_manifold_id(
5216  old_child[i]->manifold_id());
5217  new_child[i]->set_user_index(
5218  old_child[i]->user_index());
5219  if (old_child[i]->user_flag_set())
5220  new_child[i]->set_user_flag();
5221  else
5222  new_child[i]->clear_user_flag();
5223 
5224  new_child[i]->clear_children();
5225 
5226  old_child[i]->clear_user_flag();
5227  old_child[i]->clear_user_index();
5228  old_child[i]->clear_used_flag();
5229  }
5230  }
5231  // now that we cared about the lines, go on
5232  // with the quads themselves, where we might
5233  // encounter similar situations...
5234  if (aniso_quad_ref_case ==
5236  {
5237  new_line->set_children(
5238  0, quad->child(0)->line_index(1));
5239  Assert(new_line->child(1) ==
5240  quad->child(2)->line(1),
5241  ExcInternalError());
5242  // now evereything is quite
5243  // complicated. we have the children
5244  // numbered according to
5245  //
5246  // *---*---*
5247  // |n+2|n+3|
5248  // *---*---*
5249  // | n |n+1|
5250  // *---*---*
5251  //
5252  // from the original isotropic
5253  // refinement. we have to reorder them as
5254  //
5255  // *---*---*
5256  // |n+1|n+3|
5257  // *---*---*
5258  // | n |n+2|
5259  // *---*---*
5260  //
5261  // such that n and n+1 are consecutive
5262  // children of m and n+2 and n+3 are
5263  // consecutive children of m+1, where m
5264  // and m+1 are given as in
5265  //
5266  // *---*---*
5267  // | | |
5268  // | m |m+1|
5269  // | | |
5270  // *---*---*
5271  //
5272  // this is a bit ugly, of course: loop
5273  // over all cells on all levels and look
5274  // for faces n+1 (switch_1) and n+2
5275  // (switch_2).
5276  const typename Triangulation<dim, spacedim>::
5277  quad_iterator switch_1 = quad->child(1),
5278  switch_2 = quad->child(2);
5279  const int switch_1_index = switch_1->index();
5280  const int switch_2_index = switch_2->index();
5281  for (unsigned int l = 0;
5282  l < triangulation.levels.size();
5283  ++l)
5284  for (unsigned int h = 0;
5285  h <
5286  triangulation.levels[l]->cells.n_objects();
5287  ++h)
5288  for (const unsigned int q :
5290  {
5291  const int face_index =
5292  triangulation.levels[l]
5293  ->cells.get_bounding_object_indices(
5294  h)[q];
5295  if (face_index == switch_1_index)
5296  triangulation.levels[l]
5297  ->cells.get_bounding_object_indices(
5298  h)[q] = switch_2_index;
5299  else if (face_index == switch_2_index)
5300  triangulation.levels[l]
5301  ->cells.get_bounding_object_indices(
5302  h)[q] = switch_1_index;
5303  }
5304  // now we have to copy all information of
5305  // the two quads
5306  const unsigned int switch_1_lines[4] = {
5307  switch_1->line_index(0),
5308  switch_1->line_index(1),
5309  switch_1->line_index(2),
5310  switch_1->line_index(3)};
5311  const bool switch_1_line_orientations[4] = {
5312  switch_1->line_orientation(0),
5313  switch_1->line_orientation(1),
5314  switch_1->line_orientation(2),
5315  switch_1->line_orientation(3)};
5316  const types::boundary_id switch_1_boundary_id =
5317  switch_1->boundary_id();
5318  const unsigned int switch_1_user_index =
5319  switch_1->user_index();
5320  const bool switch_1_user_flag =
5321  switch_1->user_flag_set();
5322  const RefinementCase<dim - 1>
5323  switch_1_refinement_case =
5324  switch_1->refinement_case();
5325  const int switch_1_first_child_pair =
5326  (switch_1_refinement_case ?
5327  switch_1->child_index(0) :
5328  -1);
5329  const int switch_1_second_child_pair =
5330  (switch_1_refinement_case ==
5331  RefinementCase<dim - 1>::cut_xy ?
5332  switch_1->child_index(2) :
5333  -1);
5334 
5335  switch_1->set_bounding_object_indices(
5336  {switch_2->line_index(0),
5337  switch_2->line_index(1),
5338  switch_2->line_index(2),
5339  switch_2->line_index(3)});
5340  switch_1->set_line_orientation(
5341  0, switch_2->line_orientation(0));
5342  switch_1->set_line_orientation(
5343  1, switch_2->line_orientation(1));
5344  switch_1->set_line_orientation(
5345  2, switch_2->line_orientation(2));
5346  switch_1->set_line_orientation(
5347  3, switch_2->line_orientation(3));
5348  switch_1->set_boundary_id_internal(
5349  switch_2->boundary_id());
5350  switch_1->set_manifold_id(switch_2->manifold_id());
5351  switch_1->set_user_index(switch_2->user_index());
5352  if (switch_2->user_flag_set())
5353  switch_1->set_user_flag();
5354  else
5355  switch_1->clear_user_flag();
5356  switch_1->clear_refinement_case();
5357  switch_1->set_refinement_case(
5358  switch_2->refinement_case());
5359  switch_1->clear_children();
5360  if (switch_2->refinement_case())
5361  switch_1->set_children(0,
5362  switch_2->child_index(0));
5363  if (switch_2->refinement_case() ==
5364  RefinementCase<dim - 1>::cut_xy)
5365  switch_1->set_children(2,
5366  switch_2->child_index(2));
5367 
5368  switch_2->set_bounding_object_indices(
5369  {switch_1_lines[0],
5370  switch_1_lines[1],
5371  switch_1_lines[2],
5372  switch_1_lines[3]});
5373  switch_2->set_line_orientation(
5374  0, switch_1_line_orientations[0]);
5375  switch_2->set_line_orientation(
5376  1, switch_1_line_orientations[1]);
5377  switch_2->set_line_orientation(
5378  2, switch_1_line_orientations[2]);
5379  switch_2->set_line_orientation(
5380  3, switch_1_line_orientations[3]);
5381  switch_2->set_boundary_id_internal(
5382  switch_1_boundary_id);
5383  switch_2->set_manifold_id(switch_1->manifold_id());
5384  switch_2->set_user_index(switch_1_user_index);
5385  if (switch_1_user_flag)
5386  switch_2->set_user_flag();
5387  else
5388  switch_2->clear_user_flag();
5389  switch_2->clear_refinement_case();
5390  switch_2->set_refinement_case(
5391  switch_1_refinement_case);
5392  switch_2->clear_children();
5393  switch_2->set_children(0,
5394  switch_1_first_child_pair);
5395  switch_2->set_children(2,
5396  switch_1_second_child_pair);
5397 
5398  new_quads[0]->set_refinement_case(
5400  new_quads[0]->set_children(0, quad->child_index(0));
5401  new_quads[1]->set_refinement_case(
5403  new_quads[1]->set_children(0, quad->child_index(2));
5404  }
5405  else
5406  {
5407  new_quads[0]->set_refinement_case(
5409  new_quads[0]->set_children(0, quad->child_index(0));
5410  new_quads[1]->set_refinement_case(
5412  new_quads[1]->set_children(0, quad->child_index(2));
5413  new_line->set_children(
5414  0, quad->child(0)->line_index(3));
5415  Assert(new_line->child(1) ==
5416  quad->child(1)->line(3),
5417  ExcInternalError());
5418  }
5419  quad->clear_children();
5420  }
5421 
5422  // note these quads as children to the present one
5423  quad->set_children(0, new_quads[0]->index());
5424 
5425  quad->set_refinement_case(aniso_quad_ref_case);
5426 
5427  // finally clear flag indicating the need for
5428  // refinement
5429  quad->clear_user_data();
5430  } // if (anisotropic refinement)
5431 
5432  if (quad->user_flag_set())
5433  {
5434  // this quad needs to be refined isotropically
5435 
5436  // first of all: we only get here in the first run
5437  // of the loop
5438  Assert(loop == 0, ExcInternalError());
5439 
5440  // find the next unused vertex. we'll need this in
5441  // any case
5442  while (triangulation.vertices_used[next_unused_vertex] ==
5443  true)
5444  ++next_unused_vertex;
5445  Assert(
5446  next_unused_vertex < triangulation.vertices.size(),
5447  ExcMessage(
5448  "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."));
5449 
5450  // now: if the quad is refined anisotropically
5451  // already, set the anisotropic refinement flag
5452  // for both children. Additionally, we have to
5453  // refine the inner line, as it is an outer line
5454  // of the two (anisotropic) children
5455  const RefinementCase<dim - 1> quad_ref_case =
5456  quad->refinement_case();
5457 
5458  if (quad_ref_case == RefinementCase<dim - 1>::cut_x ||
5459  quad_ref_case == RefinementCase<dim - 1>::cut_y)
5460  {
5461  // set the 'opposite' refine case for children
5462  quad->child(0)->set_user_index(
5463  RefinementCase<dim - 1>::cut_xy - quad_ref_case);
5464  quad->child(1)->set_user_index(
5465  RefinementCase<dim - 1>::cut_xy - quad_ref_case);
5466  // refine the inner line
5468  middle_line;
5469  if (quad_ref_case == RefinementCase<dim - 1>::cut_x)
5470  middle_line = quad->child(0)->line(1);
5471  else
5472  middle_line = quad->child(0)->line(3);
5473 
5474  // if the face has been refined
5475  // anisotropically in the last refinement step
5476  // it might be, that it is flagged already and
5477  // that the middle line is thus refined
5478  // already. if not create children.
5479  if (!middle_line->has_children())
5480  {
5481  // set the middle vertex
5482  // appropriately. double refinement of
5483  // quads can only happen in the interior
5484  // of the domain, so we need not care
5485  // about boundary quads here
5486  triangulation.vertices[next_unused_vertex] =
5487  middle_line->center(true);
5488  triangulation.vertices_used[next_unused_vertex] =
5489  true;
5490 
5491  // now search a slot for the two
5492  // child lines
5493  next_unused_line =
5494  triangulation.faces->lines
5495  .template next_free_pair_object<1>(
5496  triangulation);
5497 
5498  // set the child pointer of the present
5499  // line
5500  middle_line->set_children(
5501  0, next_unused_line->index());
5502 
5503  // set the two new lines
5504  const typename Triangulation<dim, spacedim>::
5505  raw_line_iterator children[2] = {
5506  next_unused_line, ++next_unused_line};
5507 
5508  // some tests; if any of the iterators
5509  // should be invalid, then already
5510  // dereferencing will fail
5511  Assert(
5512  children[0]->used() == false,
5513  ExcMessage(
5514  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5515  Assert(
5516  children[1]->used() == false,
5517  ExcMessage(
5518  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5519 
5520  children[0]->set_bounding_object_indices(
5521  {middle_line->vertex_index(0),
5522  next_unused_vertex});
5523  children[1]->set_bounding_object_indices(
5524  {next_unused_vertex,
5525  middle_line->vertex_index(1)});
5526 
5527  children[0]->set_used_flag();
5528  children[1]->set_used_flag();
5529  children[0]->clear_children();
5530  children[1]->clear_children();
5531  children[0]->clear_user_data();
5532  children[1]->clear_user_data();
5533  children[0]->clear_user_flag();
5534  children[1]->clear_user_flag();
5535 
5536  children[0]->set_boundary_id_internal(
5537  middle_line->boundary_id());
5538  children[1]->set_boundary_id_internal(
5539  middle_line->boundary_id());
5540 
5541  children[0]->set_manifold_id(
5542  middle_line->manifold_id());
5543  children[1]->set_manifold_id(
5544  middle_line->manifold_id());
5545  }
5546  // now remove the flag from the quad and go to
5547  // the next quad, the actual refinement of the
5548  // quad takes place later on in this pass of
5549  // the loop or in the next one
5550  quad->clear_user_flag();
5551  continue;
5552  } // if (several refinement cases)
5553 
5554  // if we got here, we have an unrefined quad and
5555  // have to do the usual work like in an purely
5556  // isotropic refinement
5557  Assert(quad_ref_case ==
5559  ExcInternalError());
5560 
5561  // set the middle vertex appropriately: it might be that
5562  // the quad itself is not at the boundary, but that one of
5563  // its lines actually is. in this case, the newly created
5564  // vertices at the centers of the lines are not
5565  // necessarily the mean values of the adjacent vertices,
5566  // so do not compute the new vertex as the mean value of
5567  // the 4 vertices of the face, but rather as a weighted
5568  // mean value of the 8 vertices which we already have (the
5569  // four old ones, and the four ones inserted as middle
5570  // points for the four lines). summing up some more points
5571  // is generally cheaper than first asking whether one of
5572  // the lines is at the boundary
5573  //
5574  // note that the exact weights are chosen such as to
5575  // minimize the distortion of the four new quads from the
5576  // optimal shape. their description uses the formulas
5577  // underlying the TransfiniteInterpolationManifold
5578  // implementation
5579  triangulation.vertices[next_unused_vertex] =
5580  quad->center(true, true);
5581  triangulation.vertices_used[next_unused_vertex] = true;
5582 
5583  // now that we created the right point, make up
5584  // the four lines interior to the quad (++ takes
5585  // care of the end of the vector)
5587  new_lines[4];
5588 
5589  for (unsigned int i = 0; i < 4; ++i)
5590  {
5591  if (i % 2 == 0)
5592  // search a free pair of lines for 0. and
5593  // 2. line, so that two of them end up
5594  // together, which is necessary if later on
5595  // we want to refine the quad
5596  // anisotropically and the two lines end up
5597  // as children of new line
5598  next_unused_line =
5599  triangulation.faces->lines
5600  .template next_free_pair_object<1>(triangulation);
5601 
5602  new_lines[i] = next_unused_line;
5603  ++next_unused_line;
5604 
5605  Assert(
5606  new_lines[i]->used() == false,
5607  ExcMessage(
5608  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5609  }
5610 
5611  // set the data of the four lines. first collect
5612  // the indices of the five vertices:
5613  //
5614  // *--3--*
5615  // | | |
5616  // 0--4--1
5617  // | | |
5618  // *--2--*
5619  //
5620  // the lines are numbered as follows:
5621  //
5622  // *--*--*
5623  // | 1 |
5624  // *2-*-3*
5625  // | 0 |
5626  // *--*--*
5627 
5628  const unsigned int vertex_indices[5] = {
5629  quad->line(0)->child(0)->vertex_index(1),
5630  quad->line(1)->child(0)->vertex_index(1),
5631  quad->line(2)->child(0)->vertex_index(1),
5632  quad->line(3)->child(0)->vertex_index(1),
5633  next_unused_vertex};
5634 
5635  new_lines[0]->set_bounding_object_indices(
5636  {vertex_indices[2], vertex_indices[4]});
5637  new_lines[1]->set_bounding_object_indices(
5638  {vertex_indices[4], vertex_indices[3]});
5639  new_lines[2]->set_bounding_object_indices(
5640  {vertex_indices[0], vertex_indices[4]});
5641  new_lines[3]->set_bounding_object_indices(
5642  {vertex_indices[4], vertex_indices[1]});
5643 
5644  for (const auto &new_line : new_lines)
5645  {
5646  new_line->set_used_flag();
5647  new_line->clear_user_flag();
5648  new_line->clear_user_data();
5649  new_line->clear_children();
5650  new_line->set_boundary_id_internal(quad->boundary_id());
5651  new_line->set_manifold_id(quad->manifold_id());
5652  }
5653 
5654  // now for the quads. again, first collect some
5655  // data about the indices of the lines, with the
5656  // following numbering:
5657  //
5658  // .-6-.-7-.
5659  // 1 9 3
5660  // .-10.11-.
5661  // 0 8 2
5662  // .-4-.-5-.
5663 
5664  // child 0 and 1 of a line are switched if the
5665  // line orientation is false. set up a miniature
5666  // table, indicating which child to take for line
5667  // orientations false and true. first index: child
5668  // index in standard orientation, second index:
5669  // line orientation
5670  const unsigned int index[2][2] = {
5671  {1, 0}, // child 0, line_orientation=false and true
5672  {0, 1}}; // child 1, line_orientation=false and true
5673 
5674  const int line_indices[12] = {
5675  quad->line(0)
5676  ->child(index[0][quad->line_orientation(0)])
5677  ->index(),
5678  quad->line(0)
5679  ->child(index[1][quad->line_orientation(0)])
5680  ->index(),
5681  quad->line(1)
5682  ->child(index[0][quad->line_orientation(1)])
5683  ->index(),
5684  quad->line(1)
5685  ->child(index[1][quad->line_orientation(1)])
5686  ->index(),
5687  quad->line(2)
5688  ->child(index[0][quad->line_orientation(2)])
5689  ->index(),
5690  quad->line(2)
5691  ->child(index[1][quad->line_orientation(2)])
5692  ->index(),
5693  quad->line(3)
5694  ->child(index[0][quad->line_orientation(3)])
5695  ->index(),
5696  quad->line(3)
5697  ->child(index[1][quad->line_orientation(3)])
5698  ->index(),
5699  new_lines[0]->index(),
5700  new_lines[1]->index(),
5701  new_lines[2]->index(),
5702  new_lines[3]->index()};
5703 
5704  // find some space (consecutive)
5705  // for the first two newly to be
5706  // created quads.
5708  new_quads[4];
5709 
5710  next_unused_quad =
5711  triangulation.faces->quads
5712  .template next_free_pair_object<2>(triangulation);
5713 
5714  new_quads[0] = next_unused_quad;
5715  Assert(
5716  new_quads[0]->used() == false,
5717  ExcMessage(
5718  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5719 
5720  ++next_unused_quad;
5721  new_quads[1] = next_unused_quad;
5722  Assert(
5723  new_quads[1]->used() == false,
5724  ExcMessage(
5725  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5726 
5727  next_unused_quad =
5728  triangulation.faces->quads
5729  .template next_free_pair_object<2>(triangulation);
5730  new_quads[2] = next_unused_quad;
5731  Assert(
5732  new_quads[2]->used() == false,
5733  ExcMessage(
5734  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5735 
5736  ++next_unused_quad;
5737  new_quads[3] = next_unused_quad;
5738  Assert(
5739  new_quads[3]->used() == false,
5740  ExcMessage(
5741  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5742 
5743  // note these quads as children to the present one
5744  quad->set_children(0, new_quads[0]->index());
5745  quad->set_children(2, new_quads[2]->index());
5746  quad->set_refinement_case(RefinementCase<2>::cut_xy);
5747 
5748  new_quads[0]->set_bounding_object_indices(
5749  {line_indices[0],
5750  line_indices[8],
5751  line_indices[4],
5752  line_indices[10]});
5753  new_quads[1]->set_bounding_object_indices(
5754  {line_indices[8],
5755  line_indices[2],
5756  line_indices[5],
5757  line_indices[11]});
5758  new_quads[2]->set_bounding_object_indices(
5759  {line_indices[1],
5760  line_indices[9],
5761  line_indices[10],
5762  line_indices[6]});
5763  new_quads[3]->set_bounding_object_indices(
5764  {line_indices[9],
5765  line_indices[3],
5766  line_indices[11],
5767  line_indices[7]});
5768  for (const auto &new_quad : new_quads)
5769  {
5770  new_quad->set_used_flag();
5771  new_quad->clear_user_flag();
5772  new_quad->clear_user_data();
5773  new_quad->clear_children();
5774  new_quad->set_boundary_id_internal(quad->boundary_id());
5775  new_quad->set_manifold_id(quad->manifold_id());
5776  // set all line orientations to true, change
5777  // this after the loop, as we have to consider
5778  // different lines for each child
5779  for (unsigned int j = 0;
5780  j < GeometryInfo<dim>::lines_per_face;
5781  ++j)
5782  new_quad->set_line_orientation(j, true);
5783  }
5784  // now set the line orientation of children of
5785  // outer lines correctly, the lines in the
5786  // interior of the refined quad are automatically
5787  // oriented conforming to the standard
5788  new_quads[0]->set_line_orientation(
5789  0, quad->line_orientation(0));
5790  new_quads[0]->set_line_orientation(
5791  2, quad->line_orientation(2));
5792  new_quads[1]->set_line_orientation(
5793  1, quad->line_orientation(1));
5794  new_quads[1]->set_line_orientation(
5795  2, quad->line_orientation(2));
5796  new_quads[2]->set_line_orientation(
5797  0, quad->line_orientation(0));
5798  new_quads[2]->set_line_orientation(
5799  3, quad->line_orientation(3));
5800  new_quads[3]->set_line_orientation(
5801  1, quad->line_orientation(1));
5802  new_quads[3]->set_line_orientation(
5803  3, quad->line_orientation(3));
5804 
5805  // finally clear flag indicating the need for
5806  // refinement
5807  quad->clear_user_flag();
5808  } // if (isotropic refinement)
5809  } // for all quads
5810  } // looped two times over all quads, all quads refined now
5811 
5813  // Now, finally, set up the new
5814  // cells
5816 
5818  cells_with_distorted_children;
5819 
5820  for (unsigned int level = 0; level != triangulation.levels.size() - 1;
5821  ++level)
5822  {
5823  // only active objects can be refined further; remember
5824  // that we won't operate on the finest level, so
5825  // triangulation.begin_*(level+1) is allowed
5827  hex = triangulation.begin_active_hex(level),
5828  endh = triangulation.begin_active_hex(level + 1);
5830  next_unused_hex = triangulation.begin_raw_hex(level + 1);
5831 
5832  for (; hex != endh; ++hex)
5833  if (hex->refine_flag_set())
5834  {
5835  // this hex needs to be refined
5836 
5837  // clear flag indicating the need for refinement. do
5838  // it here already, since we can't do it anymore
5839  // once the cell has children
5840  const RefinementCase<dim> ref_case = hex->refine_flag_set();
5841  hex->clear_refine_flag();
5842  hex->set_refinement_case(ref_case);
5843 
5844  // depending on the refine case we might have to
5845  // create additional vertices, lines and quads
5846  // interior of the hex before the actual children
5847  // can be set up.
5848 
5849  // in a first step: reserve the needed space for
5850  // lines, quads and hexes and initialize them
5851  // correctly
5852 
5853  unsigned int n_new_lines = 0;
5854  unsigned int n_new_quads = 0;
5855  unsigned int n_new_hexes = 0;
5856  switch (ref_case)
5857  {
5861  n_new_lines = 0;
5862  n_new_quads = 1;
5863  n_new_hexes = 2;
5864  break;
5868  n_new_lines = 1;
5869  n_new_quads = 4;
5870  n_new_hexes = 4;
5871  break;
5873  n_new_lines = 6;
5874  n_new_quads = 12;
5875  n_new_hexes = 8;
5876  break;
5877  default:
5878  Assert(false, ExcInternalError());
5879  break;
5880  }
5881 
5882  // find some space for the newly to be created
5883  // interior lines and initialize them.
5884  std::vector<
5886  new_lines(n_new_lines);
5887  for (unsigned int i = 0; i < n_new_lines; ++i)
5888  {
5889  new_lines[i] =
5890  triangulation.faces->lines
5891  .template next_free_single_object<1>(triangulation);
5892 
5893  Assert(
5894  new_lines[i]->used() == false,
5895  ExcMessage(
5896  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5897  new_lines[i]->set_used_flag();
5898  new_lines[i]->clear_user_flag();
5899  new_lines[i]->clear_user_data();
5900  new_lines[i]->clear_children();
5901  // interior line
5902  new_lines[i]->set_boundary_id_internal(
5904  // they inherit geometry description of the hex they
5905  // belong to
5906  new_lines[i]->set_manifold_id(hex->manifold_id());
5907  }
5908 
5909  // find some space for the newly to be created
5910  // interior quads and initialize them.
5911  std::vector<
5913  new_quads(n_new_quads);
5914  for (unsigned int i = 0; i < n_new_quads; ++i)
5915  {
5916  new_quads[i] =
5917  triangulation.faces->quads
5918  .template next_free_single_object<2>(triangulation);
5919 
5920  Assert(
5921  new_quads[i]->used() == false,
5922  ExcMessage(
5923  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5924  new_quads[i]->set_used_flag();
5925  new_quads[i]->clear_user_flag();
5926  new_quads[i]->clear_user_data();
5927  new_quads[i]->clear_children();
5928  // interior quad
5929  new_quads[i]->set_boundary_id_internal(
5931  // they inherit geometry description of the hex they
5932  // belong to
5933  new_quads[i]->set_manifold_id(hex->manifold_id());
5934  // set all line orientation flags to true by
5935  // default, change this afterwards, if necessary
5936  for (unsigned int j = 0;
5937  j < GeometryInfo<dim>::lines_per_face;
5938  ++j)
5939  new_quads[i]->set_line_orientation(j, true);
5940  }
5941 
5942  types::subdomain_id subdomainid = hex->subdomain_id();
5943 
5944  // find some space for the newly to be created hexes
5945  // and initialize them.
5946  std::vector<
5948  new_hexes(n_new_hexes);
5949  for (unsigned int i = 0; i < n_new_hexes; ++i)
5950  {
5951  if (i % 2 == 0)
5952  next_unused_hex =
5953  triangulation.levels[level + 1]->cells.next_free_hex(
5954  triangulation, level + 1);
5955  else
5956  ++next_unused_hex;
5957 
5958  new_hexes[i] = next_unused_hex;
5959 
5960  Assert(
5961  new_hexes[i]->used() == false,
5962  ExcMessage(
5963  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5964  new_hexes[i]->set_used_flag();
5965  new_hexes[i]->clear_user_flag();
5966  new_hexes[i]->clear_user_data();
5967  new_hexes[i]->clear_children();
5968  // inherit material
5969  // properties
5970  new_hexes[i]->set_material_id(hex->material_id());
5971  new_hexes[i]->set_manifold_id(hex->manifold_id());
5972  new_hexes[i]->set_subdomain_id(subdomainid);
5973 
5974  if (i % 2)
5975  new_hexes[i]->set_parent(hex->index());
5976  // set the face_orientation flag to true for all
5977  // faces initially, as this is the default value
5978  // which is true for all faces interior to the
5979  // hex. later on go the other way round and
5980  // reset faces that are at the boundary of the
5981  // mother cube
5982  //
5983  // the same is true for the face_flip and
5984  // face_rotation flags. however, the latter two
5985  // are set to false by default as this is the
5986  // standard value
5987  for (const unsigned int f :
5989  {
5990  new_hexes[i]->set_face_orientation(f, true);
5991  new_hexes[i]->set_face_flip(f, false);
5992  new_hexes[i]->set_face_rotation(f, false);
5993  }
5994  }
5995  // note these hexes as children to the present cell
5996  for (unsigned int i = 0; i < n_new_hexes / 2; ++i)
5997  hex->set_children(2 * i, new_hexes[2 * i]->index());
5998 
5999  // we have to take into account whether the
6000  // different faces are oriented correctly or in the
6001  // opposite direction, so store that up front
6002 
6003  // face_orientation
6004  const bool f_or[6] = {hex->face_orientation(0),
6005  hex->face_orientation(1),
6006  hex->face_orientation(2),
6007  hex->face_orientation(3),
6008  hex->face_orientation(4),
6009  hex->face_orientation(5)};
6010 
6011  // face_flip
6012  const bool f_fl[6] = {hex->face_flip(0),
6013  hex->face_flip(1),
6014  hex->face_flip(2),
6015  hex->face_flip(3),
6016  hex->face_flip(4),
6017  hex->face_flip(5)};
6018 
6019  // face_rotation
6020  const bool f_ro[6] = {hex->face_rotation(0),
6021  hex->face_rotation(1),
6022  hex->face_rotation(2),
6023  hex->face_rotation(3),
6024  hex->face_rotation(4),
6025  hex->face_rotation(5)};
6026 
6027  // little helper table, indicating, whether the
6028  // child with index 0 or with index 1 can be found
6029  // at the standard origin of an anisotropically
6030  // refined quads in real orientation index 1:
6031  // (RefineCase - 1) index 2: face_flip
6032 
6033  // index 3: face rotation
6034  // note: face orientation has no influence
6035  const unsigned int child_at_origin[2][2][2] = {
6036  {{0, 0}, // RefinementCase<dim>::cut_x, face_flip=false,
6037  // face_rotation=false and true
6038  {1, 1}}, // RefinementCase<dim>::cut_x, face_flip=true,
6039  // face_rotation=false and true
6040  {{0, 1}, // RefinementCase<dim>::cut_y, face_flip=false,
6041  // face_rotation=false and true
6042  {1, 0}}}; // RefinementCase<dim>::cut_y, face_flip=true,
6043  // face_rotation=false and true
6044 
6046  //
6047  // in the following we will do the same thing for
6048  // each refinement case: create a new vertex (if
6049  // needed), create new interior lines (if needed),
6050  // create new interior quads and afterwards build
6051  // the children hexes out of these and the existing
6052  // subfaces of the outer quads (which have been
6053  // created above). However, even if the steps are
6054  // quite similar, the actual work strongly depends
6055  // on the actual refinement case. therefore, we use
6056  // separate blocks of code for each of these cases,
6057  // which hopefully increases the readability to some
6058  // extend.
6059 
6060  switch (ref_case)
6061  {
6063  {
6065  //
6066  // RefinementCase<dim>::cut_x
6067  //
6068  // the refined cube will look
6069  // like this:
6070  //
6071  // *----*----*
6072  // / / /|
6073  // / / / |
6074  // / / / |
6075  // *----*----* |
6076  // | | | |
6077  // | | | *
6078  // | | | /
6079  // | | | /
6080  // | | |/
6081  // *----*----*
6082  //
6083  // again, first collect some data about the
6084  // indices of the lines, with the following
6085  // numbering:
6086 
6087  // face 2: front plane
6088  // (note: x,y exchanged)
6089  // *---*---*
6090  // | | |
6091  // | 0 |
6092  // | | |
6093  // *---*---*
6094  // m0
6095  // face 3: back plane
6096  // (note: x,y exchanged)
6097  // m1
6098  // *---*---*
6099  // | | |
6100  // | 1 |
6101  // | | |
6102  // *---*---*
6103  // face 4: bottom plane
6104  // *---*---*
6105  // / / /
6106  // / 2 /
6107  // / / /
6108  // *---*---*
6109  // m0
6110  // face 5: top plane
6111  // m1
6112  // *---*---*
6113  // / / /
6114  // / 3 /
6115  // / / /
6116  // *---*---*
6117 
6118  // set up a list of line iterators first. from
6119  // this, construct lists of line_indices and
6120  // line orientations later on
6121  const typename Triangulation<dim, spacedim>::
6122  raw_line_iterator lines[4] = {
6123  hex->face(2)->child(0)->line(
6124  (hex->face(2)->refinement_case() ==
6126  1 :
6127  3), // 0
6128  hex->face(3)->child(0)->line(
6129  (hex->face(3)->refinement_case() ==
6131  1 :
6132  3), // 1
6133  hex->face(4)->child(0)->line(
6134  (hex->face(4)->refinement_case() ==
6136  1 :
6137  3), // 2
6138  hex->face(5)->child(0)->line(
6139  (hex->face(5)->refinement_case() ==
6141  1 :
6142  3) // 3
6143  };
6144 
6145  unsigned int line_indices[4];
6146  for (unsigned int i = 0; i < 4; ++i)
6147  line_indices[i] = lines[i]->index();
6148 
6149  // the orientation of lines for the inner quads
6150  // is quite tricky. as these lines are newly
6151  // created ones and thus have no parents, they
6152  // cannot inherit this property. set up an array
6153  // and fill it with the respective values
6154  bool line_orientation[4];
6155 
6156  // the middle vertex marked as m0 above is the
6157  // start vertex for lines 0 and 2 in standard
6158  // orientation, whereas m1 is the end vertex of
6159  // lines 1 and 3 in standard orientation
6160  const unsigned int middle_vertices[2] = {
6161  hex->line(2)->child(0)->vertex_index(1),
6162  hex->line(7)->child(0)->vertex_index(1)};
6163 
6164  for (unsigned int i = 0; i < 4; ++i)
6165  if (lines[i]->vertex_index(i % 2) ==
6166  middle_vertices[i % 2])
6167  line_orientation[i] = true;
6168  else
6169  {
6170  // it must be the other
6171  // way round then
6172  Assert(lines[i]->vertex_index((i + 1) % 2) ==
6173  middle_vertices[i % 2],
6174  ExcInternalError());
6175  line_orientation[i] = false;
6176  }
6177 
6178  // set up the new quad, line numbering is as
6179  // indicated above
6180  new_quads[0]->set_bounding_object_indices(
6181  {line_indices[0],
6182  line_indices[1],
6183  line_indices[2],
6184  line_indices[3]});
6185 
6186  new_quads[0]->set_line_orientation(
6187  0, line_orientation[0]);
6188  new_quads[0]->set_line_orientation(
6189  1, line_orientation[1]);
6190  new_quads[0]->set_line_orientation(
6191  2, line_orientation[2]);
6192  new_quads[0]->set_line_orientation(
6193  3, line_orientation[3]);
6194 
6195  // the quads are numbered as follows:
6196  //
6197  // planes in the interior of the old hex:
6198  //
6199  // *
6200  // /|
6201  // / | x
6202  // / | *-------* *---------*
6203  // * | | | / /
6204  // | 0 | | | / /
6205  // | * | | / /
6206  // | / *-------*y *---------*x
6207  // | /
6208  // |/
6209  // *
6210  //
6211  // children of the faces of the old hex
6212  //
6213  // *---*---* *---*---*
6214  // /| | | / / /|
6215  // / | | | / 9 / 10/ |
6216  // / | 5 | 6 | / / / |
6217  // * | | | *---*---* |
6218  // | 1 *---*---* | | | 2 *
6219  // | / / / | | | /
6220  // | / 7 / 8 / | 3 | 4 | /
6221  // |/ / / | | |/
6222  // *---*---* *---*---*
6223  //
6224  // note that we have to take care of the
6225  // orientation of faces.
6226  const int quad_indices[11] = {
6227  new_quads[0]->index(), // 0
6228 
6229  hex->face(0)->index(), // 1
6230 
6231  hex->face(1)->index(), // 2
6232 
6233  hex->face(2)->child_index(
6234  child_at_origin[hex->face(2)->refinement_case() -
6235  1][f_fl[2]][f_ro[2]]), // 3
6236  hex->face(2)->child_index(
6237  1 -
6238  child_at_origin[hex->face(2)->refinement_case() -
6239  1][f_fl[2]][f_ro[2]]),
6240 
6241  hex->face(3)->child_index(
6242  child_at_origin[hex->face(3)->refinement_case() -
6243  1][f_fl[3]][f_ro[3]]), // 5
6244  hex->face(3)->child_index(
6245  1 -
6246  child_at_origin[hex->face(3)->refinement_case() -
6247  1][f_fl[3]][f_ro[3]]),
6248 
6249  hex->face(4)->child_index(
6250  child_at_origin[hex->face(4)->refinement_case() -
6251  1][f_fl[4]][f_ro[4]]), // 7
6252  hex->face(4)->child_index(
6253  1 -
6254  child_at_origin[hex->face(4)->refinement_case() -
6255  1][f_fl[4]][f_ro[4]]),
6256 
6257  hex->face(5)->child_index(
6258  child_at_origin[hex->face(5)->refinement_case() -
6259  1][f_fl[5]][f_ro[5]]), // 9
6260  hex->face(5)->child_index(
6261  1 -
6262  child_at_origin[hex->face(5)->refinement_case() -
6263  1][f_fl[5]][f_ro[5]])
6264 
6265  };
6266 
6267  new_hexes[0]->set_bounding_object_indices(
6268  {quad_indices[1],
6269  quad_indices[0],
6270  quad_indices[3],
6271  quad_indices[5],
6272  quad_indices[7],
6273  quad_indices[9]});
6274  new_hexes[1]->set_bounding_object_indices(
6275  {quad_indices[0],
6276  quad_indices[2],
6277  quad_indices[4],
6278  quad_indices[6],
6279  quad_indices[8],
6280  quad_indices[10]});
6281  break;
6282  }
6283 
6285  {
6287  //
6288  // RefinementCase<dim>::cut_y
6289  //
6290  // the refined cube will look like this:
6291  //
6292  // *---------*
6293  // / /|
6294  // *---------* |
6295  // / /| |
6296  // *---------* | |
6297  // | | | |
6298  // | | | *
6299  // | | |/
6300  // | | *
6301  // | |/
6302  // *---------*
6303  //
6304  // again, first collect some data about the
6305  // indices of the lines, with the following
6306  // numbering:
6307 
6308  // face 0: left plane
6309  // *
6310  // /|
6311  // * |
6312  // /| |
6313  // * | |
6314  // | 0 |
6315  // | | *
6316  // | |/
6317  // | *m0
6318  // |/
6319  // *
6320  // face 1: right plane
6321  // *
6322  // /|
6323  // m1* |
6324  // /| |
6325  // * | |
6326  // | 1 |
6327  // | | *
6328  // | |/
6329  // | *
6330  // |/
6331  // *
6332  // face 4: bottom plane
6333  // *-------*
6334  // / /
6335  // m0*---2---*
6336  // / /
6337  // *-------*
6338  // face 5: top plane
6339  // *-------*
6340  // / /
6341  // *---3---*m1
6342  // / /
6343  // *-------*
6344 
6345  // set up a list of line iterators first. from
6346  // this, construct lists of line_indices and
6347  // line orientations later on
6348  const typename Triangulation<dim, spacedim>::
6349  raw_line_iterator lines[4] = {
6350  hex->face(0)->child(0)->line(
6351  (hex->face(0)->refinement_case() ==
6353  1 :
6354  3), // 0
6355  hex->face(1)->child(0)->line(
6356  (hex->face(1)->refinement_case() ==
6358  1 :
6359  3), // 1
6360  hex->face(4)->child(0)->line(
6361  (hex->face(4)->refinement_case() ==
6363  1 :
6364  3), // 2
6365  hex->face(5)->child(0)->line(
6366  (hex->face(5)->refinement_case() ==
6368  1 :
6369  3) // 3
6370  };
6371 
6372  unsigned int line_indices[4];
6373  for (unsigned int i = 0; i < 4; ++i)
6374  line_indices[i] = lines[i]->index();
6375 
6376  // the orientation of lines for the inner quads
6377  // is quite tricky. as these lines are newly
6378  // created ones and thus have no parents, they
6379  // cannot inherit this property. set up an array
6380  // and fill it with the respective values
6381  bool line_orientation[4];
6382 
6383  // the middle vertex marked as m0 above is the
6384  // start vertex for lines 0 and 2 in standard
6385  // orientation, whereas m1 is the end vertex of
6386  // lines 1 and 3 in standard orientation
6387  const unsigned int middle_vertices[2] = {
6388  hex->line(0)->child(0)->vertex_index(1),
6389  hex->line(5)->child(0)->vertex_index(1)};
6390 
6391  for (unsigned int i = 0; i < 4; ++i)
6392  if (lines[i]->vertex_index(i % 2) ==
6393  middle_vertices[i % 2])
6394  line_orientation[i] = true;
6395  else
6396  {
6397  // it must be the other way round then
6398  Assert(lines[i]->vertex_index((i + 1) % 2) ==
6399  middle_vertices[i % 2],
6400  ExcInternalError());
6401  line_orientation[i] = false;
6402  }
6403 
6404  // set up the new quad, line numbering is as
6405  // indicated above
6406  new_quads[0]->set_bounding_object_indices(
6407  {line_indices[2],
6408  line_indices[3],
6409  line_indices[0],
6410  line_indices[1]});
6411 
6412  new_quads[0]->set_line_orientation(
6413  0, line_orientation[2]);
6414  new_quads[0]->set_line_orientation(
6415  1, line_orientation[3]);
6416  new_quads[0]->set_line_orientation(
6417  2, line_orientation[0]);
6418  new_quads[0]->set_line_orientation(
6419  3, line_orientation[1]);
6420 
6421  // the quads are numbered as follows:
6422  //
6423  // planes in the interior of the old hex:
6424  //
6425  // *
6426  // /|
6427  // / | x
6428  // / | *-------* *---------*
6429  // * | | | / /
6430  // | | | 0 | / /
6431  // | * | | / /
6432  // | / *-------*y *---------*x
6433  // | /
6434  // |/
6435  // *
6436  //
6437  // children of the faces of the old hex
6438  //
6439  // *-------* *-------*
6440  // /| | / 10 /|
6441  // * | | *-------* |
6442  // /| | 6 | / 9 /| |
6443  // * |2| | *-------* |4|
6444  // | | *-------* | | | *
6445  // |1|/ 8 / | |3|/
6446  // | *-------* | 5 | *
6447  // |/ 7 / | |/
6448  // *-------* *-------*
6449  //
6450  // note that we have to take care of the
6451  // orientation of faces.
6452  const int quad_indices[11] = {
6453  new_quads[0]->index(), // 0
6454 
6455  hex->face(0)->child_index(
6456  child_at_origin[hex->face(0)->refinement_case() -
6457  1][f_fl[0]][f_ro[0]]), // 1
6458  hex->face(0)->child_index(
6459  1 -
6460  child_at_origin[hex->face(0)->refinement_case() -
6461  1][f_fl[0]][f_ro[0]]),
6462 
6463  hex->face(1)->child_index(
6464  child_at_origin[hex->face(1)->refinement_case() -
6465  1][f_fl[1]][f_ro[1]]), // 3
6466  hex->face(1)->child_index(
6467  1 -
6468  child_at_origin[hex->face(1)->refinement_case() -
6469  1][f_fl[1]][f_ro[1]]),
6470 
6471  hex->face(2)->index(), // 5
6472 
6473  hex->face(3)->index(), // 6
6474 
6475  hex->face(4)->child_index(
6476  child_at_origin[hex->face(4)->refinement_case() -
6477  1][f_fl[4]][f_ro[4]]), // 7
6478  hex->face(4)->child_index(
6479  1 -
6480  child_at_origin[hex->face(4)->refinement_case() -
6481  1][f_fl[4]][f_ro[4]]),
6482 
6483  hex->face(5)->child_index(
6484  child_at_origin[hex->face(5)->refinement_case() -
6485  1][f_fl[5]][f_ro[5]]), // 9
6486  hex->face(5)->child_index(
6487  1 -
6488  child_at_origin[hex->face(5)->refinement_case() -
6489  1][f_fl[5]][f_ro[5]])
6490 
6491  };
6492 
6493  new_hexes[0]->set_bounding_object_indices(
6494  {quad_indices[1],
6495  quad_indices[3],
6496  quad_indices[5],
6497  quad_indices[0],
6498  quad_indices[7],
6499  quad_indices[9]});
6500  new_hexes[1]->set_bounding_object_indices(
6501  {quad_indices[2],
6502  quad_indices[4],
6503  quad_indices[0],
6504  quad_indices[6],
6505  quad_indices[8],
6506  quad_indices[10]});
6507  break;
6508  }
6509 
6511  {
6513  //
6514  // RefinementCase<dim>::cut_z
6515  //
6516  // the refined cube will look like this:
6517  //
6518  // *---------*
6519  // / /|
6520  // / / |
6521  // / / *
6522  // *---------* /|
6523  // | | / |
6524  // | |/ *
6525  // *---------* /
6526  // | | /
6527  // | |/
6528  // *---------*
6529  //
6530  // again, first collect some data about the
6531  // indices of the lines, with the following
6532  // numbering:
6533 
6534  // face 0: left plane
6535  // *
6536  // /|
6537  // / |
6538  // / *
6539  // * /|
6540  // | 0 |
6541  // |/ *
6542  // m0* /
6543  // | /
6544  // |/
6545  // *
6546  // face 1: right plane
6547  // *
6548  // /|
6549  // / |
6550  // / *m1
6551  // * /|
6552  // | 1 |
6553  // |/ *
6554  // * /
6555  // | /
6556  // |/
6557  // *
6558  // face 2: front plane
6559  // (note: x,y exchanged)
6560  // *-------*
6561  // | |
6562  // m0*---2---*
6563  // | |
6564  // *-------*
6565  // face 3: back plane
6566  // (note: x,y exchanged)
6567  // *-------*
6568  // | |
6569  // *---3---*m1
6570  // | |
6571  // *-------*
6572 
6573  // set up a list of line iterators first. from
6574  // this, construct lists of line_indices and
6575  // line orientations later on
6576  const typename Triangulation<dim, spacedim>::
6577  raw_line_iterator lines[4] = {
6578  hex->face(0)->child(0)->line(
6579  (hex->face(0)->refinement_case() ==
6581  1 :
6582  3), // 0
6583  hex->face(1)->child(0)->line(
6584  (hex->face(1)->refinement_case() ==
6586  1 :
6587  3), // 1
6588  hex->face(2)->child(0)->line(
6589  (hex->face(2)->refinement_case() ==
6591  1 :
6592  3), // 2
6593  hex->face(3)->child(0)->line(
6594  (hex->face(3)->refinement_case() ==
6596  1 :
6597  3) // 3
6598  };
6599 
6600  unsigned int line_indices[4];
6601  for (unsigned int i = 0; i < 4; ++i)
6602  line_indices[i] = lines[i]->index();
6603 
6604  // the orientation of lines for the inner quads
6605  // is quite tricky. as these lines are newly
6606  // created ones and thus have no parents, they
6607  // cannot inherit this property. set up an array
6608  // and fill it with the respective values
6609  bool line_orientation[4];
6610 
6611  // the middle vertex marked as m0 above is the
6612  // start vertex for lines 0 and 2 in standard
6613  // orientation, whereas m1 is the end vertex of
6614  // lines 1 and 3 in standard orientation
6615  const unsigned int middle_vertices[2] = {
6616  middle_vertex_index<dim, spacedim>(hex->line(8)),
6617  middle_vertex_index<dim, spacedim>(hex->line(11))};
6618 
6619  for (unsigned int i = 0; i < 4; ++i)
6620  if (lines[i]->vertex_index(i % 2) ==
6621  middle_vertices[i % 2])
6622  line_orientation[i] = true;
6623  else
6624  {
6625  // it must be the other way round then
6626  Assert(lines[i]->vertex_index((i + 1) % 2) ==
6627  middle_vertices[i % 2],
<