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