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