Reference documentation for deal.II version Git 9297d75edf 2020-11-26 18:52:14 +0100
\(\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. It 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  tria_level.global_active_cell_indices.reserve(total_cells);
1436  tria_level.global_active_cell_indices.insert(
1437  tria_level.global_active_cell_indices.end(),
1438  total_cells - tria_level.global_active_cell_indices.size(),
1440 
1441  tria_level.global_level_cell_indices.reserve(total_cells);
1442  tria_level.global_level_cell_indices.insert(
1443  tria_level.global_level_cell_indices.end(),
1444  total_cells - tria_level.global_level_cell_indices.size(),
1446 
1447  if (dimension < space_dimension)
1448  {
1449  tria_level.direction_flags.reserve(total_cells);
1450  tria_level.direction_flags.insert(
1451  tria_level.direction_flags.end(),
1452  total_cells - tria_level.direction_flags.size(),
1453  true);
1454  }
1455  else
1456  tria_level.direction_flags.clear();
1457 
1458  tria_level.parents.reserve((total_cells + 1) / 2);
1459  tria_level.parents.insert(tria_level.parents.end(),
1460  (total_cells + 1) / 2 -
1461  tria_level.parents.size(),
1462  -1);
1463 
1464  tria_level.neighbors.reserve(total_cells * (2 * dimension));
1465  tria_level.neighbors.insert(tria_level.neighbors.end(),
1466  total_cells * (2 * dimension) -
1467  tria_level.neighbors.size(),
1468  std::make_pair(-1, -1));
1469 
1470 
1471  if (tria_level.dim == 3)
1472  {
1473  tria_level.face_orientations.reserve(
1474  total_cells * GeometryInfo<3>::faces_per_cell);
1475  tria_level.face_orientations.insert(
1476  tria_level.face_orientations.end(),
1477  total_cells * GeometryInfo<3>::faces_per_cell -
1478  tria_level.face_orientations.size(),
1479  true);
1480  }
1481 
1482  if (tria_level.dim == 2 || tria_level.dim == 3)
1483  {
1484  tria_level.reference_cell_type.reserve(total_cells);
1485  tria_level.reference_cell_type.insert(
1486  tria_level.reference_cell_type.end(),
1487  total_cells - tria_level.reference_cell_type.size(),
1488  tria_level.dim == 2 ? ReferenceCell::Type::Quad :
1490  }
1491  }
1492  }
1493 
1494 
1495 
1500  int,
1501  int,
1502  << "The containers have sizes " << arg1 << " and " << arg2
1503  << ", which is not as expected.");
1504 
1510  void
1511  monitor_memory(const TriaLevel & tria_level,
1512  const unsigned int true_dimension)
1513  {
1514  (void)tria_level;
1515  (void)true_dimension;
1516  Assert(2 * true_dimension * tria_level.refine_flags.size() ==
1517  tria_level.neighbors.size(),
1518  ExcMemoryInexact(tria_level.refine_flags.size(),
1519  tria_level.neighbors.size()));
1520  Assert(2 * true_dimension * tria_level.coarsen_flags.size() ==
1521  tria_level.neighbors.size(),
1522  ExcMemoryInexact(tria_level.coarsen_flags.size(),
1523  tria_level.neighbors.size()));
1524  }
1525 
1526 
1527 
1540  void
1541  reserve_space(TriaObjects & tria_objects,
1542  const unsigned int new_objects_in_pairs,
1543  const unsigned int new_objects_single = 0)
1544  {
1545  if (tria_objects.structdim <= 2)
1546  {
1547  Assert(new_objects_in_pairs % 2 == 0, ExcInternalError());
1548 
1549  tria_objects.next_free_single = 0;
1550  tria_objects.next_free_pair = 0;
1551  tria_objects.reverse_order_next_free_single = false;
1552 
1553  // count the number of objects, of unused single objects and of
1554  // unused pairs of objects
1555  unsigned int n_objects = 0;
1556  unsigned int n_unused_pairs = 0;
1557  unsigned int n_unused_singles = 0;
1558  for (unsigned int i = 0; i < tria_objects.used.size(); ++i)
1559  {
1560  if (tria_objects.used[i])
1561  ++n_objects;
1562  else if (i + 1 < tria_objects.used.size())
1563  {
1564  if (tria_objects.used[i + 1])
1565  {
1566  ++n_unused_singles;
1567  if (tria_objects.next_free_single == 0)
1568  tria_objects.next_free_single = i;
1569  }
1570  else
1571  {
1572  ++n_unused_pairs;
1573  if (tria_objects.next_free_pair == 0)
1574  tria_objects.next_free_pair = i;
1575  ++i;
1576  }
1577  }
1578  else
1579  ++n_unused_singles;
1580  }
1581  Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
1582  tria_objects.used.size(),
1583  ExcInternalError());
1584 
1585  // how many single objects are needed in addition to
1586  // n_unused_objects?
1587  const int additional_single_objects =
1588  new_objects_single - n_unused_singles;
1589 
1590  unsigned int new_size = tria_objects.used.size() +
1591  new_objects_in_pairs - 2 * n_unused_pairs;
1592  if (additional_single_objects > 0)
1593  new_size += additional_single_objects;
1594 
1595  // only allocate space if necessary
1596  if (new_size > tria_objects.n_objects())
1597  {
1598  unsigned int faces_per_cell = 1;
1599  unsigned int max_children_per_cell = 1;
1600 
1601  if (tria_objects.structdim == 1)
1602  faces_per_cell = GeometryInfo<1>::faces_per_cell;
1603  else if (tria_objects.structdim == 2)
1604  faces_per_cell = GeometryInfo<2>::faces_per_cell;
1605  else if (tria_objects.structdim == 3)
1606  faces_per_cell = GeometryInfo<3>::faces_per_cell;
1607  else
1608  AssertThrow(false, ExcNotImplemented());
1609 
1610  if (tria_objects.structdim == 1)
1611  max_children_per_cell = GeometryInfo<1>::max_children_per_cell;
1612  else if (tria_objects.structdim == 2)
1613  max_children_per_cell = GeometryInfo<2>::max_children_per_cell;
1614  else if (tria_objects.structdim == 3)
1615  max_children_per_cell = GeometryInfo<3>::max_children_per_cell;
1616  else
1617  AssertThrow(false, ExcNotImplemented());
1618 
1619  tria_objects.cells.reserve(new_size * faces_per_cell);
1620  tria_objects.cells.insert(tria_objects.cells.end(),
1621  (new_size - tria_objects.n_objects()) *
1622  faces_per_cell,
1623  -1);
1624 
1625  tria_objects.used.reserve(new_size);
1626  tria_objects.used.insert(tria_objects.used.end(),
1627  new_size - tria_objects.used.size(),
1628  false);
1629 
1630  tria_objects.user_flags.reserve(new_size);
1631  tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1632  new_size -
1633  tria_objects.user_flags.size(),
1634  false);
1635 
1636  const unsigned int factor = max_children_per_cell / 2;
1637  tria_objects.children.reserve(factor * new_size);
1638  tria_objects.children.insert(tria_objects.children.end(),
1639  factor * new_size -
1640  tria_objects.children.size(),
1641  -1);
1642 
1643  if (tria_objects.structdim > 1)
1644  {
1645  tria_objects.refinement_cases.reserve(new_size);
1646  tria_objects.refinement_cases.insert(
1647  tria_objects.refinement_cases.end(),
1648  new_size - tria_objects.refinement_cases.size(),
1649  /*RefinementCase::no_refinement=*/0);
1650  }
1651 
1652  // first reserve, then resize. Otherwise the std library can
1653  // decide to allocate more entries.
1654  tria_objects.boundary_or_material_id.reserve(new_size);
1655  tria_objects.boundary_or_material_id.resize(new_size);
1656 
1657  tria_objects.user_data.reserve(new_size);
1658  tria_objects.user_data.resize(new_size);
1659 
1660  tria_objects.manifold_id.reserve(new_size);
1661  tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1662  new_size -
1663  tria_objects.manifold_id.size(),
1665  }
1666 
1667  if (n_unused_singles == 0)
1668  {
1669  tria_objects.next_free_single = new_size - 1;
1670  tria_objects.reverse_order_next_free_single = true;
1671  }
1672  }
1673  else
1674  {
1675  const unsigned int new_hexes = new_objects_in_pairs;
1676 
1677  const unsigned int new_size =
1678  new_hexes + std::count(tria_objects.used.begin(),
1679  tria_objects.used.end(),
1680  true);
1681 
1682  // see above...
1683  if (new_size > tria_objects.n_objects())
1684  {
1685  unsigned int faces_per_cell = 1;
1686 
1687  if (tria_objects.structdim == 1)
1688  faces_per_cell = GeometryInfo<1>::faces_per_cell;
1689  else if (tria_objects.structdim == 2)
1690  faces_per_cell = GeometryInfo<2>::faces_per_cell;
1691  else if (tria_objects.structdim == 3)
1692  faces_per_cell = GeometryInfo<3>::faces_per_cell;
1693  else
1694  AssertThrow(false, ExcNotImplemented());
1695 
1696  tria_objects.cells.reserve(new_size * faces_per_cell);
1697  tria_objects.cells.insert(tria_objects.cells.end(),
1698  (new_size - tria_objects.n_objects()) *
1699  faces_per_cell,
1700  -1);
1701 
1702  tria_objects.used.reserve(new_size);
1703  tria_objects.used.insert(tria_objects.used.end(),
1704  new_size - tria_objects.used.size(),
1705  false);
1706 
1707  tria_objects.user_flags.reserve(new_size);
1708  tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1709  new_size -
1710  tria_objects.user_flags.size(),
1711  false);
1712 
1713  tria_objects.children.reserve(4 * new_size);
1714  tria_objects.children.insert(tria_objects.children.end(),
1715  4 * new_size -
1716  tria_objects.children.size(),
1717  -1);
1718 
1719  // for the following fields, we know exactly how many elements
1720  // we need, so first reserve then resize (resize itself, at least
1721  // with some compiler libraries, appears to round up the size it
1722  // actually reserves)
1723  tria_objects.boundary_or_material_id.reserve(new_size);
1724  tria_objects.boundary_or_material_id.resize(new_size);
1725 
1726  tria_objects.manifold_id.reserve(new_size);
1727  tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1728  new_size -
1729  tria_objects.manifold_id.size(),
1731 
1732  tria_objects.user_data.reserve(new_size);
1733  tria_objects.user_data.resize(new_size);
1734 
1735  tria_objects.refinement_cases.reserve(new_size);
1736  tria_objects.refinement_cases.insert(
1737  tria_objects.refinement_cases.end(),
1738  new_size - tria_objects.refinement_cases.size(),
1739  /*RefinementCase::no_refinement=*/0);
1740  }
1741  tria_objects.next_free_single = tria_objects.next_free_pair = 0;
1742  }
1743  }
1744 
1745 
1746 
1752  void
1753  monitor_memory(const TriaObjects &tria_object, const unsigned int)
1754  {
1755  Assert(tria_object.n_objects() == tria_object.used.size(),
1756  ExcMemoryInexact(tria_object.n_objects(),
1757  tria_object.used.size()));
1758  Assert(tria_object.n_objects() == tria_object.user_flags.size(),
1759  ExcMemoryInexact(tria_object.n_objects(),
1760  tria_object.user_flags.size()));
1761  Assert(tria_object.n_objects() ==
1762  tria_object.boundary_or_material_id.size(),
1763  ExcMemoryInexact(tria_object.n_objects(),
1764  tria_object.boundary_or_material_id.size()));
1765  Assert(tria_object.n_objects() == tria_object.manifold_id.size(),
1766  ExcMemoryInexact(tria_object.n_objects(),
1767  tria_object.manifold_id.size()));
1768  Assert(tria_object.n_objects() == tria_object.user_data.size(),
1769  ExcMemoryInexact(tria_object.n_objects(),
1770  tria_object.user_data.size()));
1771 
1772  if (tria_object.structdim == 1)
1773  {
1774  Assert(1 * tria_object.n_objects() == tria_object.children.size(),
1775  ExcMemoryInexact(tria_object.n_objects(),
1776  tria_object.children.size()));
1777  }
1778  else if (tria_object.structdim == 2)
1779  {
1780  Assert(2 * tria_object.n_objects() == tria_object.children.size(),
1781  ExcMemoryInexact(tria_object.n_objects(),
1782  tria_object.children.size()));
1783  }
1784  else if (tria_object.structdim == 3)
1785  {
1786  Assert(4 * tria_object.n_objects() == tria_object.children.size(),
1787  ExcMemoryInexact(tria_object.n_objects(),
1788  tria_object.children.size()));
1789  }
1790  }
1791 
1792 
1793 
1798  template <int dim, int spacedim>
1799  class Policy
1800  {
1801  public:
1805  virtual ~Policy() = default;
1806 
1810  virtual void
1811  delete_children(
1814  std::vector<unsigned int> & line_cell_count,
1815  std::vector<unsigned int> &quad_cell_count) = 0;
1816 
1821  execute_refinement(Triangulation<dim, spacedim> &triangulation,
1822  const bool check_for_distorted_cells) = 0;
1823 
1827  virtual void
1828  prevent_distorted_boundary_cells(
1829  Triangulation<dim, spacedim> &triangulation) = 0;
1830 
1834  virtual void
1835  prepare_refinement_dim_dependent(
1836  Triangulation<dim, spacedim> &triangulation) = 0;
1837 
1841  virtual bool
1842  coarsening_allowed(
1843  const typename Triangulation<dim, spacedim>::cell_iterator &cell) = 0;
1844 
1851  virtual std::unique_ptr<Policy<dim, spacedim>>
1852  clone() = 0;
1853  };
1854 
1855 
1856 
1862  template <int dim, int spacedim, typename T>
1863  class PolicyWrapper : public Policy<dim, spacedim>
1864  {
1865  public:
1866  void
1870  std::vector<unsigned int> & line_cell_count,
1871  std::vector<unsigned int> &quad_cell_count) override
1872  {
1873  T::delete_children(tria, cell, line_cell_count, quad_cell_count);
1874  }
1875 
1878  const bool check_for_distorted_cells) override
1879  {
1880  return T::execute_refinement(triangulation, check_for_distorted_cells);
1881  }
1882 
1883  void
1886  {
1887  T::prevent_distorted_boundary_cells(triangulation);
1888  }
1889 
1890  void
1893  {
1894  T::prepare_refinement_dim_dependent(triangulation);
1895  }
1896 
1897  bool
1899  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
1900  override
1901  {
1902  return T::template coarsening_allowed<dim, spacedim>(cell);
1903  }
1904 
1905  std::unique_ptr<Policy<dim, spacedim>>
1906  clone() override
1907  {
1908  return std::unique_ptr<Policy<dim, spacedim>>(
1910  }
1911  };
1912 
2009  {
2021  template <int dim, int spacedim>
2022  static void
2025  const unsigned int level_objects,
2027  {
2028  using line_iterator =
2030 
2031  number_cache.n_levels = 0;
2032  if (level_objects > 0)
2033  // find the last level on which there are used cells
2034  for (unsigned int level = 0; level < level_objects; ++level)
2035  if (triangulation.begin(level) != triangulation.end(level))
2036  number_cache.n_levels = level + 1;
2037 
2038  // no cells at all?
2039  Assert(number_cache.n_levels > 0, ExcInternalError());
2040 
2042  // update the number of lines on the different levels in the
2043  // cache
2044  number_cache.n_lines = 0;
2045  number_cache.n_active_lines = 0;
2046 
2047  // for 1d, lines have levels so take count the objects per
2048  // level and globally
2049  if (dim == 1)
2050  {
2051  number_cache.n_lines_level.resize(number_cache.n_levels);
2052  number_cache.n_active_lines_level.resize(number_cache.n_levels);
2053 
2054  for (unsigned int level = 0; level < number_cache.n_levels; ++level)
2055  {
2056  // count lines on this level
2057  number_cache.n_lines_level[level] = 0;
2058  number_cache.n_active_lines_level[level] = 0;
2059 
2060  line_iterator line = triangulation.begin_line(level),
2061  endc =
2062  (level == number_cache.n_levels - 1 ?
2063  line_iterator(triangulation.end_line()) :
2064  triangulation.begin_line(level + 1));
2065  for (; line != endc; ++line)
2066  {
2067  ++number_cache.n_lines_level[level];
2068  if (line->has_children() == false)
2069  ++number_cache.n_active_lines_level[level];
2070  }
2071 
2072  // update total number of lines
2073  number_cache.n_lines += number_cache.n_lines_level[level];
2074  number_cache.n_active_lines +=
2075  number_cache.n_active_lines_level[level];
2076  }
2077  }
2078  else
2079  {
2080  // for dim>1, there are no levels for lines
2081  number_cache.n_lines_level.clear();
2082  number_cache.n_active_lines_level.clear();
2083 
2084  line_iterator line = triangulation.begin_line(),
2085  endc = triangulation.end_line();
2086  for (; line != endc; ++line)
2087  {
2088  ++number_cache.n_lines;
2089  if (line->has_children() == false)
2090  ++number_cache.n_active_lines;
2091  }
2092  }
2093  }
2094 
2109  template <int dim, int spacedim>
2110  static void
2113  const unsigned int level_objects,
2115  {
2116  // update lines and n_levels in number_cache. since we don't
2117  // access any of these numbers, we can do this in the
2118  // background
2119  Threads::Task<void> update_lines = Threads::new_task(
2120  static_cast<
2121  void (*)(const Triangulation<dim, spacedim> &,
2122  const unsigned int,
2124  &compute_number_cache<dim, spacedim>),
2125  triangulation,
2126  level_objects,
2128  number_cache));
2129 
2130  using quad_iterator =
2132 
2134  // update the number of quads on the different levels in the
2135  // cache
2136  number_cache.n_quads = 0;
2137  number_cache.n_active_quads = 0;
2138 
2139  // for 2d, quads have levels so take count the objects per
2140  // level and globally
2141  if (dim == 2)
2142  {
2143  // count the number of levels; the function we called above
2144  // on a separate Task for lines also does this and puts it into
2145  // number_cache.n_levels, but this datum may not yet be
2146  // available as we call the function on a separate task
2147  unsigned int n_levels = 0;
2148  if (level_objects > 0)
2149  // find the last level on which there are used cells
2150  for (unsigned int level = 0; level < level_objects; ++level)
2151  if (triangulation.begin(level) != triangulation.end(level))
2152  n_levels = level + 1;
2153 
2154  number_cache.n_quads_level.resize(n_levels);
2155  number_cache.n_active_quads_level.resize(n_levels);
2156 
2157  for (unsigned int level = 0; level < n_levels; ++level)
2158  {
2159  // count quads on this level
2160  number_cache.n_quads_level[level] = 0;
2161  number_cache.n_active_quads_level[level] = 0;
2162 
2163  quad_iterator quad = triangulation.begin_quad(level),
2164  endc =
2165  (level == n_levels - 1 ?
2166  quad_iterator(triangulation.end_quad()) :
2167  triangulation.begin_quad(level + 1));
2168  for (; quad != endc; ++quad)
2169  {
2170  ++number_cache.n_quads_level[level];
2171  if (quad->has_children() == false)
2172  ++number_cache.n_active_quads_level[level];
2173  }
2174 
2175  // update total number of quads
2176  number_cache.n_quads += number_cache.n_quads_level[level];
2177  number_cache.n_active_quads +=
2178  number_cache.n_active_quads_level[level];
2179  }
2180  }
2181  else
2182  {
2183  // for dim>2, there are no levels for quads
2184  number_cache.n_quads_level.clear();
2185  number_cache.n_active_quads_level.clear();
2186 
2187  quad_iterator quad = triangulation.begin_quad(),
2188  endc = triangulation.end_quad();
2189  for (; quad != endc; ++quad)
2190  {
2191  ++number_cache.n_quads;
2192  if (quad->has_children() == false)
2193  ++number_cache.n_active_quads;
2194  }
2195  }
2196 
2197  // wait for the background computation for lines
2198  update_lines.join();
2199  }
2200 
2216  template <int dim, int spacedim>
2217  static void
2220  const unsigned int level_objects,
2222  {
2223  // update quads, lines and n_levels in number_cache. since we
2224  // don't access any of these numbers, we can do this in the
2225  // background
2226  Threads::Task<void> update_quads_and_lines = Threads::new_task(
2227  static_cast<
2228  void (*)(const Triangulation<dim, spacedim> &,
2229  const unsigned int,
2231  &compute_number_cache<dim, spacedim>),
2232  triangulation,
2233  level_objects,
2235  number_cache));
2236 
2237  using hex_iterator =
2239 
2241  // update the number of hexes on the different levels in the
2242  // cache
2243  number_cache.n_hexes = 0;
2244  number_cache.n_active_hexes = 0;
2245 
2246  // for 3d, hexes have levels so take count the objects per
2247  // level and globally
2248  if (dim == 3)
2249  {
2250  // count the number of levels; the function we called
2251  // above on a separate Task for quads (recursively, via
2252  // the lines function) also does this and puts it into
2253  // number_cache.n_levels, but this datum may not yet be
2254  // available as we call the function on a separate task
2255  unsigned int n_levels = 0;
2256  if (level_objects > 0)
2257  // find the last level on which there are used cells
2258  for (unsigned int level = 0; level < level_objects; ++level)
2259  if (triangulation.begin(level) != triangulation.end(level))
2260  n_levels = level + 1;
2261 
2262  number_cache.n_hexes_level.resize(n_levels);
2263  number_cache.n_active_hexes_level.resize(n_levels);
2264 
2265  for (unsigned int level = 0; level < n_levels; ++level)
2266  {
2267  // count hexes on this level
2268  number_cache.n_hexes_level[level] = 0;
2269  number_cache.n_active_hexes_level[level] = 0;
2270 
2271  hex_iterator hex = triangulation.begin_hex(level),
2272  endc = (level == n_levels - 1 ?
2273  hex_iterator(triangulation.end_hex()) :
2274  triangulation.begin_hex(level + 1));
2275  for (; hex != endc; ++hex)
2276  {
2277  ++number_cache.n_hexes_level[level];
2278  if (hex->has_children() == false)
2279  ++number_cache.n_active_hexes_level[level];
2280  }
2281 
2282  // update total number of hexes
2283  number_cache.n_hexes += number_cache.n_hexes_level[level];
2284  number_cache.n_active_hexes +=
2285  number_cache.n_active_hexes_level[level];
2286  }
2287  }
2288  else
2289  {
2290  // for dim>3, there are no levels for hexes
2291  number_cache.n_hexes_level.clear();
2292  number_cache.n_active_hexes_level.clear();
2293 
2294  hex_iterator hex = triangulation.begin_hex(),
2295  endc = triangulation.end_hex();
2296  for (; hex != endc; ++hex)
2297  {
2298  ++number_cache.n_hexes;
2299  if (hex->has_children() == false)
2300  ++number_cache.n_active_hexes;
2301  }
2302  }
2303 
2304  // wait for the background computation for quads
2305  update_quads_and_lines.join();
2306  }
2307 
2308 
2312  template <int dim, int spacedim>
2313  static void
2315  const std::vector<CellData<dim>> & cells,
2316  const SubCellData & subcelldata,
2318  {
2319  AssertThrow(vertices.size() > 0, ExcMessage("No vertices given"));
2320  AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
2321 
2322  // Check that all cells have positive volume.
2323 #ifndef _MSC_VER
2324  // TODO: The following code does not compile with MSVC. Find a way
2325  // around it
2326  if (dim == spacedim)
2327  for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2328  {
2329  // If we should check for distorted cells, then we permit them
2330  // to exist. If a cell has negative measure, then it must be
2331  // distorted (the converse is not necessarily true); hence
2332  // throw an exception if no such cells should exist.
2333  if (tria.check_for_distorted_cells)
2334  {
2335  const double cell_measure = GridTools::cell_measure<spacedim>(
2336  vertices,
2337  ArrayView<const unsigned int>(cells[cell_no].vertices));
2338  AssertThrow(cell_measure > 0, ExcGridHasInvalidCell(cell_no));
2339  }
2340  }
2341 #endif
2342 
2343  // clear old content
2344  tria.levels.clear();
2345  tria.levels.push_back(
2346  std::make_unique<
2348 
2349  if (dim > 1)
2350  tria.faces = std::make_unique<
2352 
2353  // copy vertices
2354  tria.vertices = vertices;
2355  tria.vertices_used.assign(vertices.size(), true);
2356 
2357  // compute connectivity
2358  const auto connectivity = build_connectivity<unsigned int>(cells);
2359  const unsigned int n_cell = cells.size();
2360 
2361  // TriaObjects: lines
2362  if (dim >= 2)
2363  {
2364  auto &lines_0 = tria.faces->lines; // data structure to be filled
2365 
2366  // get connectivity between quads and lines
2367  const auto & crs = connectivity.entity_to_entities(1, 0);
2368  const unsigned int n_lines = crs.ptr.size() - 1;
2369 
2370  // allocate memory
2371  reserve_space_(lines_0, n_lines);
2372 
2373  // loop over lines
2374  for (unsigned int line = 0; line < n_lines; ++line)
2375  for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
2376  ++i, ++j)
2377  lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
2378  crs.col[i]; // set vertex indices
2379  }
2380 
2381  // TriaObjects: quads
2382  if (dim == 3)
2383  {
2384  auto &quads_0 = tria.faces->quads; // data structures to be filled
2385  auto &faces = *tria.faces;
2386 
2387  // get connectivity between quads and lines
2388  const auto & crs = connectivity.entity_to_entities(2, 1);
2389  const unsigned int n_quads = crs.ptr.size() - 1;
2390 
2391  // allocate memory
2392  reserve_space_(quads_0, n_quads);
2393  reserve_space_(faces, 2 /*structdim*/, n_quads);
2394 
2395  // loop over all quads -> entity type, line indices/orientations
2396  for (unsigned int q = 0, k = 0; q < n_quads; ++q)
2397  {
2398  // set entity type of quads
2399  faces.quad_reference_cell_type[q] =
2400  connectivity.entity_types(2)[q];
2401 
2402  // loop over all its lines
2403  for (unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
2404  ++i, ++j, ++k)
2405  {
2406  // set line index
2407  quads_0.cells[q * GeometryInfo<2>::faces_per_cell + j] =
2408  crs.col[i];
2409 
2410  // set line orientations
2411  faces.quads_line_orientations
2413  connectivity.entity_orientations(1)[k];
2414  }
2415  }
2416  }
2417 
2418  // TriaObjects/TriaLevel: cell
2419  {
2420  auto &cells_0 = tria.levels[0]->cells; // data structure to be filled
2421  auto &level = *tria.levels[0];
2422 
2423  // get connectivity between cells/faces and cells/cells
2424  const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2425  const auto &nei = connectivity.entity_to_entities(dim, dim);
2426 
2427  // in 2D optional: since in in pure QUAD meshes same line
2428  // orientations can be guaranteed
2429  const bool orientation_needed =
2430  dim == 3 ||
2431  (dim == 2 &&
2432  std::any_of(connectivity.entity_orientations(1).begin(),
2433  connectivity.entity_orientations(1).end(),
2434  [](const auto &i) { return i == 0; }));
2435 
2436  // allocate memory
2437  reserve_space_(cells_0, n_cell);
2438  reserve_space_(level, spacedim, n_cell, orientation_needed);
2439 
2440  // loop over all cells
2441  for (unsigned int cell = 0; cell < n_cell; ++cell)
2442  {
2443  // set material ids
2444  cells_0.boundary_or_material_id[cell].material_id =
2445  cells[cell].material_id;
2446 
2447  // set manifold ids
2448  cells_0.manifold_id[cell] = cells[cell].manifold_id;
2449 
2450  // set entity types
2451  level.reference_cell_type[cell] =
2452  connectivity.entity_types(dim)[cell];
2453 
2454  // loop over faces
2455  for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2456  ++i, ++j)
2457  {
2458  // set neighbor if not at boundary
2459  if (nei.col[i] != static_cast<unsigned int>(-1))
2460  level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
2461  j] = {0, nei.col[i]};
2462 
2463  // set face indices
2464  cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
2465  crs.col[i];
2466 
2467  // set face orientation if needed
2468  if (orientation_needed)
2469  level.face_orientations
2470  [cell * GeometryInfo<dim>::faces_per_cell + j] =
2471  connectivity.entity_orientations(dim - 1)[i];
2472  }
2473  }
2474  }
2475 
2476  // TriaFaces: boundary id of boundary faces
2477  if (dim > 1)
2478  {
2479  auto &bids_face = dim == 3 ?
2480  tria.faces->quads.boundary_or_material_id :
2481  tria.faces->lines.boundary_or_material_id;
2482 
2483  // count number of cells a face is belonging to
2484  std::vector<unsigned int> count(bids_face.size(), 0);
2485 
2486  // get connectivity between cells/faces
2487  const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2488 
2489  // count how many cells are adjacent to the same face
2490  for (unsigned int cell = 0; cell < cells.size(); ++cell)
2491  for (unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
2492  count[crs.col[i]]++;
2493 
2494  // loop over all faces
2495  for (unsigned int face = 0; face < count.size(); ++face)
2496  {
2497  if (count[face] != 1) // inner face
2498  continue;
2499 
2500  // boundary faces ...
2501  bids_face[face].boundary_id = 0;
2502 
2503  if (dim != 3)
2504  continue;
2505 
2506  // ... and the lines of quads in 3D
2507  const auto &crs = connectivity.entity_to_entities(2, 1);
2508  for (unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
2509  tria.faces->lines.boundary_or_material_id[crs.col[i]]
2510  .boundary_id = 0;
2511  }
2512  }
2513  else // 1D
2514  {
2515  static const unsigned int t_tba = static_cast<unsigned int>(-1);
2516  static const unsigned int t_inner = static_cast<unsigned int>(-2);
2517 
2518  std::vector<unsigned int> type(vertices.size(), t_tba);
2519 
2520  const auto &crs = connectivity.entity_to_entities(1, 0);
2521 
2522  for (unsigned int cell = 0; cell < cells.size(); ++cell)
2523  for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2524  ++i, ++j)
2525  if (type[crs.col[i]] != t_inner)
2526  type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
2527 
2528  for (unsigned int face = 0; face < type.size(); ++face)
2529  {
2530  // note: we also treat manifolds here!?
2531  (*tria.vertex_to_manifold_id_map_1d)[face] =
2533  if (type[face] != t_inner && type[face] != t_tba)
2534  (*tria.vertex_to_boundary_id_map_1d)[face] = type[face];
2535  }
2536  }
2537 
2538  // SubCellData: line
2539  if (dim >= 2)
2540  process_subcelldata(connectivity.entity_to_entities(1, 0),
2541  tria.faces->lines,
2542  subcelldata.boundary_lines);
2543 
2544  // SubCellData: quad
2545  if (dim == 3)
2546  process_subcelldata(connectivity.entity_to_entities(2, 0),
2547  tria.faces->quads,
2548  subcelldata.boundary_quads);
2549  }
2550 
2551 
2552  template <int structdim, typename T>
2553  static void
2555  const CRS<T> & crs,
2556  TriaObjects & obj,
2557  const std::vector<CellData<structdim>> &boundary_objects_in)
2558  {
2559  AssertDimension(obj.structdim, structdim);
2560 
2561  if (boundary_objects_in.size() == 0)
2562  return; // empty subcelldata -> nothing to do
2563 
2564  // pre-sort subcelldata
2565  auto boundary_objects = boundary_objects_in;
2566 
2567  // ... sort vertices
2568  for (auto &boundary_object : boundary_objects)
2569  std::sort(boundary_object.vertices.begin(),
2570  boundary_object.vertices.end());
2571 
2572  // ... sort cells
2573  std::sort(boundary_objects.begin(),
2574  boundary_objects.end(),
2575  [](const auto &a, const auto &b) {
2576  return a.vertices < b.vertices;
2577  });
2578 
2579  unsigned int counter = 0;
2580 
2581  std::vector<unsigned int> key;
2583 
2584  for (unsigned int o = 0; o < obj.n_objects(); ++o)
2585  {
2586  auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
2587  auto &manifold_id = obj.manifold_id[o];
2588 
2589  // assert that object has not been visited yet and its value
2590  // has not been modified yet
2591  AssertThrow(boundary_id == 0 ||
2593  ExcNotImplemented());
2595  ExcNotImplemented());
2596 
2597  // create key
2598  key.assign(crs.col.data() + crs.ptr[o],
2599  crs.col.data() + crs.ptr[o + 1]);
2600  std::sort(key.begin(), key.end());
2601 
2602  // is subcelldata provided? -> binary search
2603  const auto subcell_object =
2604  std::lower_bound(boundary_objects.begin(),
2605  boundary_objects.end(),
2606  key,
2607  [&](const auto &cell, const auto &key) {
2608  return cell.vertices < key;
2609  });
2610 
2611  // no subcelldata provided for this object
2612  if (subcell_object == boundary_objects.end() ||
2613  subcell_object->vertices != key)
2614  continue;
2615 
2616  counter++;
2617 
2618  // set manifold id
2619  manifold_id = subcell_object->manifold_id;
2620 
2621  // set boundary id
2622  if (subcell_object->boundary_id !=
2624  {
2626  ExcNotImplemented());
2627  boundary_id = subcell_object->boundary_id;
2628  }
2629  }
2630 
2631  // make sure that all subcelldata entries have been processed
2632  // TODO: this is not guaranteed, why?
2633  // AssertDimension(counter, boundary_objects_in.size());
2634  }
2635 
2636 
2637 
2638  static void
2640  const unsigned structdim,
2641  const unsigned int size)
2642  {
2643  const unsigned int dim = faces.dim;
2644 
2645  const unsigned int faces_per_cell =
2646  structdim == 1 ? GeometryInfo<1>::faces_per_cell :
2647  (structdim == 2 ? GeometryInfo<2>::faces_per_cell :
2649 
2650  if (dim == 3 && structdim == 2)
2651  {
2652  // quad entity types
2653  faces.quad_reference_cell_type.assign(size,
2655 
2656  // quad line orientations
2657  faces.quads_line_orientations.assign(size * faces_per_cell, -1);
2658  }
2659  }
2660 
2661 
2662 
2663  static void
2665  const unsigned int spacedim,
2666  const unsigned int size,
2667  const bool orientation_needed)
2668  {
2669  const unsigned int dim = level.dim;
2670 
2671  const unsigned int faces_per_cell =
2672  dim == 1 ? GeometryInfo<1>::faces_per_cell :
2673  (dim == 2 ? GeometryInfo<2>::faces_per_cell :
2675 
2676  level.active_cell_indices.assign(size, -1);
2677  level.subdomain_ids.assign(size, 0);
2678  level.level_subdomain_ids.assign(size, 0);
2679 
2680  level.refine_flags.assign(size, false);
2681  level.coarsen_flags.assign(size, false);
2682 
2683  level.parents.assign((size + 1) / 2, -1);
2684 
2685  if (dim < spacedim)
2686  level.direction_flags.assign(size, true);
2687 
2688  level.neighbors.assign(size * faces_per_cell, {-1, -1});
2689 
2691 
2692  if (orientation_needed)
2693  level.face_orientations.assign(size * faces_per_cell, -1);
2694 
2695  level.global_active_cell_indices.assign(size,
2697  level.global_level_cell_indices.assign(size,
2699  }
2700 
2701 
2702 
2703  static void
2704  reserve_space_(TriaObjects &obj, const unsigned int size)
2705  {
2706  const unsigned int structdim = obj.structdim;
2707 
2708  const unsigned int max_children_per_cell =
2709  structdim == 1 ?
2711  (structdim == 2 ? GeometryInfo<2>::max_children_per_cell :
2713  const unsigned int faces_per_cell =
2714  structdim == 1 ? GeometryInfo<1>::faces_per_cell :
2715  (structdim == 2 ? GeometryInfo<2>::faces_per_cell :
2717 
2718  obj.used.assign(size, true);
2719  obj.boundary_or_material_id.assign(
2720  size,
2722  BoundaryOrMaterialId());
2723  obj.manifold_id.assign(size, -1);
2724  obj.user_flags.assign(size, false);
2725  obj.user_data.resize(size);
2726 
2727  if (structdim > 1) // TODO: why?
2728  obj.refinement_cases.assign(size, 0);
2729 
2730  obj.children.assign(max_children_per_cell / 2 * size, -1);
2731 
2732  obj.cells.assign(faces_per_cell * size, -1);
2733 
2734  if (structdim <= 2)
2735  {
2736  obj.next_free_single = size - 1;
2737  obj.next_free_pair = 0;
2738  obj.reverse_order_next_free_single = true;
2739  }
2740  else
2741  {
2742  obj.next_free_single = obj.next_free_pair = 0;
2743  }
2744  }
2745 
2746 
2762  template <int spacedim>
2763  static void delete_children(
2766  std::vector<unsigned int> &,
2767  std::vector<unsigned int> &)
2768  {
2769  const unsigned int dim = 1;
2770 
2771  // first we need to reset the
2772  // neighbor pointers of the
2773  // neighbors of this cell's
2774  // children to this cell. This is
2775  // different for one dimension,
2776  // since there neighbors can have a
2777  // refinement level differing from
2778  // that of this cell's children by
2779  // more than one level.
2780 
2781  Assert(!cell->child(0)->has_children() &&
2782  !cell->child(1)->has_children(),
2783  ExcInternalError());
2784 
2785  // first do it for the cells to the
2786  // left
2787  if (cell->neighbor(0).state() == IteratorState::valid)
2788  if (cell->neighbor(0)->has_children())
2789  {
2791  cell->neighbor(0);
2792  Assert(neighbor->level() == cell->level(), ExcInternalError());
2793 
2794  // right child
2795  neighbor = neighbor->child(1);
2796  while (true)
2797  {
2798  Assert(neighbor->neighbor(1) == cell->child(0),
2799  ExcInternalError());
2800  neighbor->set_neighbor(1, cell);
2801 
2802  // move on to further
2803  // children on the
2804  // boundary between this
2805  // cell and its neighbor
2806  if (neighbor->has_children())
2807  neighbor = neighbor->child(1);
2808  else
2809  break;
2810  }
2811  }
2812 
2813  // now do it for the cells to the
2814  // left
2815  if (cell->neighbor(1).state() == IteratorState::valid)
2816  if (cell->neighbor(1)->has_children())
2817  {
2819  cell->neighbor(1);
2820  Assert(neighbor->level() == cell->level(), ExcInternalError());
2821 
2822  // left child
2823  neighbor = neighbor->child(0);
2824  while (true)
2825  {
2826  Assert(neighbor->neighbor(0) == cell->child(1),
2827  ExcInternalError());
2828  neighbor->set_neighbor(0, cell);
2829 
2830  // move on to further
2831  // children on the
2832  // boundary between this
2833  // cell and its neighbor
2834  if (neighbor->has_children())
2835  neighbor = neighbor->child(0);
2836  else
2837  break;
2838  }
2839  }
2840 
2841 
2842  // delete the vertex which will not
2843  // be needed anymore. This vertex
2844  // is the second of the first child
2845  triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
2846 
2847  // invalidate children. clear user
2848  // pointers, to avoid that they may
2849  // appear at unwanted places later
2850  // on...
2851  for (unsigned int child = 0; child < cell->n_children(); ++child)
2852  {
2853  cell->child(child)->clear_user_data();
2854  cell->child(child)->clear_user_flag();
2855  cell->child(child)->clear_used_flag();
2856  }
2857 
2858 
2859  // delete pointer to children
2860  cell->clear_children();
2861  cell->clear_user_flag();
2862  }
2863 
2864 
2865 
2866  template <int spacedim>
2867  static void delete_children(
2870  std::vector<unsigned int> & line_cell_count,
2871  std::vector<unsigned int> &)
2872  {
2873  const unsigned int dim = 2;
2874  const RefinementCase<dim> ref_case = cell->refinement_case();
2875 
2876  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2877  ExcInternalError());
2878 
2879  // vectors to hold all lines which
2880  // may be deleted
2881  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2882  lines_to_delete(0);
2883 
2884  lines_to_delete.reserve(4 * 2 + 4);
2885 
2886  // now we decrease the counters for
2887  // lines contained in the child
2888  // cells
2889  for (unsigned int c = 0; c < cell->n_children(); ++c)
2890  {
2892  cell->child(c);
2893  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2894  --line_cell_count[child->line_index(l)];
2895  }
2896 
2897 
2898  // delete the vertex which will not
2899  // be needed anymore. This vertex
2900  // is the second of the second line
2901  // of the first child, if the cell
2902  // is refined with cut_xy, else there
2903  // is no inner vertex.
2904  // additionally delete unneeded inner
2905  // lines
2906  if (ref_case == RefinementCase<dim>::cut_xy)
2907  {
2908  triangulation
2909  .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
2910 
2911  lines_to_delete.push_back(cell->child(0)->line(1));
2912  lines_to_delete.push_back(cell->child(0)->line(3));
2913  lines_to_delete.push_back(cell->child(3)->line(0));
2914  lines_to_delete.push_back(cell->child(3)->line(2));
2915  }
2916  else
2917  {
2918  unsigned int inner_face_no =
2919  ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
2920 
2921  // the inner line will not be
2922  // used any more
2923  lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2924  }
2925 
2926  // invalidate children
2927  for (unsigned int child = 0; child < cell->n_children(); ++child)
2928  {
2929  cell->child(child)->clear_user_data();
2930  cell->child(child)->clear_user_flag();
2931  cell->child(child)->clear_used_flag();
2932  }
2933 
2934 
2935  // delete pointer to children
2936  cell->clear_children();
2937  cell->clear_refinement_case();
2938  cell->clear_user_flag();
2939 
2940  // look at the refinement of outer
2941  // lines. if nobody needs those
2942  // anymore we can add them to the
2943  // list of lines to be deleted.
2944  for (unsigned int line_no = 0;
2945  line_no < GeometryInfo<dim>::lines_per_cell;
2946  ++line_no)
2947  {
2949  cell->line(line_no);
2950 
2951  if (line->has_children())
2952  {
2953  // if one of the cell counters is
2954  // zero, the other has to be as well
2955 
2956  Assert((line_cell_count[line->child_index(0)] == 0 &&
2957  line_cell_count[line->child_index(1)] == 0) ||
2958  (line_cell_count[line->child_index(0)] > 0 &&
2959  line_cell_count[line->child_index(1)] > 0),
2960  ExcInternalError());
2961 
2962  if (line_cell_count[line->child_index(0)] == 0)
2963  {
2964  for (unsigned int c = 0; c < 2; ++c)
2965  Assert(!line->child(c)->has_children(),
2966  ExcInternalError());
2967 
2968  // we may delete the line's
2969  // children and the middle vertex
2970  // as no cell references them
2971  // anymore
2972  triangulation
2973  .vertices_used[line->child(0)->vertex_index(1)] = false;
2974 
2975  lines_to_delete.push_back(line->child(0));
2976  lines_to_delete.push_back(line->child(1));
2977 
2978  line->clear_children();
2979  }
2980  }
2981  }
2982 
2983  // finally, delete unneeded lines
2984 
2985  // clear user pointers, to avoid that
2986  // they may appear at unwanted places
2987  // later on...
2988  // same for user flags, then finally
2989  // delete the lines
2990  typename std::vector<
2992  line = lines_to_delete.begin(),
2993  endline = lines_to_delete.end();
2994  for (; line != endline; ++line)
2995  {
2996  (*line)->clear_user_data();
2997  (*line)->clear_user_flag();
2998  (*line)->clear_used_flag();
2999  }
3000  }
3001 
3002 
3003 
3004  template <int spacedim>
3005  static void delete_children(
3008  std::vector<unsigned int> & line_cell_count,
3009  std::vector<unsigned int> & quad_cell_count)
3010  {
3011  const unsigned int dim = 3;
3012 
3013  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
3014  ExcInternalError());
3015  Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
3016  ExcInternalError());
3017 
3018  // first of all, we store the RefineCase of
3019  // this cell
3020  const RefinementCase<dim> ref_case = cell->refinement_case();
3021  // vectors to hold all lines and quads which
3022  // may be deleted
3023  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3024  lines_to_delete(0);
3025  std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3026  quads_to_delete(0);
3027 
3028  lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3029  quads_to_delete.reserve(6 * 4 + 12);
3030 
3031  // now we decrease the counters for lines and
3032  // quads contained in the child cells
3033  for (unsigned int c = 0; c < cell->n_children(); ++c)
3034  {
3036  cell->child(c);
3037  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3038  --line_cell_count[child->line_index(l)];
3039  for (auto f : GeometryInfo<dim>::face_indices())
3040  --quad_cell_count[child->quad_index(f)];
3041  }
3042 
3044  // delete interior quads and lines and the
3045  // interior vertex, depending on the
3046  // refinement case of the cell
3047  //
3048  // for append quads and lines: only append
3049  // them to the list of objects to be deleted
3050 
3051  switch (ref_case)
3052  {
3054  quads_to_delete.push_back(cell->child(0)->face(1));
3055  break;
3057  quads_to_delete.push_back(cell->child(0)->face(3));
3058  break;
3060  quads_to_delete.push_back(cell->child(0)->face(5));
3061  break;
3063  quads_to_delete.push_back(cell->child(0)->face(1));
3064  quads_to_delete.push_back(cell->child(0)->face(3));
3065  quads_to_delete.push_back(cell->child(3)->face(0));
3066  quads_to_delete.push_back(cell->child(3)->face(2));
3067 
3068  lines_to_delete.push_back(cell->child(0)->line(11));
3069  break;
3071  quads_to_delete.push_back(cell->child(0)->face(1));
3072  quads_to_delete.push_back(cell->child(0)->face(5));
3073  quads_to_delete.push_back(cell->child(3)->face(0));
3074  quads_to_delete.push_back(cell->child(3)->face(4));
3075 
3076  lines_to_delete.push_back(cell->child(0)->line(5));
3077  break;
3079  quads_to_delete.push_back(cell->child(0)->face(3));
3080  quads_to_delete.push_back(cell->child(0)->face(5));
3081  quads_to_delete.push_back(cell->child(3)->face(2));
3082  quads_to_delete.push_back(cell->child(3)->face(4));
3083 
3084  lines_to_delete.push_back(cell->child(0)->line(7));
3085  break;
3087  quads_to_delete.push_back(cell->child(0)->face(1));
3088  quads_to_delete.push_back(cell->child(2)->face(1));
3089  quads_to_delete.push_back(cell->child(4)->face(1));
3090  quads_to_delete.push_back(cell->child(6)->face(1));
3091 
3092  quads_to_delete.push_back(cell->child(0)->face(3));
3093  quads_to_delete.push_back(cell->child(1)->face(3));
3094  quads_to_delete.push_back(cell->child(4)->face(3));
3095  quads_to_delete.push_back(cell->child(5)->face(3));
3096 
3097  quads_to_delete.push_back(cell->child(0)->face(5));
3098  quads_to_delete.push_back(cell->child(1)->face(5));
3099  quads_to_delete.push_back(cell->child(2)->face(5));
3100  quads_to_delete.push_back(cell->child(3)->face(5));
3101 
3102  lines_to_delete.push_back(cell->child(0)->line(5));
3103  lines_to_delete.push_back(cell->child(0)->line(7));
3104  lines_to_delete.push_back(cell->child(0)->line(11));
3105  lines_to_delete.push_back(cell->child(7)->line(0));
3106  lines_to_delete.push_back(cell->child(7)->line(2));
3107  lines_to_delete.push_back(cell->child(7)->line(8));
3108  // delete the vertex which will not
3109  // be needed anymore. This vertex
3110  // is the vertex at the heart of
3111  // this cell, which is the sixth of
3112  // the first child
3113  triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3114  false;
3115  break;
3116  default:
3117  // only remaining case is
3118  // no_refinement, thus an error
3119  Assert(false, ExcInternalError());
3120  break;
3121  }
3122 
3123 
3124  // invalidate children
3125  for (unsigned int child = 0; child < cell->n_children(); ++child)
3126  {
3127  cell->child(child)->clear_user_data();
3128  cell->child(child)->clear_user_flag();
3129 
3130  for (auto f : GeometryInfo<dim>::face_indices())
3131  {
3132  // set flags denoting deviations from
3133  // standard orientation of faces back
3134  // to initialization values
3135  cell->child(child)->set_face_orientation(f, true);
3136  cell->child(child)->set_face_flip(f, false);
3137  cell->child(child)->set_face_rotation(f, false);
3138  }
3139 
3140  cell->child(child)->clear_used_flag();
3141  }
3142 
3143 
3144  // delete pointer to children
3145  cell->clear_children();
3146  cell->clear_refinement_case();
3147  cell->clear_user_flag();
3148 
3149  // so far we only looked at inner quads,
3150  // lines and vertices. Now we have to
3151  // consider outer ones as well. here, we have
3152  // to check, whether there are other cells
3153  // still needing these objects. otherwise we
3154  // can delete them. first for quads (and
3155  // their inner lines).
3156 
3157  for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
3158  {
3160  cell->face(quad_no);
3161 
3162  Assert(
3163  (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
3164  quad->has_children()) ||
3165  GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
3166  RefinementCase<dim - 1>::no_refinement,
3167  ExcInternalError());
3168 
3169  switch (quad->refinement_case())
3170  {
3171  case RefinementCase<dim - 1>::no_refinement:
3172  // nothing to do as the quad
3173  // is not refined
3174  break;
3175  case RefinementCase<dim - 1>::cut_x:
3176  case RefinementCase<dim - 1>::cut_y:
3177  {
3178  // if one of the cell counters is
3179  // zero, the other has to be as
3180  // well
3181  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3182  quad_cell_count[quad->child_index(1)] == 0) ||
3183  (quad_cell_count[quad->child_index(0)] > 0 &&
3184  quad_cell_count[quad->child_index(1)] > 0),
3185  ExcInternalError());
3186  // it might be, that the quad is
3187  // refined twice anisotropically,
3188  // first check, whether we may
3189  // delete possible grand_children
3190  unsigned int deleted_grandchildren = 0;
3191  unsigned int number_of_child_refinements = 0;
3192 
3193  for (unsigned int c = 0; c < 2; ++c)
3194  if (quad->child(c)->has_children())
3195  {
3196  ++number_of_child_refinements;
3197  // if one of the cell counters is
3198  // zero, the other has to be as
3199  // well
3200  Assert(
3201  (quad_cell_count[quad->child(c)->child_index(0)] ==
3202  0 &&
3203  quad_cell_count[quad->child(c)->child_index(1)] ==
3204  0) ||
3205  (quad_cell_count[quad->child(c)->child_index(0)] >
3206  0 &&
3207  quad_cell_count[quad->child(c)->child_index(1)] >
3208  0),
3209  ExcInternalError());
3210  if (quad_cell_count[quad->child(c)->child_index(0)] ==
3211  0)
3212  {
3213  // Assert, that the two
3214  // anisotropic
3215  // refinements add up to
3216  // isotropic refinement
3217  Assert(quad->refinement_case() +
3218  quad->child(c)->refinement_case() ==
3220  ExcInternalError());
3221  // we may delete the
3222  // quad's children and
3223  // the inner line as no
3224  // cell references them
3225  // anymore
3226  quads_to_delete.push_back(
3227  quad->child(c)->child(0));
3228  quads_to_delete.push_back(
3229  quad->child(c)->child(1));
3230  if (quad->child(c)->refinement_case() ==
3232  lines_to_delete.push_back(
3233  quad->child(c)->child(0)->line(1));
3234  else
3235  lines_to_delete.push_back(
3236  quad->child(c)->child(0)->line(3));
3237  quad->child(c)->clear_children();
3238  quad->child(c)->clear_refinement_case();
3239  ++deleted_grandchildren;
3240  }
3241  }
3242  // if no grandchildren are left, we
3243  // may as well delete the
3244  // refinement of the inner line
3245  // between our children and the
3246  // corresponding vertex
3247  if (number_of_child_refinements > 0 &&
3248  deleted_grandchildren == number_of_child_refinements)
3249  {
3251  middle_line;
3252  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3253  middle_line = quad->child(0)->line(1);
3254  else
3255  middle_line = quad->child(0)->line(3);
3256 
3257  lines_to_delete.push_back(middle_line->child(0));
3258  lines_to_delete.push_back(middle_line->child(1));
3259  triangulation
3260  .vertices_used[middle_vertex_index<dim, spacedim>(
3261  middle_line)] = false;
3262  middle_line->clear_children();
3263  }
3264 
3265  // now consider the direct children
3266  // of the given quad
3267  if (quad_cell_count[quad->child_index(0)] == 0)
3268  {
3269  // we may delete the quad's
3270  // children and the inner line
3271  // as no cell references them
3272  // anymore
3273  quads_to_delete.push_back(quad->child(0));
3274  quads_to_delete.push_back(quad->child(1));
3275  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3276  lines_to_delete.push_back(quad->child(0)->line(1));
3277  else
3278  lines_to_delete.push_back(quad->child(0)->line(3));
3279 
3280  // if the counters just dropped
3281  // to zero, otherwise the
3282  // children would have been
3283  // deleted earlier, then this
3284  // cell's children must have
3285  // contained the anisotropic
3286  // quad children. thus, if
3287  // those have again anisotropic
3288  // children, which are in
3289  // effect isotropic children of
3290  // the original quad, those are
3291  // still needed by a
3292  // neighboring cell and we
3293  // cannot delete them. instead,
3294  // we have to reset this quad's
3295  // refine case to isotropic and
3296  // set the children
3297  // accordingly.
3298  if (quad->child(0)->has_children())
3299  if (quad->refinement_case() ==
3301  {
3302  // now evereything is
3303  // quite complicated. we
3304  // have the children
3305  // numbered according to
3306  //
3307  // *---*---*
3308  // |n+1|m+1|
3309  // *---*---*
3310  // | n | m |
3311  // *---*---*
3312  //
3313  // from the original
3314  // anisotropic
3315  // refinement. we have to
3316  // reorder them as
3317  //
3318  // *---*---*
3319  // | m |m+1|
3320  // *---*---*
3321  // | n |n+1|
3322  // *---*---*
3323  //
3324  // for isotropic refinement.
3325  //
3326  // this is a bit ugly, of
3327  // course: loop over all
3328  // cells on all levels
3329  // and look for faces n+1
3330  // (switch_1) and m
3331  // (switch_2).
3332  const typename Triangulation<dim, spacedim>::
3333  quad_iterator switch_1 =
3334  quad->child(0)->child(1),
3335  switch_2 =
3336  quad->child(1)->child(0);
3337 
3338  Assert(!switch_1->has_children(),
3339  ExcInternalError());
3340  Assert(!switch_2->has_children(),
3341  ExcInternalError());
3342 
3343  const int switch_1_index = switch_1->index();
3344  const int switch_2_index = switch_2->index();
3345  for (unsigned int l = 0;
3346  l < triangulation.levels.size();
3347  ++l)
3348  for (unsigned int h = 0;
3349  h <
3350  triangulation.levels[l]->cells.n_objects();
3351  ++h)
3352  for (const unsigned int q :
3354  {
3355  const int index =
3356  triangulation.levels[l]
3357  ->cells.get_bounding_object_indices(
3358  h)[q];
3359  if (index == switch_1_index)
3360  triangulation.levels[l]
3361  ->cells.get_bounding_object_indices(
3362  h)[q] = switch_2_index;
3363  else if (index == switch_2_index)
3364  triangulation.levels[l]
3365  ->cells.get_bounding_object_indices(
3366  h)[q] = switch_1_index;
3367  }
3368  // now we have to copy
3369  // all information of the
3370  // two quads
3371  const int switch_1_lines[4] = {
3372  static_cast<signed int>(
3373  switch_1->line_index(0)),
3374  static_cast<signed int>(
3375  switch_1->line_index(1)),
3376  static_cast<signed int>(
3377  switch_1->line_index(2)),
3378  static_cast<signed int>(
3379  switch_1->line_index(3))};
3380  const bool switch_1_line_orientations[4] = {
3381  switch_1->line_orientation(0),
3382  switch_1->line_orientation(1),
3383  switch_1->line_orientation(2),
3384  switch_1->line_orientation(3)};
3385  const types::boundary_id switch_1_boundary_id =
3386  switch_1->boundary_id();
3387  const unsigned int switch_1_user_index =
3388  switch_1->user_index();
3389  const bool switch_1_user_flag =
3390  switch_1->user_flag_set();
3391 
3392  switch_1->set_bounding_object_indices(
3393  {switch_2->line_index(0),
3394  switch_2->line_index(1),
3395  switch_2->line_index(2),
3396  switch_2->line_index(3)});
3397  switch_1->set_line_orientation(
3398  0, switch_2->line_orientation(0));
3399  switch_1->set_line_orientation(
3400  1, switch_2->line_orientation(1));
3401  switch_1->set_line_orientation(
3402  2, switch_2->line_orientation(2));
3403  switch_1->set_line_orientation(
3404  3, switch_2->line_orientation(3));
3405  switch_1->set_boundary_id_internal(
3406  switch_2->boundary_id());
3407  switch_1->set_manifold_id(
3408  switch_2->manifold_id());
3409  switch_1->set_user_index(switch_2->user_index());
3410  if (switch_2->user_flag_set())
3411  switch_1->set_user_flag();
3412  else
3413  switch_1->clear_user_flag();
3414 
3415  switch_2->set_bounding_object_indices(
3416  {switch_1_lines[0],
3417  switch_1_lines[1],
3418  switch_1_lines[2],
3419  switch_1_lines[3]});
3420  switch_2->set_line_orientation(
3421  0, switch_1_line_orientations[0]);
3422  switch_2->set_line_orientation(
3423  1, switch_1_line_orientations[1]);
3424  switch_2->set_line_orientation(
3425  2, switch_1_line_orientations[2]);
3426  switch_2->set_line_orientation(
3427  3, switch_1_line_orientations[3]);
3428  switch_2->set_boundary_id_internal(
3429  switch_1_boundary_id);
3430  switch_2->set_manifold_id(
3431  switch_1->manifold_id());
3432  switch_2->set_user_index(switch_1_user_index);
3433  if (switch_1_user_flag)
3434  switch_2->set_user_flag();
3435  else
3436  switch_2->clear_user_flag();
3437 
3438  const unsigned int child_0 =
3439  quad->child(0)->child_index(0);
3440  const unsigned int child_2 =
3441  quad->child(1)->child_index(0);
3442  quad->clear_children();
3443  quad->clear_refinement_case();
3444  quad->set_refinement_case(
3446  quad->set_children(0, child_0);
3447  quad->set_children(2, child_2);
3448  std::swap(quad_cell_count[child_0 + 1],
3449  quad_cell_count[child_2]);
3450  }
3451  else
3452  {
3453  // the face was refined
3454  // with cut_y, thus the
3455  // children are already
3456  // in correct order. we
3457  // only have to set them
3458  // correctly, deleting
3459  // the indirection of two
3460  // anisotropic refinement
3461  // and going directly
3462  // from the quad to
3463  // isotropic children
3464  const unsigned int child_0 =
3465  quad->child(0)->child_index(0);
3466  const unsigned int child_2 =
3467  quad->child(1)->child_index(0);
3468  quad->clear_children();
3469  quad->clear_refinement_case();
3470  quad->set_refinement_case(
3472  quad->set_children(0, child_0);
3473  quad->set_children(2, child_2);
3474  }
3475  else
3476  {
3477  quad->clear_children();
3478  quad->clear_refinement_case();
3479  }
3480  }
3481  break;
3482  }
3483  case RefinementCase<dim - 1>::cut_xy:
3484  {
3485  // if one of the cell counters is
3486  // zero, the others have to be as
3487  // well
3488 
3489  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3490  quad_cell_count[quad->child_index(1)] == 0 &&
3491  quad_cell_count[quad->child_index(2)] == 0 &&
3492  quad_cell_count[quad->child_index(3)] == 0) ||
3493  (quad_cell_count[quad->child_index(0)] > 0 &&
3494  quad_cell_count[quad->child_index(1)] > 0 &&
3495  quad_cell_count[quad->child_index(2)] > 0 &&
3496  quad_cell_count[quad->child_index(3)] > 0),
3497  ExcInternalError());
3498 
3499  if (quad_cell_count[quad->child_index(0)] == 0)
3500  {
3501  // we may delete the quad's
3502  // children, the inner lines
3503  // and the middle vertex as no
3504  // cell references them anymore
3505  lines_to_delete.push_back(quad->child(0)->line(1));
3506  lines_to_delete.push_back(quad->child(3)->line(0));
3507  lines_to_delete.push_back(quad->child(0)->line(3));
3508  lines_to_delete.push_back(quad->child(3)->line(2));
3509 
3510  for (unsigned int child = 0; child < quad->n_children();
3511  ++child)
3512  quads_to_delete.push_back(quad->child(child));
3513 
3514  triangulation
3515  .vertices_used[quad->child(0)->vertex_index(3)] =
3516  false;
3517 
3518  quad->clear_children();
3519  quad->clear_refinement_case();
3520  }
3521  }
3522  break;
3523 
3524  default:
3525  Assert(false, ExcInternalError());
3526  break;
3527  }
3528  }
3529 
3530  // now we repeat a similar procedure
3531  // for the outer lines of this cell.
3532 
3533  // if in debug mode: check that each
3534  // of the lines for which we consider
3535  // deleting the children in fact has
3536  // children (the bits/coarsening_3d
3537  // test tripped over this initially)
3538  for (unsigned int line_no = 0;
3539  line_no < GeometryInfo<dim>::lines_per_cell;
3540  ++line_no)
3541  {
3543  cell->line(line_no);
3544 
3545  Assert(
3546  (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
3547  line->has_children()) ||
3548  GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
3550  ExcInternalError());
3551 
3552  if (line->has_children())
3553  {
3554  // if one of the cell counters is
3555  // zero, the other has to be as well
3556 
3557  Assert((line_cell_count[line->child_index(0)] == 0 &&
3558  line_cell_count[line->child_index(1)] == 0) ||
3559  (line_cell_count[line->child_index(0)] > 0 &&
3560  line_cell_count[line->child_index(1)] > 0),
3561  ExcInternalError());
3562 
3563  if (line_cell_count[line->child_index(0)] == 0)
3564  {
3565  for (unsigned int c = 0; c < 2; ++c)
3566  Assert(!line->child(c)->has_children(),
3567  ExcInternalError());
3568 
3569  // we may delete the line's
3570  // children and the middle vertex
3571  // as no cell references them
3572  // anymore
3573  triangulation
3574  .vertices_used[line->child(0)->vertex_index(1)] = false;
3575 
3576  lines_to_delete.push_back(line->child(0));
3577  lines_to_delete.push_back(line->child(1));
3578 
3579  line->clear_children();
3580  }
3581  }
3582  }
3583 
3584  // finally, delete unneeded quads and lines
3585 
3586  // clear user pointers, to avoid that
3587  // they may appear at unwanted places
3588  // later on...
3589  // same for user flags, then finally
3590  // delete the quads and lines
3591  typename std::vector<
3593  line = lines_to_delete.begin(),
3594  endline = lines_to_delete.end();
3595  for (; line != endline; ++line)
3596  {
3597  (*line)->clear_user_data();
3598  (*line)->clear_user_flag();
3599  (*line)->clear_used_flag();
3600  }
3601 
3602  typename std::vector<
3604  quad = quads_to_delete.begin(),
3605  endquad = quads_to_delete.end();
3606  for (; quad != endquad; ++quad)
3607  {
3608  (*quad)->clear_user_data();
3609  (*quad)->clear_children();
3610  (*quad)->clear_refinement_case();
3611  (*quad)->clear_user_flag();
3612  (*quad)->clear_used_flag();
3613  }
3614  }
3615 
3616 
3634  template <int spacedim>
3635  static void create_children(
3637  unsigned int & next_unused_vertex,
3639  &next_unused_line,
3641  &next_unused_cell,
3642  const typename Triangulation<2, spacedim>::cell_iterator &cell)
3643  {
3644  const unsigned int dim = 2;
3645  // clear refinement flag
3646  const RefinementCase<dim> ref_case = cell->refine_flag_set();
3647  cell->clear_refine_flag();
3648 
3649  /* For the refinement process: since we go the levels up from the
3650  lowest, there are (unlike above) only two possibilities: a neighbor
3651  cell is on the same level or one level up (in both cases, it may or
3652  may not be refined later on, but we don't care here).
3653 
3654  First:
3655  Set up an array of the 3x3 vertices, which are distributed on the
3656  cell (the array consists of indices into the @p{vertices} std::vector
3657 
3658  2--7--3
3659  | | |
3660  4--8--5
3661  | | |
3662  0--6--1
3663 
3664  note: in case of cut_x or cut_y not all these vertices are needed for
3665  the new cells
3666 
3667  Second:
3668  Set up an array of the new lines (the array consists of iterator
3669  pointers into the lines arrays)
3670 
3671  .-6-.-7-. The directions are: .->-.->-.
3672  1 9 3 ^ ^ ^
3673  .-10.11-. .->-.->-.
3674  0 8 2 ^ ^ ^
3675  .-4-.-5-. .->-.->-.
3676 
3677  cut_x:
3678  .-4-.-5-.
3679  | | |
3680  0 6 1
3681  | | |
3682  .-2-.-3-.
3683 
3684  cut_y:
3685  .---5---.
3686  1 3
3687  .---6---.
3688  0 2
3689  .---4---.
3690 
3691 
3692  Third:
3693  Set up an array of neighbors:
3694 
3695  6 7
3696  .--.--.
3697  1| | |3
3698  .--.--.
3699  0| | |2
3700  .--.--.
3701  4 5
3702 
3703  We need this array for two reasons: first to get the lines which will
3704  bound the four subcells (if the neighboring cell is refined, these
3705  lines already exist), and second to update neighborship information.
3706  Since if a neighbor is not refined, its neighborship record only
3707  points to the present, unrefined, cell rather than the children we
3708  are presently creating, we only need the neighborship information
3709  if the neighbor cells are refined. In all other cases, we store
3710  the unrefined neighbor address
3711 
3712  We also need for every neighbor (if refined) which number among its
3713  neighbors the present (unrefined) cell has, since that number is to
3714  be replaced and because that also is the number of the subline which
3715  will be the interface between that neighbor and the to be created
3716  cell. We will store this number (between 0 and 3) in the field
3717  @p{neighbors_neighbor}.
3718 
3719  It would be sufficient to use the children of the common line to the
3720  neighbor, if we only wanted to get the new sublines and the new
3721  vertex, but because we need to update the neighborship information of
3722  the two refined subcells of the neighbor, we need to search these
3723  anyway.
3724 
3725  Convention:
3726  The created children are numbered like this:
3727 
3728  .--.--.
3729  |2 . 3|
3730  .--.--.
3731  |0 | 1|
3732  .--.--.
3733  */
3734  // collect the
3735  // indices of the
3736  // eight
3737  // surrounding
3738  // vertices
3739  // 2--7--3
3740  // | | |
3741  // 4--9--5
3742  // | | |
3743  // 0--6--1
3744  int new_vertices[9];
3745  for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3746  new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3747  for (unsigned int line_no = 0; line_no < 4; ++line_no)
3748  if (cell->line(line_no)->has_children())
3749  new_vertices[4 + line_no] =
3750  cell->line(line_no)->child(0)->vertex_index(1);
3751 
3752  if (ref_case == RefinementCase<dim>::cut_xy)
3753  {
3754  // find the next
3755  // unused vertex and
3756  // allocate it for
3757  // the new vertex we
3758  // need here
3759  while (triangulation.vertices_used[next_unused_vertex] == true)
3760  ++next_unused_vertex;
3761  Assert(
3762  next_unused_vertex < triangulation.vertices.size(),
3763  ExcMessage(
3764  "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."));
3765  triangulation.vertices_used[next_unused_vertex] = true;
3766 
3767  new_vertices[8] = next_unused_vertex;
3768 
3769  // if this quad lives
3770  // in 2d, then we can
3771  // compute the new
3772  // central vertex
3773  // location just from
3774  // the surrounding
3775  // ones. If this is
3776  // not the case, then
3777  // we need to ask a
3778  // boundary object
3779  if (dim == spacedim)
3780  {
3781  // triangulation.vertices[next_unused_vertex] = new_point;
3782  triangulation.vertices[next_unused_vertex] = cell->center(true);
3783 
3784  // if the user_flag is set, i.e. if the cell is at the
3785  // boundary, use a different calculation of the middle vertex
3786  // here. this is of advantage if the boundary is strongly
3787  // curved (whereas the cell is not) and the cell has a high
3788  // aspect ratio.
3789  if (cell->user_flag_set())
3790  {
3791  // first reset the user_flag and then refine
3792  cell->clear_user_flag();
3793  triangulation.vertices[next_unused_vertex] =
3794  cell->center(true, true);
3795  }
3796  }
3797  else
3798  {
3799  // if this quad lives in a higher dimensional space
3800  // then we don't need to worry if it is at the
3801  // boundary of the manifold -- we always have to use
3802  // the boundary object anyway; so ignore whether the
3803  // user flag is set or not
3804  cell->clear_user_flag();
3805 
3806  // determine middle vertex by transfinite interpolation to be
3807  // consistent with what happens to quads in a Triangulation<3,
3808  // 3> when they are refined
3809  triangulation.vertices[next_unused_vertex] =
3810  cell->center(true, true);
3811  }
3812  }
3813 
3814 
3815  // Now the lines:
3816  typename Triangulation<dim, spacedim>::raw_line_iterator new_lines[12];
3817  unsigned int lmin = 8;
3818  unsigned int lmax = 12;
3819  if (ref_case != RefinementCase<dim>::cut_xy)
3820  {
3821  lmin = 6;
3822  lmax = 7;
3823  }
3824 
3825  for (unsigned int l = lmin; l < lmax; ++l)
3826  {
3827  while (next_unused_line->used() == true)
3828  ++next_unused_line;
3829  new_lines[l] = next_unused_line;
3830  ++next_unused_line;
3831 
3832  Assert(
3833  new_lines[l]->used() == false,
3834  ExcMessage(
3835  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3836  }
3837 
3838  if (ref_case == RefinementCase<dim>::cut_xy)
3839  {
3840  // .-6-.-7-.
3841  // 1 9 3
3842  // .-10.11-.
3843  // 0 8 2
3844  // .-4-.-5-.
3845 
3846  // lines 0-7 already exist, create only the four interior
3847  // lines 8-11
3848  unsigned int l = 0;
3849  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
3850  for (unsigned int c = 0; c < 2; ++c, ++l)
3851  new_lines[l] = cell->line(face_no)->child(c);
3852  Assert(l == 8, ExcInternalError());
3853 
3854  new_lines[8]->set_bounding_object_indices(
3855  {new_vertices[6], new_vertices[8]});
3856  new_lines[9]->set_bounding_object_indices(
3857  {new_vertices[8], new_vertices[7]});
3858  new_lines[10]->set_bounding_object_indices(
3859  {new_vertices[4], new_vertices[8]});
3860  new_lines[11]->set_bounding_object_indices(
3861  {new_vertices[8], new_vertices[5]});
3862  }
3863  else if (ref_case == RefinementCase<dim>::cut_x)
3864  {
3865  // .-4-.-5-.
3866  // | | |
3867  // 0 6 1
3868  // | | |
3869  // .-2-.-3-.
3870  new_lines[0] = cell->line(0);
3871  new_lines[1] = cell->line(1);
3872  new_lines[2] = cell->line(2)->child(0);
3873  new_lines[3] = cell->line(2)->child(1);
3874  new_lines[4] = cell->line(3)->child(0);
3875  new_lines[5] = cell->line(3)->child(1);
3876  new_lines[6]->set_bounding_object_indices(
3877  {new_vertices[6], new_vertices[7]});
3878  }
3879  else
3880  {
3882  // .---5---.
3883  // 1 3
3884  // .---6---.
3885  // 0 2
3886  // .---4---.
3887  new_lines[0] = cell->line(0)->child(0);
3888  new_lines[1] = cell->line(0)->child(1);
3889  new_lines[2] = cell->line(1)->child(0);
3890  new_lines[3] = cell->line(1)->child(1);
3891  new_lines[4] = cell->line(2);
3892  new_lines[5] = cell->line(3);
3893  new_lines[6]->set_bounding_object_indices(
3894  {new_vertices[4], new_vertices[5]});
3895  }
3896 
3897  for (unsigned int l = lmin; l < lmax; ++l)
3898  {
3899  new_lines[l]->set_used_flag();
3900  new_lines[l]->clear_user_flag();
3901  new_lines[l]->clear_user_data();
3902  new_lines[l]->clear_children();
3903  // interior line
3904  new_lines[l]->set_boundary_id_internal(
3906  new_lines[l]->set_manifold_id(cell->manifold_id());
3907  }
3908 
3909  // Now add the four (two)
3910  // new cells!
3913  while (next_unused_cell->used() == true)
3914  ++next_unused_cell;
3915 
3916  const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
3917  for (unsigned int i = 0; i < n_children; ++i)
3918  {
3919  Assert(
3920  next_unused_cell->used() == false,
3921  ExcMessage(
3922  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3923  subcells[i] = next_unused_cell;
3924  ++next_unused_cell;
3925  if (i % 2 == 1 && i < n_children - 1)
3926  while (next_unused_cell->used() == true)
3927  ++next_unused_cell;
3928  }
3929 
3930  if (ref_case == RefinementCase<dim>::cut_xy)
3931  {
3932  // children:
3933  // .--.--.
3934  // |2 . 3|
3935  // .--.--.
3936  // |0 | 1|
3937  // .--.--.
3938  // lines:
3939  // .-6-.-7-.
3940  // 1 9 3
3941  // .-10.11-.
3942  // 0 8 2
3943  // .-4-.-5-.
3944  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3945  new_lines[8]->index(),
3946  new_lines[4]->index(),
3947  new_lines[10]->index()});
3948  subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3949  new_lines[2]->index(),
3950  new_lines[5]->index(),
3951  new_lines[11]->index()});
3952  subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3953  new_lines[9]->index(),
3954  new_lines[10]->index(),
3955  new_lines[6]->index()});
3956  subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3957  new_lines[3]->index(),
3958  new_lines[11]->index(),
3959  new_lines[7]->index()});
3960  }
3961  else if (ref_case == RefinementCase<dim>::cut_x)
3962  {
3963  // children:
3964  // .--.--.
3965  // | . |
3966  // .0 . 1.
3967  // | | |
3968  // .--.--.
3969  // lines:
3970  // .-4-.-5-.
3971  // | | |
3972  // 0 6 1
3973  // | | |
3974  // .-2-.-3-.
3975  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3976  new_lines[6]->index(),
3977  new_lines[2]->index(),
3978  new_lines[4]->index()});
3979  subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3980  new_lines[1]->index(),
3981  new_lines[3]->index(),
3982  new_lines[5]->index()});
3983  }
3984  else
3985  {
3987  // children:
3988  // .-----.
3989  // | 1 |
3990  // .-----.
3991  // | 0 |
3992  // .-----.
3993  // lines:
3994  // .---5---.
3995  // 1 3
3996  // .---6---.
3997  // 0 2
3998  // .---4---.
3999  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4000  new_lines[2]->index(),
4001  new_lines[4]->index(),
4002  new_lines[6]->index()});
4003  subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4004  new_lines[3]->index(),
4005  new_lines[6]->index(),
4006  new_lines[5]->index()});
4007  }
4008 
4009  types::subdomain_id subdomainid = cell->subdomain_id();
4010 
4011  for (unsigned int i = 0; i < n_children; ++i)
4012  {
4013  subcells[i]->set_used_flag();
4014  subcells[i]->clear_refine_flag();
4015  subcells[i]->clear_user_flag();
4016  subcells[i]->clear_user_data();
4017  subcells[i]->clear_children();
4018  // inherit material
4019  // properties
4020  subcells[i]->set_material_id(cell->material_id());
4021  subcells[i]->set_manifold_id(cell->manifold_id());
4022  subcells[i]->set_subdomain_id(subdomainid);
4023 
4024  if (i % 2 == 0)
4025  subcells[i]->set_parent(cell->index());
4026  }
4027 
4028 
4029 
4030  // set child index for
4031  // even children children
4032  // i=0,2 (0)
4033  for (unsigned int i = 0; i < n_children / 2; ++i)
4034  cell->set_children(2 * i, subcells[2 * i]->index());
4035  // set the refine case
4036  cell->set_refinement_case(ref_case);
4037 
4038  // note that the
4039  // refinement flag was
4040  // already cleared at the
4041  // beginning of this function
4042 
4043  if (dim < spacedim)
4044  for (unsigned int c = 0; c < n_children; ++c)
4045  cell->child(c)->set_direction_flag(cell->direction_flag());
4046  }
4047 
4048 
4049 
4054  template <int spacedim>
4057  const bool /*check_for_distorted_cells*/)
4058  {
4059  const unsigned int dim = 1;
4060 
4061  // check whether a new level is needed we have to check for
4062  // this on the highest level only (on this, all used cells are
4063  // also active, so we only have to check for this)
4064  {
4066  cell = triangulation.begin_active(triangulation.levels.size() - 1),
4067  endc = triangulation.end();
4068  for (; cell != endc; ++cell)
4069  if (cell->used())
4070  if (cell->refine_flag_set())
4071  {
4072  triangulation.levels.push_back(
4073  std::make_unique<
4075  break;
4076  }
4077  }
4078 
4079 
4080  // check how much space is needed on every level we need not
4081  // check the highest level since either - on the highest level
4082  // no cells are flagged for refinement - there are, but
4083  // prepare_refinement added another empty level
4084  unsigned int needed_vertices = 0;
4085  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4086  {
4087  // count number of flagged
4088  // cells on this level
4089  unsigned int flagged_cells = 0;
4090 
4091  for (const auto &acell :
4092  triangulation.active_cell_iterators_on_level(level))
4093  if (acell->refine_flag_set())
4094  ++flagged_cells;
4095 
4096  // count number of used cells
4097  // on the next higher level
4098  const unsigned int used_cells =
4099  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4100  triangulation.levels[level + 1]->cells.used.end(),
4101  true);
4102 
4103  // reserve space for the used_cells cells already existing
4104  // on the next higher level as well as for the
4105  // 2*flagged_cells that will be created on that level
4106  reserve_space(*triangulation.levels[level + 1],
4108  flagged_cells,
4109  1,
4110  spacedim);
4111  // reserve space for 2*flagged_cells new lines on the next
4112  // higher level
4113  reserve_space(triangulation.levels[level + 1]->cells,
4115  flagged_cells,
4116  0);
4117 
4118  needed_vertices += flagged_cells;
4119  }
4120 
4121  // add to needed vertices how many
4122  // vertices are already in use
4123  needed_vertices += std::count(triangulation.vertices_used.begin(),
4124  triangulation.vertices_used.end(),
4125  true);
4126  // if we need more vertices: create them, if not: leave the
4127  // array as is, since shrinking is not really possible because
4128  // some of the vertices at the end may be in use
4129  if (needed_vertices > triangulation.vertices.size())
4130  {
4131  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4132  triangulation.vertices_used.resize(needed_vertices, false);
4133  }
4134 
4135 
4136  // Do REFINEMENT on every level; exclude highest level as
4137  // above
4138 
4139  // index of next unused vertex
4140  unsigned int next_unused_vertex = 0;
4141 
4142  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4143  {
4145  next_unused_cell = triangulation.begin_raw(level + 1);
4146 
4147  for (const auto &cell :
4148  triangulation.active_cell_iterators_on_level(level))
4149  if (cell->refine_flag_set())
4150  {
4151  // clear refinement flag
4152  cell->clear_refine_flag();
4153 
4154  // search for next unused
4155  // vertex
4156  while (triangulation.vertices_used[next_unused_vertex] ==
4157  true)
4158  ++next_unused_vertex;
4159  Assert(
4160  next_unused_vertex < triangulation.vertices.size(),
4161  ExcMessage(
4162  "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."));
4163 
4164  // Now we always ask the cell itself where to put
4165  // the new point. The cell in turn will query the
4166  // manifold object internally.
4167  triangulation.vertices[next_unused_vertex] =
4168  cell->center(true);
4169 
4170  triangulation.vertices_used[next_unused_vertex] = true;
4171 
4172  // search for next two unused cell (++ takes care of
4173  // the end of the vector)
4175  first_child,
4176  second_child;
4177  while (next_unused_cell->used() == true)
4178  ++next_unused_cell;
4179  first_child = next_unused_cell;
4180  first_child->set_used_flag();
4181  first_child->clear_user_data();
4182  ++next_unused_cell;
4183  Assert(
4184  next_unused_cell->used() == false,
4185  ExcMessage(
4186  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4187  second_child = next_unused_cell;
4188  second_child->set_used_flag();
4189  second_child->clear_user_data();
4190 
4191  types::subdomain_id subdomainid = cell->subdomain_id();
4192 
4193  // insert first child
4194  cell->set_children(0, first_child->index());
4195  first_child->clear_children();
4196  first_child->set_bounding_object_indices(
4197  {cell->vertex_index(0), next_unused_vertex});
4198  first_child->set_material_id(cell->material_id());
4199  first_child->set_manifold_id(cell->manifold_id());
4200  first_child->set_subdomain_id(subdomainid);
4201  first_child->set_direction_flag(cell->direction_flag());
4202 
4203  first_child->set_parent(cell->index());
4204 
4205  // Set manifold id of the right face. Only do this
4206  // on the first child.
4207  first_child->face(1)->set_manifold_id(cell->manifold_id());
4208 
4209  // reset neighborship info (refer to
4210  // internal::TriangulationImplementation::TriaLevel<0> for
4211  // details)
4212  first_child->set_neighbor(1, second_child);
4213  if (cell->neighbor(0).state() != IteratorState::valid)
4214  first_child->set_neighbor(0, cell->neighbor(0));
4215  else if (cell->neighbor(0)->is_active())
4216  {
4217  // since the neighbors level is always <=level,
4218  // if the cell is active, then there are no
4219  // cells to the left which may want to know
4220  // about this new child cell.
4221  Assert(cell->neighbor(0)->level() <= cell->level(),
4222  ExcInternalError());
4223  first_child->set_neighbor(0, cell->neighbor(0));
4224  }
4225  else
4226  // left neighbor is refined
4227  {
4228  // set neighbor to cell on same level
4229  const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4230  first_child->set_neighbor(0,
4231  cell->neighbor(0)->child(nbnb));
4232 
4233  // reset neighbor info of all right descendant
4234  // of the left neighbor of cell
4236  left_neighbor = cell->neighbor(0);
4237  while (left_neighbor->has_children())
4238  {
4239  left_neighbor = left_neighbor->child(nbnb);
4240  left_neighbor->set_neighbor(nbnb, first_child);
4241  }
4242  }
4243 
4244  // insert second child
4245  second_child->clear_children();
4246  second_child->set_bounding_object_indices(
4247  {next_unused_vertex, cell->vertex_index(1)});
4248  second_child->set_neighbor(0, first_child);
4249  second_child->set_material_id(cell->material_id());
4250  second_child->set_manifold_id(cell->manifold_id());
4251  second_child->set_subdomain_id(subdomainid);
4252  second_child->set_direction_flag(cell->direction_flag());
4253 
4254  if (cell->neighbor(1).state() != IteratorState::valid)
4255  second_child->set_neighbor(1, cell->neighbor(1));
4256  else if (cell->neighbor(1)->is_active())
4257  {
4258  Assert(cell->neighbor(1)->level() <= cell->level(),
4259  ExcInternalError());
4260  second_child->set_neighbor(1, cell->neighbor(1));
4261  }
4262  else
4263  // right neighbor is refined same as above
4264  {
4265  const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4266  second_child->set_neighbor(
4267  1, cell->neighbor(1)->child(nbnb));
4268 
4270  right_neighbor = cell->neighbor(1);
4271  while (right_neighbor->has_children())
4272  {
4273  right_neighbor = right_neighbor->child(nbnb);
4274  right_neighbor->set_neighbor(nbnb, second_child);
4275  }
4276  }
4277  // inform all listeners that cell refinement is done
4278  triangulation.signals.post_refinement_on_cell(cell);
4279  }
4280  }
4281 
4282  // in 1d, we can not have distorted children unless the parent
4283  // was already distorted (that is because we don't use
4284  // boundary information for 1d triangulations). so return an
4285  // empty list
4287  }
4288 
4289 
4294  template <int spacedim>
4297  const bool check_for_distorted_cells)
4298  {
4299  const unsigned int dim = 2;
4300 
4301  // check whether a new level is needed we have to check for
4302  // this on the highest level only (on this, all used cells are
4303  // also active, so we only have to check for this)
4304  {
4306  cell = triangulation.begin_active(triangulation.levels.size() - 1),
4307  endc = triangulation.end();
4308  for (; cell != endc; ++cell)
4309  if (cell->used())
4310  if (cell->refine_flag_set())
4311  {
4312  triangulation.levels.push_back(
4313  std::make_unique<
4315  break;
4316  }
4317  }
4318 
4319  // TODO[WB]: we clear user flags and pointers of lines; we're going
4320  // to use them to flag which lines need refinement
4321  for (typename Triangulation<dim, spacedim>::line_iterator line =
4322  triangulation.begin_line();
4323  line != triangulation.end_line();
4324  ++line)
4325  {
4326  line->clear_user_flag();
4327  line->clear_user_data();
4328  }
4329  // running over all cells and lines count the number
4330  // n_single_lines of lines which can be stored as single
4331  // lines, e.g. inner lines
4332  unsigned int n_single_lines = 0;
4333 
4334  // New lines to be created: number lines which are stored in
4335  // pairs (the children of lines must be stored in pairs)
4336  unsigned int n_lines_in_pairs = 0;
4337 
4338  // check how much space is needed on every level we need not
4339  // check the highest level since either - on the highest level
4340  // no cells are flagged for refinement - there are, but
4341  // prepare_refinement added another empty level
4342  unsigned int needed_vertices = 0;
4343  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4344  {
4345  // count number of flagged cells on this level and compute
4346  // how many new vertices and new lines will be needed
4347  unsigned int needed_cells = 0;
4348 
4349  for (const auto &cell :
4350  triangulation.active_cell_iterators_on_level(level))
4351  if (cell->refine_flag_set())
4352  {
4353  if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
4354  {
4355  needed_cells += 4;
4356 
4357  // new vertex at center of cell is needed in any
4358  // case
4359  ++needed_vertices;
4360 
4361  // the four inner lines can be stored as singles
4362  n_single_lines += 4;
4363  }
4364  else // cut_x || cut_y
4365  {
4366  // set the flag showing that anisotropic
4367  // refinement is used for at least one cell
4368  triangulation.anisotropic_refinement = true;
4369 
4370  needed_cells += 2;
4371  // no vertex at center
4372 
4373  // the inner line can be stored as single
4374  n_single_lines += 1;
4375  }
4376 
4377  // mark all faces (lines) for refinement; checking
4378  // locally whether the neighbor would also like to
4379  // refine them is rather difficult for lines so we
4380  // only flag them and after visiting all cells, we
4381  // decide which lines need refinement;
4382  for (const unsigned int line_no :
4384  {
4386  cell->refine_flag_set(), line_no) ==
4388  {
4390  line = cell->line(line_no);
4391  if (line->has_children() == false)
4392  line->set_user_flag();
4393  }
4394  }
4395  }
4396 
4397 
4398  // count number of used cells on the next higher level
4399  const unsigned int used_cells =
4400  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4401  triangulation.levels[level + 1]->cells.used.end(),
4402  true);
4403 
4404 
4405  // reserve space for the used_cells cells already existing
4406  // on the next higher level as well as for the
4407  // needed_cells that will be created on that level
4408  reserve_space(*triangulation.levels[level + 1],
4409  used_cells + needed_cells,
4410  2,
4411  spacedim);
4412 
4413  // reserve space for needed_cells new quads on the next
4414  // higher level
4415  reserve_space(triangulation.levels[level + 1]->cells,
4416  needed_cells,
4417  0);
4418  }
4419 
4420  // now count the lines which were flagged for refinement
4421  for (typename Triangulation<dim, spacedim>::line_iterator line =
4422  triangulation.begin_line();
4423  line != triangulation.end_line();
4424  ++line)
4425  if (line->user_flag_set())
4426  {
4427  Assert(line->has_children() == false, ExcInternalError());
4428  n_lines_in_pairs += 2;
4429  needed_vertices += 1;
4430  }
4431  // reserve space for n_lines_in_pairs new lines. note, that
4432  // we can't reserve space for the single lines here as well,
4433  // as all the space reserved for lines in pairs would be
4434  // counted as unused and we would end up with too little space
4435  // to store all lines. memory reservation for n_single_lines
4436  // can only be done AFTER we refined the lines of the current
4437  // cells
4438  reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4439 
4440  // add to needed vertices how many vertices are already in use
4441  needed_vertices += std::count(triangulation.vertices_used.begin(),
4442  triangulation.vertices_used.end(),
4443  true);
4444  // if we need more vertices: create them, if not: leave the
4445  // array as is, since shrinking is not really possible because
4446  // some of the vertices at the end may be in use
4447  if (needed_vertices > triangulation.vertices.size())
4448  {
4449  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4450  triangulation.vertices_used.resize(needed_vertices, false);
4451  }
4452 
4453 
4454  // Do REFINEMENT on every level; exclude highest level as
4455  // above
4456 
4457  // index of next unused vertex
4458  unsigned int next_unused_vertex = 0;
4459 
4460  // first the refinement of lines. children are stored
4461  // pairwise
4462  {
4463  // only active objects can be refined further
4465  line = triangulation.begin_active_line(),
4466  endl = triangulation.end_line();
4468  next_unused_line = triangulation.begin_raw_line();
4469 
4470  for (; line != endl; ++line)
4471  if (line->user_flag_set())
4472  {
4473  // this line needs to be refined
4474 
4475  // find the next unused vertex and set it
4476  // appropriately
4477  while (triangulation.vertices_used[next_unused_vertex] == true)
4478  ++next_unused_vertex;
4479  Assert(
4480  next_unused_vertex < triangulation.vertices.size(),
4481  ExcMessage(
4482  "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."));
4483  triangulation.vertices_used[next_unused_vertex] = true;
4484 
4485  triangulation.vertices[next_unused_vertex] = line->center(true);
4486 
4487  // now that we created the right point, make up the
4488  // two child lines. To this end, find a pair of
4489  // unused lines
4490  bool pair_found = false;
4491  (void)pair_found;
4492  for (; next_unused_line != endl; ++next_unused_line)
4493  if (!next_unused_line->used() &&
4494  !(++next_unused_line)->used())
4495  {
4496  // go back to the first of the two unused
4497  // lines
4498  --next_unused_line;
4499  pair_found = true;
4500  break;
4501  }
4502  Assert(pair_found, ExcInternalError());
4503 
4504  // there are now two consecutive unused lines, such
4505  // that the children of a line will be consecutive.
4506  // then set the child pointer of the present line
4507  line->set_children(0, next_unused_line->index());
4508 
4509  // set the two new lines
4511  children[2] = {next_unused_line, ++next_unused_line};
4512  // some tests; if any of the iterators should be
4513  // invalid, then already dereferencing will fail
4514  Assert(
4515  children[0]->used() == false,
4516  ExcMessage(
4517  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4518  Assert(
4519  children[1]->used() == false,
4520  ExcMessage(
4521  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4522 
4523  children[0]->set_bounding_object_indices(
4524  {line->vertex_index(0), next_unused_vertex});
4525  children[1]->set_bounding_object_indices(
4526  {next_unused_vertex, line->vertex_index(1)});
4527 
4528  children[0]->set_used_flag();
4529  children[1]->set_used_flag();
4530  children[0]->clear_children();
4531  children[1]->clear_children();
4532  children[0]->clear_user_data();
4533  children[1]->clear_user_data();
4534  children[0]->clear_user_flag();
4535  children[1]->clear_user_flag();
4536 
4537 
4538  children[0]->set_boundary_id_internal(line->boundary_id());
4539  children[1]->set_boundary_id_internal(line->boundary_id());
4540 
4541  children[0]->set_manifold_id(line->manifold_id());
4542  children[1]->set_manifold_id(line->manifold_id());
4543 
4544  // finally clear flag indicating the need for
4545  // refinement
4546  line->clear_user_flag();
4547  }
4548  }
4549 
4550 
4551  // Now set up the new cells
4552 
4553  // reserve space for inner lines (can be stored as single
4554  // lines)
4555  reserve_space(triangulation.faces->lines, 0, n_single_lines);
4556 
4558  cells_with_distorted_children;
4559 
4560  // reset next_unused_line, as now also single empty places in
4561  // the vector can be used
4563  next_unused_line = triangulation.begin_raw_line();
4564 
4565  for (int level = 0;
4566  level < static_cast<int>(triangulation.levels.size()) - 1;
4567  ++level)
4568  {
4570  next_unused_cell = triangulation.begin_raw(level + 1);
4571 
4572  for (const auto &cell :
4573  triangulation.active_cell_iterators_on_level(level))
4574  if (cell->refine_flag_set())
4575  {
4576  // set the user flag to indicate, that at least one
4577  // line is at the boundary
4578 
4579  // TODO[Tobias Leicht] find a better place to set
4580  // this flag, so that we do not need so much time to
4581  // check each cell here
4582  if (cell->at_boundary())
4583  cell->set_user_flag();
4584 
4585  // actually set up the children and update neighbor
4586  // information
4587  create_children(triangulation,
4588  next_unused_vertex,
4589  next_unused_line,
4590  next_unused_cell,
4591  cell);
4592 
4593  if (check_for_distorted_cells &&
4594  has_distorted_children(
4595  cell,
4596  std::integral_constant<int, dim>(),
4597  std::integral_constant<int, spacedim>()))
4598  cells_with_distorted_children.distorted_cells.push_back(
4599  cell);
4600  // inform all listeners that cell refinement is done
4601  triangulation.signals.post_refinement_on_cell(cell);
4602  }
4603  }
4604 
4605  return cells_with_distorted_children;
4606  }
4607 
4608 
4613  template <int spacedim>
4616  const bool check_for_distorted_cells)
4617  {
4618  const unsigned int dim = 3;
4619 
4620  // this function probably also works for spacedim>3 but it
4621  // isn't tested. it will probably be necessary to pull new
4622  // vertices onto the manifold just as we do for the other
4623  // functions above.
4624  Assert(spacedim == 3, ExcNotImplemented());
4625 
4626  // check whether a new level is needed we have to check for
4627  // this on the highest level only (on this, all used cells are
4628  // also active, so we only have to check for this)
4629  {
4631  cell = triangulation.begin_active(triangulation.levels.size() - 1),
4632  endc = triangulation.end();
4633  for (; cell != endc; ++cell)
4634  if (cell->used())
4635  if (cell->refine_flag_set())
4636  {
4637  triangulation.levels.push_back(
4638  std::make_unique<
4640  break;
4641  }
4642  }
4643 
4644 
4645  // first clear user flags for quads and lines; we're going to
4646  // use them to flag which lines and quads need refinement
4647  triangulation.faces->quads.clear_user_data();
4648 
4649  for (typename Triangulation<dim, spacedim>::line_iterator line =
4650  triangulation.begin_line();
4651  line != triangulation.end_line();
4652  ++line)
4653  line->clear_user_flag();
4654  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
4655  triangulation.begin_quad();
4656  quad != triangulation.end_quad();
4657  ++quad)
4658  quad->clear_user_flag();
4659 
4660  // create an array of face refine cases. User indices of faces
4661  // will be set to values corresponding with indices in this
4662  // array.
4663  const RefinementCase<dim - 1> face_refinement_cases[4] = {
4664  RefinementCase<dim - 1>::no_refinement,
4665  RefinementCase<dim - 1>::cut_x,
4666  RefinementCase<dim - 1>::cut_y,
4667  RefinementCase<dim - 1>::cut_xy};
4668 
4669  // check how much space is needed on every level we need not
4670  // check the highest level since either
4671  // - on the highest level no cells are flagged for refinement
4672  // - there are, but prepare_refinement added another empty
4673  // level which then is the highest level
4674 
4675  // variables to hold the number of newly to be created
4676  // vertices, lines and quads. as these are stored globally,
4677  // declare them outside the loop over al levels. we need lines
4678  // and quads in pairs for refinement of old ones and lines and
4679  // quads, that can be stored as single ones, as they are newly
4680  // created in the inside of an existing cell
4681  unsigned int needed_vertices = 0;
4682  unsigned int needed_lines_single = 0;
4683  unsigned int needed_quads_single = 0;
4684  unsigned int needed_lines_pair = 0;
4685  unsigned int needed_quads_pair = 0;
4686  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4687  {
4688  // count number of flagged cells on this level and compute
4689  // how many new vertices and new lines will be needed
4690  unsigned int new_cells = 0;
4691 
4692  for (const auto &acell :
4693  triangulation.active_cell_iterators_on_level(level))
4694  if (acell->refine_flag_set())
4695  {
4696  RefinementCase<dim> ref_case = acell->refine_flag_set();
4697 
4698  // now for interior vertices, lines and quads, which
4699  // are needed in any case
4700  if (ref_case == RefinementCase<dim>::cut_x ||
4701  ref_case == RefinementCase<dim>::cut_y ||
4702  ref_case == RefinementCase<dim>::cut_z)
4703  {
4704  ++needed_quads_single;
4705  new_cells += 2;
4706  triangulation.anisotropic_refinement = true;
4707  }
4708  else if (ref_case == RefinementCase<dim>::cut_xy ||
4709  ref_case == RefinementCase<dim>::cut_xz ||
4710  ref_case == RefinementCase<dim>::cut_yz)
4711  {
4712  ++needed_lines_single;
4713  needed_quads_single += 4;
4714  new_cells += 4;
4715  triangulation.anisotropic_refinement = true;
4716  }
4717  else if (ref_case == RefinementCase<dim>::cut_xyz)
4718  {
4719  ++needed_vertices;
4720  needed_lines_single += 6;
4721  needed_quads_single += 12;
4722  new_cells += 8;
4723  }
4724  else
4725  {
4726  // we should never get here
4727  Assert(false, ExcInternalError());
4728  }
4729 
4730  // mark all faces for refinement; checking locally
4731  // if and how the neighbor would like to refine
4732  // these is difficult so we only flag them and after
4733  // visiting all cells, we decide which faces need
4734  // which refinement;
4735  for (const unsigned int face :
4737  {
4739  aface = acell->face(face);
4740  // get the RefineCase this faces has for the
4741  // given RefineCase of the cell
4742  RefinementCase<dim - 1> face_ref_case =
4744  ref_case,
4745  face,
4746  acell->face_orientation(face),
4747  acell->face_flip(face),
4748  acell->face_rotation(face));
4749  // only do something, if this face has to be
4750  // refined
4751  if (face_ref_case)
4752  {
4753  if (face_ref_case ==
4755  {
4756  if (aface->number_of_children() < 4)
4757  // we use user_flags to denote needed
4758  // isotropic refinement
4759  aface->set_user_flag();
4760  }
4761  else if (aface->refinement_case() != face_ref_case)
4762  // we use user_indices to denote needed
4763  // anisotropic refinement. note, that we
4764  // can have at most one anisotropic
4765  // refinement case for this face, as
4766  // otherwise prepare_refinement() would
4767  // have changed one of the cells to yield
4768  // isotropic refinement at this
4769  // face. therefore we set the user_index
4770  // uniquely
4771  {
4772  Assert(aface->refinement_case() ==
4774  dim - 1>::isotropic_refinement ||
4775  aface->refinement_case() ==
4776  RefinementCase<dim - 1>::no_refinement,
4777  ExcInternalError());
4778  aface->set_user_index(face_ref_case);
4779  }
4780  }
4781  } // for all faces
4782 
4783  // flag all lines, that have to be refined
4784  for (unsigned int line = 0;
4785  line < GeometryInfo<dim>::lines_per_cell;
4786  ++line)
4788  line) &&
4789  !acell->line(line)->has_children())
4790  acell->line(line)->set_user_flag();
4791 
4792  } // if refine_flag set and for all cells on this level
4793 
4794 
4795  // count number of used cells on the next higher level
4796  const unsigned int used_cells =
4797  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4798  triangulation.levels[level + 1]->cells.used.end(),
4799  true);
4800 
4801 
4802  // reserve space for the used_cells cells already existing
4803  // on the next higher level as well as for the
4804  // 8*flagged_cells that will be created on that level
4805  reserve_space(*triangulation.levels[level + 1],
4806  used_cells + new_cells,
4807  3,
4808  spacedim);
4809  // reserve space for 8*flagged_cells new hexes on the next
4810  // higher level
4811  reserve_space(triangulation.levels[level + 1]->cells, new_cells);
4812  } // for all levels
4813  // now count the quads and lines which were flagged for
4814  // refinement
4815  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
4816  triangulation.begin_quad();
4817  quad != triangulation.end_quad();
4818  ++quad)
4819  {
4820  if (quad->user_flag_set())
4821  {
4822  // isotropic refinement: 1 interior vertex, 4 quads
4823  // and 4 interior lines. we store the interior lines
4824  // in pairs in case the face is already or will be
4825  // refined anisotropically
4826  needed_quads_pair += 4;
4827  needed_lines_pair += 4;
4828  needed_vertices += 1;
4829  }
4830  if (quad->user_index())
4831  {
4832  // anisotropic refinement: 1 interior
4833  // line and two quads
4834  needed_quads_pair += 2;
4835  needed_lines_single += 1;
4836  // there is a kind of complicated situation here which
4837  // requires our attention. if the quad is refined
4838  // isotropcally, two of the interior lines will get a
4839  // new mother line - the interior line of our
4840  // anisotropically refined quad. if those two lines
4841  // are not consecutive, we cannot do so and have to
4842  // replace them by two lines that are consecutive. we
4843  // try to avoid that situation, but it may happen
4844  // nevertheless through repeated refinement and
4845  // coarsening. thus we have to check here, as we will
4846  // need some additional space to store those new lines
4847  // in case we need them...
4848  if (quad->has_children())
4849  {
4850  Assert(quad->refinement_case() ==
4851  RefinementCase<dim - 1>::isotropic_refinement,
4852  ExcInternalError());
4853  if ((face_refinement_cases[quad->user_index()] ==
4854  RefinementCase<dim - 1>::cut_x &&
4855  (quad->child(0)->line_index(1) + 1 !=
4856  quad->child(2)->line_index(1))) ||
4857  (face_refinement_cases[quad->user_index()] ==
4858  RefinementCase<dim - 1>::cut_y &&
4859  (quad->child(0)->line_index(3) + 1 !=
4860  quad->child(1)->line_index(3))))
4861  needed_lines_pair += 2;
4862  }
4863  }
4864  }
4865 
4866  for (typename Triangulation<dim, spacedim>::line_iterator line =
4867  triangulation.begin_line();
4868  line != triangulation.end_line();
4869  ++line)
4870  if (line->user_flag_set())
4871  {
4872  needed_lines_pair += 2;
4873  needed_vertices += 1;
4874  }
4875 
4876  // reserve space for needed_lines new lines stored in pairs
4877  reserve_space(triangulation.faces->lines,
4878  needed_lines_pair,
4879  needed_lines_single);
4880  // reserve space for needed_quads new quads stored in pairs
4881  reserve_space(*triangulation.faces,
4882  needed_quads_pair,
4883  needed_quads_single);
4884  reserve_space(triangulation.faces->quads,
4885  needed_quads_pair,
4886  needed_quads_single);
4887 
4888 
4889  // add to needed vertices how many vertices are already in use
4890  needed_vertices += std::count(triangulation.vertices_used.begin(),
4891  triangulation.vertices_used.end(),
4892  true);
4893  // if we need more vertices: create them, if not: leave the
4894  // array as is, since shrinking is not really possible because
4895  // some of the vertices at the end may be in use
4896  if (needed_vertices > triangulation.vertices.size())
4897  {
4898  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4899  triangulation.vertices_used.resize(needed_vertices, false);
4900  }
4901 
4902 
4904  // Before we start with the actual refinement, we do some
4905  // sanity checks if in debug mode. especially, we try to catch
4906  // the notorious problem with lines being twice refined,
4907  // i.e. there are cells adjacent at one line ("around the
4908  // edge", but not at a face), with two cells differing by more
4909  // than one refinement level
4910  //
4911  // this check is very simple to implement here, since we have
4912  // all lines flagged if they shall be refined
4913 #ifdef DEBUG
4914  for (const auto &cell : triangulation.active_cell_iterators())
4915  if (!cell->refine_flag_set())
4916  for (unsigned int line = 0;
4917  line < GeometryInfo<dim>::lines_per_cell;
4918  ++line)
4919  if (cell->line(line)->has_children())
4920  for (unsigned int c = 0; c < 2; ++c)
4921  Assert(cell->line(line)->child(c)->user_flag_set() == false,
4922  ExcInternalError());
4923 #endif
4924 
4926  // Do refinement on every level
4927  //
4928  // To make life a bit easier, we first refine those lines and
4929  // quads that were flagged for refinement and then compose the
4930  // newly to be created cells.
4931  //
4932  // index of next unused vertex
4933  unsigned int next_unused_vertex = 0;
4934 
4935  // first for lines
4936  {
4937  // only active objects can be refined further
4939  line = triangulation.begin_active_line(),
4940  endl = triangulation.end_line();
4942  next_unused_line = triangulation.begin_raw_line();
4943 
4944  for (; line != endl; ++line)
4945  if (line->user_flag_set())
4946  {
4947  // this line needs to be refined
4948 
4949  // find the next unused vertex and set it
4950  // appropriately
4951  while (triangulation.vertices_used[next_unused_vertex] == true)
4952  ++next_unused_vertex;
4953  Assert(
4954  next_unused_vertex < triangulation.vertices.size(),
4955  ExcMessage(
4956  "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."));
4957  triangulation.vertices_used[next_unused_vertex] = true;
4958 
4959  triangulation.vertices[next_unused_vertex] = line->center(true);
4960 
4961  // now that we created the right point, make up the
4962  // two child lines (++ takes care of the end of the
4963  // vector)
4964  next_unused_line =
4965  triangulation.faces->lines.template next_free_pair_object<1>(
4966  triangulation);
4967  Assert(next_unused_line.state() == IteratorState::valid,
4968  ExcInternalError());
4969 
4970  // now we found two consecutive unused lines, such
4971  // that the children of a line will be consecutive.
4972  // then set the child pointer of the present line
4973  line->set_children(0, next_unused_line->index());
4974 
4975  // set the two new lines
4977  children[2] = {next_unused_line, ++next_unused_line};
4978 
4979  // some tests; if any of the iterators should be
4980  // invalid, then already dereferencing will fail
4981  Assert(
4982  children[0]->used() == false,
4983  ExcMessage(
4984  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4985  Assert(
4986  children[1]->used() == false,
4987  ExcMessage(
4988  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4989 
4990  children[0]->set_bounding_object_indices(
4991  {line->vertex_index(0), next_unused_vertex});
4992  children[1]->set_bounding_object_indices(
4993  {next_unused_vertex, line->vertex_index(1)});
4994 
4995  children[0]->set_used_flag();
4996  children[1]->set_used_flag();
4997  children[0]->clear_children();
4998  children[1]->clear_children();
4999  children[0]->clear_user_data();
5000  children[1]->clear_user_data();
5001  children[0]->clear_user_flag();
5002  children[1]->clear_user_flag();
5003 
5004  children[0]->set_boundary_id_internal(line->boundary_id());
5005  children[1]->set_boundary_id_internal(line->boundary_id());
5006 
5007  children[0]->set_manifold_id(line->manifold_id());
5008  children[1]->set_manifold_id(line->manifold_id());
5009 
5010  // finally clear flag
5011  // indicating the need
5012  // for refinement
5013  line->clear_user_flag();
5014  }
5015  }
5016 
5017 
5019  // now refine marked quads
5021 
5022  // here we encounter several cases:
5023 
5024  // a) the quad is unrefined and shall be refined isotropically
5025 
5026  // b) the quad is unrefined and shall be refined
5027  // anisotropically
5028 
5029  // c) the quad is unrefined and shall be refined both
5030  // anisotropically and isotropically (this is reduced to case
5031  // b) and then case b) for the children again)
5032 
5033  // d) the quad is refined anisotropically and shall be refined
5034  // isotropically (this is reduced to case b) for the
5035  // anisotropic children)
5036 
5037  // e) the quad is refined isotropically and shall be refined
5038  // anisotropically (this is transformed to case c), however we
5039  // might have to renumber/rename children...)
5040 
5041  // we need a loop in cases c) and d), as the anisotropic
5042  // children migt have a lower index than the mother quad
5043  for (unsigned int loop = 0; loop < 2; ++loop)
5044  {
5045  // usually, only active objects can be refined
5046  // further. however, in cases d) and e) that is not true,
5047  // so we have to use 'normal' iterators here
5049  quad = triangulation.begin_quad(),
5050  endq = triangulation.end_quad();
5052  next_unused_line = triangulation.begin_raw_line();
5054  next_unused_quad = triangulation.begin_raw_quad();
5055 
5056  for (; quad != endq; ++quad)
5057  {
5058  if (quad->user_index())
5059  {
5060  RefinementCase<dim - 1> aniso_quad_ref_case =
5061  face_refinement_cases[quad->user_index()];
5062  // there is one unlikely event here, where we
5063  // already have refind the face: if the face was
5064  // refined anisotropically and we want to refine
5065  // it isotropically, both children are flagged for
5066  // anisotropic refinement. however, if those
5067  // children were already flagged for anisotropic
5068  // refinement, they might already be processed and
5069  // refined.
5070  if (aniso_quad_ref_case == quad->refinement_case())
5071  continue;
5072 
5073  Assert(quad->refinement_case() ==
5074  RefinementCase<dim - 1>::cut_xy ||
5075  quad->refinement_case() ==
5076  RefinementCase<dim - 1>::no_refinement,
5077  ExcInternalError());
5078 
5079  // this quad needs to be refined anisotropically
5080  Assert(quad->user_index() ==
5081  RefinementCase<dim - 1>::cut_x ||
5082  quad->user_index() ==
5083  RefinementCase<dim - 1>::cut_y,
5084  ExcInternalError());
5085 
5086  // make the new line interior to the quad
5088  new_line;
5089 
5090  new_line =
5091  triangulation.faces->lines
5092  .template next_free_single_object<1>(triangulation);
5093  Assert(
5094  new_line->used() == false,
5095  ExcMessage(
5096  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5097 
5098  // first collect the
5099  // indices of the vertices:
5100  // *--1--*
5101  // | | |
5102  // | | | cut_x
5103  // | | |
5104  // *--0--*
5105  //
5106  // *-----*
5107  // | |
5108  // 0-----1 cut_y
5109  // | |
5110  // *-----*
5111  unsigned int vertex_indices[2];
5112  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5113  {
5114  vertex_indices[0] =
5115  quad->line(2)->child(0)->vertex_index(1);
5116  vertex_indices[1] =
5117  quad->line(3)->child(0)->vertex_index(1);
5118  }
5119  else
5120  {
5121  vertex_indices[0] =
5122  quad->line(0)->child(0)->vertex_index(1);
5123  vertex_indices[1] =
5124  quad->line(1)->child(0)->vertex_index(1);
5125  }
5126 
5127  new_line->set_bounding_object_indices(
5128  {vertex_indices[0], vertex_indices[1]});
5129  new_line->set_used_flag();
5130  new_line->clear_user_flag();
5131  new_line->clear_user_data();
5132  new_line->clear_children();
5133  new_line->set_boundary_id_internal(quad->boundary_id());
5134  new_line->set_manifold_id(quad->manifold_id());
5135 
5136  // child 0 and 1 of a line are switched if the
5137  // line orientation is false. set up a miniature
5138  // table, indicating which child to take for line
5139  // orientations false and true. first index: child
5140  // index in standard orientation, second index:
5141  // line orientation
5142  const unsigned int index[2][2] = {
5143  {1, 0}, // child 0, line_orientation=false and true
5144  {0, 1}}; // child 1, line_orientation=false and true
5145 
5146  // find some space (consecutive) for the two newly
5147  // to be created quads.
5149  new_quads[2];
5150 
5151  next_unused_quad =
5152  triangulation.faces->quads
5153  .template next_free_pair_object<2>(triangulation);
5154  new_quads[0] = next_unused_quad;
5155  Assert(
5156  new_quads[0]->used() == false,
5157  ExcMessage(
5158  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5159 
5160  ++next_unused_quad;
5161  new_quads[1] = next_unused_quad;
5162  Assert(
5163  new_quads[1]->used() == false,
5164  ExcMessage(
5165  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5166 
5167 
5168  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5169  {
5170  new_quads[0]->set_bounding_object_indices(
5171  {static_cast<int>(quad->line_index(0)),
5172  new_line->index(),
5173  quad->line(2)
5174  ->child(index[0][quad->line_orientation(2)])
5175  ->index(),
5176  quad->line(3)
5177  ->child(index[0][quad->line_orientation(3)])
5178  ->index()});
5179  new_quads[1]->set_bounding_object_indices(
5180  {new_line->index(),
5181  static_cast<int>(quad->line_index(1)),
5182  quad->line(2)
5183  ->child(index[1][quad->line_orientation(2)])
5184  ->index(),
5185  quad->line(3)
5186  ->child(index[1][quad->line_orientation(3)])
5187  ->index()});
5188  }
5189  else
5190  {
5191  new_quads[0]->set_bounding_object_indices(
5192  {quad->line(0)
5193  ->child(index[0][quad->line_orientation(0)])
5194  ->index(),
5195  quad->line(1)
5196  ->child(index[0][quad->line_orientation(1)])
5197  ->index(),
5198  static_cast<int>(quad->line_index(2)),
5199  new_line->index()});
5200  new_quads[1]->set_bounding_object_indices(
5201  {quad->line(0)
5202  ->child(index[1][quad->line_orientation(0)])
5203  ->index(),
5204  quad->line(1)
5205  ->child(index[1][quad->line_orientation(1)])
5206  ->index(),
5207  new_line->index(),
5208  static_cast<int>(quad->line_index(3))});
5209  }
5210 
5211  for (const auto &new_quad : new_quads)
5212  {
5213  new_quad->set_used_flag();
5214  new_quad->clear_user_flag();
5215  new_quad->clear_user_data();
5216  new_quad->clear_children();
5217  new_quad->set_boundary_id_internal(quad->boundary_id());
5218  new_quad->set_manifold_id(quad->manifold_id());
5219  // set all line orientations to true, change
5220  // this after the loop, as we have to consider
5221  // different lines for each child
5222  for (unsigned int j = 0;
5223  j < GeometryInfo<dim>::lines_per_face;
5224  ++j)
5225  new_quad->set_line_orientation(j, true);
5226  }
5227  // now set the line orientation of children of
5228  // outer lines correctly, the lines in the
5229  // interior of the refined quad are automatically
5230  // oriented conforming to the standard
5231  new_quads[0]->set_line_orientation(
5232  0, quad->line_orientation(0));
5233  new_quads[0]->set_line_orientation(
5234  2, quad->line_orientation(2));
5235  new_quads[1]->set_line_orientation(
5236  1, quad->line_orientation(1));
5237  new_quads[1]->set_line_orientation(
5238  3, quad->line_orientation(3));
5239  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5240  {
5241  new_quads[0]->set_line_orientation(
5242  3, quad->line_orientation(3));
5243  new_quads[1]->set_line_orientation(
5244  2, quad->line_orientation(2));
5245  }
5246  else
5247  {
5248  new_quads[0]->set_line_orientation(
5249  1, quad->line_orientation(1));
5250  new_quads[1]->set_line_orientation(
5251  0, quad->line_orientation(0));
5252  }
5253 
5254  // test, whether this face is refined
5255  // isotropically already. if so, set the correct
5256  // children pointers.
5257  if (quad->refinement_case() ==
5258  RefinementCase<dim - 1>::cut_xy)
5259  {
5260  // we will put a new refinemnt level of
5261  // anisotropic refinement between the
5262  // unrefined and isotropically refined quad
5263  // ending up with the same fine quads but
5264  // introducing anisotropically refined ones as
5265  // children of the unrefined quad and mother
5266  // cells of the original fine ones.
5267 
5268  // this process includes the creation of a new
5269  // middle line which we will assign as the
5270  // mother line of two of the existing inner
5271  // lines. If those inner lines are not
5272  // consecutive in memory, we won't find them
5273  // later on, so we have to create new ones
5274  // instead and replace all occurrences of the
5275  // old ones with those new ones. As this is
5276  // kind of ugly, we hope we don't have to do
5277  // it often...
5279  old_child[2];
5280  if (aniso_quad_ref_case ==
5282  {
5283  old_child[0] = quad->child(0)->line(1);
5284  old_child[1] = quad->child(2)->line(1);
5285  }
5286  else
5287  {
5288  Assert(aniso_quad_ref_case ==
5290  ExcInternalError());
5291 
5292  old_child[0] = quad->child(0)->line(3);
5293  old_child[1] = quad->child(1)->line(3);
5294  }
5295 
5296  if (old_child[0]->index() + 1 != old_child[1]->index())
5297  {
5298  // this is exactly the ugly case we taked
5299  // about. so, no coimplaining, lets get
5300  // two new lines and copy all info
5301  typename Triangulation<dim,
5302  spacedim>::raw_line_iterator
5303  new_child[2];
5304 
5305  new_child[0] = new_child[1] =
5306  triangulation.faces->lines
5307  .template next_free_pair_object<1>(
5308  triangulation);
5309  ++new_child[1];
5310 
5311  new_child[0]->set_used_flag();
5312  new_child[1]->set_used_flag();
5313 
5314  const int old_index_0 = old_child[0]->index(),
5315  old_index_1 = old_child[1]->index(),
5316  new_index_0 = new_child[0]->index(),
5317  new_index_1 = new_child[1]->index();
5318 
5319  // loop over all quads and replace the old
5320  // lines
5321  for (unsigned int q = 0;
5322  q < triangulation.faces->quads.n_objects();
5323  ++q)
5324  for (unsigned int l = 0;
5325  l < GeometryInfo<dim>::lines_per_face;
5326  ++l)
5327  {
5328  const int this_index =
5329  triangulation.faces->quads
5330  .get_bounding_object_indices(q)[l];
5331  if (this_index == old_index_0)
5332  triangulation.faces->quads
5333  .get_bounding_object_indices(q)[l] =
5334  new_index_0;
5335  else if (this_index == old_index_1)
5336  triangulation.faces->quads
5337  .get_bounding_object_indices(q)[l] =
5338  new_index_1;
5339  }
5340  // now we have to copy all information of
5341  // the two lines
5342  for (unsigned int i = 0; i < 2; ++i)
5343  {
5344  Assert(!old_child[i]->has_children(),
5345  ExcInternalError());
5346 
5347  new_child[i]->set_bounding_object_indices(
5348  {old_child[i]->vertex_index(0),
5349  old_child[i]->vertex_index(1)});
5350  new_child[i]->set_boundary_id_internal(
5351  old_child[i]->boundary_id());
5352  new_child[i]->set_manifold_id(
5353  old_child[i]->manifold_id());
5354  new_child[i]->set_user_index(
5355  old_child[i]->user_index());
5356  if (old_child[i]->user_flag_set())
5357  new_child[i]->set_user_flag();
5358  else
5359  new_child[i]->clear_user_flag();
5360 
5361  new_child[i]->clear_children();
5362 
5363  old_child[i]->clear_user_flag();
5364  old_child[i]->clear_user_index();
5365  old_child[i]->clear_used_flag();
5366  }
5367  }
5368  // now that we cared about the lines, go on
5369  // with the quads themselves, where we might
5370  // encounter similar situations...
5371  if (aniso_quad_ref_case ==
5373  {
5374  new_line->set_children(
5375  0, quad->child(0)->line_index(1));
5376  Assert(new_line->child(1) ==
5377  quad->child(2)->line(1),
5378  ExcInternalError());
5379  // now evereything is quite
5380  // complicated. we have the children
5381  // numbered according to
5382  //
5383  // *---*---*
5384  // |n+2|n+3|
5385  // *---*---*
5386  // | n |n+1|
5387  // *---*---*
5388  //
5389  // from the original isotropic
5390  // refinement. we have to reorder them as
5391  //
5392  // *---*---*
5393  // |n+1|n+3|
5394  // *---*---*
5395  // | n |n+2|
5396  // *---*---*
5397  //
5398  // such that n and n+1 are consecutive
5399  // children of m and n+2 and n+3 are
5400  // consecutive children of m+1, where m
5401  // and m+1 are given as in
5402  //
5403  // *---*---*
5404  // | | |
5405  // | m |m+1|
5406  // | | |
5407  // *---*---*
5408  //
5409  // this is a bit ugly, of course: loop
5410  // over all cells on all levels and look
5411  // for faces n+1 (switch_1) and n+2
5412  // (switch_2).
5413  const typename Triangulation<dim, spacedim>::
5414  quad_iterator switch_1 = quad->child(1),
5415  switch_2 = quad->child(2);
5416  const int switch_1_index = switch_1->index();
5417  const int switch_2_index = switch_2->index();
5418  for (unsigned int l = 0;
5419  l < triangulation.levels.size();
5420  ++l)
5421  for (unsigned int h = 0;
5422  h <
5423  triangulation.levels[l]->cells.n_objects();
5424  ++h)
5425  for (const unsigned int q :
5427  {
5428  const int face_index =
5429  triangulation.levels[l]
5430  ->cells.get_bounding_object_indices(
5431  h)[q];
5432  if (face_index == switch_1_index)
5433  triangulation.levels[l]
5434  ->cells.get_bounding_object_indices(
5435  h)[q] = switch_2_index;
5436  else if (face_index == switch_2_index)
5437  triangulation.levels[l]
5438  ->cells.get_bounding_object_indices(
5439  h)[q] = switch_1_index;
5440  }
5441  // now we have to copy all information of
5442  // the two quads
5443  const unsigned int switch_1_lines[4] = {
5444  switch_1->line_index(0),
5445  switch_1->line_index(1),
5446  switch_1->line_index(2),
5447  switch_1->line_index(3)};
5448  const bool switch_1_line_orientations[4] = {
5449  switch_1->line_orientation(0),
5450  switch_1->line_orientation(1),
5451  switch_1->line_orientation(2),
5452  switch_1->line_orientation(3)};
5453  const types::boundary_id switch_1_boundary_id =
5454  switch_1->boundary_id();
5455  const unsigned int switch_1_user_index =
5456  switch_1->user_index();
5457  const bool switch_1_user_flag =
5458  switch_1->user_flag_set();
5459  const RefinementCase<dim - 1>
5460  switch_1_refinement_case =
5461  switch_1->refinement_case();
5462  const int switch_1_first_child_pair =
5463  (switch_1_refinement_case ?
5464  switch_1->child_index(0) :
5465  -1);
5466  const int switch_1_second_child_pair =
5467  (switch_1_refinement_case ==
5468  RefinementCase<dim - 1>::cut_xy ?
5469  switch_1->child_index(2) :
5470  -1);
5471 
5472  switch_1->set_bounding_object_indices(
5473  {switch_2->line_index(0),
5474  switch_2->line_index(1),
5475  switch_2->line_index(2),
5476  switch_2->line_index(3)});
5477  switch_1->set_line_orientation(
5478  0, switch_2->line_orientation(0));
5479  switch_1->set_line_orientation(
5480  1, switch_2->line_orientation(1));
5481  switch_1->set_line_orientation(
5482  2, switch_2->line_orientation(2));
5483  switch_1->set_line_orientation(
5484  3, switch_2->line_orientation(3));
5485  switch_1->set_boundary_id_internal(
5486  switch_2->boundary_id());
5487  switch_1->set_manifold_id(switch_2->manifold_id());
5488  switch_1->set_user_index(switch_2->user_index());
5489  if (switch_2->user_flag_set())
5490  switch_1->set_user_flag();
5491  else
5492  switch_1->clear_user_flag();
5493  switch_1->clear_refinement_case();
5494  switch_1->set_refinement_case(
5495  switch_2->refinement_case());
5496  switch_1->clear_children();
5497  if (switch_2->refinement_case())
5498  switch_1->set_children(0,
5499  switch_2->child_index(0));
5500  if (switch_2->refinement_case() ==
5501  RefinementCase<dim - 1>::cut_xy)
5502  switch_1->set_children(2,
5503  switch_2->child_index(2));
5504 
5505  switch_2->set_bounding_object_indices(
5506  {switch_1_lines[0],
5507  switch_1_lines[1],
5508  switch_1_lines[2],
5509  switch_1_lines[3]});
5510  switch_2->set_line_orientation(
5511  0, switch_1_line_orientations[0]);
5512  switch_2->set_line_orientation(
5513  1, switch_1_line_orientations[1]);
5514  switch_2->set_line_orientation(
5515  2, switch_1_line_orientations[2]);
5516  switch_2->set_line_orientation(
5517  3, switch_1_line_orientations[3]);
5518  switch_2->set_boundary_id_internal(
5519  switch_1_boundary_id);
5520  switch_2->set_manifold_id(switch_1->manifold_id());
5521  switch_2->set_user_index(switch_1_user_index);
5522  if (switch_1_user_flag)
5523  switch_2->set_user_flag();
5524  else
5525  switch_2->clear_user_flag();
5526  switch_2->clear_refinement_case();
5527  switch_2->set_refinement_case(
5528  switch_1_refinement_case);
5529  switch_2->clear_children();
5530  switch_2->set_children(0,
5531  switch_1_first_child_pair);
5532  switch_2->set_children(2,
5533  switch_1_second_child_pair);
5534 
5535  new_quads[0]->set_refinement_case(
5537  new_quads[0]->set_children(0, quad->child_index(0));
5538  new_quads[1]->set_refinement_case(
5540  new_quads[1]->set_children(0, quad->child_index(2));
5541  }
5542  else
5543  {
5544  new_quads[0]->set_refinement_case(
5546  new_quads[0]->set_children(0, quad->child_index(0));
5547  new_quads[1]->set_refinement_case(
5549  new_quads[1]->set_children(0, quad->child_index(2));
5550  new_line->set_children(
5551  0, quad->child(0)->line_index(3));
5552  Assert(new_line->child(1) ==
5553  quad->child(1)->line(3),
5554  ExcInternalError());
5555  }
5556  quad->clear_children();
5557  }
5558 
5559  // note these quads as children to the present one
5560  quad->set_children(0, new_quads[0]->index());
5561 
5562  quad->set_refinement_case(aniso_quad_ref_case);
5563 
5564  // finally clear flag indicating the need for
5565  // refinement
5566  quad->clear_user_data();
5567  } // if (anisotropic refinement)
5568 
5569  if (quad->user_flag_set())
5570  {
5571  // this quad needs to be refined isotropically
5572 
5573  // first of all: we only get here in the first run
5574  // of the loop
5575  Assert(loop == 0, ExcInternalError());
5576 
5577  // find the next unused vertex. we'll need this in
5578  // any case
5579  while (triangulation.vertices_used[next_unused_vertex] ==
5580  true)
5581  ++next_unused_vertex;
5582  Assert(
5583  next_unused_vertex < triangulation.vertices.size(),
5584  ExcMessage(
5585  "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."));
5586 
5587  // now: if the quad is refined anisotropically
5588  // already, set the anisotropic refinement flag
5589  // for both children. Additionally, we have to
5590  // refine the inner line, as it is an outer line
5591  // of the two (anisotropic) children
5592  const RefinementCase<dim - 1> quad_ref_case =
5593  quad->refinement_case();
5594 
5595  if (quad_ref_case == RefinementCase<dim - 1>::cut_x ||
5596  quad_ref_case == RefinementCase<dim - 1>::cut_y)
5597  {
5598  // set the 'opposite' refine case for children
5599  quad->child(0)->set_user_index(
5600  RefinementCase<dim - 1>::cut_xy - quad_ref_case);
5601  quad->child(1)->set_user_index(
5602  RefinementCase<dim - 1>::cut_xy - quad_ref_case);
5603  // refine the inner line
5605  middle_line;
5606  if (quad_ref_case == RefinementCase<dim - 1>::cut_x)
5607  middle_line = quad->child(0)->line(1);
5608  else
5609  middle_line = quad->child(0)->line(3);
5610 
5611  // if the face has been refined
5612  // anisotropically in the last refinement step
5613  // it might be, that it is flagged already and
5614  // that the middle line is thus refined
5615  // already. if not create children.
5616  if (!middle_line->has_children())
5617  {
5618  // set the middle vertex
5619  // appropriately. double refinement of
5620  // quads can only happen in the interior
5621  // of the domain, so we need not care
5622  // about boundary quads here
5623  triangulation.vertices[next_unused_vertex] =
5624  middle_line->center(true);
5625  triangulation.vertices_used[next_unused_vertex] =
5626  true;
5627 
5628  // now search a slot for the two
5629  // child lines
5630  next_unused_line =
5631  triangulation.faces->lines
5632  .template next_free_pair_object<1>(
5633  triangulation);
5634 
5635  // set the child pointer of the present
5636  // line
5637  middle_line->set_children(
5638  0, next_unused_line->index());
5639 
5640  // set the two new lines
5641  const typename Triangulation<dim, spacedim>::
5642  raw_line_iterator children[2] = {
5643  next_unused_line, ++next_unused_line};
5644 
5645  // some tests; if any of the iterators
5646  // should be invalid, then already
5647  // dereferencing will fail
5648  Assert(
5649  children[0]->used() == false,
5650  ExcMessage(
5651  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5652  Assert(
5653  children[1]->used() == false,
5654  ExcMessage(
5655  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5656 
5657  children[0]->set_bounding_object_indices(
5658  {middle_line->vertex_index(0),
5659  next_unused_vertex});
5660  children[1]->set_bounding_object_indices(
5661  {next_unused_vertex,
5662  middle_line->vertex_index(1)});
5663 
5664  children[0]->set_used_flag();
5665  children[1]->set_used_flag();
5666  children[0]->clear_children();
5667  children[1]->clear_children();
5668  children[0]->clear_user_data();
5669  children[1]->clear_user_data();
5670  children[0]->clear_user_flag();
5671  children[1]->clear_user_flag();
5672 
5673  children[0]->set_boundary_id_internal(
5674  middle_line->boundary_id());
5675  children[1]->set_boundary_id_internal(
5676  middle_line->boundary_id());
5677 
5678  children[0]->set_manifold_id(
5679  middle_line->manifold_id());
5680  children[1]->set_manifold_id(
5681  middle_line->manifold_id());
5682  }
5683  // now remove the flag from the quad and go to
5684  // the next quad, the actual refinement of the
5685  // quad takes place later on in this pass of
5686  // the loop or in the next one
5687  quad->clear_user_flag();
5688  continue;
5689  } // if (several refinement cases)
5690 
5691  // if we got here, we have an unrefined quad and
5692  // have to do the usual work like in an purely
5693  // isotropic refinement
5694  Assert(quad_ref_case ==
5696  ExcInternalError());
5697 
5698  // set the middle vertex appropriately: it might be that
5699  // the quad itself is not at the boundary, but that one of
5700  // its lines actually is. in this case, the newly created
5701  // vertices at the centers of the lines are not
5702  // necessarily the mean values of the adjacent vertices,
5703  // so do not compute the new vertex as the mean value of
5704  // the 4 vertices of the face, but rather as a weighted
5705  // mean value of the 8 vertices which we already have (the
5706  // four old ones, and the four ones inserted as middle
5707  // points for the four lines). summing up some more points
5708  // is generally cheaper than first asking whether one of
5709  // the lines is at the boundary
5710  //
5711  // note that the exact weights are chosen such as to
5712  // minimize the distortion of the four new quads from the
5713  // optimal shape. their description uses the formulas
5714  // underlying the TransfiniteInterpolationManifold
5715  // implementation
5716  triangulation.vertices[next_unused_vertex] =
5717  quad->center(true, true);
5718  triangulation.vertices_used[next_unused_vertex] = true;
5719 
5720  // now that we created the right point, make up
5721  // the four lines interior to the quad (++ takes
5722  // care of the end of the vector)
5724  new_lines[4];
5725 
5726  for (unsigned int i = 0; i < 4; ++i)
5727  {
5728  if (i % 2 == 0)
5729  // search a free pair of lines for 0. and
5730  // 2. line, so that two of them end up
5731  // together, which is necessary if later on
5732  // we want to refine the quad
5733  // anisotropically and the two lines end up
5734  // as children of new line
5735  next_unused_line =
5736  triangulation.faces->lines
5737  .template next_free_pair_object<1>(triangulation);
5738 
5739  new_lines[i] = next_unused_line;
5740  ++next_unused_line;
5741 
5742  Assert(
5743  new_lines[i]->used() == false,
5744  ExcMessage(
5745  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5746  }
5747 
5748  // set the data of the four lines. first collect
5749  // the indices of the five vertices:
5750  //
5751  // *--3--*
5752  // | | |
5753  // 0--4--1
5754  // | | |
5755  // *--2--*
5756  //
5757  // the lines are numbered as follows:
5758  //
5759  // *--*--*
5760  // | 1 |
5761  // *2-*-3*
5762  // | 0 |
5763  // *--*--*
5764 
5765  const unsigned int vertex_indices[5] = {
5766  quad->line(0)->child(0)->vertex_index(1),
5767  quad->line(1)->child(0)->vertex_index(1),
5768  quad->line(2)->child(0)->vertex_index(1),
5769  quad->line(3)->child(0)->vertex_index(1),
5770  next_unused_vertex};
5771 
5772  new_lines[0]->set_bounding_object_indices(
5773  {vertex_indices[2], vertex_indices[4]});
5774  new_lines[1]->set_bounding_object_indices(
5775  {vertex_indices[4], vertex_indices[3]});
5776  new_lines[2]->set_bounding_object_indices(
5777  {vertex_indices[0], vertex_indices[4]});
5778  new_lines[3]->set_bounding_object_indices(
5779  {vertex_indices[4], vertex_indices[1]});
5780 
5781  for (const auto &new_line : new_lines)
5782  {
5783  new_line->set_used_flag();
5784  new_line->clear_user_flag();
5785  new_line->clear_user_data();
5786  new_line->clear_children();
5787  new_line->set_boundary_id_internal(quad->boundary_id());
5788  new_line->set_manifold_id(quad->manifold_id());
5789  }
5790 
5791  // now for the quads. again, first collect some
5792  // data about the indices of the lines, with the
5793  // following numbering:
5794  //
5795  // .-6-.-7-.
5796  // 1 9 3
5797  // .-10.11-.
5798  // 0 8 2
5799  // .-4-.-5-.
5800 
5801  // child 0 and 1 of a line are switched if the
5802  // line orientation is false. set up a miniature
5803  // table, indicating which child to take for line
5804  // orientations false and true. first index: child
5805  // index in standard orientation, second index:
5806  // line orientation
5807  const unsigned int index[2][2] = {
5808  {1, 0}, // child 0, line_orientation=false and true
5809  {0, 1}}; // child 1, line_orientation=false and true
5810 
5811  const int line_indices[12] = {
5812  quad->line(0)
5813  ->child(index[0][quad->line_orientation(0)])
5814  ->index(),
5815  quad->line(0)
5816  ->child(index[1][quad->line_orientation(0)])
5817  ->index(),
5818  quad->line(1)
5819  ->child(index[0][quad->line_orientation(1)])
5820  ->index(),
5821  quad->line(1)
5822  ->child(index[1][quad->line_orientation(1)])
5823  ->index(),
5824  quad->line(2)
5825  ->child(index[0][quad->line_orientation(2)])
5826  ->index(),
5827  quad->line(2)
5828  ->child(index[1][quad->line_orientation(2)])
5829  ->index(),
5830  quad->line(3)
5831  ->child(index[0][quad->line_orientation(3)])
5832  ->index(),
5833  quad->line(3)
5834  ->child(index[1][quad->line_orientation(3)])
5835  ->index(),
5836  new_lines[0]->index(),
5837  new_lines[1]->index(),
5838  new_lines[2]->index(),
5839  new_lines[3]->index()};
5840 
5841  // find some space (consecutive)
5842  // for the first two newly to be
5843  // created quads.
5845  new_quads[4];
5846 
5847  next_unused_quad =
5848  triangulation.faces->quads
5849  .template next_free_pair_object<2>(triangulation);
5850 
5851  new_quads[0] = next_unused_quad;
5852  Assert(
5853  new_quads[0]->used() == false,
5854  ExcMessage(
5855  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5856 
5857  ++next_unused_quad;
5858  new_quads[1] = next_unused_quad;
5859  Assert(
5860  new_quads[1]->used() == false,
5861  ExcMessage(
5862  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5863 
5864  next_unused_quad =
5865  triangulation.faces->quads
5866  .template next_free_pair_object<2>(triangulation);
5867  new_quads[2] = next_unused_quad;
5868  Assert(
5869  new_quads[2]->used() == false,
5870  ExcMessage(
5871  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5872 
5873  ++next_unused_quad;
5874  new_quads[3] = next_unused_quad;
5875  Assert(
5876  new_quads[3]->used() == false,
5877  ExcMessage(
5878  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5879 
5880  // note these quads as children to the present one
5881  quad->set_children(0, new_quads[0]->index());
5882  quad->set_children(2, new_quads[2]->index());
5883  quad->set_refinement_case(RefinementCase<2>::cut_xy);
5884 
5885  new_quads[0]->set_bounding_object_indices(
5886  {line_indices[0],
5887  line_indices[8],
5888  line_indices[4],
5889  line_indices[10]});
5890  new_quads[1]->set_bounding_object_indices(
5891  {line_indices[8],
5892  line_indices[2],
5893  line_indices[5],
5894  line_indices[11]});
5895  new_quads[2]->set_bounding_object_indices(
5896  {line_indices[1],
5897  line_indices[9],
5898  line_indices[10],
5899  line_indices[6]});
5900  new_quads[3]->set_bounding_object_indices(
5901  {line_indices[9],
5902  line_indices[3],
5903  line_indices[11],
5904  line_indices[7]});
5905  for (const auto &new_quad : new_quads)
5906  {
5907  new_quad->set_used_flag();
5908  new_quad->clear_user_flag();
5909  new_quad->clear_user_data();
5910  new_quad->clear_children();
5911  new_quad->set_boundary_id_internal(quad->boundary_id());
5912  new_quad->set_manifold_id(quad->manifold_id());
5913  // set all line orientations to true, change
5914  // this after the loop, as we have to consider
5915  // different lines for each child
5916  for (unsigned int j = 0;
5917  j < GeometryInfo<dim>::lines_per_face;
5918  ++j)
5919  new_quad->set_line_orientation(j, true);
5920  }
5921  // now set the line orientation of children of
5922  // outer lines correctly, the lines in the
5923  // interior of the refined quad are automatically
5924  // oriented conforming to the standard
5925  new_quads[0]->set_line_orientation(
5926  0, quad->line_orientation(0));
5927  new_quads[0]->set_line_orientation(
5928  2, quad->line_orientation(2));
5929  new_quads[1]->set_line_orientation(
5930  1, quad->line_orientation(1));
5931  new_quads[1]->set_line_orientation(
5932  2, quad->line_orientation(2));
5933  new_quads[2]->set_line_orientation(
5934  0, quad->line_orientation(0));
5935  new_quads[2]->set_line_orientation(
5936  3, quad->line_orientation(3));
5937  new_quads[3]->set_line_orientation(
5938  1, quad->line_orientation(1));
5939  new_quads[3]->set_line_orientation(
5940  3, quad->line_orientation(3));
5941 
5942  // finally clear flag indicating the need for
5943  // refinement
5944  quad->clear_user_flag();
5945  } // if (isotropic refinement)
5946  } // for all quads
5947  } // looped two times over all quads, all quads refined now
5948 
5950  // Now, finally, set up the new
5951  // cells
5953 
5955  cells_with_distorted_children;
5956 
5957  for (unsigned int level = 0; level != triangulation.levels.size() - 1;
5958  ++level)
5959  {
5960  // only active objects can be refined further; remember
5961  // that we won't operate on the finest level, so
5962  // triangulation.begin_*(level+1) is allowed
5964  hex = triangulation.begin_active_hex(level),
5965  endh = triangulation.begin_active_hex(level + 1);
5967  next_unused_hex = triangulation.begin_raw_hex(level + 1);
5968 
5969  for (; hex != endh; ++hex)
5970  if (hex->refine_flag_set())
5971  {
5972  // this hex needs to be refined
5973 
5974  // clear flag indicating the need for refinement. do
5975  // it here already, since we can't do it anymore
5976  // once the cell has children
5977  const RefinementCase<dim> ref_case = hex->refine_flag_set();
5978  hex->clear_refine_flag();
5979  hex->set_refinement_case(ref_case);
5980 
5981  // depending on the refine case we might have to
5982  // create additional vertices, lines and quads
5983  // interior of the hex before the actual children
5984  // can be set up.
5985 
5986  // in a first step: reserve the needed space for
5987  // lines, quads and hexes and initialize them
5988  // correctly
5989 
5990  unsigned int n_new_lines = 0;
5991  unsigned int n_new_quads = 0;
5992  unsigned int n_new_hexes = 0;
5993  switch (ref_case)
5994  {
5998  n_new_lines = 0;
5999  n_new_quads = 1;
6000  n_new_hexes = 2;
6001  break;
6005  n_new_lines = 1;
6006  n_new_quads = 4;
6007  n_new_hexes = 4;
6008  break;
6010  n_new_lines = 6;
6011  n_new_quads = 12;
6012  n_new_hexes = 8;
6013  break;
6014  default:
6015  Assert(false, ExcInternalError());
6016  break;
6017  }
6018 
6019  // find some space for the newly to be created
6020  // interior lines and initialize them.
6021  std::vector<
6023  new_lines(n_new_lines);
6024  for (unsigned int i = 0; i < n_new_lines; ++i)
6025  {
6026  new_lines[i] =
6027  triangulation.faces->lines
6028  .template next_free_single_object<1>(triangulation);
6029 
6030  Assert(
6031  new_lines[i]->used() == false,
6032  ExcMessage(
6033  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6034  new_lines[i]->set_used_flag();
6035  new_lines[i]->clear_user_flag();
6036  new_lines[i]->clear_user_data();
6037  new_lines[i]->clear_children();
6038  // interior line
6039  new_lines[i]->set_boundary_id_internal(
6041  // they inherit geometry description of the hex they
6042  // belong to
6043  new_lines[i]->set_manifold_id(hex->manifold_id());
6044  }
6045 
6046  // find some space for the newly to be created
6047  // interior quads and initialize them.
6048  std::vector<
6050  new_quads(n_new_quads);
6051  for (unsigned int i = 0; i < n_new_quads; ++i)
6052  {
6053  new_quads[i] =
6054  triangulation.faces->quads
6055  .template next_free_single_object<2>(triangulation);
6056 
6057  Assert(
6058  new_quads[i]->used() == false,
6059  ExcMessage(
6060  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6061  new_quads[i]->set_used_flag();
6062  new_quads[i]->clear_user_flag();
6063  new_quads[i]->clear_user_data();
6064  new_quads[i]->clear_children();
6065  // interior quad
6066  new_quads[i]->set_boundary_id_internal(
6068  // they inherit geometry description of the hex they
6069  // belong to
6070  new_quads[i]->set_manifold_id(hex->manifold_id());
6071  // set all line orientation flags to true by
6072  // default, change this afterwards, if necessary
6073  for (unsigned int j = 0;
6074  j < GeometryInfo<dim>::lines_per_face;
6075  ++j)
6076  new_quads[i]->set_line_orientation(j, true);
6077  }
6078 
6079  types::subdomain_id subdomainid = hex->subdomain_id();
6080 
6081  // find some space for the newly to be created hexes
6082  // and initialize them.
6083  std::vector<
6085  new_hexes(n_new_hexes);
6086  for (unsigned int i = 0; i < n_new_hexes; ++i)
6087  {
6088  if (i % 2 == 0)
6089  next_unused_hex =
6090  triangulation.levels[level + 1]->cells.next_free_hex(
6091  triangulation, level + 1);
6092  else
6093  ++next_unused_hex;
6094 
6095  new_hexes[i] = next_unused_hex;
6096 
6097  Assert(
6098  new_hexes[i]->used() == false,
6099  ExcMessage(
6100  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6101  new_hexes[i]->set_used_flag();
6102  new_hexes[i]->clear_user_flag();
6103  new_hexes[i]->clear_user_data();
6104  new_hexes[i]->clear_children();
6105  // inherit material
6106  // properties
6107  new_hexes[i]->set_material_id(hex->material_id());
6108  new_hexes[i]->set_manifold_id(hex->manifold_id());
6109  new_hexes[i]->set_subdomain_id(subdomainid);
6110 
6111  if (i % 2)
6112  new_hexes[i]->set_parent(hex->index());
6113  // set the face_orientation flag to true for all
6114  // faces initially, as this is the default value
6115  // which is true for all faces interior to the
6116  // hex. later on go the other way round and
6117  // reset faces that are at the boundary of the
6118  // mother cube
6119  //
6120  // the same is true for the face_flip and
6121  // face_rotation flags. however, the latter two
6122  // are set to false by default as this is the
6123  // standard value
6124  for (const unsigned int f :
6126  {
6127  new_hexes[i]->set_face_orientation(f, true);
6128  new_hexes[i]->set_face_flip(f, false);
6129  new_hexes[i]->set_face_rotation(f, false);
6130  }
6131  }
6132  // note these hexes as children to the present cell
6133  for (unsigned int i = 0; i < n_new_hexes / 2; ++i)
6134  hex->set_children(2 * i, new_hexes[2 * i]->index());
6135 
6136  // we have to take into account whether the
6137  // different faces are oriented correctly or in the
6138  // opposite direction, so store that up front
6139 
6140  // face_orientation
6141  const bool f_or[6] = {hex->face_orientation(0),
6142  hex->face_orientation(1),
6143  hex->face_orientation(2),
6144  hex->face_orientation(3),
6145  hex->face_orientation(4),
6146  hex->face_orientation(5)};
6147 
6148  // face_flip
6149  const bool f_fl[6] = {hex->face_flip(0),
6150  hex->face_flip(1),
6151  hex->face_flip(2),
6152  hex->face_flip(3),
6153  hex->face_flip(4),
6154  hex->face_flip(5)};
6155 
6156  // face_rotation
6157  const bool f_ro[6] = {hex->face_rotation(0),
6158  hex->face_rotation(1),
6159  hex->face_rotation(2),
6160  hex->face_rotation(3),
6161  hex->face_rotation(4),
6162  hex->face_rotation(5)};
6163 
6164  // little helper table, indicating, whether the
6165  // child with index 0 or with index 1 can be found
6166  // at the standard origin of an anisotropically
6167  // refined quads in real orientation index 1:
6168  // (RefineCase - 1) index 2: face_flip
6169 
6170  // index 3: face rotation
6171  // note: face orientation has no influence
6172  const unsigned int child_at_origin[2][2][2] = {
6173  {{0, 0}, // RefinementCase<dim>::cut_x, face_flip=false,
6174  // face_rotation=false and true
6175  {1, 1}}, // RefinementCase<dim>::cut_x, face_flip=true,
6176  // face_rotation=false and true
6177  {{0, 1}, // RefinementCase<dim>::cut_y, face_flip=false,
6178  // face_rotation=false and true
6179  {1, 0}}}; // RefinementCase<dim>::cut_y, face_flip=true,
6180  // face_rotation=false and true
6181 
6183  //
6184  // in the following we will do the same thing for
6185  // each refinement case: create a new vertex (if
6186  // needed), create new interior lines (if needed),
6187  // create new interior quads and afterwards build
6188  // the children hexes out of these and the existing
6189  // subfaces of the outer quads (which have been
6190  // created above). However, even if the steps are
6191  // quite similar, the actual work strongly depends
6192  // on the actual refinement case. therefore, we use
6193  // separate blocks of code for each of these cases,
6194  // which hopefully increases the readability to some
6195  // extend.
6196 
6197  switch (ref_case)
6198  {
6200  {
6202  //
6203  // RefinementCase<dim>::cut_x
6204  //
6205  // the refined cube will look
6206  // like this:
6207  //
6208  // *----*----*
6209  // / / /|
6210  // / / / |
6211  // / / / |
6212  // *----*----* |
6213  // | | | |
6214  // | | | *
6215  // | | | /
6216  // | | | /
6217  // | | |/
6218  // *----*----*
6219  //
6220  // again, first collect some data about the
6221  // indices of the lines, with the following
6222  // numbering:
6223 
6224  // face 2: front plane
6225  // (note: x,y exchanged)
6226  // *---*---*
6227  // | | |
6228  // | 0 |
6229  // | | |
6230  // *---*---*
6231  // m0
6232  // face 3: back plane
6233  // (note: x,y exchanged)
6234  // m1
6235  // *---*---*
6236  // | | |
6237  // | 1 |
6238  // | | |
6239  // *---*---*
6240  // face 4: bottom plane
6241  // *---*---*
6242  // / / /
6243  // / 2 /
6244  // / / /
6245  // *---*---*
6246  // m0
6247  // face 5: top plane
6248  // m1
6249  // *---*---*
6250  // / / /
6251  // / 3 /
6252  // / / /
6253  // *---*---*
6254 
6255  // set up a list of line iterators first. from
6256  // this, construct lists of line_indices and
6257  // line orientations later on
6258  const typename Triangulation<dim, spacedim>::
6259  raw_line_iterator lines[4] = {
6260  hex->face(2)->child(0)->line(
6261  (hex->face(2)->refinement_case() ==
6263  1 :
6264  3), // 0
6265  hex->face(3)->child(0)->line(
6266  (hex->face(3)->refinement_case() ==
6268  1 :
6269  3), // 1
6270  hex->face(4)->child(0)->line(
6271  (hex->face(4)->refinement_case() ==
6273  1 :
6274  3), // 2
6275  hex->face(5)->child(0)->line(
6276  (hex->face(5)->refinement_case() ==
6278  1 :
6279  3) // 3
6280  };
6281 
6282  unsigned int line_indices[4];
6283  for (unsigned int i = 0; i < 4; ++i)
6284  line_indices[i] = lines[i]->index();
6285 
6286  // the orientation of lines for the inner quads
6287  // is quite tricky. as these lines are newly
6288  // created ones and thus have no parents, they
6289  // cannot inherit this property. set up an array
6290  // and fill it with the respective values
6291  bool line_orientation[4];
6292 
6293  // the middle vertex marked as m0 above is the
6294  // start vertex for lines 0 and 2 in standard
6295  // orientation, whereas m1 is the end vertex of
6296  // lines 1 and 3 in standard orientation
6297  const unsigned int middle_vertices[2] = {
6298  hex->line(2)->child(0)->vertex_index(1),
6299  hex->line(7)->child(0)->vertex_index(1)};
6300 
6301  for (unsigned int i = 0; i < 4; ++i)
6302  if (lines[i]->vertex_index(i % 2) ==
6303  middle_vertices[i % 2])
6304  line_orientation[i] = true;
6305  else
6306  {
6307  // it must be the other
6308  // way round then
6309  Assert(lines[i]->vertex_index((i + 1) % 2) ==
6310  middle_vertices[i % 2],
6311  ExcInternalError());
6312  line_orientation[i] = false;
6313  }
6314 
6315  // set up the new quad, line numbering is as
6316  // indicated above
6317  new_quads[0]->set_bounding_object_indices(
6318  {line_indices[0],
6319  line_indices[1],
6320  line_indices[2],
6321  line_indices[3]});
6322 
6323  new_quads[0]->set_line_orientation(
6324  0, line_orientation[0]);
6325  new_quads[0]->set_line_orientation(
6326  1, line_orientation[1]);
6327  new_quads[0]->set_line_orientation(
6328  2, line_orientation[2]);
6329  new_quads[0]->set_line_orientation(
6330  3, line_orientation[3]);
6331 
6332  // the quads are numbered as follows:
6333  //
6334  // planes in the interior of the old hex:
6335  //
6336  // *
6337  // /|
6338  // / | x
6339  // / | *-------* *---------*
6340  // * | | | / /
6341  // | 0 | | | / /
6342  // | * | | / /
6343  // | / *-------*y *---------*x
6344  // | /
6345  // |/
6346  // *
6347  //
6348  // children of the faces of the old hex
6349  //
6350  // *---*---* *---*---*
6351  // /| | | / / /|
6352  // / | | | / 9 / 10/ |
6353  // / | 5 | 6 | / / / |
6354  // * | | | *---*---* |
6355  // | 1 *---*---* | | | 2 *
6356  // | / / / | | | /
6357  // | / 7 / 8 / | 3 | 4 | /
6358  // |/ / / | | |/
6359  // *---*---* *---*---*
6360  //
6361  // note that we have to take care of the
6362  // orientation of faces.
6363  const int quad_indices[11] = {
6364  new_quads[0]->index(), // 0
6365 
6366  hex->face(0)->index(), // 1
6367 
6368  hex->face(1)->index(), // 2
6369 
6370  hex->face(2)->child_index(
6371  child_at_origin[hex->face(2)->refinement_case() -
6372  1][f_fl[2]][f_ro[2]]), // 3
6373  hex->face(2)->child_index(
6374  1 -
6375  child_at_origin[hex->face(2)->refinement_case() -
6376  1][f_fl[2]][f_ro[2]]),
6377 
6378  hex->face(3)->child_index(
6379  child_at_origin[hex->face(3)->refinement_case() -
6380  1][f_fl[3]][f_ro[3]]), // 5
6381  hex->face(3)->child_index(
6382  1 -
6383  child_at_origin[hex->face(3)->refinement_case() -
6384  1][f_fl[3]][f_ro[3]]),
6385 
6386  hex->face(4)->child_index(
6387  child_at_origin[hex->face(4)->refinement_case() -
6388  1][f_fl[4]][f_ro[4]]), // 7
6389  hex->face(4)->child_index(
6390  1 -
6391  child_at_origin[hex->face(4)->refinement_case() -
6392  1][f_fl[4]][f_ro[4]]),
6393 
6394  hex->face(5)->child_index(
6395  child_at_origin[hex->face(5)->refinement_case() -
6396  1][f_fl[5]][f_ro[5]]), // 9
6397  hex->face(5)->child_index(
6398  1 -
6399  child_at_origin[hex->face(5)->refinement_case() -
6400  1][f_fl[5]][f_ro[5]])
6401 
6402  };
6403 
6404  new_hexes[0]->set_bounding_object_indices(
6405  {quad_indices[1],
6406  quad_indices[0],
6407  quad_indices[3],
6408  quad_indices[5],
6409  quad_indices[7],
6410  quad_indices[9]});
6411  new_hexes[1]->set_bounding_object_indices(
6412  {quad_indices[0],
6413  quad_indices[2],
6414  quad_indices[4],
6415  quad_indices[6],
6416  quad_indices[8],
6417  quad_indices[10]});
6418  break;
6419  }
6420 
6422  {
6424  //
6425  // RefinementCase<dim>::cut_y
6426  //
6427  // the refined cube will look like this:
6428  //
6429  // *---------*
6430  // / /|
6431  // *---------* |
6432  // / /| |
6433  // *---------* | |
6434  // | | | |
6435  // | | | *
6436  // | | |/
6437  // | | *
6438  // | |/
6439  // *---------*
6440  //
6441  // again, first collect some data about the
6442  // indices of the lines, with the following
6443  // numbering:
6444 
6445  // face 0: left plane
6446  // *
6447  // /|
6448  // * |
6449  // /| |
6450  // * | |
6451  // | 0 |
6452  // | | *
6453  // | |/
6454  // | *m0
6455  // |/
6456  // *
6457  // face 1: right plane
6458  // *
6459  // /|
6460  // m1* |
6461  // /| |
6462  // * | |
6463  // | 1 |
6464  // | | *
6465  // | |/
6466  // | *
6467  // |/
6468  // *
6469  // face 4: bottom plane
6470  // *-------*
6471  // / /
6472  // m0*---2---*
6473  // / /
6474  // *-------*
6475  // face 5: top plane
6476  // *-------*
6477  // / /
6478  // *---3---*m1
6479  // / /
6480  // *-------*
6481 
6482  // set up a list of line iterators first. from
6483  // this, construct lists of line_indices and
6484  // line orientations later on
6485  const typename Triangulation<dim, spacedim>::
6486  raw_line_iterator lines[4] = {
6487  hex->face(0)->child(0)->line(
6488  (hex->face(0)->refinement_case() ==
6490  1 :
6491  3), // 0
6492  hex->face(1)->child(0)->line(
6493  (hex->face(1)->refinement_case() ==
6495  1 :
6496  3), // 1
6497  hex->face(4)->child(0)->line(
6498  (hex->face(4)->refinement_case() ==
6500  1 :
6501  3), // 2
6502  hex->face(5)->child(0)->line(
6503  (hex->face(5)->refinement_case() ==
6505  1 :
6506  3) // 3
6507  };
6508 
6509  unsigned int line_indices[4];
6510  for (unsigned int i = 0; i < 4; ++i)
6511  line_indices[i] = lines[i]->index();
6512 
6513  // the orientation of lines for the inner quads
6514  // is quite tricky. as these lines are newly
6515  // created ones and thus have no parents, they
6516  // cannot inherit this property. set up an array
6517  // and fill it with the respective values
6518  bool line_orientation[4];
6519 
6520  // the middle vertex marked as m0 above is the
6521  // start vertex for lines 0 and 2 in standard
6522  // orientation, whereas m1 is the end vertex of
6523  // lines 1 and 3 in standard orientation
6524  const unsigned int middle_vertices[2] = {
6525  hex->line(0)->child(0)->vertex_index(1),
6526  hex->line(5)->child(0)->vertex_index(1)};
6527 
6528  for (unsigned int i = 0; i < 4; ++i)
6529  if (lines[i]->vertex_index(i % 2) ==
6530  middle_vertices[i % 2])
6531  line_orientation[i] = true;
6532  else
6533  {
6534  // it must be the other way round then
6535  Assert(lines[i]->vertex_index((i + 1) % 2) ==
6536  middle_vertices[i % 2],
6537  ExcInternalError());
6538  line_orientation[i] = false;
6539  }
6540 
6541  // set up the new quad, line numbering is as
6542  // indicated above
6543  new_quads[0]->set_bounding_object_indices(
6544  {line_indices[2],
6545  line_indices[3],
6546  line_indices[0],
6547  line_indices[1]});
6548 
6549  new_quads[0]->set_line_orientation(
6550  0, line_orientation[2]);
6551  new_quads[0]->set_line_orientation(
6552  1, line_orientation[3]);
6553  new_quads[0]->set_line_orientation(
6554  2, line_orientation[0]);
6555  new_quads[0]->set_line_orientation(
6556  3, line_orientation[1]);
6557 
6558  // the quads are numbered as follows:
6559  //
6560  // planes in the interior of the old hex:
6561  //
6562  // *
6563  // /|
6564  // / | x
6565  // / | *-------* *---------*
6566  // * | | | / /
6567  // | | | 0 | / /
6568  // | * | | / /
6569  // | / *-------*y *---------*x
6570  // | /
6571  // |/
6572  // *
6573  //
6574  // children of the faces of the old hex
6575  //
6576  // *-------* *-------*
6577  // /| | / 10 /|
6578  // * | | *-------* |
6579  // /| | 6 | / 9 /| |
6580  // * |2| | *-------* |4|
6581  // | | *-------* | | | *
6582  // |1|/ 8 / | |3|/
6583  // | *-------* | 5 | *
6584  // |/ 7 / | |/
6585  // *-------* *-------*
6586  //
6587  // note that we have to take care of the
6588  // orientation of faces.
6589  const int quad_indices[11] = {
6590  new_quads[0]->index(), // 0
6591 
6592  hex->face(0)->child_index(
6593  child_at_origin[hex->face(0)->refinement_case() -
6594  1][f_fl[0]][f_ro[0]]), // 1
6595  hex->face(0)->child_index(
6596  1 -
6597  child_at_origin[hex->face(0)->refinement_case() -
6598  1][f_fl[0]][f_ro[0]]),
6599 
6600  hex->face(1)->child_index(
6601  child_at_origin[hex->face(1)->refinement_case() -
6602  1][f_fl[1]][f_ro[1]]), // 3
6603  hex->face(1)->child_index(
6604  1 -
6605  child_at_origin[hex->face(1)->refinement_case() -
6606  1][f_fl[1]][f_ro[1]]),
6607 
6608  hex->face(2)->index(), // 5
6609 
6610  hex->face(3)->index(), // 6
6611 
6612  hex->face(4)->child_index(
6613  child_at_origin[hex->face(4)->refinement_case() -
6614  1][f_fl[4]][f_ro[4]]), // 7
6615  hex->face(4)->child_index(
6616  1 -
6617  child_at_origin[hex->face(4)->refinement_case() -
6618  1][f_fl[4]][f_ro[4]]),
6619 
6620  hex->face(5)->child_index(
6621  child_at_origin[hex->face(5)->refinement_case() -
6622  1][f_fl[5]][f_ro[5]]), // 9
6623  hex->face(5)->child_index(
6624  1 -
6625  child_at_origin[hex->face(5)->refinement_case() -
6626  1][f_fl[5]][f_ro[5]])
6627 
6628  };
6629 
6630  new_hexes[0]->set_bounding_object_indices(
6631  {quad_indices[1],
6632  quad_indices[3],
6633  quad_indices[5],
6634  quad_indices[0],
6635  quad_indices[7],
6636  quad_indices[9]});
6637  new_hexes[1]->set_bounding_object_indices(
6638  {quad_indices[2],
6639  quad_indices[4],
6640  quad_indices[0],
6641  quad_indices[6],
6642  quad_indices[8],
6643  quad_indices[10]});
66