Reference documentation for deal.II version GIT ffb4c3937f 2023-03-31 14:25:02+00:00
\(\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 - 2022 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 
20 
24 #include <deal.II/grid/manifold.h>
25 #include <deal.II/grid/tria.h>
30 
31 #include <algorithm>
32 #include <array>
33 #include <cmath>
34 #include <functional>
35 #include <list>
36 #include <map>
37 #include <memory>
38 #include <numeric>
39 
40 
42 
43 
44 namespace internal
45 {
46  namespace TriangulationImplementation
47  {
49  : n_levels(0)
50  , n_lines(0)
51  , n_active_lines(0)
52  // all other fields are
53  // default constructed
54  {}
55 
56 
57 
58  std::size_t
60  {
61  std::size_t mem =
66  MemoryConsumption::memory_consumption(n_active_lines_level);
67 
68  if (active_cell_index_partitioner)
69  mem += active_cell_index_partitioner->memory_consumption();
70 
71  for (const auto &partitioner : level_cell_index_partitioners)
72  if (partitioner)
73  mem += partitioner->memory_consumption();
74 
75  return mem;
76  }
77 
78 
80  : n_quads(0)
81  , n_active_quads(0)
82  // all other fields are
83  // default constructed
84  {}
85 
86 
87 
88  std::size_t
90  {
95  MemoryConsumption::memory_consumption(n_active_quads_level));
96  }
97 
98 
99 
101  : n_hexes(0)
102  , n_active_hexes(0)
103  // all other fields are
104  // default constructed
105  {}
106 
107 
108 
109  std::size_t
111  {
115  MemoryConsumption::memory_consumption(n_active_hexes) +
116  MemoryConsumption::memory_consumption(n_active_hexes_level));
117  }
118  } // namespace TriangulationImplementation
119 } // namespace internal
120 
121 // anonymous namespace for internal helper functions
122 namespace
123 {
124  // return whether the given cell is
125  // patch_level_1, i.e. determine
126  // whether either all or none of
127  // its children are further
128  // refined. this function can only
129  // be called for non-active cells.
130  template <int dim, int spacedim>
131  bool
132  cell_is_patch_level_1(
134  {
135  Assert(cell->is_active() == false, ExcInternalError());
136 
137  unsigned int n_active_children = 0;
138  for (unsigned int i = 0; i < cell->n_children(); ++i)
139  if (cell->child(i)->is_active())
140  ++n_active_children;
141 
142  return (n_active_children == 0) ||
143  (n_active_children == cell->n_children());
144  }
145 
146 
147 
148  // return, whether a given @p cell will be
149  // coarsened, which is the case if all
150  // children are active and have their coarsen
151  // flag set. In case only part of the coarsen
152  // flags are set, remove them.
153  template <int dim, int spacedim>
154  bool
155  cell_will_be_coarsened(
157  {
158  // only cells with children should be
159  // considered for coarsening
160 
161  if (cell->has_children())
162  {
163  unsigned int children_to_coarsen = 0;
164  const unsigned int n_children = cell->n_children();
165 
166  for (unsigned int c = 0; c < n_children; ++c)
167  if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
168  ++children_to_coarsen;
169  if (children_to_coarsen == n_children)
170  return true;
171  else
172  for (unsigned int c = 0; c < n_children; ++c)
173  if (cell->child(c)->is_active())
174  cell->child(c)->clear_coarsen_flag();
175  }
176  // no children, so no coarsening
177  // possible. however, no children also
178  // means that this cell will be in the same
179  // state as if it had children and was
180  // coarsened. So, what should we return -
181  // false or true?
182  // make sure we do not have to do this at
183  // all...
184  Assert(cell->has_children(), ExcInternalError());
185  // ... and then simply return false
186  return false;
187  }
188 
189 
190  // return, whether the face @p face_no of the
191  // given @p cell will be refined after the
192  // current refinement step, considering
193  // refine and coarsen flags and considering
194  // only those refinemnts that will be caused
195  // by the neighboring cell.
196 
197  // this function is used on both active cells
198  // and cells with children. on cells with
199  // children it also of interest to know 'how'
200  // the face will be refined. thus there is an
201  // additional third argument @p
202  // expected_face_ref_case returning just
203  // that. be aware, that this variable will
204  // only contain useful information if this
205  // function is called for an active cell.
206  //
207  // thus, this is an internal function, users
208  // should call one of the two alternatives
209  // following below.
210  template <int dim, int spacedim>
211  bool
212  face_will_be_refined_by_neighbor_internal(
214  const unsigned int face_no,
215  RefinementCase<dim - 1> &expected_face_ref_case)
216  {
217  // first of all: set the default value for
218  // expected_face_ref_case, which is no
219  // refinement at all
220  expected_face_ref_case = RefinementCase<dim - 1>::no_refinement;
221 
222  const typename Triangulation<dim, spacedim>::cell_iterator neighbor =
223  cell->neighbor(face_no);
224 
225  // If we are at the boundary, there is no
226  // neighbor which could refine the face
227  if (neighbor.state() != IteratorState::valid)
228  return false;
229 
230  if (neighbor->has_children())
231  {
232  // if the neighbor is refined, it may be
233  // coarsened. if so, then it won't refine
234  // the face, no matter what else happens
235  if (cell_will_be_coarsened(neighbor))
236  return false;
237  else
238  // if the neighbor is refined, then it
239  // is also refined at our current
240  // face. It will stay so without
241  // coarsening, so return true in that
242  // case.
243  {
244  expected_face_ref_case = cell->face(face_no)->refinement_case();
245  return true;
246  }
247  }
248 
249  // now, the neighbor is not refined, but
250  // perhaps it will be
251  const RefinementCase<dim> nb_ref_flag = neighbor->refine_flag_set();
252  if (nb_ref_flag != RefinementCase<dim>::no_refinement)
253  {
254  // now we need to know, which of the
255  // neighbors faces points towards us
256  const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
257  // check, whether the cell will be
258  // refined in a way that refines our
259  // face
260  const RefinementCase<dim - 1> face_ref_case =
262  nb_ref_flag,
263  neighbor_neighbor,
264  neighbor->face_orientation(neighbor_neighbor),
265  neighbor->face_flip(neighbor_neighbor),
266  neighbor->face_rotation(neighbor_neighbor));
267  if (face_ref_case != RefinementCase<dim - 1>::no_refinement)
268  {
270  neighbor_face = neighbor->face(neighbor_neighbor);
271  const int this_face_index = cell->face_index(face_no);
272 
273  // there are still two basic
274  // possibilities here: the neighbor
275  // might be coarser or as coarse
276  // as we are
277  if (neighbor_face->index() == this_face_index)
278  // the neighbor is as coarse as
279  // we are and will be refined at
280  // the face of consideration, so
281  // return true
282  {
283  expected_face_ref_case = face_ref_case;
284  return true;
285  }
286  else
287  {
288  // the neighbor is coarser.
289  // this is the most complicated
290  // case. It might be, that the
291  // neighbor's face will be
292  // refined, but that we will
293  // not see this, as we are
294  // refined in a similar way.
295 
296  // so, the neighbor's face must
297  // have children. check, if our
298  // cell's face is one of these
299  // (it could also be a
300  // grand_child)
301  for (unsigned int c = 0; c < neighbor_face->n_children(); ++c)
302  if (neighbor_face->child_index(c) == this_face_index)
303  {
304  // if the flagged refine
305  // case of the face is a
306  // subset or the same as
307  // the current refine case,
308  // then the face, as seen
309  // from our cell, won't be
310  // refined by the neighbor
311  if ((neighbor_face->refinement_case() | face_ref_case) ==
312  neighbor_face->refinement_case())
313  return false;
314  else
315  {
316  // if we are active, we
317  // must be an
318  // anisotropic child
319  // and the coming
320  // face_ref_case is
321  // isotropic. Thus,
322  // from our cell we
323  // will see exactly the
324  // opposite refine case
325  // that the face has
326  // now...
327  Assert(
328  face_ref_case ==
330  ExcInternalError());
331  expected_face_ref_case =
332  ~neighbor_face->refinement_case();
333  return true;
334  }
335  }
336 
337  // so, obviously we were not
338  // one of the children, but a
339  // grandchild. This is only
340  // possible in 3d.
341  Assert(dim == 3, ExcInternalError());
342  // In that case, however, no
343  // matter what the neighbor
344  // does, it won't be finer
345  // after the next refinement
346  // step.
347  return false;
348  }
349  } // if face will be refined
350  } // if neighbor is flagged for refinement
351 
352  // no cases left, so the neighbor will not
353  // refine the face
354  return false;
355  }
356 
357  // version of above function for both active
358  // and non-active cells
359  template <int dim, int spacedim>
360  bool
361  face_will_be_refined_by_neighbor(
363  const unsigned int face_no)
364  {
365  RefinementCase<dim - 1> dummy = RefinementCase<dim - 1>::no_refinement;
366  return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
367  }
368 
369  // version of above function for active cells
370  // only. Additionally returning the refine
371  // case (to come) of the face under
372  // consideration
373  template <int dim, int spacedim>
374  bool
375  face_will_be_refined_by_neighbor(
377  const unsigned int face_no,
378  RefinementCase<dim - 1> &expected_face_ref_case)
379  {
380  return face_will_be_refined_by_neighbor_internal(cell,
381  face_no,
382  expected_face_ref_case);
383  }
384 
385 
386 
387  template <int dim, int spacedim>
388  bool
389  satisfies_level1_at_vertex_rule(
391  {
392  std::vector<unsigned int> min_adjacent_cell_level(
393  triangulation.n_vertices(), triangulation.n_levels());
394  std::vector<unsigned int> max_adjacent_cell_level(
395  triangulation.n_vertices(), 0);
396 
397  for (const auto &cell : triangulation.active_cell_iterators())
398  for (const unsigned int v : cell->vertex_indices())
399  {
400  min_adjacent_cell_level[cell->vertex_index(v)] =
401  std::min<unsigned int>(
402  min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
403  max_adjacent_cell_level[cell->vertex_index(v)] =
404  std::max<unsigned int>(
405  min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
406  }
407 
408  for (unsigned int k = 0; k < triangulation.n_vertices(); ++k)
409  if (triangulation.vertex_used(k))
410  if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
411  return false;
412  return true;
413  }
414 
415 
416 
428  void
429  reorder_compatibility(const std::vector<CellData<1>> &, const SubCellData &)
430  {
431  // nothing to do here: the format
432  // hasn't changed for 1d
433  }
434 
435 
436  void
437  reorder_compatibility(std::vector<CellData<2>> &cells, const SubCellData &)
438  {
439  for (auto &cell : cells)
440  if (cell.vertices.size() == GeometryInfo<2>::vertices_per_cell)
441  std::swap(cell.vertices[2], cell.vertices[3]);
442  }
443 
444 
445  void
446  reorder_compatibility(std::vector<CellData<3>> &cells,
447  SubCellData & subcelldata)
448  {
449  unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
450  static constexpr std::array<unsigned int,
452  local_vertex_numbering{{0, 1, 5, 4, 2, 3, 7, 6}};
453  for (auto &cell : cells)
454  if (cell.vertices.size() == GeometryInfo<3>::vertices_per_cell)
455  {
456  for (const unsigned int i : GeometryInfo<3>::vertex_indices())
457  tmp[i] = cell.vertices[i];
458  for (const unsigned int i : GeometryInfo<3>::vertex_indices())
459  cell.vertices[local_vertex_numbering[i]] = tmp[i];
460  }
461 
462  // now points in boundary quads
463  for (auto &boundary_quad : subcelldata.boundary_quads)
464  if (boundary_quad.vertices.size() == GeometryInfo<2>::vertices_per_cell)
465  std::swap(boundary_quad.vertices[2], boundary_quad.vertices[3]);
466  }
467 
468 
469 
487  template <int dim, int spacedim>
488  unsigned int
489  middle_vertex_index(
490  const typename Triangulation<dim, spacedim>::line_iterator &line)
491  {
492  if (line->has_children())
493  return line->child(0)->vertex_index(1);
495  }
496 
497 
498  template <int dim, int spacedim>
499  unsigned int
500  middle_vertex_index(
501  const typename Triangulation<dim, spacedim>::quad_iterator &quad)
502  {
503  switch (static_cast<unsigned char>(quad->refinement_case()))
504  {
506  return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
507  break;
509  return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
510  break;
512  return quad->child(0)->vertex_index(3);
513  break;
514  default:
515  break;
516  }
518  }
519 
520 
521  template <int dim, int spacedim>
522  unsigned int
523  middle_vertex_index(
524  const typename Triangulation<dim, spacedim>::hex_iterator &hex)
525  {
526  switch (static_cast<unsigned char>(hex->refinement_case()))
527  {
529  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
530  break;
532  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
533  break;
535  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
536  break;
538  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
539  break;
541  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
542  break;
544  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
545  break;
547  return hex->child(0)->vertex_index(7);
548  break;
549  default:
550  break;
551  }
553  }
554 
555 
568  template <class TRIANGULATION>
569  inline typename TRIANGULATION::DistortedCellList
570  collect_distorted_coarse_cells(const TRIANGULATION &)
571  {
572  return typename TRIANGULATION::DistortedCellList();
573  }
574 
575 
576 
585  template <int dim>
587  collect_distorted_coarse_cells(const Triangulation<dim, dim> &triangulation)
588  {
589  typename Triangulation<dim, dim>::DistortedCellList distorted_cells;
590  for (const auto &cell : triangulation.cell_iterators_on_level(0))
591  {
593  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
594  vertices[i] = cell->vertex(i);
595 
598 
599  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
600  if (determinants[i] <= 1e-9 * std::pow(cell->diameter(), 1. * dim))
601  {
602  distorted_cells.distorted_cells.push_back(cell);
603  break;
604  }
605  }
606 
607  return distorted_cells;
608  }
609 
610 
617  template <int dim>
618  bool
619  has_distorted_children(
620  const typename Triangulation<dim, dim>::cell_iterator &cell)
621  {
622  Assert(cell->has_children(), ExcInternalError());
623 
624  for (unsigned int c = 0; c < cell->n_children(); ++c)
625  {
627  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
628  vertices[i] = cell->child(c)->vertex(i);
629 
632 
633  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
634  if (determinants[i] <=
635  1e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
636  return true;
637  }
638 
639  return false;
640  }
641 
642 
650  template <int dim, int spacedim>
651  bool
652  has_distorted_children(
654  {
655  return false;
656  }
657 
658 
659  template <int dim, int spacedim>
660  void
661  update_periodic_face_map_recursively(
662  const typename Triangulation<dim, spacedim>::cell_iterator &cell_1,
663  const typename Triangulation<dim, spacedim>::cell_iterator &cell_2,
664  unsigned int n_face_1,
665  unsigned int n_face_2,
666  const std::bitset<3> & orientation,
667  typename std::map<
669  unsigned int>,
670  std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
671  unsigned int>,
672  std::bitset<3>>> &periodic_face_map)
673  {
674  using FaceIterator = typename Triangulation<dim, spacedim>::face_iterator;
675  const FaceIterator face_1 = cell_1->face(n_face_1);
676  const FaceIterator face_2 = cell_2->face(n_face_2);
677 
678  const bool face_orientation = orientation[0];
679  const bool face_flip = orientation[1];
680  const bool face_rotation = orientation[2];
681 
682  Assert((dim != 1) || (face_orientation == true && face_flip == false &&
683  face_rotation == false),
684  ExcMessage("The supplied orientation "
685  "(face_orientation, face_flip, face_rotation) "
686  "is invalid for 1d"));
687 
688  Assert((dim != 2) || (face_orientation == true && face_rotation == false),
689  ExcMessage("The supplied orientation "
690  "(face_orientation, face_flip, face_rotation) "
691  "is invalid for 2d"));
692 
693  Assert(face_1 != face_2, ExcMessage("face_1 and face_2 are equal!"));
694 
695  Assert(face_1->at_boundary() && face_2->at_boundary(),
696  ExcMessage("Periodic faces must be on the boundary"));
697 
698  // Check if the requirement that each edge can only have at most one hanging
699  // node, and as a consequence neighboring cells can differ by at most
700  // one refinement level is enforced. In 1d, there are no hanging nodes and
701  // so neighboring cells can differ by more than one refinement level.
702  Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
703  ExcInternalError());
704 
705  // insert periodic face pair for both cells
706  using CellFace =
707  std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
708  unsigned int>;
709  const CellFace cell_face_1(cell_1, n_face_1);
710  const CellFace cell_face_2(cell_2, n_face_2);
711  const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
712  cell_face_2, orientation);
713 
714  const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
715  periodic_faces(cell_face_1, cell_face_orientation_2);
716 
717  // Only one periodic neighbor is allowed
718  Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
719  periodic_face_map.insert(periodic_faces);
720 
721  if (dim == 1)
722  {
723  if (cell_1->has_children())
724  {
725  if (cell_2->has_children())
726  {
727  update_periodic_face_map_recursively<dim, spacedim>(
728  cell_1->child(n_face_1),
729  cell_2->child(n_face_2),
730  n_face_1,
731  n_face_2,
732  orientation,
733  periodic_face_map);
734  }
735  else // only face_1 has children
736  {
737  update_periodic_face_map_recursively<dim, spacedim>(
738  cell_1->child(n_face_1),
739  cell_2,
740  n_face_1,
741  n_face_2,
742  orientation,
743  periodic_face_map);
744  }
745  }
746  }
747  else // dim == 2 || dim == 3
748  {
749  // A lookup table on how to go through the child cells depending on the
750  // orientation:
751  // see Documentation of GeometryInfo for details
752 
753  static const int lookup_table_2d[2][2] =
754  // flip:
755  {
756  {0, 1}, // false
757  {1, 0} // true
758  };
759 
760  static const int lookup_table_3d[2][2][2][4] =
761  // orientation flip rotation
762  {{{
763  {0, 2, 1, 3}, // false false false
764  {2, 3, 0, 1} // false false true
765  },
766  {
767  {3, 1, 2, 0}, // false true false
768  {1, 0, 3, 2} // false true true
769  }},
770  {{
771  {0, 1, 2, 3}, // true false false
772  {1, 3, 0, 2} // true false true
773  },
774  {
775  {3, 2, 1, 0}, // true true false
776  {2, 0, 3, 1} // true true true
777  }}};
778 
779  if (cell_1->has_children())
780  {
781  if (cell_2->has_children())
782  {
783  // In the case that both faces have children, we loop over all
784  // children and apply update_periodic_face_map_recursively
785  // recursively:
786 
787  Assert(face_1->n_children() ==
789  face_2->n_children() ==
792 
793  for (unsigned int i = 0;
794  i < GeometryInfo<dim>::max_children_per_face;
795  ++i)
796  {
797  // Lookup the index for the second face
798  unsigned int j = 0;
799  switch (dim)
800  {
801  case 2:
802  j = lookup_table_2d[face_flip][i];
803  break;
804  case 3:
805  j = lookup_table_3d[face_orientation][face_flip]
806  [face_rotation][i];
807  break;
808  default:
809  AssertThrow(false, ExcNotImplemented());
810  }
811 
812  // find subcell ids that belong to the subface indices
813  unsigned int child_cell_1 =
815  cell_1->refinement_case(),
816  n_face_1,
817  i,
818  cell_1->face_orientation(n_face_1),
819  cell_1->face_flip(n_face_1),
820  cell_1->face_rotation(n_face_1),
821  face_1->refinement_case());
822  unsigned int child_cell_2 =
824  cell_2->refinement_case(),
825  n_face_2,
826  j,
827  cell_2->face_orientation(n_face_2),
828  cell_2->face_flip(n_face_2),
829  cell_2->face_rotation(n_face_2),
830  face_2->refinement_case());
831 
832  Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
833  face_1->child(i),
834  ExcInternalError());
835  Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
836  face_2->child(j),
837  ExcInternalError());
838 
839  // precondition: subcell has the same orientation as cell
840  // (so that the face numbers coincide) recursive call
841  update_periodic_face_map_recursively<dim, spacedim>(
842  cell_1->child(child_cell_1),
843  cell_2->child(child_cell_2),
844  n_face_1,
845  n_face_2,
846  orientation,
847  periodic_face_map);
848  }
849  }
850  else // only face_1 has children
851  {
852  for (unsigned int i = 0;
853  i < GeometryInfo<dim>::max_children_per_face;
854  ++i)
855  {
856  // find subcell ids that belong to the subface indices
857  unsigned int child_cell_1 =
859  cell_1->refinement_case(),
860  n_face_1,
861  i,
862  cell_1->face_orientation(n_face_1),
863  cell_1->face_flip(n_face_1),
864  cell_1->face_rotation(n_face_1),
865  face_1->refinement_case());
866 
867  // recursive call
868  update_periodic_face_map_recursively<dim, spacedim>(
869  cell_1->child(child_cell_1),
870  cell_2,
871  n_face_1,
872  n_face_2,
873  orientation,
874  periodic_face_map);
875  }
876  }
877  }
878  }
879  }
880 
881 
882 } // end of anonymous namespace
883 
884 
885 namespace internal
886 {
887  namespace TriangulationImplementation
888  {
889  // make sure that if in the following we
890  // write Triangulation<dim,spacedim>
891  // we mean the *class*
892  // ::Triangulation, not the
893  // enclosing namespace
894  // internal::TriangulationImplementation
895  using ::Triangulation;
896 
902  int,
903  << "Something went wrong upon construction of cell "
904  << arg1);
915  int,
916  << "Cell " << arg1
917  << " has negative measure. This typically "
918  << "indicates some distortion in the cell, or a mistakenly "
919  << "swapped pair of vertices in the input to "
920  << "Triangulation::create_triangulation().");
929  int,
930  int,
931  int,
932  << "Error while creating cell " << arg1
933  << ": the vertex index " << arg2 << " must be between 0 and "
934  << arg3 << '.');
940  int,
941  int,
942  << "While trying to assign a boundary indicator to a line: "
943  << "the line with end vertices " << arg1 << " and " << arg2
944  << " does not exist.");
950  int,
951  int,
952  int,
953  int,
954  << "While trying to assign a boundary indicator to a quad: "
955  << "the quad with bounding lines " << arg1 << ", " << arg2
956  << ", " << arg3 << ", " << arg4 << " does not exist.");
963  int,
964  int,
966  << "The input data for creating a triangulation contained "
967  << "information about a line with indices " << arg1 << " and " << arg2
968  << " that is described to have boundary indicator "
969  << static_cast<int>(arg3)
970  << ". However, this is an internal line not located on the "
971  << "boundary. You cannot assign a boundary indicator to it." << std::endl
972  << std::endl
973  << "If this happened at a place where you call "
974  << "Triangulation::create_triangulation() yourself, you need "
975  << "to check the SubCellData object you pass to this function."
976  << std::endl
977  << std::endl
978  << "If this happened in a place where you are reading a mesh "
979  << "from a file, then you need to investigate why such a line "
980  << "ended up in the input file. A typical case is a geometry "
981  << "that consisted of multiple parts and for which the mesh "
982  << "generator program assumes that the interface between "
983  << "two parts is a boundary when that isn't supposed to be "
984  << "the case, or where the mesh generator simply assigns "
985  << "'geometry indicators' to lines at the perimeter of "
986  << "a part that are not supposed to be interpreted as "
987  << "'boundary indicators'.");
994  int,
995  int,
996  int,
997  int,
999  << "The input data for creating a triangulation contained "
1000  << "information about a quad with indices " << arg1 << ", " << arg2
1001  << ", " << arg3 << ", and " << arg4
1002  << " that is described to have boundary indicator "
1003  << static_cast<int>(arg5)
1004  << ". However, this is an internal quad not located on the "
1005  << "boundary. You cannot assign a boundary indicator to it." << std::endl
1006  << std::endl
1007  << "If this happened at a place where you call "
1008  << "Triangulation::create_triangulation() yourself, you need "
1009  << "to check the SubCellData object you pass to this function."
1010  << std::endl
1011  << std::endl
1012  << "If this happened in a place where you are reading a mesh "
1013  << "from a file, then you need to investigate why such a quad "
1014  << "ended up in the input file. A typical case is a geometry "
1015  << "that consisted of multiple parts and for which the mesh "
1016  << "generator program assumes that the interface between "
1017  << "two parts is a boundary when that isn't supposed to be "
1018  << "the case, or where the mesh generator simply assigns "
1019  << "'geometry indicators' to quads at the surface of "
1020  << "a part that are not supposed to be interpreted as "
1021  << "'boundary indicators'.");
1028  int,
1029  int,
1030  << "In SubCellData the line info of the line with vertex indices " << arg1
1031  << " and " << arg2 << " appears more than once. "
1032  << "This is not allowed.");
1039  int,
1040  int,
1041  std::string,
1042  << "In SubCellData the line info of the line with vertex indices " << arg1
1043  << " and " << arg2 << " appears multiple times with different (valid) "
1044  << arg3 << ". This is not allowed.");
1051  int,
1052  int,
1053  int,
1054  int,
1055  std::string,
1056  << "In SubCellData the quad info of the quad with line indices " << arg1
1057  << ", " << arg2 << ", " << arg3 << " and " << arg4
1058  << " appears multiple times with different (valid) " << arg5
1059  << ". This is not allowed.");
1060 
1061  /*
1062  * Reserve space for TriaFaces. Details:
1063  *
1064  * Reserve space for line_orientations.
1065  *
1066  * @note Used only for dim=3.
1067  */
1068  void
1069  reserve_space(TriaFaces & tria_faces,
1070  const unsigned int new_quads_in_pairs,
1071  const unsigned int new_quads_single)
1072  {
1073  AssertDimension(tria_faces.dim, 3);
1074 
1075  Assert(new_quads_in_pairs % 2 == 0, ExcInternalError());
1076 
1077  unsigned int next_free_single = 0;
1078  unsigned int next_free_pair = 0;
1079 
1080  // count the number of objects, of unused single objects and of
1081  // unused pairs of objects
1082  unsigned int n_quads = 0;
1083  unsigned int n_unused_pairs = 0;
1084  unsigned int n_unused_singles = 0;
1085  for (unsigned int i = 0; i < tria_faces.quads.used.size(); ++i)
1086  {
1087  if (tria_faces.quads.used[i])
1088  ++n_quads;
1089  else if (i + 1 < tria_faces.quads.used.size())
1090  {
1091  if (tria_faces.quads.used[i + 1])
1092  {
1093  ++n_unused_singles;
1094  if (next_free_single == 0)
1095  next_free_single = i;
1096  }
1097  else
1098  {
1099  ++n_unused_pairs;
1100  if (next_free_pair == 0)
1101  next_free_pair = i;
1102  ++i;
1103  }
1104  }
1105  else
1106  ++n_unused_singles;
1107  }
1108  Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
1109  tria_faces.quads.used.size(),
1110  ExcInternalError());
1111  (void)n_quads;
1112 
1113  // how many single quads are needed in addition to n_unused_quads?
1114  const int additional_single_quads = new_quads_single - n_unused_singles;
1115 
1116  unsigned int new_size =
1117  tria_faces.quads.used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
1118  if (additional_single_quads > 0)
1119  new_size += additional_single_quads;
1120 
1121  // see above...
1122  if (new_size > tria_faces.quads.n_objects())
1123  {
1124  // reserve the field of the derived class
1125  tria_faces.quads_line_orientations.resize(
1126  new_size * GeometryInfo<3>::lines_per_face, true);
1127 
1128  auto &q_is_q = tria_faces.quad_is_quadrilateral;
1129  q_is_q.reserve(new_size);
1130  q_is_q.insert(q_is_q.end(), new_size - q_is_q.size(), true);
1131  }
1132  }
1133 
1134 
1135 
1149  void
1150  reserve_space(TriaLevel & tria_level,
1151  const unsigned int total_cells,
1152  const unsigned int dimension,
1153  const unsigned int space_dimension)
1154  {
1155  // we need space for total_cells cells. Maybe we have more already
1156  // with those cells which are unused, so only allocate new space if
1157  // needed.
1158  //
1159  // note that all arrays should have equal sizes (checked by
1160  // @p{monitor_memory}
1161  if (total_cells > tria_level.refine_flags.size())
1162  {
1163  tria_level.refine_flags.reserve(total_cells);
1164  tria_level.refine_flags.insert(tria_level.refine_flags.end(),
1165  total_cells -
1166  tria_level.refine_flags.size(),
1167  /*RefinementCase::no_refinement=*/0);
1168 
1169  tria_level.coarsen_flags.reserve(total_cells);
1170  tria_level.coarsen_flags.insert(tria_level.coarsen_flags.end(),
1171  total_cells -
1172  tria_level.coarsen_flags.size(),
1173  false);
1174 
1175  tria_level.active_cell_indices.reserve(total_cells);
1176  tria_level.active_cell_indices.insert(
1177  tria_level.active_cell_indices.end(),
1178  total_cells - tria_level.active_cell_indices.size(),
1180 
1181  tria_level.subdomain_ids.reserve(total_cells);
1182  tria_level.subdomain_ids.insert(tria_level.subdomain_ids.end(),
1183  total_cells -
1184  tria_level.subdomain_ids.size(),
1185  0);
1186 
1187  tria_level.level_subdomain_ids.reserve(total_cells);
1188  tria_level.level_subdomain_ids.insert(
1189  tria_level.level_subdomain_ids.end(),
1190  total_cells - tria_level.level_subdomain_ids.size(),
1191  0);
1192 
1193  tria_level.global_active_cell_indices.reserve(total_cells);
1194  tria_level.global_active_cell_indices.insert(
1195  tria_level.global_active_cell_indices.end(),
1196  total_cells - tria_level.global_active_cell_indices.size(),
1198 
1199  tria_level.global_level_cell_indices.reserve(total_cells);
1200  tria_level.global_level_cell_indices.insert(
1201  tria_level.global_level_cell_indices.end(),
1202  total_cells - tria_level.global_level_cell_indices.size(),
1204 
1205  if (dimension < space_dimension)
1206  {
1207  tria_level.direction_flags.reserve(total_cells);
1208  tria_level.direction_flags.insert(
1209  tria_level.direction_flags.end(),
1210  total_cells - tria_level.direction_flags.size(),
1211  true);
1212  }
1213  else
1214  tria_level.direction_flags.clear();
1215 
1216  tria_level.parents.reserve((total_cells + 1) / 2);
1217  tria_level.parents.insert(tria_level.parents.end(),
1218  (total_cells + 1) / 2 -
1219  tria_level.parents.size(),
1220  -1);
1221 
1222  tria_level.neighbors.reserve(total_cells * (2 * dimension));
1223  tria_level.neighbors.insert(tria_level.neighbors.end(),
1224  total_cells * (2 * dimension) -
1225  tria_level.neighbors.size(),
1226  std::make_pair(-1, -1));
1227 
1228  if (tria_level.dim == 2 || tria_level.dim == 3)
1229  {
1230  const unsigned int max_faces_per_cell = 2 * dimension;
1231  tria_level.face_orientations.resize(total_cells *
1232  max_faces_per_cell);
1233 
1234  tria_level.reference_cell.reserve(total_cells);
1235  tria_level.reference_cell.insert(
1236  tria_level.reference_cell.end(),
1237  total_cells - tria_level.reference_cell.size(),
1238  tria_level.dim == 2 ? ReferenceCells::Quadrilateral :
1240  }
1241  }
1242  }
1243 
1244 
1245 
1250  int,
1251  int,
1252  << "The containers have sizes " << arg1 << " and " << arg2
1253  << ", which is not as expected.");
1254 
1260  void
1261  monitor_memory(const TriaLevel & tria_level,
1262  const unsigned int true_dimension)
1263  {
1264  (void)tria_level;
1265  (void)true_dimension;
1266  Assert(2 * true_dimension * tria_level.refine_flags.size() ==
1267  tria_level.neighbors.size(),
1268  ExcMemoryInexact(tria_level.refine_flags.size(),
1269  tria_level.neighbors.size()));
1270  Assert(2 * true_dimension * tria_level.coarsen_flags.size() ==
1271  tria_level.neighbors.size(),
1272  ExcMemoryInexact(tria_level.coarsen_flags.size(),
1273  tria_level.neighbors.size()));
1274  }
1275 
1276 
1277 
1290  void
1291  reserve_space(TriaObjects & tria_objects,
1292  const unsigned int new_objects_in_pairs,
1293  const unsigned int new_objects_single = 0)
1294  {
1295  if (tria_objects.structdim <= 2)
1296  {
1297  Assert(new_objects_in_pairs % 2 == 0, ExcInternalError());
1298 
1299  tria_objects.next_free_single = 0;
1300  tria_objects.next_free_pair = 0;
1301  tria_objects.reverse_order_next_free_single = false;
1302 
1303  // count the number of objects, of unused single objects and of
1304  // unused pairs of objects
1305  unsigned int n_objects = 0;
1306  unsigned int n_unused_pairs = 0;
1307  unsigned int n_unused_singles = 0;
1308  for (unsigned int i = 0; i < tria_objects.used.size(); ++i)
1309  {
1310  if (tria_objects.used[i])
1311  ++n_objects;
1312  else if (i + 1 < tria_objects.used.size())
1313  {
1314  if (tria_objects.used[i + 1])
1315  {
1316  ++n_unused_singles;
1317  if (tria_objects.next_free_single == 0)
1318  tria_objects.next_free_single = i;
1319  }
1320  else
1321  {
1322  ++n_unused_pairs;
1323  if (tria_objects.next_free_pair == 0)
1324  tria_objects.next_free_pair = i;
1325  ++i;
1326  }
1327  }
1328  else
1329  ++n_unused_singles;
1330  }
1331  Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
1332  tria_objects.used.size(),
1333  ExcInternalError());
1334  (void)n_objects;
1335 
1336  // how many single objects are needed in addition to
1337  // n_unused_objects?
1338  const int additional_single_objects =
1339  new_objects_single - n_unused_singles;
1340 
1341  unsigned int new_size = tria_objects.used.size() +
1342  new_objects_in_pairs - 2 * n_unused_pairs;
1343  if (additional_single_objects > 0)
1344  new_size += additional_single_objects;
1345 
1346  // only allocate space if necessary
1347  if (new_size > tria_objects.n_objects())
1348  {
1349  const unsigned int max_faces_per_cell =
1350  2 * tria_objects.structdim;
1351  const unsigned int max_children_per_cell =
1352  1 << tria_objects.structdim;
1353 
1354  tria_objects.cells.reserve(new_size * max_faces_per_cell);
1355  tria_objects.cells.insert(tria_objects.cells.end(),
1356  (new_size - tria_objects.n_objects()) *
1357  max_faces_per_cell,
1358  -1);
1359 
1360  tria_objects.used.reserve(new_size);
1361  tria_objects.used.insert(tria_objects.used.end(),
1362  new_size - tria_objects.used.size(),
1363  false);
1364 
1365  tria_objects.user_flags.reserve(new_size);
1366  tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1367  new_size -
1368  tria_objects.user_flags.size(),
1369  false);
1370 
1371  const unsigned int factor = max_children_per_cell / 2;
1372  tria_objects.children.reserve(factor * new_size);
1373  tria_objects.children.insert(tria_objects.children.end(),
1374  factor * new_size -
1375  tria_objects.children.size(),
1376  -1);
1377 
1378  if (tria_objects.structdim > 1)
1379  {
1380  tria_objects.refinement_cases.reserve(new_size);
1381  tria_objects.refinement_cases.insert(
1382  tria_objects.refinement_cases.end(),
1383  new_size - tria_objects.refinement_cases.size(),
1384  /*RefinementCase::no_refinement=*/0);
1385  }
1386 
1387  // first reserve, then resize. Otherwise the std library can
1388  // decide to allocate more entries.
1389  tria_objects.boundary_or_material_id.reserve(new_size);
1390  tria_objects.boundary_or_material_id.resize(new_size);
1391 
1392  tria_objects.user_data.reserve(new_size);
1393  tria_objects.user_data.resize(new_size);
1394 
1395  tria_objects.manifold_id.reserve(new_size);
1396  tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1397  new_size -
1398  tria_objects.manifold_id.size(),
1400  }
1401 
1402  if (n_unused_singles == 0)
1403  {
1404  tria_objects.next_free_single = new_size - 1;
1405  tria_objects.reverse_order_next_free_single = true;
1406  }
1407  }
1408  else
1409  {
1410  const unsigned int new_hexes = new_objects_in_pairs;
1411 
1412  const unsigned int new_size =
1413  new_hexes + std::count(tria_objects.used.begin(),
1414  tria_objects.used.end(),
1415  true);
1416 
1417  // see above...
1418  if (new_size > tria_objects.n_objects())
1419  {
1420  const unsigned int max_faces_per_cell =
1421  2 * tria_objects.structdim;
1422 
1423  tria_objects.cells.reserve(new_size * max_faces_per_cell);
1424  tria_objects.cells.insert(tria_objects.cells.end(),
1425  (new_size - tria_objects.n_objects()) *
1426  max_faces_per_cell,
1427  -1);
1428 
1429  tria_objects.used.reserve(new_size);
1430  tria_objects.used.insert(tria_objects.used.end(),
1431  new_size - tria_objects.used.size(),
1432  false);
1433 
1434  tria_objects.user_flags.reserve(new_size);
1435  tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1436  new_size -
1437  tria_objects.user_flags.size(),
1438  false);
1439 
1440  tria_objects.children.reserve(4 * new_size);
1441  tria_objects.children.insert(tria_objects.children.end(),
1442  4 * new_size -
1443  tria_objects.children.size(),
1444  -1);
1445 
1446  // for the following fields, we know exactly how many elements
1447  // we need, so first reserve then resize (resize itself, at least
1448  // with some compiler libraries, appears to round up the size it
1449  // actually reserves)
1450  tria_objects.boundary_or_material_id.reserve(new_size);
1451  tria_objects.boundary_or_material_id.resize(new_size);
1452 
1453  tria_objects.manifold_id.reserve(new_size);
1454  tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1455  new_size -
1456  tria_objects.manifold_id.size(),
1458 
1459  tria_objects.user_data.reserve(new_size);
1460  tria_objects.user_data.resize(new_size);
1461 
1462  tria_objects.refinement_cases.reserve(new_size);
1463  tria_objects.refinement_cases.insert(
1464  tria_objects.refinement_cases.end(),
1465  new_size - tria_objects.refinement_cases.size(),
1466  /*RefinementCase::no_refinement=*/0);
1467  }
1468  tria_objects.next_free_single = tria_objects.next_free_pair = 0;
1469  }
1470  }
1471 
1472 
1473 
1479  void
1480  monitor_memory(const TriaObjects &tria_object, const unsigned int)
1481  {
1482  Assert(tria_object.n_objects() == tria_object.used.size(),
1483  ExcMemoryInexact(tria_object.n_objects(),
1484  tria_object.used.size()));
1485  Assert(tria_object.n_objects() == tria_object.user_flags.size(),
1486  ExcMemoryInexact(tria_object.n_objects(),
1487  tria_object.user_flags.size()));
1488  Assert(tria_object.n_objects() ==
1489  tria_object.boundary_or_material_id.size(),
1490  ExcMemoryInexact(tria_object.n_objects(),
1491  tria_object.boundary_or_material_id.size()));
1492  Assert(tria_object.n_objects() == tria_object.manifold_id.size(),
1493  ExcMemoryInexact(tria_object.n_objects(),
1494  tria_object.manifold_id.size()));
1495  Assert(tria_object.n_objects() == tria_object.user_data.size(),
1496  ExcMemoryInexact(tria_object.n_objects(),
1497  tria_object.user_data.size()));
1498 
1499  if (tria_object.structdim == 1)
1500  {
1501  Assert(1 * tria_object.n_objects() == tria_object.children.size(),
1502  ExcMemoryInexact(tria_object.n_objects(),
1503  tria_object.children.size()));
1504  }
1505  else if (tria_object.structdim == 2)
1506  {
1507  Assert(2 * tria_object.n_objects() == tria_object.children.size(),
1508  ExcMemoryInexact(tria_object.n_objects(),
1509  tria_object.children.size()));
1510  }
1511  else if (tria_object.structdim == 3)
1512  {
1513  Assert(4 * tria_object.n_objects() == tria_object.children.size(),
1514  ExcMemoryInexact(tria_object.n_objects(),
1515  tria_object.children.size()));
1516  }
1517  }
1518 
1519 
1520 
1525  template <int dim, int spacedim>
1526  class Policy
1527  {
1528  public:
1532  virtual ~Policy() = default;
1533 
1537  virtual void
1539 
1543  virtual void
1547  std::vector<unsigned int> & line_cell_count,
1548  std::vector<unsigned int> &quad_cell_count) = 0;
1549 
1555  const bool check_for_distorted_cells) = 0;
1556 
1560  virtual void
1563 
1567  virtual void
1570 
1574  virtual bool
1576  const typename Triangulation<dim, spacedim>::cell_iterator &cell) = 0;
1577 
1584  virtual std::unique_ptr<Policy<dim, spacedim>>
1585  clone() = 0;
1586  };
1587 
1588 
1589 
1595  template <int dim, int spacedim, typename T>
1596  class PolicyWrapper : public Policy<dim, spacedim>
1597  {
1598  public:
1599  void
1601  {
1602  T::update_neighbors(tria);
1603  }
1604 
1605  void
1609  std::vector<unsigned int> & line_cell_count,
1610  std::vector<unsigned int> &quad_cell_count) override
1611  {
1612  T::delete_children(tria, cell, line_cell_count, quad_cell_count);
1613  }
1614 
1617  const bool check_for_distorted_cells) override
1618  {
1619  return T::execute_refinement(triangulation, check_for_distorted_cells);
1620  }
1621 
1622  void
1625  {
1626  T::prevent_distorted_boundary_cells(triangulation);
1627  }
1628 
1629  void
1632  {
1633  T::prepare_refinement_dim_dependent(triangulation);
1634  }
1635 
1636  bool
1638  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
1639  override
1640  {
1641  return T::template coarsening_allowed<dim, spacedim>(cell);
1642  }
1643 
1644  std::unique_ptr<Policy<dim, spacedim>>
1645  clone() override
1646  {
1647  return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
1648  }
1649  };
1650 
1651 
1652 
1749  {
1761  template <int dim, int spacedim>
1762  static void
1765  const unsigned int level_objects,
1767  {
1768  using line_iterator =
1770 
1771  number_cache.n_levels = 0;
1772  if (level_objects > 0)
1773  // find the last level on which there are used cells
1774  for (unsigned int level = 0; level < level_objects; ++level)
1775  if (triangulation.begin(level) != triangulation.end(level))
1776  number_cache.n_levels = level + 1;
1777 
1778  // no cells at all?
1779  Assert(number_cache.n_levels > 0, ExcInternalError());
1780 
1781  //---------------------------------
1782  // update the number of lines on the different levels in the
1783  // cache
1784  number_cache.n_lines = 0;
1785  number_cache.n_active_lines = 0;
1786 
1787  // for 1d, lines have levels so take count the objects per
1788  // level and globally
1789  if (dim == 1)
1790  {
1791  number_cache.n_lines_level.resize(number_cache.n_levels);
1792  number_cache.n_active_lines_level.resize(number_cache.n_levels);
1793 
1794  for (unsigned int level = 0; level < number_cache.n_levels; ++level)
1795  {
1796  // count lines on this level
1797  number_cache.n_lines_level[level] = 0;
1798  number_cache.n_active_lines_level[level] = 0;
1799 
1800  line_iterator line = triangulation.begin_line(level),
1801  endc =
1802  (level == number_cache.n_levels - 1 ?
1803  line_iterator(triangulation.end_line()) :
1804  triangulation.begin_line(level + 1));
1805  for (; line != endc; ++line)
1806  {
1807  ++number_cache.n_lines_level[level];
1808  if (line->has_children() == false)
1809  ++number_cache.n_active_lines_level[level];
1810  }
1811 
1812  // update total number of lines
1813  number_cache.n_lines += number_cache.n_lines_level[level];
1814  number_cache.n_active_lines +=
1815  number_cache.n_active_lines_level[level];
1816  }
1817  }
1818  else
1819  {
1820  // for dim>1, there are no levels for lines
1821  number_cache.n_lines_level.clear();
1822  number_cache.n_active_lines_level.clear();
1823 
1824  line_iterator line = triangulation.begin_line(),
1825  endc = triangulation.end_line();
1826  for (; line != endc; ++line)
1827  {
1828  ++number_cache.n_lines;
1829  if (line->has_children() == false)
1830  ++number_cache.n_active_lines;
1831  }
1832  }
1833  }
1834 
1849  template <int dim, int spacedim>
1850  static void
1853  const unsigned int level_objects,
1855  {
1856  // update lines and n_levels in number_cache. since we don't
1857  // access any of these numbers, we can do this in the
1858  // background
1859  Threads::Task<void> update_lines = Threads::new_task(
1860  static_cast<
1861  void (*)(const Triangulation<dim, spacedim> &,
1862  const unsigned int,
1864  &compute_number_cache_dim<dim, spacedim>),
1865  triangulation,
1866  level_objects,
1868  number_cache));
1869 
1870  using quad_iterator =
1872 
1873  //---------------------------------
1874  // update the number of quads on the different levels in the
1875  // cache
1876  number_cache.n_quads = 0;
1877  number_cache.n_active_quads = 0;
1878 
1879  // for 2d, quads have levels so take count the objects per
1880  // level and globally
1881  if (dim == 2)
1882  {
1883  // count the number of levels; the function we called above
1884  // on a separate Task for lines also does this and puts it into
1885  // number_cache.n_levels, but this datum may not yet be
1886  // available as we call the function on a separate task
1887  unsigned int n_levels = 0;
1888  if (level_objects > 0)
1889  // find the last level on which there are used cells
1890  for (unsigned int level = 0; level < level_objects; ++level)
1891  if (triangulation.begin(level) != triangulation.end(level))
1892  n_levels = level + 1;
1893 
1894  number_cache.n_quads_level.resize(n_levels);
1895  number_cache.n_active_quads_level.resize(n_levels);
1896 
1897  for (unsigned int level = 0; level < n_levels; ++level)
1898  {
1899  // count quads on this level
1900  number_cache.n_quads_level[level] = 0;
1901  number_cache.n_active_quads_level[level] = 0;
1902 
1903  quad_iterator quad = triangulation.begin_quad(level),
1904  endc =
1905  (level == n_levels - 1 ?
1906  quad_iterator(triangulation.end_quad()) :
1907  triangulation.begin_quad(level + 1));
1908  for (; quad != endc; ++quad)
1909  {
1910  ++number_cache.n_quads_level[level];
1911  if (quad->has_children() == false)
1912  ++number_cache.n_active_quads_level[level];
1913  }
1914 
1915  // update total number of quads
1916  number_cache.n_quads += number_cache.n_quads_level[level];
1917  number_cache.n_active_quads +=
1918  number_cache.n_active_quads_level[level];
1919  }
1920  }
1921  else
1922  {
1923  // for dim>2, there are no levels for quads
1924  number_cache.n_quads_level.clear();
1925  number_cache.n_active_quads_level.clear();
1926 
1927  quad_iterator quad = triangulation.begin_quad(),
1928  endc = triangulation.end_quad();
1929  for (; quad != endc; ++quad)
1930  {
1931  ++number_cache.n_quads;
1932  if (quad->has_children() == false)
1933  ++number_cache.n_active_quads;
1934  }
1935  }
1936 
1937  // wait for the background computation for lines
1938  update_lines.join();
1939  }
1940 
1956  template <int dim, int spacedim>
1957  static void
1960  const unsigned int level_objects,
1962  {
1963  // update quads, lines and n_levels in number_cache. since we
1964  // don't access any of these numbers, we can do this in the
1965  // background
1966  Threads::Task<void> update_quads_and_lines = Threads::new_task(
1967  static_cast<
1968  void (*)(const Triangulation<dim, spacedim> &,
1969  const unsigned int,
1971  &compute_number_cache_dim<dim, spacedim>),
1972  triangulation,
1973  level_objects,
1975  number_cache));
1976 
1977  using hex_iterator =
1979 
1980  //---------------------------------
1981  // update the number of hexes on the different levels in the
1982  // cache
1983  number_cache.n_hexes = 0;
1984  number_cache.n_active_hexes = 0;
1985 
1986  // for 3d, hexes have levels so take count the objects per
1987  // level and globally
1988  if (dim == 3)
1989  {
1990  // count the number of levels; the function we called
1991  // above on a separate Task for quads (recursively, via
1992  // the lines function) also does this and puts it into
1993  // number_cache.n_levels, but this datum may not yet be
1994  // available as we call the function on a separate task
1995  unsigned int n_levels = 0;
1996  if (level_objects > 0)
1997  // find the last level on which there are used cells
1998  for (unsigned int level = 0; level < level_objects; ++level)
1999  if (triangulation.begin(level) != triangulation.end(level))
2000  n_levels = level + 1;
2001 
2002  number_cache.n_hexes_level.resize(n_levels);
2003  number_cache.n_active_hexes_level.resize(n_levels);
2004 
2005  for (unsigned int level = 0; level < n_levels; ++level)
2006  {
2007  // count hexes on this level
2008  number_cache.n_hexes_level[level] = 0;
2009  number_cache.n_active_hexes_level[level] = 0;
2010 
2011  hex_iterator hex = triangulation.begin_hex(level),
2012  endc = (level == n_levels - 1 ?
2013  hex_iterator(triangulation.end_hex()) :
2014  triangulation.begin_hex(level + 1));
2015  for (; hex != endc; ++hex)
2016  {
2017  ++number_cache.n_hexes_level[level];
2018  if (hex->has_children() == false)
2019  ++number_cache.n_active_hexes_level[level];
2020  }
2021 
2022  // update total number of hexes
2023  number_cache.n_hexes += number_cache.n_hexes_level[level];
2024  number_cache.n_active_hexes +=
2025  number_cache.n_active_hexes_level[level];
2026  }
2027  }
2028  else
2029  {
2030  // for dim>3, there are no levels for hexes
2031  number_cache.n_hexes_level.clear();
2032  number_cache.n_active_hexes_level.clear();
2033 
2034  hex_iterator hex = triangulation.begin_hex(),
2035  endc = triangulation.end_hex();
2036  for (; hex != endc; ++hex)
2037  {
2038  ++number_cache.n_hexes;
2039  if (hex->has_children() == false)
2040  ++number_cache.n_active_hexes;
2041  }
2042  }
2043 
2044  // wait for the background computation for quads
2045  update_quads_and_lines.join();
2046  }
2047 
2048 
2049  template <int dim, int spacedim>
2050  static void
2053  const unsigned int level_objects,
2055  {
2056  compute_number_cache_dim(triangulation, level_objects, number_cache);
2057 
2058  number_cache.active_cell_index_partitioner =
2059  std::make_shared<const Utilities::MPI::Partitioner>(
2060  triangulation.n_active_cells());
2061 
2062  number_cache.level_cell_index_partitioners.resize(
2063  triangulation.n_levels());
2064  for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
2065  number_cache.level_cell_index_partitioners[level] =
2066  std::make_shared<const Utilities::MPI::Partitioner>(
2067  triangulation.n_cells(level));
2068  }
2069 
2070 
2071  template <int spacedim>
2072  static void
2074  {}
2075 
2076 
2077  template <int dim, int spacedim>
2078  static void
2080  {
2081  // each face can be neighbored on two sides
2082  // by cells. according to the face's
2083  // intrinsic normal we define the left
2084  // neighbor as the one for which the face
2085  // normal points outward, and store that
2086  // one first; the second one is then
2087  // the right neighbor for which the
2088  // face normal points inward. This
2089  // information depends on the type of cell
2090  // and local number of face for the
2091  // 'standard ordering and orientation' of
2092  // faces and then on the face_orientation
2093  // information for the real mesh. Set up a
2094  // table to have fast access to those
2095  // offsets (0 for left and 1 for
2096  // right). Some of the values are invalid
2097  // as they reference too large face
2098  // numbers, but we just leave them at a
2099  // zero value.
2100  //
2101  // Note, that in 2d for lines as faces the
2102  // normal direction given in the
2103  // GeometryInfo class is not consistent. We
2104  // thus define here that the normal for a
2105  // line points to the right if the line
2106  // points upwards.
2107  //
2108  // There is one more point to
2109  // consider, however: if we have
2110  // dim<spacedim, then we may have
2111  // cases where cells are
2112  // inverted. In effect, both
2113  // cells think they are the left
2114  // neighbor of an edge, for
2115  // example, which leads us to
2116  // forget neighborship
2117  // information (a case that shows
2118  // this is
2119  // codim_one/hanging_nodes_02). We
2120  // store whether a cell is
2121  // inverted using the
2122  // direction_flag, so if a cell
2123  // has a false direction_flag,
2124  // then we need to invert our
2125  // selection whether we are a
2126  // left or right neighbor in all
2127  // following computations.
2128  //
2129  // first index: dimension (minus 2)
2130  // second index: local face index
2131  // third index: face_orientation (false and true)
2132  static const unsigned int left_right_offset[2][6][2] = {
2133  // quadrilateral
2134  {{0, 1}, // face 0, face_orientation = false and true
2135  {1, 0}, // face 1, face_orientation = false and true
2136  {1, 0}, // face 2, face_orientation = false and true
2137  {0, 1}, // face 3, face_orientation = false and true
2138  {0, 0}, // face 4, invalid face
2139  {0, 0}}, // face 5, invalid face
2140  // hexahedron
2141  {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
2142 
2143  // now create a vector of the two active
2144  // neighbors (left and right) for each face
2145  // and fill it by looping over all cells. For
2146  // cases with anisotropic refinement and more
2147  // then one cell neighboring at a given side
2148  // of the face we will automatically get the
2149  // active one on the highest level as we loop
2150  // over cells from lower levels first.
2151  const typename Triangulation<dim, spacedim>::cell_iterator dummy;
2152  std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
2153  adjacent_cells(2 * triangulation.n_raw_faces(), dummy);
2154 
2155  for (const auto &cell : triangulation.cell_iterators())
2156  for (auto f : cell->face_indices())
2157  {
2158  const typename Triangulation<dim, spacedim>::face_iterator face =
2159  cell->face(f);
2160 
2161  const unsigned int offset =
2162  (cell->direction_flag() ?
2163  left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2164  1 -
2165  left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2166 
2167  adjacent_cells[2 * face->index() + offset] = cell;
2168 
2169  // if this cell is not refined, but the
2170  // face is, then we'll have to set our
2171  // cell as neighbor for the child faces
2172  // as well. Fortunately the normal
2173  // orientation of children will be just
2174  // the same.
2175  if (dim == 2)
2176  {
2177  if (cell->is_active() && face->has_children())
2178  {
2179  adjacent_cells[2 * face->child(0)->index() + offset] =
2180  cell;
2181  adjacent_cells[2 * face->child(1)->index() + offset] =
2182  cell;
2183  }
2184  }
2185  else // -> dim == 3
2186  {
2187  // We need the same as in 2d
2188  // here. Furthermore, if the face is
2189  // refined with cut_x or cut_y then
2190  // those children again in the other
2191  // direction, and if this cell is
2192  // refined isotropically (along the
2193  // face) then the neighbor will
2194  // (probably) be refined as cut_x or
2195  // cut_y along the face. For those
2196  // neighboring children cells, their
2197  // neighbor will be the current,
2198  // inactive cell, as our children are
2199  // too fine to be neighbors. Catch that
2200  // case by also acting on inactive
2201  // cells with isotropic refinement
2202  // along the face. If the situation
2203  // described is not present, the data
2204  // will be overwritten later on when we
2205  // visit cells on finer levels, so no
2206  // harm will be done.
2207  if (face->has_children() &&
2208  (cell->is_active() ||
2210  cell->refinement_case(), f) ==
2211  RefinementCase<dim - 1>::isotropic_refinement))
2212  {
2213  for (unsigned int c = 0; c < face->n_children(); ++c)
2214  adjacent_cells[2 * face->child(c)->index() + offset] =
2215  cell;
2216  if (face->child(0)->has_children())
2217  {
2218  adjacent_cells[2 * face->child(0)->child(0)->index() +
2219  offset] = cell;
2220  adjacent_cells[2 * face->child(0)->child(1)->index() +
2221  offset] = cell;
2222  }
2223  if (face->child(1)->has_children())
2224  {
2225  adjacent_cells[2 * face->child(1)->child(0)->index() +
2226  offset] = cell;
2227  adjacent_cells[2 * face->child(1)->child(1)->index() +
2228  offset] = cell;
2229  }
2230  } // if cell active and face refined
2231  } // else -> dim==3
2232  } // for all faces of all cells
2233 
2234  // now loop again over all cells and set the
2235  // corresponding neighbor cell. Note, that we
2236  // have to use the opposite of the
2237  // left_right_offset in this case as we want
2238  // the offset of the neighbor, not our own.
2239  for (const auto &cell : triangulation.cell_iterators())
2240  for (auto f : cell->face_indices())
2241  {
2242  const unsigned int offset =
2243  (cell->direction_flag() ?
2244  left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2245  1 -
2246  left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2247  cell->set_neighbor(
2248  f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
2249  }
2250  }
2251 
2252 
2256  template <int dim, int spacedim>
2257  static void
2259  const std::vector<CellData<dim>> & cells,
2260  const SubCellData & subcelldata,
2262  {
2263  AssertThrow(vertices.size() > 0, ExcMessage("No vertices given"));
2264  AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
2265 
2266  // Check that all cells have positive volume.
2267 #ifndef _MSC_VER
2268  // TODO: The following code does not compile with MSVC. Find a way
2269  // around it
2270  if (dim == spacedim)
2271  for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2272  {
2273  // If we should check for distorted cells, then we permit them
2274  // to exist. If a cell has negative measure, then it must be
2275  // distorted (the converse is not necessarily true); hence
2276  // throw an exception if no such cells should exist.
2278  {
2279  const double cell_measure = GridTools::cell_measure<spacedim>(
2280  vertices,
2281  ArrayView<const unsigned int>(cells[cell_no].vertices));
2283  }
2284  }
2285 #endif
2286 
2287  // clear old content
2288  tria.levels.clear();
2289  tria.levels.push_back(
2290  std::make_unique<
2292 
2293  if (dim > 1)
2294  tria.faces = std::make_unique<
2296 
2297  // copy vertices
2299  tria.vertices_used.assign(vertices.size(), true);
2300 
2301  // compute connectivity
2302  const auto connectivity = build_connectivity<unsigned int>(cells);
2303  const unsigned int n_cell = cells.size();
2304 
2305  // TriaObjects: lines
2306  if (dim >= 2)
2307  {
2308  auto &lines_0 = tria.faces->lines; // data structure to be filled
2309 
2310  // get connectivity between quads and lines
2311  const auto & crs = connectivity.entity_to_entities(1, 0);
2312  const unsigned int n_lines = crs.ptr.size() - 1;
2313 
2314  // allocate memory
2315  reserve_space_(lines_0, n_lines);
2316 
2317  // loop over lines
2318  for (unsigned int line = 0; line < n_lines; ++line)
2319  for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
2320  ++i, ++j)
2321  lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
2322  crs.col[i]; // set vertex indices
2323  }
2324 
2325  // TriaObjects: quads
2326  if (dim == 3)
2327  {
2328  auto &quads_0 = tria.faces->quads; // data structures to be filled
2329  auto &faces = *tria.faces;
2330 
2331  // get connectivity between quads and lines
2332  const auto & crs = connectivity.entity_to_entities(2, 1);
2333  const unsigned int n_quads = crs.ptr.size() - 1;
2334 
2335  // allocate memory
2336  reserve_space_(quads_0, n_quads);
2337  reserve_space_(faces, 2 /*structdim*/, n_quads);
2338 
2339  // loop over all quads -> entity type, line indices/orientations
2340  for (unsigned int q = 0, k = 0; q < n_quads; ++q)
2341  {
2342  // set entity type of quads
2343  faces.set_quad_type(q, connectivity.entity_types(2)[q]);
2344 
2345  // loop over all its lines
2346  for (unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
2347  ++i, ++j, ++k)
2348  {
2349  // set line index
2350  quads_0.cells[q * GeometryInfo<3>::lines_per_face + j] =
2351  crs.col[i];
2352 
2353  // set line orientations
2354  const unsigned char combined_orientation =
2355  connectivity.entity_orientations(1)
2356  .get_combined_orientation(k);
2357  // it doesn't make sense to set any flags except
2358  // orientation for a line
2359  Assert(
2360  combined_orientation ==
2362  combined_orientation ==
2364  ExcInternalError());
2365  faces.quads_line_orientations
2367  combined_orientation ==
2369  }
2370  }
2371  }
2372 
2373  // TriaObjects/TriaLevel: cell
2374  {
2375  auto &cells_0 = tria.levels[0]->cells; // data structure to be filled
2376  auto &level = *tria.levels[0];
2377 
2378  // get connectivity between cells/faces and cells/cells
2379  const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2380  const auto &nei = connectivity.entity_to_entities(dim, dim);
2381 
2382  // in 2d optional: since in in pure QUAD meshes same line
2383  // orientations can be guaranteed
2384  bool orientation_needed = false;
2385  if (dim == 3)
2386  orientation_needed = true;
2387  else if (dim == 2)
2388  {
2389  const auto &orientations = connectivity.entity_orientations(1);
2390  for (unsigned int i = 0; i < orientations.n_objects(); ++i)
2391  if (orientations.get_combined_orientation(i) !=
2393  {
2394  orientation_needed = true;
2395  break;
2396  }
2397  }
2398 
2399  // allocate memory
2400  reserve_space_(cells_0, n_cell);
2401  reserve_space_(level, spacedim, n_cell, orientation_needed);
2402 
2403  // loop over all cells
2404  for (unsigned int cell = 0; cell < n_cell; ++cell)
2405  {
2406  // set material ids
2407  cells_0.boundary_or_material_id[cell].material_id =
2408  cells[cell].material_id;
2409 
2410  // set manifold ids
2411  cells_0.manifold_id[cell] = cells[cell].manifold_id;
2412 
2413  // set entity types
2414  level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
2415 
2416  // loop over faces
2417  for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2418  ++i, ++j)
2419  {
2420  // set neighbor if not at boundary
2421  if (nei.col[i] != static_cast<unsigned int>(-1))
2422  level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
2423  j] = {0, nei.col[i]};
2424 
2425  // set face indices
2426  cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
2427  crs.col[i];
2428 
2429  // set face orientation if needed
2430  if (orientation_needed)
2431  {
2432  level.face_orientations.set_combined_orientation(
2434  connectivity.entity_orientations(dim - 1)
2435  .get_combined_orientation(i));
2436  }
2437  }
2438  }
2439  }
2440 
2441  // TriaFaces: boundary id of boundary faces
2442  if (dim > 1)
2443  {
2444  auto &bids_face = dim == 3 ?
2445  tria.faces->quads.boundary_or_material_id :
2446  tria.faces->lines.boundary_or_material_id;
2447 
2448  // count number of cells a face is belonging to
2449  std::vector<unsigned int> count(bids_face.size(), 0);
2450 
2451  // get connectivity between cells/faces
2452  const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2453 
2454  // count how many cells are adjacent to the same face
2455  for (unsigned int cell = 0; cell < cells.size(); ++cell)
2456  for (unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
2457  count[crs.col[i]]++;
2458 
2459  // loop over all faces
2460  for (unsigned int face = 0; face < count.size(); ++face)
2461  {
2462  if (count[face] != 1) // inner face
2463  continue;
2464 
2465  // boundary faces ...
2466  bids_face[face].boundary_id = 0;
2467 
2468  if (dim != 3)
2469  continue;
2470 
2471  // ... and the lines of quads in 3d
2472  const auto &crs = connectivity.entity_to_entities(2, 1);
2473  for (unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
2474  tria.faces->lines.boundary_or_material_id[crs.col[i]]
2475  .boundary_id = 0;
2476  }
2477  }
2478  else // 1d
2479  {
2480  static const unsigned int t_tba = static_cast<unsigned int>(-1);
2481  static const unsigned int t_inner = static_cast<unsigned int>(-2);
2482 
2483  std::vector<unsigned int> type(vertices.size(), t_tba);
2484 
2485  const auto &crs = connectivity.entity_to_entities(1, 0);
2486 
2487  for (unsigned int cell = 0; cell < cells.size(); ++cell)
2488  for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2489  ++i, ++j)
2490  if (type[crs.col[i]] != t_inner)
2491  type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
2492 
2493  for (unsigned int face = 0; face < type.size(); ++face)
2494  {
2495  // note: we also treat manifolds here!?
2498  if (type[face] != t_inner && type[face] != t_tba)
2499  (*tria.vertex_to_boundary_id_map_1d)[face] = type[face];
2500  }
2501  }
2502 
2503  // SubCellData: line
2504  if (dim >= 2)
2505  process_subcelldata(connectivity.entity_to_entities(1, 0),
2506  tria.faces->lines,
2507  subcelldata.boundary_lines,
2508  vertices);
2509 
2510  // SubCellData: quad
2511  if (dim == 3)
2512  process_subcelldata(connectivity.entity_to_entities(2, 0),
2513  tria.faces->quads,
2514  subcelldata.boundary_quads,
2515  vertices);
2516  }
2517 
2518 
2519  template <int structdim, int spacedim, typename T>
2520  static void
2522  const CRS<T> & crs,
2523  TriaObjects & obj,
2524  const std::vector<CellData<structdim>> &boundary_objects_in,
2525  const std::vector<Point<spacedim>> & vertex_locations)
2526  {
2527  AssertDimension(obj.structdim, structdim);
2528 
2529  if (boundary_objects_in.size() == 0)
2530  return; // empty subcelldata -> nothing to do
2531 
2532  // pre-sort subcelldata
2533  auto boundary_objects = boundary_objects_in;
2534 
2535  // ... sort vertices
2536  for (auto &boundary_object : boundary_objects)
2537  std::sort(boundary_object.vertices.begin(),
2538  boundary_object.vertices.end());
2539 
2540  // ... sort cells
2541  std::sort(boundary_objects.begin(),
2542  boundary_objects.end(),
2543  [](const auto &a, const auto &b) {
2544  return a.vertices < b.vertices;
2545  });
2546 
2547  unsigned int counter = 0;
2548 
2549  std::vector<unsigned int> key;
2551 
2552  for (unsigned int o = 0; o < obj.n_objects(); ++o)
2553  {
2554  auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
2555  auto &manifold_id = obj.manifold_id[o];
2556 
2557  // assert that object has not been visited yet and its value
2558  // has not been modified yet
2559  AssertThrow(boundary_id == 0 ||
2561  ExcNotImplemented());
2563  ExcNotImplemented());
2564 
2565  // create key
2566  key.assign(crs.col.data() + crs.ptr[o],
2567  crs.col.data() + crs.ptr[o + 1]);
2568  std::sort(key.begin(), key.end());
2569 
2570  // is subcelldata provided? -> binary search
2571  const auto subcell_object =
2572  std::lower_bound(boundary_objects.begin(),
2573  boundary_objects.end(),
2574  key,
2575  [&](const auto &cell, const auto &key) {
2576  return cell.vertices < key;
2577  });
2578 
2579  // no subcelldata provided for this object
2580  if (subcell_object == boundary_objects.end() ||
2581  subcell_object->vertices != key)
2582  continue;
2583 
2584  counter++;
2585 
2586  // set manifold id
2587  manifold_id = subcell_object->manifold_id;
2588 
2589  // set boundary id
2590  if (subcell_object->boundary_id !=
2592  {
2593  (void)vertex_locations;
2594  AssertThrow(
2596  ExcMessage(
2597  "The input arguments for creating a triangulation "
2598  "specified a boundary id for an internal face. This "
2599  "is not allowed."
2600  "\n\n"
2601  "The object in question has vertex indices " +
2602  [subcell_object]() {
2603  std::string s;
2604  for (const auto v : subcell_object->vertices)
2605  s += std::to_string(v) + ',';
2606  return s;
2607  }() +
2608  " which are located at positions " +
2609  [vertex_locations, subcell_object]() {
2610  std::ostringstream s;
2611  for (const auto v : subcell_object->vertices)
2612  s << '(' << vertex_locations[v] << ')';
2613  return s.str();
2614  }() +
2615  "."));
2616  boundary_id = subcell_object->boundary_id;
2617  }
2618  }
2619 
2620  // make sure that all subcelldata entries have been processed
2621  // TODO: this is not guaranteed, why?
2622  // AssertDimension(counter, boundary_objects_in.size());
2623  (void)counter;
2624  }
2625 
2626 
2627 
2628  static void
2630  const unsigned structdim,
2631  const unsigned int size)
2632  {
2633  const unsigned int dim = faces.dim;
2634 
2635  const unsigned int max_lines_per_face = 2 * structdim;
2636 
2637  if (dim == 3 && structdim == 2)
2638  {
2639  // quad entity types
2640  faces.quad_is_quadrilateral.assign(size, true);
2641 
2642  // quad line orientations
2643  faces.quads_line_orientations.assign(size * max_lines_per_face,
2644  true);
2645  }
2646  }
2647 
2648 
2649 
2650  static void
2652  const unsigned int spacedim,
2653  const unsigned int size,
2654  const bool orientation_needed)
2655  {
2656  const unsigned int dim = level.dim;
2657 
2658  const unsigned int max_faces_per_cell = 2 * dim;
2659 
2660  level.active_cell_indices.assign(size, -1);
2661  level.subdomain_ids.assign(size, 0);
2662  level.level_subdomain_ids.assign(size, 0);
2663 
2664  level.refine_flags.assign(size, 0u);
2665  level.coarsen_flags.assign(size, false);
2666 
2667  level.parents.assign((size + 1) / 2, -1);
2668 
2669  if (dim < spacedim)
2670  level.direction_flags.assign(size, true);
2671 
2672  level.neighbors.assign(size * max_faces_per_cell, {-1, -1});
2673 
2674  level.reference_cell.assign(size, ReferenceCells::Invalid);
2675 
2676  if (orientation_needed)
2677  level.face_orientations.reinit(size * max_faces_per_cell);
2678 
2679 
2680  level.global_active_cell_indices.assign(size,
2682  level.global_level_cell_indices.assign(size,
2684  }
2685 
2686 
2687 
2688  static void
2689  reserve_space_(TriaObjects &obj, const unsigned int size)
2690  {
2691  const unsigned int structdim = obj.structdim;
2692 
2693  const unsigned int max_children_per_cell = 1 << structdim;
2694  const unsigned int max_faces_per_cell = 2 * structdim;
2695 
2696  obj.used.assign(size, true);
2697  obj.boundary_or_material_id.assign(
2698  size,
2700  BoundaryOrMaterialId());
2701  obj.manifold_id.assign(size, -1);
2702  obj.user_flags.assign(size, false);
2703  obj.user_data.resize(size);
2704 
2705  if (structdim > 1) // TODO: why?
2706  obj.refinement_cases.assign(size, 0);
2707 
2708  obj.children.assign(max_children_per_cell / 2 * size, -1);
2709 
2710  obj.cells.assign(max_faces_per_cell * size, -1);
2711 
2712  if (structdim <= 2)
2713  {
2714  obj.next_free_single = size - 1;
2715  obj.next_free_pair = 0;
2716  obj.reverse_order_next_free_single = true;
2717  }
2718  else
2719  {
2720  obj.next_free_single = obj.next_free_pair = 0;
2721  }
2722  }
2723 
2724 
2740  template <int spacedim>
2741  static void
2744  std::vector<unsigned int> &,
2745  std::vector<unsigned int> &)
2746  {
2747  const unsigned int dim = 1;
2748 
2749  // first we need to reset the
2750  // neighbor pointers of the
2751  // neighbors of this cell's
2752  // children to this cell. This is
2753  // different for one dimension,
2754  // since there neighbors can have a
2755  // refinement level differing from
2756  // that of this cell's children by
2757  // more than one level.
2758 
2759  Assert(!cell->child(0)->has_children() &&
2760  !cell->child(1)->has_children(),
2761  ExcInternalError());
2762 
2763  // first do it for the cells to the
2764  // left
2765  if (cell->neighbor(0).state() == IteratorState::valid)
2766  if (cell->neighbor(0)->has_children())
2767  {
2769  cell->neighbor(0);
2770  Assert(neighbor->level() == cell->level(), ExcInternalError());
2771 
2772  // right child
2773  neighbor = neighbor->child(1);
2774  while (true)
2775  {
2776  Assert(neighbor->neighbor(1) == cell->child(0),
2777  ExcInternalError());
2778  neighbor->set_neighbor(1, cell);
2779 
2780  // move on to further
2781  // children on the
2782  // boundary between this
2783  // cell and its neighbor
2784  if (neighbor->has_children())
2785  neighbor = neighbor->child(1);
2786  else
2787  break;
2788  }
2789  }
2790 
2791  // now do it for the cells to the
2792  // left
2793  if (cell->neighbor(1).state() == IteratorState::valid)
2794  if (cell->neighbor(1)->has_children())
2795  {
2797  cell->neighbor(1);
2798  Assert(neighbor->level() == cell->level(), ExcInternalError());
2799 
2800  // left child
2801  neighbor = neighbor->child(0);
2802  while (true)
2803  {
2804  Assert(neighbor->neighbor(0) == cell->child(1),
2805  ExcInternalError());
2806  neighbor->set_neighbor(0, cell);
2807 
2808  // move on to further
2809  // children on the
2810  // boundary between this
2811  // cell and its neighbor
2812  if (neighbor->has_children())
2813  neighbor = neighbor->child(0);
2814  else
2815  break;
2816  }
2817  }
2818 
2819 
2820  // delete the vertex which will not
2821  // be needed anymore. This vertex
2822  // is the second of the first child
2823  triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
2824 
2825  // invalidate children. clear user
2826  // pointers, to avoid that they may
2827  // appear at unwanted places later
2828  // on...
2829  for (unsigned int child = 0; child < cell->n_children(); ++child)
2830  {
2831  cell->child(child)->clear_user_data();
2832  cell->child(child)->clear_user_flag();
2833  cell->child(child)->clear_used_flag();
2834  }
2835 
2836 
2837  // delete pointer to children
2838  cell->clear_children();
2839  cell->clear_user_flag();
2840  }
2841 
2842 
2843 
2844  template <int spacedim>
2845  static void
2848  std::vector<unsigned int> &line_cell_count,
2849  std::vector<unsigned int> &)
2850  {
2851  const unsigned int dim = 2;
2852  const RefinementCase<dim> ref_case = cell->refinement_case();
2853 
2854  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2855  ExcInternalError());
2856 
2857  // vectors to hold all lines which
2858  // may be deleted
2859  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2860  lines_to_delete(0);
2861 
2862  lines_to_delete.reserve(4 * 2 + 4);
2863 
2864  // now we decrease the counters for
2865  // lines contained in the child
2866  // cells
2867  for (unsigned int c = 0; c < cell->n_children(); ++c)
2868  {
2870  cell->child(c);
2871  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2872  --line_cell_count[child->line_index(l)];
2873  }
2874 
2875 
2876  // delete the vertex which will not
2877  // be needed anymore. This vertex
2878  // is the second of the second line
2879  // of the first child, if the cell
2880  // is refined with cut_xy, else there
2881  // is no inner vertex.
2882  // additionally delete unneeded inner
2883  // lines
2884  if (ref_case == RefinementCase<dim>::cut_xy)
2885  {
2887  .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
2888 
2889  lines_to_delete.push_back(cell->child(0)->line(1));
2890  lines_to_delete.push_back(cell->child(0)->line(3));
2891  lines_to_delete.push_back(cell->child(3)->line(0));
2892  lines_to_delete.push_back(cell->child(3)->line(2));
2893  }
2894  else
2895  {
2896  unsigned int inner_face_no =
2897  ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
2898 
2899  // the inner line will not be
2900  // used any more
2901  lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2902  }
2903 
2904  // invalidate children
2905  for (unsigned int child = 0; child < cell->n_children(); ++child)
2906  {
2907  cell->child(child)->clear_user_data();
2908  cell->child(child)->clear_user_flag();
2909  cell->child(child)->clear_used_flag();
2910  }
2911 
2912 
2913  // delete pointer to children
2914  cell->clear_children();
2915  cell->clear_refinement_case();
2916  cell->clear_user_flag();
2917 
2918  // look at the refinement of outer
2919  // lines. if nobody needs those
2920  // anymore we can add them to the
2921  // list of lines to be deleted.
2922  for (unsigned int line_no = 0;
2923  line_no < GeometryInfo<dim>::lines_per_cell;
2924  ++line_no)
2925  {
2927  cell->line(line_no);
2928 
2929  if (line->has_children())
2930  {
2931  // if one of the cell counters is
2932  // zero, the other has to be as well
2933 
2934  Assert((line_cell_count[line->child_index(0)] == 0 &&
2935  line_cell_count[line->child_index(1)] == 0) ||
2936  (line_cell_count[line->child_index(0)] > 0 &&
2937  line_cell_count[line->child_index(1)] > 0),
2938  ExcInternalError());
2939 
2940  if (line_cell_count[line->child_index(0)] == 0)
2941  {
2942  for (unsigned int c = 0; c < 2; ++c)
2943  Assert(!line->child(c)->has_children(),
2944  ExcInternalError());
2945 
2946  // we may delete the line's
2947  // children and the middle vertex
2948  // as no cell references them
2949  // anymore
2951  .vertices_used[line->child(0)->vertex_index(1)] = false;
2952 
2953  lines_to_delete.push_back(line->child(0));
2954  lines_to_delete.push_back(line->child(1));
2955 
2956  line->clear_children();
2957  }
2958  }
2959  }
2960 
2961  // finally, delete unneeded lines
2962 
2963  // clear user pointers, to avoid that
2964  // they may appear at unwanted places
2965  // later on...
2966  // same for user flags, then finally
2967  // delete the lines
2968  typename std::vector<
2970  line = lines_to_delete.begin(),
2971  endline = lines_to_delete.end();
2972  for (; line != endline; ++line)
2973  {
2974  (*line)->clear_user_data();
2975  (*line)->clear_user_flag();
2976  (*line)->clear_used_flag();
2977  }
2978  }
2979 
2980 
2981 
2982  template <int spacedim>
2983  static void
2986  std::vector<unsigned int> &line_cell_count,
2987  std::vector<unsigned int> &quad_cell_count)
2988  {
2989  const unsigned int dim = 3;
2990 
2991  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2992  ExcInternalError());
2993  Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
2994  ExcInternalError());
2995 
2996  // first of all, we store the RefineCase of
2997  // this cell
2998  const RefinementCase<dim> ref_case = cell->refinement_case();
2999  // vectors to hold all lines and quads which
3000  // may be deleted
3001  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3002  lines_to_delete(0);
3003  std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3004  quads_to_delete(0);
3005 
3006  lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3007  quads_to_delete.reserve(6 * 4 + 12);
3008 
3009  // now we decrease the counters for lines and
3010  // quads contained in the child cells
3011  for (unsigned int c = 0; c < cell->n_children(); ++c)
3012  {
3014  cell->child(c);
3015  const auto line_indices = TriaAccessorImplementation::
3016  Implementation::get_line_indices_of_cell(*child);
3017  for (const unsigned int l : cell->line_indices())
3018  --line_cell_count[line_indices[l]];
3019  for (auto f : GeometryInfo<dim>::face_indices())
3020  --quad_cell_count[child->quad_index(f)];
3021  }
3022 
3023  //-------------------------------------
3024  // delete interior quads and lines and the
3025  // interior vertex, depending on the
3026  // refinement case of the cell
3027  //
3028  // for append quads and lines: only append
3029  // them to the list of objects to be deleted
3030 
3031  switch (ref_case)
3032  {
3034  quads_to_delete.push_back(cell->child(0)->face(1));
3035  break;
3037  quads_to_delete.push_back(cell->child(0)->face(3));
3038  break;
3040  quads_to_delete.push_back(cell->child(0)->face(5));
3041  break;
3043  quads_to_delete.push_back(cell->child(0)->face(1));
3044  quads_to_delete.push_back(cell->child(0)->face(3));
3045  quads_to_delete.push_back(cell->child(3)->face(0));
3046  quads_to_delete.push_back(cell->child(3)->face(2));
3047 
3048  lines_to_delete.push_back(cell->child(0)->line(11));
3049  break;
3051  quads_to_delete.push_back(cell->child(0)->face(1));
3052  quads_to_delete.push_back(cell->child(0)->face(5));
3053  quads_to_delete.push_back(cell->child(3)->face(0));
3054  quads_to_delete.push_back(cell->child(3)->face(4));
3055 
3056  lines_to_delete.push_back(cell->child(0)->line(5));
3057  break;
3059  quads_to_delete.push_back(cell->child(0)->face(3));
3060  quads_to_delete.push_back(cell->child(0)->face(5));
3061  quads_to_delete.push_back(cell->child(3)->face(2));
3062  quads_to_delete.push_back(cell->child(3)->face(4));
3063 
3064  lines_to_delete.push_back(cell->child(0)->line(7));
3065  break;
3067  quads_to_delete.push_back(cell->child(0)->face(1));
3068  quads_to_delete.push_back(cell->child(2)->face(1));
3069  quads_to_delete.push_back(cell->child(4)->face(1));
3070  quads_to_delete.push_back(cell->child(6)->face(1));
3071 
3072  quads_to_delete.push_back(cell->child(0)->face(3));
3073  quads_to_delete.push_back(cell->child(1)->face(3));
3074  quads_to_delete.push_back(cell->child(4)->face(3));
3075  quads_to_delete.push_back(cell->child(5)->face(3));
3076 
3077  quads_to_delete.push_back(cell->child(0)->face(5));
3078  quads_to_delete.push_back(cell->child(1)->face(5));
3079  quads_to_delete.push_back(cell->child(2)->face(5));
3080  quads_to_delete.push_back(cell->child(3)->face(5));
3081 
3082  lines_to_delete.push_back(cell->child(0)->line(5));
3083  lines_to_delete.push_back(cell->child(0)->line(7));
3084  lines_to_delete.push_back(cell->child(0)->line(11));
3085  lines_to_delete.push_back(cell->child(7)->line(0));
3086  lines_to_delete.push_back(cell->child(7)->line(2));
3087  lines_to_delete.push_back(cell->child(7)->line(8));
3088  // delete the vertex which will not
3089  // be needed anymore. This vertex
3090  // is the vertex at the heart of
3091  // this cell, which is the sixth of
3092  // the first child
3093  triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3094  false;
3095  break;
3096  default:
3097  // only remaining case is
3098  // no_refinement, thus an error
3099  Assert(false, ExcInternalError());
3100  break;
3101  }
3102 
3103 
3104  // invalidate children
3105  for (unsigned int child = 0; child < cell->n_children(); ++child)
3106  {
3107  cell->child(child)->clear_user_data();
3108  cell->child(child)->clear_user_flag();
3109 
3110  for (auto f : GeometryInfo<dim>::face_indices())
3111  // set flags denoting deviations from standard orientation of
3112  // faces back to initialization values
3113  cell->child(child)->set_combined_face_orientation(
3115 
3116  cell->child(child)->clear_used_flag();
3117  }
3118 
3119 
3120  // delete pointer to children
3121  cell->clear_children();
3122  cell->clear_refinement_case();
3123  cell->clear_user_flag();
3124 
3125  // so far we only looked at inner quads,
3126  // lines and vertices. Now we have to
3127  // consider outer ones as well. here, we have
3128  // to check, whether there are other cells
3129  // still needing these objects. otherwise we
3130  // can delete them. first for quads (and
3131  // their inner lines).
3132 
3133  for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
3134  {
3136  cell->face(quad_no);
3137 
3138  Assert(
3139  (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
3140  quad->has_children()) ||
3141  GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
3143  ExcInternalError());
3144 
3145  switch (quad->refinement_case())
3146  {
3147  case RefinementCase<dim - 1>::no_refinement:
3148  // nothing to do as the quad
3149  // is not refined
3150  break;
3151  case RefinementCase<dim - 1>::cut_x:
3152  case RefinementCase<dim - 1>::cut_y:
3153  {
3154  // if one of the cell counters is
3155  // zero, the other has to be as
3156  // well
3157  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3158  quad_cell_count[quad->child_index(1)] == 0) ||
3159  (quad_cell_count[quad->child_index(0)] > 0 &&
3160  quad_cell_count[quad->child_index(1)] > 0),
3161  ExcInternalError());
3162  // it might be, that the quad is
3163  // refined twice anisotropically,
3164  // first check, whether we may
3165  // delete possible grand_children
3166  unsigned int deleted_grandchildren = 0;
3167  unsigned int number_of_child_refinements = 0;
3168 
3169  for (unsigned int c = 0; c < 2; ++c)
3170  if (quad->child(c)->has_children())
3171  {
3172  ++number_of_child_refinements;
3173  // if one of the cell counters is
3174  // zero, the other has to be as
3175  // well
3176  Assert(
3177  (quad_cell_count[quad->child(c)->child_index(0)] ==
3178  0 &&
3179  quad_cell_count[quad->child(c)->child_index(1)] ==
3180  0) ||
3181  (quad_cell_count[quad->child(c)->child_index(0)] >
3182  0 &&
3183  quad_cell_count[quad->child(c)->child_index(1)] >
3184  0),
3185  ExcInternalError());
3186  if (quad_cell_count[quad->child(c)->child_index(0)] ==
3187  0)
3188  {
3189  // Assert, that the two
3190  // anisotropic
3191  // refinements add up to
3192  // isotropic refinement
3193  Assert(quad->refinement_case() +
3194  quad->child(c)->refinement_case() ==
3196  ExcInternalError());
3197  // we may delete the
3198  // quad's children and
3199  // the inner line as no
3200  // cell references them
3201  // anymore
3202  quads_to_delete.push_back(
3203  quad->child(c)->child(0));
3204  quads_to_delete.push_back(
3205  quad->child(c)->child(1));
3206  if (quad->child(c)->refinement_case() ==
3208  lines_to_delete.push_back(
3209  quad->child(c)->child(0)->line(1));
3210  else
3211  lines_to_delete.push_back(
3212  quad->child(c)->child(0)->line(3));
3213  quad->child(c)->clear_children();
3214  quad->child(c)->clear_refinement_case();
3215  ++deleted_grandchildren;
3216  }
3217  }
3218  // if no grandchildren are left, we
3219  // may as well delete the
3220  // refinement of the inner line
3221  // between our children and the
3222  // corresponding vertex
3223  if (number_of_child_refinements > 0 &&
3224  deleted_grandchildren == number_of_child_refinements)
3225  {
3227  middle_line;
3228  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3229  middle_line = quad->child(0)->line(1);
3230  else
3231  middle_line = quad->child(0)->line(3);
3232 
3233  lines_to_delete.push_back(middle_line->child(0));
3234  lines_to_delete.push_back(middle_line->child(1));
3236  .vertices_used[middle_vertex_index<dim, spacedim>(
3237  middle_line)] = false;
3238  middle_line->clear_children();
3239  }
3240 
3241  // now consider the direct children
3242  // of the given quad
3243  if (quad_cell_count[quad->child_index(0)] == 0)
3244  {
3245  // we may delete the quad's
3246  // children and the inner line
3247  // as no cell references them
3248  // anymore
3249  quads_to_delete.push_back(quad->child(0));
3250  quads_to_delete.push_back(quad->child(1));
3251  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3252  lines_to_delete.push_back(quad->child(0)->line(1));
3253  else
3254  lines_to_delete.push_back(quad->child(0)->line(3));
3255 
3256  // if the counters just dropped
3257  // to zero, otherwise the
3258  // children would have been
3259  // deleted earlier, then this
3260  // cell's children must have
3261  // contained the anisotropic
3262  // quad children. thus, if
3263  // those have again anisotropic
3264  // children, which are in
3265  // effect isotropic children of
3266  // the original quad, those are
3267  // still needed by a
3268  // neighboring cell and we
3269  // cannot delete them. instead,
3270  // we have to reset this quad's
3271  // refine case to isotropic and
3272  // set the children
3273  // accordingly.
3274  if (quad->child(0)->has_children())
3275  if (quad->refinement_case() ==
3277  {
3278  // now evereything is
3279  // quite complicated. we
3280  // have the children
3281  // numbered according to
3282  //
3283  // *---*---*
3284  // |n+1|m+1|
3285  // *---*---*
3286  // | n | m |
3287  // *---*---*
3288  //
3289  // from the original
3290  // anisotropic
3291  // refinement. we have to
3292  // reorder them as
3293  //
3294  // *---*---*
3295  // | m |m+1|
3296  // *---*---*
3297  // | n |n+1|
3298  // *---*---*
3299  //
3300  // for isotropic refinement.
3301  //
3302  // this is a bit ugly, of
3303  // course: loop over all
3304  // cells on all levels
3305  // and look for faces n+1
3306  // (switch_1) and m
3307  // (switch_2).
3308  const typename Triangulation<dim, spacedim>::
3309  quad_iterator switch_1 =
3310  quad->child(0)->child(1),
3311  switch_2 =
3312  quad->child(1)->child(0);
3313 
3314  Assert(!switch_1->has_children(),
3315  ExcInternalError());
3316  Assert(!switch_2->has_children(),
3317  ExcInternalError());
3318 
3319  const int switch_1_index = switch_1->index();
3320  const int switch_2_index = switch_2->index();
3321  for (unsigned int l = 0;
3322  l < triangulation.levels.size();
3323  ++l)
3324  for (unsigned int h = 0;
3325  h <
3326  triangulation.levels[l]->cells.n_objects();
3327  ++h)
3328  for (const unsigned int q :
3330  {
3331  const int index =
3332  triangulation.levels[l]
3333  ->cells.get_bounding_object_indices(
3334  h)[q];
3335  if (index == switch_1_index)
3336  triangulation.levels[l]
3337  ->cells.get_bounding_object_indices(
3338  h)[q] = switch_2_index;
3339  else if (index == switch_2_index)
3340  triangulation.levels[l]
3341  ->cells.get_bounding_object_indices(
3342  h)[q] = switch_1_index;
3343  }
3344  // now we have to copy
3345  // all information of the
3346  // two quads
3347  const int switch_1_lines[4] = {
3348  static_cast<signed int>(
3349  switch_1->line_index(0)),
3350  static_cast<signed int>(
3351  switch_1->line_index(1)),
3352  static_cast<signed int>(
3353  switch_1->line_index(2)),
3354  static_cast<signed int>(
3355  switch_1->line_index(3))};
3356  const bool switch_1_line_orientations[4] = {
3357  switch_1->line_orientation(0),
3358  switch_1->line_orientation(1),
3359  switch_1->line_orientation(2),
3360  switch_1->line_orientation(3)};
3361  const types::boundary_id switch_1_boundary_id =
3362  switch_1->boundary_id();
3363  const unsigned int switch_1_user_index =
3364  switch_1->user_index();
3365  const bool switch_1_user_flag =
3366  switch_1->user_flag_set();
3367 
3368  switch_1->set_bounding_object_indices(
3369  {switch_2->line_index(0),
3370  switch_2->line_index(1),
3371  switch_2->line_index(2),
3372  switch_2->line_index(3)});
3373  switch_1->set_line_orientation(
3374  0, switch_2->line_orientation(0));
3375  switch_1->set_line_orientation(
3376  1, switch_2->line_orientation(1));
3377  switch_1->set_line_orientation(
3378  2, switch_2->line_orientation(2));
3379  switch_1->set_line_orientation(
3380  3, switch_2->line_orientation(3));
3381  switch_1->set_boundary_id_internal(
3382  switch_2->boundary_id());
3383  switch_1->set_manifold_id(
3384  switch_2->manifold_id());
3385  switch_1->set_user_index(switch_2->user_index());
3386  if (switch_2->user_flag_set())
3387  switch_1->set_user_flag();
3388  else
3389  switch_1->clear_user_flag();
3390 
3391  switch_2->set_bounding_object_indices(
3392  {switch_1_lines[0],
3393  switch_1_lines[1],
3394  switch_1_lines[2],
3395  switch_1_lines[3]});
3396  switch_2->set_line_orientation(
3397  0, switch_1_line_orientations[0]);
3398  switch_2->set_line_orientation(
3399  1, switch_1_line_orientations[1]);
3400  switch_2->set_line_orientation(
3401  2, switch_1_line_orientations[2]);
3402  switch_2->set_line_orientation(
3403  3, switch_1_line_orientations[3]);
3404  switch_2->set_boundary_id_internal(
3405  switch_1_boundary_id);
3406  switch_2->set_manifold_id(
3407  switch_1->manifold_id());
3408  switch_2->set_user_index(switch_1_user_index);
3409  if (switch_1_user_flag)
3410  switch_2->set_user_flag();
3411  else
3412  switch_2->clear_user_flag();
3413 
3414  const unsigned int child_0 =
3415  quad->child(0)->child_index(0);
3416  const unsigned int child_2 =
3417  quad->child(1)->child_index(0);
3418  quad->clear_children();
3419  quad->clear_refinement_case();
3420  quad->set_refinement_case(
3422  quad->set_children(0, child_0);
3423  quad->set_children(2, child_2);
3424  std::swap(quad_cell_count[child_0 + 1],
3425  quad_cell_count[child_2]);
3426  }
3427  else
3428  {
3429  // the face was refined
3430  // with cut_y, thus the
3431  // children are already
3432  // in correct order. we
3433  // only have to set them
3434  // correctly, deleting
3435  // the indirection of two
3436  // anisotropic refinement
3437  // and going directly
3438  // from the quad to
3439  // isotropic children
3440  const unsigned int child_0 =
3441  quad->child(0)->child_index(0);
3442  const unsigned int child_2 =
3443  quad->child(1)->child_index(0);
3444  quad->clear_children();
3445  quad->clear_refinement_case();
3446  quad->set_refinement_case(
3448  quad->set_children(0, child_0);
3449  quad->set_children(2, child_2);
3450  }
3451  else
3452  {
3453  quad->clear_children();
3454  quad->clear_refinement_case();
3455  }
3456  }
3457  break;
3458  }
3459  case RefinementCase<dim - 1>::cut_xy:
3460  {
3461  // if one of the cell counters is
3462  // zero, the others have to be as
3463  // well
3464 
3465  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3466  quad_cell_count[quad->child_index(1)] == 0 &&
3467  quad_cell_count[quad->child_index(2)] == 0 &&
3468  quad_cell_count[quad->child_index(3)] == 0) ||
3469  (quad_cell_count[quad->child_index(0)] > 0 &&
3470  quad_cell_count[quad->child_index(1)] > 0 &&
3471  quad_cell_count[quad->child_index(2)] > 0 &&
3472  quad_cell_count[quad->child_index(3)] > 0),
3473  ExcInternalError());
3474 
3475  if (quad_cell_count[quad->child_index(0)] == 0)
3476  {
3477  // we may delete the quad's
3478  // children, the inner lines
3479  // and the middle vertex as no
3480  // cell references them anymore
3481  lines_to_delete.push_back(quad->child(0)->line(1));
3482  lines_to_delete.push_back(quad->child(3)->line(0));
3483  lines_to_delete.push_back(quad->child(0)->line(3));
3484  lines_to_delete.push_back(quad->child(3)->line(2));
3485 
3486  for (unsigned int child = 0; child < quad->n_children();
3487  ++child)
3488  quads_to_delete.push_back(quad->child(child));
3489 
3491  .vertices_used[quad->child(0)->vertex_index(3)] =
3492  false;
3493 
3494  quad->clear_children();
3495  quad->clear_refinement_case();
3496  }
3497  }
3498  break;
3499 
3500  default:
3501  Assert(false, ExcInternalError());
3502  break;
3503  }
3504  }
3505 
3506  // now we repeat a similar procedure
3507  // for the outer lines of this cell.
3508 
3509  // if in debug mode: check that each
3510  // of the lines for which we consider
3511  // deleting the children in fact has
3512  // children (the bits/coarsening_3d
3513  // test tripped over this initially)
3514  for (unsigned int line_no = 0;
3515  line_no < GeometryInfo<dim>::lines_per_cell;
3516  ++line_no)
3517  {
3519  cell->line(line_no);
3520 
3521  Assert(
3522  (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
3523  line->has_children()) ||
3524  GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
3526  ExcInternalError());
3527 
3528  if (line->has_children())
3529  {
3530  // if one of the cell counters is
3531  // zero, the other has to be as well
3532 
3533  Assert((line_cell_count[line->child_index(0)] == 0 &&
3534  line_cell_count[line->child_index(1)] == 0) ||
3535  (line_cell_count[line->child_index(0)] > 0 &&
3536  line_cell_count[line->child_index(1)] > 0),
3537  ExcInternalError());
3538 
3539  if (line_cell_count[line->child_index(0)] == 0)
3540  {
3541  for (unsigned int c = 0; c < 2; ++c)
3542  Assert(!line->child(c)->has_children(),
3543  ExcInternalError());
3544 
3545  // we may delete the line's
3546  // children and the middle vertex
3547  // as no cell references them
3548  // anymore
3550  .vertices_used[line->child(0)->vertex_index(1)] = false;
3551 
3552  lines_to_delete.push_back(line->child(0));
3553  lines_to_delete.push_back(line->child(1));
3554 
3555  line->clear_children();
3556  }
3557  }
3558  }
3559 
3560  // finally, delete unneeded quads and lines
3561 
3562  // clear user pointers, to avoid that
3563  // they may appear at unwanted places
3564  // later on...
3565  // same for user flags, then finally
3566  // delete the quads and lines
3567  typename std::vector<
3569  line = lines_to_delete.begin(),
3570  endline = lines_to_delete.end();
3571  for (; line != endline; ++line)
3572  {
3573  (*line)->clear_user_data();
3574  (*line)->clear_user_flag();
3575  (*line)->clear_used_flag();
3576  }
3577 
3578  typename std::vector<
3580  quad = quads_to_delete.begin(),
3581  endquad = quads_to_delete.end();
3582  for (; quad != endquad; ++quad)
3583  {
3584  (*quad)->clear_user_data();
3585  (*quad)->clear_children();
3586  (*quad)->clear_refinement_case();
3587  (*quad)->clear_user_flag();
3588  (*quad)->clear_used_flag();
3589  }
3590  }
3591 
3592 
3610  template <int spacedim>
3611  static void
3614  unsigned int & next_unused_vertex,
3616  &next_unused_line,
3618  &next_unused_cell,
3619  const typename Triangulation<2, spacedim>::cell_iterator &cell)
3620  {
3621  const unsigned int dim = 2;
3622  // clear refinement flag
3623  const RefinementCase<dim> ref_case = cell->refine_flag_set();
3624  cell->clear_refine_flag();
3625 
3626  /* For the refinement process: since we go the levels up from the
3627  lowest, there are (unlike above) only two possibilities: a neighbor
3628  cell is on the same level or one level up (in both cases, it may or
3629  may not be refined later on, but we don't care here).
3630 
3631  First:
3632  Set up an array of the 3x3 vertices, which are distributed on the
3633  cell (the array consists of indices into the @p{vertices} std::vector
3634 
3635  2--7--3
3636  | | |
3637  4--8--5
3638  | | |
3639  0--6--1
3640 
3641  note: in case of cut_x or cut_y not all these vertices are needed for
3642  the new cells
3643 
3644  Second:
3645  Set up an array of the new lines (the array consists of iterator
3646  pointers into the lines arrays)
3647 
3648  .-6-.-7-. The directions are: .->-.->-.
3649  1 9 3 ^ ^ ^
3650  .-10.11-. .->-.->-.
3651  0 8 2 ^ ^ ^
3652  .-4-.-5-. .->-.->-.
3653 
3654  cut_x:
3655  .-4-.-5-.
3656  | | |
3657  0 6 1
3658  | | |
3659  .-2-.-3-.
3660 
3661  cut_y:
3662  .---5---.
3663  1 3
3664  .---6---.
3665  0 2
3666  .---4---.
3667 
3668 
3669  Third:
3670  Set up an array of neighbors:
3671 
3672  6 7
3673  .--.--.
3674  1| | |3
3675  .--.--.
3676  0| | |2
3677  .--.--.
3678  4 5
3679 
3680  We need this array for two reasons: first to get the lines which will
3681  bound the four subcells (if the neighboring cell is refined, these
3682  lines already exist), and second to update neighborship information.
3683  Since if a neighbor is not refined, its neighborship record only
3684  points to the present, unrefined, cell rather than the children we
3685  are presently creating, we only need the neighborship information
3686  if the neighbor cells are refined. In all other cases, we store
3687  the unrefined neighbor address
3688 
3689  We also need for every neighbor (if refined) which number among its
3690  neighbors the present (unrefined) cell has, since that number is to
3691  be replaced and because that also is the number of the subline which
3692  will be the interface between that neighbor and the to be created
3693  cell. We will store this number (between 0 and 3) in the field
3694  @p{neighbors_neighbor}.
3695 
3696  It would be sufficient to use the children of the common line to the
3697  neighbor, if we only wanted to get the new sublines and the new
3698  vertex, but because we need to update the neighborship information of
3699  the two refined subcells of the neighbor, we need to search these
3700  anyway.
3701 
3702  Convention:
3703  The created children are numbered like this:
3704 
3705  .--.--.
3706  |2 . 3|
3707  .--.--.
3708  |0 | 1|
3709  .--.--.
3710  */
3711  // collect the indices of the eight surrounding vertices
3712  // 2--7--3
3713  // | | |
3714  // 4--8--5
3715  // | | |
3716  // 0--6--1
3717  int new_vertices[9];
3718  for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3719  new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3720  for (unsigned int line_no = 0; line_no < 4; ++line_no)
3721  if (cell->line(line_no)->has_children())
3722  new_vertices[4 + line_no] =
3723  cell->line(line_no)->child(0)->vertex_index(1);
3724 
3725  if (ref_case == RefinementCase<dim>::cut_xy)
3726  {
3727  // find the next
3728  // unused vertex and
3729  // allocate it for
3730  // the new vertex we
3731  // need here
3732  while (triangulation.vertices_used[next_unused_vertex] == true)
3733  ++next_unused_vertex;
3734  Assert(next_unused_vertex < triangulation.vertices.size(),
3735  ExcMessage(
3736  "Internal error: During refinement, the triangulation "
3737  "wants to access an element of the 'vertices' array "
3738  "but it turns out that the array is not large enough."));
3739  triangulation.vertices_used[next_unused_vertex] = true;
3740 
3741  new_vertices[8] = next_unused_vertex;
3742 
3743  // determine middle vertex by transfinite interpolation to be
3744  // consistent with what happens to quads in a
3745  // Triangulation<3,3> when they are refined
3746  triangulation.vertices[next_unused_vertex] =
3747  cell->center(true, true);
3748  }
3749 
3750 
3751  // Now the lines:
3752  typename Triangulation<dim, spacedim>::raw_line_iterator new_lines[12];
3753  unsigned int lmin = 8;
3754  unsigned int lmax = 12;
3755  if (ref_case != RefinementCase<dim>::cut_xy)
3756  {
3757  lmin = 6;
3758  lmax = 7;
3759  }
3760 
3761  for (unsigned int l = lmin; l < lmax; ++l)
3762  {
3763  while (next_unused_line->used() == true)
3764  ++next_unused_line;
3765  new_lines[l] = next_unused_line;
3766  ++next_unused_line;
3767 
3768  AssertIsNotUsed(new_lines[l]);
3769  }
3770 
3771  if (ref_case == RefinementCase<dim>::cut_xy)
3772  {
3773  // .-6-.-7-.
3774  // 1 9 3
3775  // .-10.11-.
3776  // 0 8 2
3777  // .-4-.-5-.
3778 
3779  // lines 0-7 already exist, create only the four interior
3780  // lines 8-11
3781  unsigned int l = 0;
3782  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
3783  for (unsigned int c = 0; c < 2; ++c, ++l)
3784  new_lines[l] = cell->line(face_no)->child(c);
3785  Assert(l == 8, ExcInternalError());
3786 
3787  new_lines[8]->set_bounding_object_indices(
3788  {new_vertices[6], new_vertices[8]});
3789  new_lines[9]->set_bounding_object_indices(
3790  {new_vertices[8], new_vertices[7]});
3791  new_lines[10]->set_bounding_object_indices(
3792  {new_vertices[4], new_vertices[8]});
3793  new_lines[11]->set_bounding_object_indices(
3794  {new_vertices[8], new_vertices[5]});
3795  }
3796  else if (ref_case == RefinementCase<dim>::cut_x)
3797  {
3798  // .-4-.-5-.
3799  // | | |
3800  // 0 6 1
3801  // | | |
3802  // .-2-.-3-.
3803  new_lines[0] = cell->line(0);
3804  new_lines[1] = cell->line(1);
3805  new_lines[2] = cell->line(2)->child(0);
3806  new_lines[3] = cell->line(2)->child(1);
3807  new_lines[4] = cell->line(3)->child(0);
3808  new_lines[5] = cell->line(3)->child(1);
3809  new_lines[6]->set_bounding_object_indices(
3810  {new_vertices[6], new_vertices[7]});
3811  }
3812  else
3813  {
3815  // .---5---.
3816  // 1 3
3817  // .---6---.
3818  // 0 2
3819  // .---4---.
3820  new_lines[0] = cell->line(0)->child(0);
3821  new_lines[1] = cell->line(0)->child(1);
3822  new_lines[2] = cell->line(1)->child(0);
3823  new_lines[3] = cell->line(1)->child(1);
3824  new_lines[4] = cell->line(2);
3825  new_lines[5] = cell->line(3);
3826  new_lines[6]->set_bounding_object_indices(
3827  {new_vertices[4], new_vertices[5]});
3828  }
3829 
3830  for (unsigned int l = lmin; l < lmax; ++l)
3831  {
3832  new_lines[l]->set_used_flag();
3833  new_lines[l]->clear_user_flag();
3834  new_lines[l]->clear_user_data();
3835  new_lines[l]->clear_children();
3836  // interior line
3837  new_lines[l]->set_boundary_id_internal(
3839  new_lines[l]->set_manifold_id(cell->manifold_id());
3840  }
3841 
3842  // Now add the four (two)
3843  // new cells!
3846  while (next_unused_cell->used() == true)
3847  ++next_unused_cell;
3848 
3849  const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
3850  for (unsigned int i = 0; i < n_children; ++i)
3851  {
3852  AssertIsNotUsed(next_unused_cell);
3853  subcells[i] = next_unused_cell;
3854  ++next_unused_cell;
3855  if (i % 2 == 1 && i < n_children - 1)
3856  while (next_unused_cell->used() == true)
3857  ++next_unused_cell;
3858  }
3859 
3860  if (ref_case == RefinementCase<dim>::cut_xy)
3861  {
3862  // children:
3863  // .--.--.
3864  // |2 . 3|
3865  // .--.--.
3866  // |0 | 1|
3867  // .--.--.
3868  // lines:
3869  // .-6-.-7-.
3870  // 1 9 3
3871  // .-10.11-.
3872  // 0 8 2
3873  // .-4-.-5-.
3874  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3875  new_lines[8]->index(),
3876  new_lines[4]->index(),
3877  new_lines[10]->index()});
3878  subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3879  new_lines[2]->index(),
3880  new_lines[5]->index(),
3881  new_lines[11]->index()});
3882  subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3883  new_lines[9]->index(),
3884  new_lines[10]->index(),
3885  new_lines[6]->index()});
3886  subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3887  new_lines[3]->index(),
3888  new_lines[11]->index(),
3889  new_lines[7]->index()});
3890  }
3891  else if (ref_case == RefinementCase<dim>::cut_x)
3892  {
3893  // children:
3894  // .--.--.
3895  // | . |
3896  // .0 . 1.
3897  // | | |
3898  // .--.--.
3899  // lines:
3900  // .-4-.-5-.
3901  // | | |
3902  // 0 6 1
3903  // | | |
3904  // .-2-.-3-.
3905  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3906  new_lines[6]->index(),
3907  new_lines[2]->index(),
3908  new_lines[4]->index()});
3909  subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3910  new_lines[1]->index(),
3911  new_lines[3]->index(),
3912  new_lines[5]->index()});
3913  }
3914  else
3915  {
3917  // children:
3918  // .-----.
3919  // | 1 |
3920  // .-----.
3921  // | 0 |
3922  // .-----.
3923  // lines:
3924  // .---5---.
3925  // 1 3
3926  // .---6---.
3927  // 0 2
3928  // .---4---.
3929  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3930  new_lines[2]->index(),
3931  new_lines[4]->index(),
3932  new_lines[6]->index()});
3933  subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
3934  new_lines[3]->index(),
3935  new_lines[6]->index(),
3936  new_lines[5]->index()});
3937  }
3938 
3939  types::subdomain_id subdomainid = cell->subdomain_id();
3940 
3941  for (unsigned int i = 0; i < n_children; ++i)
3942  {
3943  subcells[i]->set_used_flag();
3944  subcells[i]->clear_refine_flag();
3945  subcells[i]->clear_user_flag();
3946  subcells[i]->clear_user_data();
3947  subcells[i]->clear_children();
3948  // inherit material properties
3949  subcells[i]->set_material_id(cell->material_id());
3950  subcells[i]->set_manifold_id(cell->manifold_id());
3951  subcells[i]->set_subdomain_id(subdomainid);
3952 
3953  if (i % 2 == 0)
3954  subcells[i]->set_parent(cell->index());
3955  }
3956 
3957 
3958 
3959  // set child index for even children i=0,2 (0)
3960  for (unsigned int i = 0; i < n_children / 2; ++i)
3961  cell->set_children(2 * i, subcells[2 * i]->index());
3962  // set the refine case
3963  cell->set_refinement_case(ref_case);
3964 
3965  // note that the
3966  // refinement flag was
3967  // already cleared at the
3968  // beginning of this function
3969 
3970  if (dim < spacedim)
3971  for (unsigned int c = 0; c < n_children; ++c)
3972  cell->child(c)->set_direction_flag(cell->direction_flag());
3973  }
3974 
3975 
3976 
3977  template <int dim, int spacedim>
3980  const bool check_for_distorted_cells)
3981  {
3982  AssertDimension(dim, 2);
3983 
3984  // Check whether a new level is needed. We have to check for
3985  // this on the highest level only
3986  for (const auto &cell : triangulation.active_cell_iterators_on_level(
3987  triangulation.levels.size() - 1))
3988  if (cell->refine_flag_set())
3989  {
3990  triangulation.levels.push_back(
3991  std::make_unique<
3993  break;
3994  }
3995 
3996  for (typename Triangulation<dim, spacedim>::line_iterator line =
3997  triangulation.begin_line();
3998  line != triangulation.end_line();
3999  ++line)
4000  {
4001  line->clear_user_flag();
4002  line->clear_user_data();
4003  }
4004 
4005  unsigned int n_single_lines = 0;
4006  unsigned int n_lines_in_pairs = 0;
4007  unsigned int needed_vertices = 0;
4008 
4009  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4010  {
4011  // count number of flagged cells on this level and compute
4012  // how many new vertices and new lines will be needed
4013  unsigned int needed_cells = 0;
4014 
4015  for (const auto &cell :
4016  triangulation.active_cell_iterators_on_level(level))
4017  if (cell->refine_flag_set())
4018  {
4019  if (cell->reference_cell() == ReferenceCells::Triangle)
4020  {
4021  needed_cells += 4;
4022  needed_vertices += 0;
4023  n_single_lines += 3;
4024  }
4025  else if (cell->reference_cell() ==
4027  {
4028  needed_cells += 4;
4029  needed_vertices += 1;
4030  n_single_lines += 4;
4031  }
4032  else
4033  {
4034  AssertThrow(false, ExcNotImplemented());
4035  }
4036 
4037  for (const auto line_no : cell->face_indices())
4038  {
4039  auto line = cell->line(line_no);
4040  if (line->has_children() == false)
4041  line->set_user_flag();
4042  }
4043  }
4044 
4045 
4046  const unsigned int used_cells =
4047  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4048  triangulation.levels[level + 1]->cells.used.end(),
4049  true);
4050 
4051 
4052  reserve_space(*triangulation.levels[level + 1],
4053  used_cells + needed_cells,
4054  2,
4055  spacedim);
4056 
4057  reserve_space(triangulation.levels[level + 1]->cells,
4058  needed_cells,
4059  0);
4060  }
4061 
4062  for (auto line = triangulation.begin_line();
4063  line != triangulation.end_line();
4064  ++line)
4065  if (line->user_flag_set())
4066  {
4067  Assert(line->has_children() == false, ExcInternalError());
4068  n_lines_in_pairs += 2;
4069  needed_vertices += 1;
4070  }
4071 
4072  reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4073 
4074  needed_vertices += std::count(triangulation.vertices_used.begin(),
4075  triangulation.vertices_used.end(),
4076  true);
4077 
4078  if (needed_vertices > triangulation.vertices.size())
4079  {
4080  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4081  triangulation.vertices_used.resize(needed_vertices, false);
4082  }
4083 
4084  unsigned int next_unused_vertex = 0;
4085 
4086  {
4088  line = triangulation.begin_active_line(),
4089  endl = triangulation.end_line();
4091  next_unused_line = triangulation.begin_raw_line();
4092 
4093  for (; line != endl; ++line)
4094  if (line->user_flag_set())
4095  {
4096  // this line needs to be refined
4097 
4098  // find the next unused vertex and set it
4099  // appropriately
4100  while (triangulation.vertices_used[next_unused_vertex] == true)
4101  ++next_unused_vertex;
4102  Assert(
4103  next_unused_vertex < triangulation.vertices.size(),
4104  ExcMessage(
4105  "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."));
4106  triangulation.vertices_used[next_unused_vertex] = true;
4107 
4108  triangulation.vertices[next_unused_vertex] = line->center(true);
4109 
4110  bool pair_found = false;
4111  (void)pair_found;
4112  for (; next_unused_line != endl; ++next_unused_line)
4113  if (!next_unused_line->used() &&
4114  !(++next_unused_line)->used())
4115  {
4116  --next_unused_line;
4117  pair_found = true;
4118  break;
4119  }
4120  Assert(pair_found, ExcInternalError());
4121 
4122  line->set_children(0, next_unused_line->index());
4123 
4125  children[2] = {next_unused_line, ++next_unused_line};
4126 
4127  AssertIsNotUsed(children[0]);
4128  AssertIsNotUsed(children[1]);
4129 
4130  children[0]->set_bounding_object_indices(
4131  {line->vertex_index(0), next_unused_vertex});
4132  children[1]->set_bounding_object_indices(
4133  {next_unused_vertex, line->vertex_index(1)});
4134 
4135  children[0]->set_used_flag();
4136  children[1]->set_used_flag();
4137  children[0]->clear_children();
4138  children[1]->clear_children();
4139  children[0]->clear_user_data();
4140  children[1]->clear_user_data();
4141  children[0]->clear_user_flag();
4142  children[1]->clear_user_flag();
4143 
4144 
4145  children[0]->set_boundary_id_internal(line->boundary_id());
4146  children[1]->set_boundary_id_internal(line->boundary_id());
4147 
4148  children[0]->set_manifold_id(line->manifold_id());
4149  children[1]->set_manifold_id(line->manifold_id());
4150 
4151  line->clear_user_flag();
4152  }
4153  }
4154 
4155  reserve_space(triangulation.faces->lines, 0, n_single_lines);
4156 
4158  cells_with_distorted_children;
4159 
4161  next_unused_line = triangulation.begin_raw_line();
4162 
4163  const auto create_children = [](auto & triangulation,
4164  unsigned int &next_unused_vertex,
4165  auto & next_unused_line,
4166  auto & next_unused_cell,
4167  const auto & cell) {
4168  const auto ref_case = cell->refine_flag_set();
4169  cell->clear_refine_flag();
4170 
4171  unsigned int n_new_vertices = 0;
4172 
4173  if (cell->reference_cell() == ReferenceCells::Triangle)
4174  n_new_vertices = 6;
4175  else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
4176  n_new_vertices = 9;
4177  else
4178  AssertThrow(false, ExcNotImplemented());
4179 
4180  std::vector<unsigned int> new_vertices(n_new_vertices,
4182  for (unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
4183  ++vertex_no)
4184  new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4185  for (unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
4186  if (cell->line(line_no)->has_children())
4187  new_vertices[cell->n_vertices() + line_no] =
4188  cell->line(line_no)->child(0)->vertex_index(1);
4189 
4190  if (cell->reference_cell() == ReferenceCells::Quadrilateral)
4191  {
4192  while (triangulation.vertices_used[next_unused_vertex] == true)
4193  ++next_unused_vertex;
4194  Assert(
4195  next_unused_vertex < triangulation.vertices.size(),
4196  ExcMessage(
4197  "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."));
4198  triangulation.vertices_used[next_unused_vertex] = true;
4199 
4200  new_vertices[8] = next_unused_vertex;
4201 
4202  triangulation.vertices[next_unused_vertex] =
4203  cell->center(true, true);
4204  }
4205 
4206  std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
4207  12>
4208  new_lines;
4209  unsigned int lmin = 0;
4210  unsigned int lmax = 0;
4211 
4212  if (cell->reference_cell() == ReferenceCells::Triangle)
4213  {
4214  lmin = 6;
4215  lmax = 9;
4216  }
4217  else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
4218  {
4219  lmin = 8;
4220  lmax = 12;
4221  }
4222  else
4223  {
4224  AssertThrow(false, ExcNotImplemented());
4225  }
4226 
4227  for (unsigned int l = lmin; l < lmax; ++l)
4228  {
4229  while (next_unused_line->used() == true)
4230  ++next_unused_line;
4231  new_lines[l] = next_unused_line;
4232  ++next_unused_line;
4233 
4234  AssertIsNotUsed(new_lines[l]);
4235  }
4236 
4237  if (cell->reference_cell() == ReferenceCells::Triangle)
4238  {
4239  // add lines in the order implied by their orientation. Here,
4240  // face_no is the cell (not subcell) face number and vertex_no is
4241  // the first vertex on that face in the standard orientation.
4242  const auto ref = [&](const unsigned int face_no,
4243  const unsigned int vertex_no) {
4244  auto l = cell->line(face_no);
4245  // if the vertex is on the first child then add the first child
4246  // first
4247  if (l->child(0)->vertex_index(0) == new_vertices[vertex_no] ||
4248  l->child(0)->vertex_index(1) == new_vertices[vertex_no])
4249  {
4250  new_lines[2 * face_no + 0] = l->child(0);
4251  new_lines[2 * face_no + 1] = l->child(1);
4252  }
4253  else
4254  {
4255  new_lines[2 * face_no + 0] = l->child(1);
4256  new_lines[2 * face_no + 1] = l->child(0);
4257  }
4258  };
4259 
4260  ref(0, 0);
4261  ref(1, 1);
4262  ref(2, 2);
4263 
4264  // set up lines which do not have parents:
4265  new_lines[6]->set_bounding_object_indices(
4266  {new_vertices[3], new_vertices[4]});
4267  new_lines[7]->set_bounding_object_indices(
4268  {new_vertices[4], new_vertices[5]});
4269  new_lines[8]->set_bounding_object_indices(
4270  {new_vertices[5], new_vertices[3]});
4271  }
4272  else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
4273  {
4274  unsigned int l = 0;
4275  for (const unsigned int face_no : cell->face_indices())
4276  for (unsigned int c = 0; c < 2; ++c, ++l)
4277  new_lines[l] = cell->line(face_no)->child(c);
4278 
4279  new_lines[8]->set_bounding_object_indices(
4280  {new_vertices[6], new_vertices[8]});
4281  new_lines[9]->set_bounding_object_indices(
4282  {new_vertices[8], new_vertices[7]});
4283  new_lines[10]->set_bounding_object_indices(
4284  {new_vertices[4], new_vertices[8]});
4285  new_lines[11]->set_bounding_object_indices(
4286  {new_vertices[8], new_vertices[5]});
4287  }
4288  else
4289  {
4290  AssertThrow(false, ExcNotImplemented());
4291  }
4292 
4293  for (unsigned int l = lmin; l < lmax; ++l)
4294  {
4295  new_lines[l]->set_used_flag();
4296  new_lines[l]->clear_user_flag();
4297  new_lines[l]->clear_user_data();
4298  new_lines[l]->clear_children();
4299  // interior line
4300  new_lines[l]->set_boundary_id_internal(
4302  new_lines[l]->set_manifold_id(cell->manifold_id());
4303  }
4304 
4307  while (next_unused_cell->used() == true)
4308  ++next_unused_cell;
4309 
4310  unsigned int n_children = 0;
4311 
4312  if (cell->reference_cell() == ReferenceCells::Triangle)
4313  n_children = 4;
4314  else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
4315  n_children = 4;
4316  else
4317  AssertThrow(false, ExcNotImplemented());
4318 
4319  for (unsigned int i = 0; i < n_children; ++i)
4320  {
4321  AssertIsNotUsed(next_unused_cell);
4322  subcells[i] = next_unused_cell;
4323  ++next_unused_cell;
4324  if (i % 2 == 1 && i < n_children - 1)
4325  while (next_unused_cell->used() == true)
4326  ++next_unused_cell;
4327  }
4328 
4329  if ((dim == 2) &&
4330  (cell->reference_cell() == ReferenceCells::Triangle))
4331  {
4332  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4333  new_lines[8]->index(),
4334  new_lines[5]->index()});
4335  subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4336  new_lines[2]->index(),
4337  new_lines[6]->index()});
4338  subcells[2]->set_bounding_object_indices({new_lines[7]->index(),
4339  new_lines[3]->index(),
4340  new_lines[4]->index()});
4341  subcells[3]->set_bounding_object_indices({new_lines[6]->index(),
4342  new_lines[7]->index(),
4343  new_lines[8]->index()});
4344 
4345  // Set subcell line orientations by checking the line's second
4346  // vertex (from the subcell's perspective) to the line's actual
4347  // second vertex.
4348  const auto fix_line_orientation =
4349  [&](const unsigned int line_no,
4350  const unsigned int vertex_no,
4351  const unsigned int subcell_no,
4352  const unsigned int subcell_line_no) {
4353  if (new_lines[line_no]->vertex_index(1) !=
4354  new_vertices[vertex_no])
4355  triangulation.levels[subcells[subcell_no]->level()]
4356  ->face_orientations.set_combined_orientation(
4357  subcells[subcell_no]->index() *
4359  subcell_line_no,
4360  0u);
4361  };
4362 
4363  fix_line_orientation(0, 3, 0, 0);
4364  fix_line_orientation(8, 5, 0, 1);
4365  fix_line_orientation(5, 0, 0, 2);
4366 
4367  fix_line_orientation(1, 1, 1, 0);
4368  fix_line_orientation(2, 4, 1, 1);
4369  fix_line_orientation(6, 3, 1, 2);
4370 
4371  fix_line_orientation(7, 4, 2, 0);
4372  fix_line_orientation(3, 2, 2, 1);
4373  fix_line_orientation(4, 5, 2, 2);
4374 
4375  // all lines of the new interior cell are oriented backwards so
4376  // that it has positive area.
4377  fix_line_orientation(6, 4, 3, 0);
4378  fix_line_orientation(7, 5, 3, 1);
4379  fix_line_orientation(8, 3, 3, 2);
4380  }
4381  else if ((dim == 2) &&
4382  (cell->reference_cell() == ReferenceCells::Quadrilateral))
4383  {
4384  subcells[0]->set_bounding_object_indices(
4385  {new_lines[0]->index(),
4386  new_lines[8]->index(),
4387  new_lines[4]->index(),
4388  new_lines[10]->index()});
4389  subcells[1]->set_bounding_object_indices(
4390  {new_lines[8]->index(),
4391  new_lines[2]->index(),
4392  new_lines[5]->index(),
4393  new_lines[11]->index()});
4394  subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
4395  new_lines[9]->index(),
4396  new_lines[10]->index(),
4397  new_lines[6]->index()});
4398  subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
4399  new_lines[3]->index(),
4400  new_lines[11]->index(),
4401  new_lines[7]->index()});
4402  }
4403  else
4404  {
4405  AssertThrow(false, ExcNotImplemented());
4406  }
4407 
4408  types::subdomain_id subdomainid = cell->subdomain_id();
4409 
4410  for (unsigned int i = 0; i < n_children; ++i)
4411  {
4412  subcells[i]->set_used_flag();
4413  subcells[i]->clear_refine_flag();
4414  subcells[i]->clear_user_flag();
4415  subcells[i]->clear_user_data();
4416  subcells[i]->clear_children();
4417  // inherit material
4418  // properties
4419  subcells[i]->set_material_id(cell->material_id());
4420  subcells[i]->set_manifold_id(cell->manifold_id());
4421  subcells[i]->set_subdomain_id(subdomainid);
4422 
4423  // TODO: here we assume that all children have the same reference
4424  // cell type as the parent! This is justified for 2d.
4425  triangulation.levels[subcells[i]->level()]
4426  ->reference_cell[subcells[i]->index()] = cell->reference_cell();
4427 
4428  if (i % 2 == 0)
4429  subcells[i]->set_parent(cell->index());
4430  }
4431 
4432  for (unsigned int i = 0; i < n_children / 2; ++i)
4433  cell->set_children(2 * i, subcells[2 * i]->index());
4434 
4435  cell->set_refinement_case(ref_case);
4436 
4437  if (dim < spacedim)
4438  for (unsigned int c = 0; c < n_children; ++c)
4439  cell->child(c)->set_direction_flag(cell->direction_flag());
4440  };
4441 
4442  for (int level = 0;
4443  level < static_cast<int>(triangulation.levels.size()) - 1;
4444  ++level)
4445  {
4447  next_unused_cell = triangulation.begin_raw(level + 1);
4448 
4449  for (const auto &cell :
4450  triangulation.active_cell_iterators_on_level(level))
4451  if (cell->refine_flag_set())
4452  {
4454  next_unused_vertex,
4455  next_unused_line,
4456  next_unused_cell,
4457  cell);
4458 
4459  if (cell->reference_cell() == ReferenceCells::Quadrilateral &&
4460  check_for_distorted_cells &&
4461  has_distorted_children<dim, spacedim>(cell))
4462  cells_with_distorted_children.distorted_cells.push_back(
4463  cell);
4464 
4465  triangulation.signals.post_refinement_on_cell(cell);
4466  }
4467  }
4468 
4469  return cells_with_distorted_children;
4470  }
4471 
4472 
4473 
4478  template <int spacedim>
4481  const bool /*check_for_distorted_cells*/)
4482  {
4483  const unsigned int dim = 1;
4484 
4485  // Check whether a new level is needed. We have to check for
4486  // this on the highest level only
4487  for (const auto &cell : triangulation.active_cell_iterators_on_level(
4488  triangulation.levels.size() - 1))
4489  if (cell->refine_flag_set())
4490  {
4491  triangulation.levels.push_back(
4492  std::make_unique<
4494  break;
4495  }
4496 
4497 
4498  // check how much space is needed on every level. We need not
4499  // check the highest level since either - on the highest level
4500  // no cells are flagged for refinement - there are, but
4501  // prepare_refinement added another empty level
4502  unsigned int needed_vertices = 0;
4503  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4504  {
4505  // count number of flagged
4506  // cells on this level
4507  unsigned int flagged_cells = 0;
4508 
4509  for (const auto &acell :
4510  triangulation.active_cell_iterators_on_level(level))
4511  if (acell->refine_flag_set())
4512  ++flagged_cells;
4513 
4514  // count number of used cells
4515  // on the next higher level
4516  const unsigned int used_cells =
4517  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4518  triangulation.levels[level + 1]->cells.used.end(),
4519  true);
4520 
4521  // reserve space for the used_cells cells already existing
4522  // on the next higher level as well as for the
4523  // 2*flagged_cells that will be created on that level
4524  reserve_space(*triangulation.levels[level + 1],
4526  flagged_cells,
4527  1,
4528  spacedim);
4529  // reserve space for 2*flagged_cells new lines on the next
4530  // higher level
4531  reserve_space(triangulation.levels[level + 1]->cells,
4533  flagged_cells,
4534  0);
4535 
4536  needed_vertices += flagged_cells;
4537  }
4538 
4539  // add to needed vertices how many
4540  // vertices are already in use
4541  needed_vertices += std::count(triangulation.vertices_used.begin(),
4542  triangulation.vertices_used.end(),
4543  true);
4544  // if we need more vertices: create them, if not: leave the
4545  // array as is, since shrinking is not really possible because
4546  // some of the vertices at the end may be in use
4547  if (needed_vertices > triangulation.vertices.size())
4548  {
4549  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4550  triangulation.vertices_used.resize(needed_vertices, false);
4551  }
4552 
4553 
4554  // Do REFINEMENT on every level; exclude highest level as
4555  // above
4556 
4557  // index of next unused vertex
4558  unsigned int next_unused_vertex = 0;
4559 
4560  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4561  {
4563  next_unused_cell = triangulation.begin_raw(level + 1);
4564 
4565  for (const auto &cell :
4566  triangulation.active_cell_iterators_on_level(level))
4567  if (cell->refine_flag_set())
4568  {
4569  // clear refinement flag
4570  cell->clear_refine_flag();
4571 
4572  // search for next unused
4573  // vertex
4574  while (triangulation.vertices_used[next_unused_vertex] ==
4575  true)
4576  ++next_unused_vertex;
4577  Assert(
4578  next_unused_vertex < triangulation.vertices.size(),
4579  ExcMessage(
4580  "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."));
4581 
4582  // Now we always ask the cell itself where to put
4583  // the new point. The cell in turn will query the
4584  // manifold object internally.
4585  triangulation.vertices[next_unused_vertex] =
4586  cell->center(true);
4587 
4588  triangulation.vertices_used[next_unused_vertex] = true;
4589 
4590  // search for next two unused cell (++ takes care of
4591  // the end of the vector)
4593  first_child,
4594  second_child;
4595  while (next_unused_cell->used() == true)
4596  ++next_unused_cell;
4597  first_child = next_unused_cell;
4598  first_child->set_used_flag();
4599  first_child->clear_user_data();
4600  ++next_unused_cell;
4601  AssertIsNotUsed(next_unused_cell);
4602  second_child = next_unused_cell;
4603  second_child->set_used_flag();
4604  second_child->clear_user_data();
4605 
4606  types::subdomain_id subdomainid = cell->subdomain_id();
4607 
4608  // insert first child
4609  cell->set_children(0, first_child->index());
4610  first_child->clear_children();
4611  first_child->set_bounding_object_indices(
4612  {cell->vertex_index(0), next_unused_vertex});
4613  first_child->set_material_id(cell->material_id());
4614  first_child->set_manifold_id(cell->manifold_id());
4615  first_child->set_subdomain_id(subdomainid);
4616  first_child->set_direction_flag(cell->direction_flag());
4617 
4618  first_child->set_parent(cell->index());
4619 
4620  // Set manifold id of the right face. Only do this
4621  // on the first child.
4622  first_child->face(1)->set_manifold_id(cell->manifold_id());
4623 
4624  // reset neighborship info (refer to
4625  // internal::TriangulationImplementation::TriaLevel<0> for
4626  // details)
4627  first_child->set_neighbor(1, second_child);
4628  if (cell->neighbor(0).state() != IteratorState::valid)
4629  first_child->set_neighbor(0, cell->neighbor(0));
4630  else if (cell->neighbor(0)->is_active())
4631  {
4632  // since the neighbors level is always <=level,
4633  // if the cell is active, then there are no
4634  // cells to the left which may want to know
4635  // about this new child cell.
4636  Assert(cell->neighbor(0)->level() <= cell->level(),
4637  ExcInternalError());
4638  first_child->set_neighbor(0, cell->neighbor(0));
4639  }
4640  else
4641  // left neighbor is refined
4642  {
4643  // set neighbor to cell on same level
4644  const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4645  first_child->set_neighbor(0,
4646  cell->neighbor(0)->child(nbnb));
4647 
4648  // reset neighbor info of all right descendant
4649  // of the left neighbor of cell
4651  left_neighbor = cell->neighbor(0);
4652  while (left_neighbor->has_children())
4653  {
4654  left_neighbor = left_neighbor->child(nbnb);
4655  left_neighbor->set_neighbor(nbnb, first_child);
4656  }
4657  }
4658 
4659  // insert second child
4660  second_child->clear_children();
4661  second_child->set_bounding_object_indices(
4662  {next_unused_vertex, cell->vertex_index(1)});
4663  second_child->set_neighbor(0, first_child);
4664  second_child->set_material_id(cell->material_id());
4665  second_child->set_manifold_id(cell->manifold_id());
4666  second_child->set_subdomain_id(subdomainid);
4667  second_child->set_direction_flag(cell->direction_flag());
4668 
4669  if (cell->neighbor(1).state() != IteratorState::valid)
4670  second_child->set_neighbor(1, cell->neighbor(1));
4671  else if (cell->neighbor(1)->is_active())
4672  {
4673  Assert(cell->neighbor(1)->level() <= cell->level(),
4674  ExcInternalError());
4675  second_child->set_neighbor(1, cell->neighbor(1));
4676  }
4677  else
4678  // right neighbor is refined same as above
4679  {
4680  const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4681  second_child->set_neighbor(
4682  1, cell->neighbor(1)->child(nbnb));
4683 
4685  right_neighbor = cell->neighbor(1);
4686  while (right_neighbor->has_children())
4687  {
4688  right_neighbor = right_neighbor->child(nbnb);
4689  right_neighbor->set_neighbor(nbnb, second_child);
4690  }
4691  }
4692  // inform all listeners that cell refinement is done
4693  triangulation.signals.post_refinement_on_cell(cell);
4694  }
4695  }
4696 
4697  // in 1d, we can not have distorted children unless the parent
4698  // was already distorted (that is because we don't use
4699  // boundary information for 1d triangulations). so return an
4700  // empty list
4702  }
4703 
4704 
4709  template <int spacedim>
4712  const bool check_for_distorted_cells)
4713  {
4714  const unsigned int dim = 2;
4715 
4716 
4717  // First check whether we can get away with isotropic refinement, or
4718  // whether we need to run through the full anisotropic algorithm
4719  {
4720  bool do_isotropic_refinement = true;
4721  for (const auto &cell : triangulation.active_cell_iterators())
4722  if (cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
4723  cell->refine_flag_set() == RefinementCase<dim>::cut_y)
4724  {
4725  do_isotropic_refinement = false;
4726  break;
4727  }
4728 
4729  if (do_isotropic_refinement)
4731  check_for_distorted_cells);
4732  }
4733 
4734  // Check whether a new level is needed. We have to check for
4735  // this on the highest level only
4736  for (const auto &cell : triangulation.active_cell_iterators_on_level(
4737  triangulation.levels.size() - 1))
4738  if (cell->refine_flag_set())
4739  {
4740  triangulation.levels.push_back(
4741  std::make_unique<
4743  break;
4744  }
4745 
4746  // TODO[WB]: we clear user flags and pointers of lines; we're going
4747  // to use them to flag which lines need refinement
4748  for (typename Triangulation<dim, spacedim>::line_iterator line =
4749  triangulation.begin_line();
4750  line != triangulation.end_line();
4751  ++line)
4752  {
4753  line->clear_user_flag();
4754  line->clear_user_data();
4755  }
4756  // running over all cells and lines count the number
4757  // n_single_lines of lines which can be stored as single
4758  // lines, e.g. inner lines
4759  unsigned int n_single_lines = 0;
4760 
4761  // New lines to be created: number lines which are stored in
4762  // pairs (the children of lines must be stored in pairs)
4763  unsigned int n_lines_in_pairs = 0;
4764 
4765  // check how much space is needed on every level. We need not
4766  // check the highest level since either - on the highest level
4767  // no cells are flagged for refinement - there are, but
4768  // prepare_refinement added another empty level
4769  unsigned int needed_vertices = 0;
4770  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4771  {
4772  // count number of flagged cells on this level and compute
4773  // how many new vertices and new lines will be needed
4774  unsigned int needed_cells = 0;
4775 
4776  for (const auto &cell :
4777  triangulation.active_cell_iterators_on_level(level))
4778  if (cell->refine_flag_set())
4779  {
4780  if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
4781  {
4782  needed_cells += 4;
4783 
4784  // new vertex at center of cell is needed in any
4785  // case
4786  ++needed_vertices;
4787 
4788  // the four inner lines can be stored as singles
4789  n_single_lines += 4;
4790  }
4791  else // cut_x || cut_y
4792  {
4793  // set the flag showing that anisotropic
4794  // refinement is used for at least one cell
4795  triangulation.anisotropic_refinement = true;
4796 
4797  needed_cells += 2;
4798  // no vertex at center
4799 
4800  // the inner line can be stored as single
4801  n_single_lines += 1;
4802  }
4803 
4804  // mark all faces (lines) for refinement; checking
4805  // locally whether the neighbor would also like to
4806  // refine them is rather difficult for lines so we
4807  // only flag them and after visiting all cells, we
4808  // decide which lines need refinement;
4809  for (const unsigned int line_no :
4811  {
4813  cell->refine_flag_set(), line_no) ==
4815  {
4817  line = cell->line(line_no);
4818  if (line->has_children() == false)
4819  line->set_user_flag();
4820  }
4821  }
4822  }
4823 
4824 
4825  // count number of used cells on the next higher level
4826  const unsigned int used_cells =
4827  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4828  triangulation.levels[level + 1]->cells.used.end(),
4829  true);
4830 
4831 
4832  // reserve space for the used_cells cells already existing
4833  // on the next higher level as well as for the
4834  // needed_cells that will be created on that level
4835  reserve_space(*triangulation.levels[level + 1],
4836  used_cells + needed_cells,
4837  2,
4838  spacedim);
4839 
4840  // reserve space for needed_cells new quads on the next
4841  // higher level
4842  reserve_space(triangulation.levels[level + 1]->cells,
4843  needed_cells,
4844  0);
4845  }
4846 
4847  // now count the lines which were flagged for refinement
4848  for (typename Triangulation<dim, spacedim>::line_iterator line =
4849  triangulation.begin_line();
4850  line != triangulation.end_line();
4851  ++line)
4852  if (line->user_flag_set())
4853  {
4854  Assert(line->has_children() == false, ExcInternalError());
4855  n_lines_in_pairs += 2;
4856  needed_vertices += 1;
4857  }
4858  // reserve space for n_lines_in_pairs new lines. note, that
4859  // we can't reserve space for the single lines here as well,
4860  // as all the space reserved for lines in pairs would be
4861  // counted as unused and we would end up with too little space
4862  // to store all lines. memory reservation for n_single_lines
4863  // can only be done AFTER we refined the lines of the current
4864  // cells
4865  reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4866 
4867  // add to needed vertices how many vertices are already in use
4868  needed_vertices += std::count(triangulation.vertices_used.begin(),
4869  triangulation.vertices_used.end(),
4870  true);
4871  // if we need more vertices: create them, if not: leave the
4872  // array as is, since shrinking is not really possible because
4873  // some of the vertices at the end may be in use
4874  if (needed_vertices > triangulation.vertices.size())
4875  {
4876  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4877  triangulation.vertices_used.resize(needed_vertices, false);
4878  }
4879 
4880 
4881  // Do REFINEMENT on every level; exclude highest level as
4882  // above
4883 
4884  // index of next unused vertex
4885  unsigned int next_unused_vertex = 0;
4886 
4887  // first the refinement of lines. children are stored
4888  // pairwise
4889  {
4890  // only active objects can be refined further
4892  line = triangulation.begin_active_line(),
4893  endl = triangulation.end_line();
4895  next_unused_line = triangulation.begin_raw_line();
4896 
4897  for (; line != endl; ++line)
4898  if (line->user_flag_set())
4899  {
4900  // this line needs to be refined
4901 
4902  // find the next unused vertex and set it
4903  // appropriately
4904  while (triangulation.vertices_used[next_unused_vertex] == true)
4905  ++next_unused_vertex;
4906  Assert(
4907  next_unused_vertex < triangulation.vertices.size(),
4908  ExcMessage(
4909  "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."));
4910  triangulation.vertices_used[next_unused_vertex] = true;
4911 
4912  triangulation.vertices[next_unused_vertex] = line->center(true);
4913 
4914  // now that we created the right point, make up the
4915  // two child lines. To this end, find a pair of
4916  // unused lines
4917  bool pair_found = false;
4918  (void)pair_found;
4919  for (; next_unused_line != endl; ++next_unused_line)
4920  if (!next_unused_line->used() &&
4921  !(++next_unused_line)->used())
4922  {
4923  // go back to the first of the two unused
4924  // lines
4925  --next_unused_line;
4926  pair_found = true;
4927  break;
4928  }
4929  Assert(pair_found, ExcInternalError());
4930 
4931  // there are now two consecutive unused lines, such
4932  // that the children of a line will be consecutive.
4933  // then set the child pointer of the present line
4934  line->set_children(0, next_unused_line->index());
4935 
4936  // set the two new lines
4938  children[2] = {next_unused_line, ++next_unused_line};
4939  // some tests; if any of the iterators should be
4940  // invalid, then already dereferencing will fail
4941  AssertIsNotUsed(children[0]);
4942  AssertIsNotUsed(children[1]);
4943 
4944  children[0]->set_bounding_object_indices(
4945  {line->vertex_index(0), next_unused_vertex});
4946  children[1]->set_bounding_object_indices(
4947  {next_unused_vertex, line->vertex_index(1)});
4948 
4949  children[0]->set_used_flag();
4950  children[1]->set_used_flag();
4951  children[0]->clear_children();
4952  children[1]->clear_children();
4953  children[0]->clear_user_data();
4954  children[1]->clear_user_data();
4955  children[0]->clear_user_flag();
4956  children[1]->clear_user_flag();
4957 
4958 
4959  children[0]->set_boundary_id_internal(line->boundary_id());
4960  children[1]->set_boundary_id_internal(line->boundary_id());
4961 
4962  children[0]->set_manifold_id(line->manifold_id());
4963  children[1]->set_manifold_id(line->manifold_id());
4964 
4965  // finally clear flag indicating the need for
4966  // refinement
4967  line->clear_user_flag();
4968  }
4969  }
4970 
4971 
4972  // Now set up the new cells
4973 
4974  // reserve space for inner lines (can be stored as single
4975  // lines)
4976  reserve_space(triangulation.faces->lines, 0, n_single_lines);
4977 
4979  cells_with_distorted_children;
4980 
4981  // reset next_unused_line, as now also single empty places in
4982  // the vector can be used
4984  next_unused_line = triangulation.begin_raw_line();
4985 
4986  for (int level = 0;
4987  level < static_cast<int>(triangulation.levels.size()) - 1;
4988  ++level)
4989  {
4991  next_unused_cell = triangulation.begin_raw(level + 1);
4992 
4993  for (const auto &cell :
4994  triangulation.active_cell_iterators_on_level(level))
4995  if (cell->refine_flag_set())
4996  {
4997  // actually set up the children and update neighbor
4998  // information
5000  next_unused_vertex,
5001  next_unused_line,
5002  next_unused_cell,
5003  cell);
5004 
5005  if (check_for_distorted_cells &&
5006  has_distorted_children<dim, spacedim>(cell))
5007  cells_with_distorted_children.distorted_cells.push_back(
5008  cell);
5009  // inform all listeners that cell refinement is done
5010  triangulation.signals.post_refinement_on_cell(cell);
5011  }
5012  }
5013 
5014  return cells_with_distorted_children;
5015  }
5016 
5017 
5018  template <int spacedim>
5021  const bool check_for_distorted_cells)
5022  {
5023  static const int dim = 3;
5024  static const unsigned int X = numbers::invalid_unsigned_int;
5025  using raw_line_iterator =
5027  using raw_quad_iterator =
5029 
5030  Assert(spacedim == 3, ExcNotImplemented());
5031 
5032  Assert(triangulation.vertices.size() ==
5033  triangulation.vertices_used.size(),
5034  ExcInternalError());
5035 
5036  // Check whether a new level is needed. We have to check for
5037  // this on the highest level only
5038  for (const auto &cell : triangulation.active_cell_iterators_on_level(
5039  triangulation.levels.size() - 1))
5040  if (cell->refine_flag_set())
5041  {
5042  triangulation.levels.push_back(
5043  std::make_unique<
5045  break;
5046  }
5047 
5048  // first clear user flags for quads and lines; we're going to
5049  // use them to flag which lines and quads need refinement
5050  triangulation.faces->quads.clear_user_data();
5051  triangulation.faces->lines.clear_user_flags();
5052  triangulation.faces->quads.clear_user_flags();
5053 
5054  // check how much space is needed on every level. We need not
5055  // check the highest level since either
5056  // - on the highest level no cells are flagged for refinement
5057  // - there are, but prepare_refinement added another empty
5058  // level which then is the highest level
5059 
5060  // variables to hold the number of newly to be created
5061  // vertices, lines and quads. as these are stored globally,
5062  // declare them outside the loop over al levels. we need lines
5063  // and quads in pairs for refinement of old ones and lines and
5064  // quads, that can be stored as single ones, as they are newly
5065  // created in the inside of an existing cell
5066  unsigned int needed_vertices = 0;
5067  unsigned int needed_lines_single = 0;
5068  unsigned int needed_quads_single = 0;
5069  unsigned int needed_lines_pair = 0;
5070  unsigned int needed_quads_pair = 0;
5071  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5072  {
5073  unsigned int new_cells = 0;
5074 
5075  for (const auto &cell :
5076  triangulation.active_cell_iterators_on_level(level))
5077  if (cell->refine_flag_set())
5078  {
5079  // Only support isotropic refinement
5080  Assert(cell->refine_flag_set() ==
5082  ExcInternalError());
5083 
5084  // Now count up how many new cells, faces, edges, and vertices
5085  // we will need to allocate to do this refinement.
5086  new_cells += cell->reference_cell().n_isotropic_children();
5087 
5088  if (cell->reference_cell() == ReferenceCells::Hexahedron)
5089  {
5090  ++needed_vertices;
5091  needed_lines_single += 6;
5092  needed_quads_single += 12;
5093  }
5094  else if (cell->reference_cell() ==
5096  {
5097  needed_lines_single += 1;
5098  needed_quads_single += 8;
5099  }
5100  else
5101  {
5102  Assert(false, ExcInternalError());
5103  }
5104 
5105  // Also check whether we have to refine any of the faces and
5106  // edges that bound this cell. They may of course already be
5107  // refined, so we only *mark* them for refinement by setting
5108  // the user flags
5109  for (const auto face : cell->face_indices())
5110  if (cell->face(face)->n_children() == 0)
5111  cell->face(face)->set_user_flag();
5112  else
5113  Assert(cell->face(face)->n_children() ==
5114  cell->reference_cell()
5115  .face_reference_cell(face)
5116  .n_isotropic_children(),
5117  ExcInternalError());
5118 
5119  for (const auto line : cell->line_indices())
5120  if (cell->line(line)->has_children() == false)
5121  cell->line(line)->set_user_flag();
5122  else
5123  Assert(cell->line(line)->n_children() == 2,
5124  ExcInternalError());
5125  }
5126 
5127  const unsigned int used_cells =
5128  std::count(triangulation.levels[level + 1]->cells.used.begin(),
5129  triangulation.levels[level + 1]->cells.used.end(),
5130  true);
5131 
5132  reserve_space(*triangulation.levels[level + 1],
5133  used_cells + new_cells,
5134  3,
5135  spacedim);
5136 
5137  reserve_space(triangulation.levels[level + 1]->cells, new_cells);
5138  }
5139 
5140  // now count the quads and lines which were flagged for
5141  // refinement
5142  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
5143  triangulation.begin_quad();
5144  quad != triangulation.end_quad();
5145  ++quad)
5146  {
5147  if (quad->user_flag_set() == false)
5148  continue;
5149 
5150  if (quad->reference_cell() == ReferenceCells::Quadrilateral)
5151  {
5152  needed_quads_pair += 4;
5153  needed_lines_pair += 4;
5154  needed_vertices += 1;
5155  }
5156  else if (quad->reference_cell() == ReferenceCells::Triangle)
5157  {
5158  needed_quads_pair += 4;
5159  needed_lines_single += 3;
5160  }
5161  else
5162  {
5163  Assert(false, ExcInternalError());
5164  }
5165  }
5166 
5167  for (typename Triangulation<dim, spacedim>::line_iterator line =
5168  triangulation.begin_line();
5169  line != triangulation.end_line();
5170  ++line)
5171  {
5172  if (line->user_flag_set() == false)
5173  continue;
5174 
5175  needed_lines_pair += 2;
5176  needed_vertices += 1;
5177  }
5178 
5179  reserve_space(triangulation.faces->lines,
5180  needed_lines_pair,
5181  needed_lines_single);
5183  needed_quads_pair,
5184  needed_quads_single);
5185  reserve_space(triangulation.faces->quads,
5186  needed_quads_pair,
5187  needed_quads_single);
5188 
5189 
5190  // add to needed vertices how many vertices are already in use
5191  needed_vertices += std::count(triangulation.vertices_used.begin(),
5192  triangulation.vertices_used.end(),
5193  true);
5194 
5195  if (needed_vertices > triangulation.vertices.size())
5196  {
5197  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5198  triangulation.vertices_used.resize(needed_vertices, false);
5199  }
5200 
5201  //-----------------------------------------
5202  // Before we start with the actual refinement, we do some
5203  // sanity checks if in debug mode. especially, we try to catch
5204  // the notorious problem with lines being twice refined,
5205  // i.e. there are cells adjacent at one line ("around the
5206  // edge", but not at a face), with two cells differing by more
5207  // than one refinement level
5208  //
5209  // this check is very simple to implement here, since we have
5210  // all lines flagged if they shall be refined
5211 #ifdef DEBUG
5212  for (const auto &cell : triangulation.active_cell_iterators())
5213  if (!cell->refine_flag_set())
5214  for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
5215  if (cell->line(line_n)->has_children())
5216  for (unsigned int c = 0; c < 2; ++c)
5217  Assert(cell->line(line_n)->child(c)->user_flag_set() == false,
5218  ExcInternalError());
5219 #endif
5220 
5221  unsigned int current_vertex = 0;
5222 
5223  // helper function - find the next available vertex number and mark it
5224  // as used.
5225  auto get_next_unused_vertex = [](const unsigned int current_vertex,
5226  std::vector<bool> &vertices_used) {
5227  unsigned int next_vertex = current_vertex;
5228  while (next_vertex < vertices_used.size() &&
5229  vertices_used[next_vertex] == true)
5230  ++next_vertex;
5231  Assert(next_vertex < vertices_used.size(), ExcInternalError());
5232  vertices_used[next_vertex] = true;
5233 
5234  return next_vertex;
5235  };
5236 
5237  // LINES
5238  {
5240  line = triangulation.begin_active_line(),
5241  endl = triangulation.end_line();
5242  raw_line_iterator next_unused_line = triangulation.begin_raw_line();
5243 
5244  for (; line != endl; ++line)
5245  {
5246  if (line->user_flag_set() == false)
5247  continue;
5248 
5249  next_unused_line =
5250  triangulation.faces->lines.template next_free_pair_object<1>(
5251  triangulation);
5252  Assert(next_unused_line.state() == IteratorState::valid,
5253  ExcInternalError());
5254 
5255  // now we found two consecutive unused lines, such
5256  // that the children of a line will be consecutive.
5257  // then set the child pointer of the present line
5258  line->set_children(0, next_unused_line->index());
5259 
5260  const std::array<raw_line_iterator, 2> children{
5261  {next_unused_line, ++next_unused_line}};
5262 
5263  AssertIsNotUsed(children[0]);
5264  AssertIsNotUsed(children[1]);
5265 
5266  current_vertex =
5267  get_next_unused_vertex(current_vertex,
5268  triangulation.vertices_used);
5269  triangulation.vertices[current_vertex] = line->center(true);
5270 
5271  children[0]->set_bounding_object_indices(
5272  {line->vertex_index(0), current_vertex});
5273  children[1]->set_bounding_object_indices(
5274  {current_vertex, line->vertex_index(1)});
5275 
5276  const auto manifold_id = line->manifold_id();
5277  const auto boundary_id = line->boundary_id();
5278  for (const auto &child : children)
5279  {
5280  child->set_used_flag();
5281  child->clear_children();
5282  child->clear_user_data();
5283  child->clear_user_flag();
5284  child->set_boundary_id_internal(boundary_id);
5285  child->set_manifold_id(manifold_id);
5286  }
5287 
5288  line->clear_user_flag();
5289  }
5290  }
5291 
5292  // QUADS
5293  {
5295  quad = triangulation.begin_quad(),
5296  endq = triangulation.end_quad();
5297 
5298  for (; quad != endq; ++quad)
5299  {
5300  if (quad->user_flag_set() == false)
5301  continue;
5302 
5303  const auto reference_face_type = quad->reference_cell();
5304 
5305  // 1) create new lines (property is set later)
5306  // maximum of 4 new lines (4 quadrilateral, 3 triangle)
5307  std::array<raw_line_iterator, 4> new_lines;
5308  if (reference_face_type == ReferenceCells::Quadrilateral)
5309  {
5310  for (unsigned int l = 0; l < 2; ++l)
5311  {
5312  auto next_unused_line =
5313  triangulation.faces->lines
5314  .template next_free_pair_object<1>(triangulation);
5315  new_lines[2 * l] = next_unused_line;
5316  new_lines[2 * l + 1] = ++next_unused_line;
5317  }
5318  }
5319  else if (reference_face_type == ReferenceCells::Triangle)
5320  {
5321  for (unsigned int l = 0; l < 3; ++l)
5322  new_lines[l] =
5323  triangulation.faces->lines
5324  .template next_free_single_object<1>(triangulation);
5325  }
5326  else
5327  {
5328  Assert(false, ExcNotImplemented());
5329  }
5330 
5331  for (const unsigned int line : quad->line_indices())
5332  {
5333  AssertIsNotUsed(new_lines[line]);
5334  (void)line;
5335  }
5336 
5337  // 2) create new quads (properties are set below). Both triangles
5338  // and quads are divided in four.
5339  std::array<raw_quad_iterator, 4> new_quads;
5340  for (unsigned int q = 0; q < 2; ++q)
5341  {
5342  auto next_unused_quad =
5343  triangulation.faces->quads
5344  .template next_free_pair_object<2>(triangulation);
5345 
5346  new_quads[2 * q] = next_unused_quad;
5347  new_quads[2 * q + 1] = ++next_unused_quad;
5348 
5349  quad->set_children(2 * q, new_quads[2 * q]->index());
5350  }
5351  quad->set_refinement_case(RefinementCase<2>::cut_xy);
5352 
5353  for (const auto &quad : new_quads)
5354  {
5355  AssertIsNotUsed(quad);
5356  (void)quad;
5357  }
5358 
5359  // 3) set vertex indices and set new vertex
5360 
5361  // Maximum of 9 vertices per refined quad (9 for Quadrilateral, 6
5362  // for Triangle)
5363  std::array<unsigned int, 9> vertex_indices = {};
5364  unsigned int k = 0;
5365  for (const auto i : quad->vertex_indices())
5366  vertex_indices[k++] = quad->vertex_index(i);
5367 
5368  for (const auto i : quad->line_indices())
5369  vertex_indices[k++] = quad->line(i)->child(0)->vertex_index(1);
5370 
5371  if (reference_face_type == ReferenceCells::Quadrilateral)
5372  {
5373  current_vertex =
5374  get_next_unused_vertex(current_vertex,
5375  triangulation.vertices_used);
5376  vertex_indices[k++] = current_vertex;
5377 
5378  triangulation.vertices[current_vertex] =
5379  quad->center(true, true);
5380  }
5381 
5382  // 4) set new lines on quads and their properties
5383  std::array<raw_line_iterator, 12> lines;
5384  unsigned int n_lines = 0;
5385  for (unsigned int l = 0; l < quad->n_lines(); ++l)
5386  for (unsigned int c = 0; c < 2; ++c)
5387  {
5388  static constexpr ::ndarray<unsigned int, 2, 2> index =
5389  {{// child 0, line_orientation=false and true
5390  {{1, 0}},
5391  // child 1, line_orientation=false and true
5392  {{0, 1}}}};
5393 
5394  lines[n_lines++] =
5395  quad->line(l)->child(index[c][quad->line_orientation(l)]);
5396  }
5397 
5398  for (unsigned int l = 0; l < quad->n_lines(); ++l)
5399  lines[n_lines++] = new_lines[l];
5400 
5401  std::array<int, 12> line_indices;
5402  for (unsigned int i = 0; i < n_lines; ++i)
5403  line_indices[i] = lines[i]->index();
5404 
5405  static constexpr ::ndarray<unsigned int, 12, 2>
5406  line_vertices_quad{{{{0, 4}},
5407  {{4, 2}},
5408  {{1, 5}},
5409  {{5, 3}},
5410  {{0, 6}},
5411  {{6, 1}},
5412  {{2, 7}},
5413  {{7, 3}},
5414  {{6, 8}},
5415  {{8, 7}},
5416  {{4, 8}},
5417  {{8, 5}}}};
5418 
5419  static constexpr ::ndarray<unsigned int, 4, 4>
5420  quad_lines_quad{{{{0, 8, 4, 10}},
5421  {{8, 2, 5, 11}},
5422  {{1, 9, 10, 6}},
5423  {{9, 3, 11, 7}}}};
5424 
5425  static constexpr ::ndarray<unsigned int, 12, 2>
5426  line_vertices_tri{{{{0, 3}},
5427  {{3, 1}},
5428  {{1, 4}},
5429  {{4, 2}},
5430  {{2, 5}},
5431  {{5, 0}},
5432  {{3, 4}},
5433  {{4, 5}},
5434  {{3, 5}},
5435  {{X, X}},
5436  {{X, X}},
5437  {{X, X}}}};
5438 
5439  static constexpr ::ndarray<unsigned int, 4, 4>
5440  quad_lines_tri{{{{0, 8, 5, X}},
5441  {{1, 2, 6, X}},
5442  {{7, 3, 4, X}},
5443  {{6, 7, 8, X}}}};
5444 
5445  static constexpr ::ndarray<unsigned int, 4, 4, 2>
5446  quad_line_vertices_tri{
5447  {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
5448  {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
5449  {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
5450  {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
5451 
5452  const auto &line_vertices =
5453  (reference_face_type == ReferenceCells::Quadrilateral) ?
5454  line_vertices_quad :
5455  line_vertices_tri;
5456  const auto &quad_lines =
5457  (reference_face_type == ReferenceCells::Quadrilateral) ?
5458  quad_lines_quad :
5459  quad_lines_tri;
5460 
5461  for (unsigned int i = 0, j = 2 * quad->n_lines();
5462  i < quad->n_lines();
5463  ++i, ++j)
5464  {
5465  auto &new_line = new_lines[i];
5466  new_line->set_bounding_object_indices(
5467  {vertex_indices[line_vertices[j][0]],
5468  vertex_indices[line_vertices[j][1]]});
5469  new_line->set_used_flag();
5470  new_line->clear_user_flag();
5471  new_line->clear_user_data();
5472  new_line->clear_children();
5473  new_line->set_boundary_id_internal(quad->boundary_id());
5474  new_line->set_manifold_id(quad->manifold_id());
5475  }
5476 
5477  // 5) set properties of quads
5478  for (unsigned int i = 0; i < new_quads.size(); ++i)
5479  {
5480  auto &new_quad = new_quads[i];
5481 
5482  // TODO: we assume here that all children have the same type
5483  // as the parent
5484  triangulation.faces->set_quad_type(new_quad->index(),
5485  reference_face_type);
5486 
5487  if (reference_face_type == ReferenceCells::Triangle)
5488  new_quad->set_bounding_object_indices(
5489  {line_indices[quad_lines[i][0]],
5490  line_indices[quad_lines[i][1]],
5491  line_indices[quad_lines[i][2]]});
5492  else if (reference_face_type == ReferenceCells::Quadrilateral)
5493  new_quad->set_bounding_object_indices(
5494  {line_indices[quad_lines[i][0]],
5495  line_indices[quad_lines[i][1]],
5496  line_indices[quad_lines[i][2]],
5497  line_indices[quad_lines[i][3]]});
5498  else
5499  Assert(false, ExcNotImplemented());
5500 
5501  new_quad->set_used_flag();
5502  new_quad->clear_user_flag();
5503  new_quad->clear_user_data();
5504  new_quad->clear_children();
5505  new_quad->set_boundary_id_internal(quad->boundary_id());
5506  new_quad->set_manifold_id(quad->manifold_id());
5507 
5508 #ifdef DEBUG
5509  std::set<unsigned int> s;
5510 #endif
5511 
5512  // ... and fix orientation of lines of face for triangles,
5513  // using an expensive algorithm, quadrilaterals are treated
5514  // a few lines below by a cheaper algorithm
5515  if (reference_face_type == ReferenceCells::Triangle)
5516  {
5517  for (const auto f : new_quad->line_indices())
5518  {
5519  const std::array<unsigned int, 2> vertices_0 = {
5520  {lines[quad_lines[i][f]]->vertex_index(0),
5521  lines[quad_lines[i][f]]->vertex_index(1)}};
5522 
5523  const std::array<unsigned int, 2> vertices_1 = {
5524  {vertex_indices[quad_line_vertices_tri[i][f][0]],
5525  vertex_indices[quad_line_vertices_tri[i][f][1]]}};
5526 
5527  const auto orientation =
5529  make_array_view(vertices_0),
5530  make_array_view(vertices_1));
5531 
5532 #ifdef DEBUG
5533  for (const auto i : vertices_0)
5534  s.insert(i);
5535  for (const auto i : vertices_1)
5536  s.insert(i);
5537 #endif
5538 
5539  new_quad->set_line_orientation(f, orientation);
5540  }
5541 #ifdef DEBUG
5542  AssertDimension(s.size(), 3);
5543 #endif
5544  }
5545  }
5546 
5547  // fix orientation of lines of faces for quadrilaterals with
5548  // cheap algorithm
5549  if (reference_face_type == ReferenceCells::Quadrilateral)
5550  {
5551  static constexpr ::ndarray<unsigned int, 4, 2>
5552  quad_child_boundary_lines{
5553  {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
5554 
5555  for (unsigned int i = 0; i < 4; ++i)
5556  for (unsigned int j = 0; j < 2; ++j)
5557  new_quads[quad_child_boundary_lines[i][j]]
5558  ->set_line_orientation(i, quad->line_orientation(i));
5559  }
5560 
5561  quad->clear_user_flag();
5562  }
5563  }
5564 
5566  cells_with_distorted_children;
5567 
5569  triangulation.begin_active_hex(0);
5570  for (unsigned int level = 0; level != triangulation.levels.size() - 1;
5571  ++level)
5572  {
5574  next_unused_hex = triangulation.begin_raw_hex(level + 1);
5575  Assert(hex == triangulation.end() ||
5576  hex->level() >= static_cast<int>(level),
5577  ExcInternalError());
5578 
5579  for (; hex != triangulation.end() &&
5580  hex->level() == static_cast<int>(level);
5581  ++hex)
5582  {
5583  if (hex->refine_flag_set() ==
5585  continue;
5586 
5587  const auto &reference_cell_type = hex->reference_cell();
5588 
5589  const RefinementCase<dim> ref_case = hex->refine_flag_set();
5590  hex->clear_refine_flag();
5591  hex->set_refinement_case(ref_case);
5592 
5593  unsigned int n_new_lines = 0;
5594  unsigned int n_new_quads = 0;
5595  unsigned int n_new_hexes = 0;
5596 
5597  if (reference_cell_type == ReferenceCells::Hexahedron)
5598  {
5599  n_new_lines = 6;
5600  n_new_quads = 12;
5601  n_new_hexes = 8;
5602  }
5603  else if (reference_cell_type == ReferenceCells::Tetrahedron)
5604  {
5605  n_new_lines = 1;
5606  n_new_quads = 8;
5607  n_new_hexes = 8;
5608  }
5609  else
5610  Assert(false, ExcNotImplemented());
5611 
5612  std::array<raw_line_iterator, 6> new_lines;
5613  for (unsigned int i = 0; i < n_new_lines; ++i)
5614  {
5615  new_lines[i] =
5616  triangulation.faces->lines
5617  .template next_free_single_object<1>(triangulation);
5618 
5619  AssertIsNotUsed(new_lines[i]);
5620  new_lines[i]->set_used_flag();
5621  new_lines[i]->clear_user_flag();
5622  new_lines[i]->clear_user_data();
5623  new_lines[i]->clear_children();
5624  new_lines[i]->set_boundary_id_internal(
5626  new_lines[i]->set_manifold_id(hex->manifold_id());
5627  }
5628 
5629  std::array<raw_quad_iterator, 12> new_quads;
5630  for (unsigned int i = 0; i < n_new_quads; ++i)
5631  {
5632  new_quads[i] =
5633  triangulation.faces->quads
5634  .template next_free_single_object<2>(triangulation);
5635 
5636  auto &new_quad = new_quads[i];
5637 
5638  // TODO: faces of children have the same type as the faces
5639  // of the parent
5640  triangulation.faces->set_quad_type(
5641  new_quad->index(),
5642  reference_cell_type.face_reference_cell(0));
5643 
5644  AssertIsNotUsed(new_quad);
5645  new_quad->set_used_flag();
5646  new_quad->clear_user_flag();
5647  new_quad->clear_user_data();
5648  new_quad->clear_children();
5649  new_quad->set_boundary_id_internal(
5651  new_quad->set_manifold_id(hex->manifold_id());
5652  for (const auto j : new_quads[i]->line_indices())
5653  new_quad->set_line_orientation(j, true);
5654  }
5655 
5656  // we always get 8 children per refined cell
5657  std::array<
5659  8>
5660  new_hexes;
5661  {
5662  for (unsigned int i = 0; i < n_new_hexes; ++i)
5663  {
5664  if (i % 2 == 0)
5665  next_unused_hex =
5666  triangulation.levels[level + 1]->cells.next_free_hex(
5667  triangulation, level + 1);
5668  else
5669  ++next_unused_hex;
5670 
5671  new_hexes[i] = next_unused_hex;
5672 
5673  auto &new_hex = new_hexes[i];
5674 
5675  // children have the same type as the parent
5676  triangulation.levels[new_hex->level()]
5677  ->reference_cell[new_hex->index()] =
5678  reference_cell_type;
5679 
5680  AssertIsNotUsed(new_hex);
5681  new_hex->set_used_flag();
5682  new_hex->clear_user_flag();
5683  new_hex->clear_user_data();
5684  new_hex->clear_children();
5685  new_hex->set_material_id(hex->material_id());
5686  new_hex->set_manifold_id(hex->manifold_id());
5687  new_hex->set_subdomain_id(hex->subdomain_id());
5688 
5689  if (i % 2)
5690  new_hex->set_parent(hex->index());
5691 
5692  // set the orientation flag to its default state for all
5693  // faces initially. later on go the other way round and
5694  // reset faces that are at the boundary of the mother cube
5695  for (const auto f : new_hex->face_indices())
5696  new_hex->set_combined_face_orientation(
5697  f,
5699  }
5700  for (unsigned int i = 0; i < n_new_hexes / 2; ++i)
5701  hex->set_children(2 * i, new_hexes[2 * i]->index());
5702  }
5703 
5704  {
5705  // load vertex indices
5706  std::array<unsigned int, 27> vertex_indices = {};
5707 
5708  {
5709  unsigned int k = 0;
5710 
5711  for (const unsigned int i : hex->vertex_indices())
5712  vertex_indices[k++] = hex->vertex_index(i);
5713 
5714  const std::array<unsigned int, 12> line_indices =
5715  TriaAccessorImplementation::Implementation::
5716  get_line_indices_of_cell(*hex);
5717  // avoid a compiler warning by fixing the max number of
5718  // loop iterations to 12
5719  const unsigned int n_lines = std::min(hex->n_lines(), 12u);
5720  for (unsigned int l = 0; l < n_lines; ++l)
5721  {
5722  raw_line_iterator line(&triangulation,
5723  0,
5724  line_indices[l]);
5725  vertex_indices[k++] = line->child(0)->vertex_index(1);
5726  }
5727 
5728  if (reference_cell_type == ReferenceCells::Hexahedron)
5729  {
5730  for (const unsigned int i : hex->face_indices())
5731  vertex_indices[k++] =
5732  hex->face(i)->child(0)->vertex_index(3);
5733 
5734  // Set single new vertex in the center
5735  current_vertex =
5736  get_next_unused_vertex(current_vertex,
5737  triangulation.vertices_used);
5738  vertex_indices[k++] = current_vertex;
5739 
5740  triangulation.vertices[current_vertex] =
5741  hex->center(true, true);
5742  }
5743  }
5744 
5745  // set up new lines
5746  {
5747  static constexpr ::ndarray<unsigned int, 6, 2>
5748  new_line_vertices_hex = {{{{22, 26}},
5749  {{26, 23}},
5750  {{20, 26}},
5751  {{26, 21}},
5752  {{24, 26}},
5753  {{26, 25}}}};
5754 
5755  static constexpr ::ndarray<unsigned int, 6, 2>
5756  new_line_vertices_tet = {{{{6, 8}},
5757  {{X, X}},
5758  {{X, X}},
5759  {{X, X}},
5760  {{X, X}},
5761  {{X, X}}}};
5762 
5763  const auto &new_line_vertices =
5764  (reference_cell_type == ReferenceCells::Hexahedron) ?
5765  new_line_vertices_hex :
5766  new_line_vertices_tet;
5767 
5768  for (unsigned int i = 0; i < n_new_lines; ++i)
5769  new_lines[i]->set_bounding_object_indices(
5770  {vertex_indices[new_line_vertices[i][0]],
5771  vertex_indices[new_line_vertices[i][1]]});
5772  }
5773 
5774  // set up new quads
5775  {
5776  boost::container::small_vector<raw_line_iterator, 30>
5777  relevant_lines;
5778 
5779  if (reference_cell_type == ReferenceCells::Hexahedron)
5780  {
5781  relevant_lines.resize(30);
5782  for (unsigned int f = 0, k = 0; f < 6; ++f)
5783  for (unsigned int c = 0; c < 4; ++c, ++k)
5784  {
5785  static constexpr ::
5786  ndarray<unsigned int, 4, 2>
5787  temp = {
5788  {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}};
5789 
5790  relevant_lines[k] =
5791  hex->face(f)
5792  ->isotropic_child(
5794  standard_to_real_face_vertex(
5795  temp[c][0],
5796  hex->face_orientation(f),
5797  hex->face_flip(f),
5798  hex->face_rotation(f)))
5799  ->line(GeometryInfo<dim>::
5800  standard_to_real_face_line(
5801  temp[c][1],
5802  hex->face_orientation(f),
5803  hex->face_flip(f),
5804  hex->face_rotation(f)));
5805  }
5806 
5807  for (unsigned int i = 0, k = 24; i < 6; ++i, ++k)
5808  relevant_lines[k] = new_lines[i];
5809  }
5810  else if (reference_cell_type == ReferenceCells::Tetrahedron)
5811  {
5812  relevant_lines.resize(13);
5813 
5814  unsigned int k = 0;
5815  for (unsigned int f = 0; f < 4; ++f)
5816  for (unsigned int l = 0; l < 3; ++l, ++k)
5817  {
5818  // TODO: add comment
5819  static const std::
5820  array<std::array<unsigned int, 3>, 6>
5821  table = {{{{1, 0, 2}}, // 0
5822  {{0, 1, 2}},
5823  {{0, 2, 1}}, // 2
5824  {{1, 2, 0}},
5825  {{2, 1, 0}}, // 4
5826  {{2, 0, 1}}}};
5827 
5828  relevant_lines[k] =
5829  hex->face(f)
5830  ->child(3 /*center triangle*/)
5831  ->line(
5832  table[triangulation.levels[hex->level()]
5833  ->face_orientations
5834  .get_combined_orientation(
5835  hex->index() * GeometryInfo<dim>::
5836  faces_per_cell +
5837  f)][l]);
5838  }
5839 
5840  relevant_lines[k++] = new_lines[0];
5841 
5842  AssertDimension(k, 13);
5843  }
5844  else
5845  Assert(false, ExcNotImplemented());
5846 
5847  boost::container::small_vector<unsigned int, 30>
5848  relevant_line_indices(relevant_lines.size());
5849  for (unsigned int i = 0; i < relevant_line_indices.size();
5850  ++i)
5851  relevant_line_indices[i] = relevant_lines[i]->index();
5852 
5853  static constexpr ::ndarray<unsigned int, 12, 4>
5854  new_quad_lines_hex = {{{{10, 28, 16, 24}},
5855  {{28, 14, 17, 25}},
5856  {{11, 29, 24, 20}},
5857  {{29, 15, 25, 21}},
5858  {{18, 26, 0, 28}},
5859  {{26, 22, 1, 29}},
5860  {{19, 27, 28, 4}},
5861  {{27, 23, 29, 5}},
5862  {{2, 24, 8, 26}},
5863  {{24, 6, 9, 27}},
5864  {{3, 25, 26, 12}},
5865  {{25, 7, 27, 13}}}};
5866 
5867  static constexpr ::ndarray<unsigned int, 12, 4>
5868  new_quad_lines_tet = {{{{2, 3, 8, X}},
5869  {{0, 9, 5, X}},
5870  {{1, 6, 11, X}},
5871  {{4, 10, 7, X}},
5872  {{2, 12, 5, X}},
5873  {{1, 9, 12, X}},
5874  {{4, 8, 12, X}},
5875  {{6, 12, 10, X}},
5876  {{X, X, X, X}},
5877  {{X, X, X, X}},
5878  {{X, X, X, X}},
5879  {{X, X, X, X}}}};
5880 
5881  static constexpr ::ndarray<unsigned int, 12, 4, 2>
5882  table_hex = {
5883  {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
5884  {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
5885  {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
5886  {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
5887  {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
5888  {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
5889  {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
5890  {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
5891  {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
5892  {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
5893  {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
5894  {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
5895 
5896  static constexpr ::ndarray<unsigned int, 12, 4, 2>
5897  table_tet = {
5898  {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
5899  {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
5900  {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
5901  {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
5902  {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
5903  {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
5904  {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
5905  {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
5906  {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5907  {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5908  {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5909  {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
5910 
5911  const auto &new_quad_lines =
5912  (reference_cell_type == ReferenceCells::Hexahedron) ?
5913  new_quad_lines_hex :
5914  new_quad_lines_tet;
5915 
5916  const auto &table =
5917  (reference_cell_type == ReferenceCells::Hexahedron) ?
5918  table_hex :
5919  table_tet;
5920 
5921  static constexpr ::ndarray<unsigned int, 4, 2>
5922  representative_lines{
5923  {{{0, 2}}, {{2, 0}}, {{3, 3}}, {{1, 1}}}};
5924 
5925  for (unsigned int q = 0; q < n_new_quads; ++q)
5926  {
5927  auto &new_quad = new_quads[q];
5928 
5929  if (new_quad->n_lines() == 3)
5930  new_quad->set_bounding_object_indices(
5931  {relevant_line_indices[new_quad_lines[q][0]],
5932  relevant_line_indices[new_quad_lines[q][1]],
5933  relevant_line_indices[new_quad_lines[q][2]]});
5934  else if (new_quad->n_lines() == 4)
5935  new_quad->set_bounding_object_indices(
5936  {relevant_line_indices[new_quad_lines[q][0]],
5937  relevant_line_indices[new_quad_lines[q][1]],
5938  relevant_line_indices[new_quad_lines[q][2]],
5939  relevant_line_indices[new_quad_lines[q][3]]});
5940  else
5941  Assert(false, ExcNotImplemented());
5942 
5943  // On hexes, we must only determine a single line
5944  // according to the representative_lines array above
5945  // (this saves expensive operations), for tets we do
5946  // all lines manually
5947  const unsigned int n_compute_lines =
5948  reference_cell_type == ReferenceCells::Hexahedron ?
5949  1 :
5950  new_quad->n_lines();
5951  for (unsigned int line = 0; line < n_compute_lines;
5952  ++line)
5953  {
5954  const unsigned int l =
5955  (reference_cell_type ==
5957  representative_lines[q % 4][0] :
5958  line;
5959 
5960  const std::array<unsigned int, 2> vertices_0 = {
5961  {relevant_lines[new_quad_lines[q][l]]
5962  ->vertex_index(0),
5963  relevant_lines[new_quad_lines[q][l]]
5964  ->vertex_index(1)}};
5965 
5966  const std::array<unsigned int, 2> vertices_1 = {
5967  {vertex_indices[table[q][l][0]],
5968  vertex_indices[table[q][l][1]]}};
5969 
5970  const auto orientation =
5972  make_array_view(vertices_0),
5973  make_array_view(vertices_1));
5974 
5975  new_quad->set_line_orientation(l, orientation);
5976 
5977  // on a hex, inject the status of the current line
5978  // also to the line on the other quad along the
5979  // same direction
5980  if (reference_cell_type ==
5982  new_quads[representative_lines[q % 4][1] + q -
5983  (q % 4)]
5984  ->set_line_orientation(l, orientation);
5985  }
5986  }
5987  }
5988 
5989  // set up new hex
5990  {
5991  std::array<int, 36> quad_indices;
5992 
5993  if (reference_cell_type == ReferenceCells::Hexahedron)
5994  {
5995  for (unsigned int i = 0; i < n_new_quads; ++i)
5996  quad_indices[i] = new_quads[i]->index();
5997 
5998  for (unsigned int f = 0, k = n_new_quads; f < 6; ++f)
5999  for (unsigned int c = 0; c < 4; ++c, ++k)
6000  quad_indices[k] =
6001  hex->face(f)->isotropic_child_index(
6003  c,
6004  hex->face_orientation(f),
6005  hex->face_flip(f),
6006  hex->face_rotation(f)));
6007  }
6008  else if (reference_cell_type == ReferenceCells::Tetrahedron)
6009  {
6010  for (unsigned int i = 0; i < n_new_quads; ++i)
6011  quad_indices[i] = new_quads[i]->index();
6012 
6013  for (unsigned int f = 0, k = n_new_quads; f < 4; ++f)
6014  for (unsigned int c = 0; c < 4; ++c, ++k)
6015  {
6016  quad_indices[k] = hex->face(f)->child_index(
6017  (c == 3) ?
6018  3 :
6019  reference_cell_type
6020  .standard_to_real_face_vertex(
6021  c,
6022  f,
6023  triangulation.levels[hex->level()]
6024  ->face_orientations
6025  .get_combined_orientation(
6026  hex->index() *
6028  f)));
6029  }
6030  }
6031  else
6032  {
6033  Assert(false, ExcNotImplemented());
6034  }
6035 
6036  static constexpr ::ndarray<unsigned int, 8, 6>
6037  cell_quads_hex = {{
6038  {{12, 0, 20, 4, 28, 8}}, // bottom children
6039  {{0, 16, 22, 6, 29, 9}}, //
6040  {{13, 1, 4, 24, 30, 10}}, //
6041  {{1, 17, 6, 26, 31, 11}}, //
6042  {{14, 2, 21, 5, 8, 32}}, // top children
6043  {{2, 18, 23, 7, 9, 33}}, //
6044  {{15, 3, 5, 25, 10, 34}}, //
6045  {{3, 19, 7, 27, 11, 35}} //
6046  }};
6047 
6048  static constexpr ::ndarray<unsigned int, 8, 6>
6049  cell_quads_tet{{{{8, 13, 16, 0, X, X}},
6050  {{9, 12, 1, 21, X, X}},
6051  {{10, 2, 17, 20, X, X}},
6052  {{3, 14, 18, 22, X, X}},
6053  {{11, 1, 4, 5, X, X}},
6054  {{15, 0, 4, 6, X, X}},
6055  {{19, 7, 6, 3, X, X}},
6056  {{23, 5, 2, 7, X, X}}}};
6057 
6058  static constexpr ::ndarray<unsigned int, 8, 6, 4>
6059  cell_face_vertices_tet{{{{{{0, 4, 6, X}},
6060  {{4, 0, 7, X}},
6061  {{0, 6, 7, X}},
6062  {{6, 4, 7, X}},
6063  {{X, X, X, X}},
6064  {{X, X, X, X}}}},
6065  {{{{4, 1, 5, X}},
6066  {{1, 4, 8, X}},
6067  {{4, 5, 8, X}},
6068  {{5, 1, 8, X}},
6069  {{X, X, X, X}},
6070  {{X, X, X, X}}}},
6071  {{{{6, 5, 2, X}},
6072  {{5, 6, 9, X}},
6073  {{6, 2, 9, X}},
6074  {{2, 5, 9, X}},
6075  {{X, X, X, X}},
6076  {{X, X, X, X}}}},
6077  {{{{7, 8, 9, X}},
6078  {{8, 7, 3, X}},
6079  {{7, 9, 3, X}},
6080  {{9, 8, 3, X}},
6081  {{X, X, X, X}},
6082  {{X, X, X, X}}}},
6083  {{{{4, 5, 6, X}},
6084  {{5, 4, 8, X}},
6085  {{4, 6, 8, X}},
6086  {{6, 5, 8, X}},
6087  {{X, X, X, X}},
6088  {{X, X, X, X}}}},
6089  {{{{4, 7, 8, X}},
6090  {{7, 4, 6, X}},
6091  {{4, 8, 6, X}},
6092  {{8, 7, 6, X}},
6093  {{X, X, X, X}},
6094  {{X, X, X, X}}}},
6095  {{{{6, 9, 7, X}},
6096  {{9, 6, 8, X}},
6097  {{6, 7, 8, X}},
6098  {{7, 9, 8, X}},
6099  {{X, X, X, X}},
6100  {{X, X, X, X}}}},
6101  {{{{5, 8, 9, X}},
6102  {{8, 5, 6, X}},
6103  {{5, 9, 6, X}},
6104  {{9, 8, 6, X}},
6105  {{X, X, X, X}},
6106  {{X, X, X, X}}}}}};
6107 
6108  const auto &cell_quads =
6109  (reference_cell_type == ReferenceCells::Hexahedron) ?
6110  cell_quads_hex :
6111  cell_quads_tet;
6112 
6113  for (unsigned int c = 0;
6114  c < GeometryInfo<dim>::max_children_per_cell;
6115  ++c)
6116  {
6117  auto &new_hex = new_hexes[c];
6118 
6119  if (new_hex->n_faces() == 4)
6120  {
6121  new_hex->set_bounding_object_indices(
6122  {quad_indices[cell_quads[c][0]],
6123  quad_indices[cell_quads[c][1]],
6124  quad_indices[cell_quads[c][2]],
6125  quad_indices[cell_quads[c][3]]});
6126 
6127  // for tets, we need to go through the faces and
6128  // figure the orientation out the hard way
6129  for (const auto f : new_hex->face_indices())
6130  {
6131  const auto &face = new_hex->face(f);
6132 
6133  Assert(face->n_vertices() == 3,
6134  ExcInternalError());
6135 
6136  const std::array<unsigned int, 3> vertices_0 = {
6137  {face->vertex_index(0),
6138  face->vertex_index(1),
6139  face->vertex_index(2)}};
6140 
6141  const std::array<unsigned int, 3> vertices_1 = {
6142  {
6143  vertex_indices[cell_face_vertices_tet[c][f]
6144  [0]],
6145  vertex_indices[cell_face_vertices_tet[c][f]
6146  [1]],
6147  vertex_indices[cell_face_vertices_tet[c][f]
6148  [2]],
6149  }};
6150 
6151  new_hex->set_combined_face_orientation(
6152  f,
6153  face->reference_cell()
6154  .get_combined_orientation(
6155  make_array_view(vertices_1),
6156  make_array_view(vertices_0)));
6157  }
6158  }
6159  else if (new_hex->n_faces() == 6)
6160  new_hex->set_bounding_object_indices(
6161  {quad_indices[cell_quads[c][0]],
6162  quad_indices[cell_quads[c][1]],
6163  quad_indices[cell_quads[c][2]],
6164  quad_indices[cell_quads[c][3]],
6165  quad_indices[cell_quads[c][4]],
6166  quad_indices[cell_quads[c][5]]});
6167  else
6168  Assert(false, ExcNotImplemented());
6169  }
6170 
6171  // for hexes, we can simply inherit the orientation values
6172  // from the parent on the outer faces; the inner faces can
6173  // be skipped as their orientation is always the default
6174  // one set above
6175  static constexpr ::ndarray<unsigned int, 6, 4>
6176  face_to_child_indices_hex{{{{0, 2, 4, 6}},
6177  {{1, 3, 5, 7}},
6178  {{0, 1, 4, 5}},
6179  {{2, 3, 6, 7}},
6180  {{0, 1, 2, 3}},
6181  {{4, 5, 6, 7}}}};
6182  if (hex->n_faces() == 6)
6183  for (const auto f : hex->face_indices())
6184  {
6185  const unsigned char combined_orientation =
6186  hex->combined_face_orientation(f);
6187  for (unsigned int c = 0; c < 4; ++c)
6188  new_hexes[face_to_child_indices_hex[f][c]]
6189  ->set_combined_face_orientation(
6190  f, combined_orientation);
6191  }
6192  }
6193  }
6194 
6195  if (check_for_distorted_cells &&
6196  has_distorted_children<dim, spacedim>(hex))
6197  cells_with_distorted_children.distorted_cells.push_back(hex);
6198 
6199  triangulation.signals.post_refinement_on_cell(hex);
6200  }
6201  }
6202 
6203  triangulation.faces->quads.clear_user_data();
6204 
6205  return cells_with_distorted_children;
6206  }
6207 
6212  template <int spacedim>
6215  const bool check_for_distorted_cells)
6216  {
6217  const unsigned int dim = 3;
6218 
6219  {
6220  bool flag_isotropic_mesh = true;
6222  cell = triangulation.begin(),
6223  endc = triangulation.end();
6224  for (; cell != endc; ++cell)
6225  if (cell->used())
6226  if (triangulation.get_anisotropic_refinement_flag() ||
6227  cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
6228  cell->refine_flag_set() == RefinementCase<dim>::cut_y ||
6229  cell->refine_flag_set() == RefinementCase<dim>::cut_z ||
6230  cell->refine_flag_set() == RefinementCase<dim>::cut_xy ||
6231  cell->refine_flag_set() == RefinementCase<dim>::cut_xz ||
6232  cell->refine_flag_set() == RefinementCase<dim>::cut_yz)
6233  {
6234  flag_isotropic_mesh = false;
6235  break;
6236  }
6237 
6238  if (flag_isotropic_mesh)
6239  return execute_refinement_isotropic(triangulation,
6240  check_for_distorted_cells);
6241  }
6242 
6243  // this function probably also works for spacedim>3 but it
6244  // isn't tested. it will probably be necessary to pull new
6245  // vertices onto the manifold just as we do for the other
6246  // functions above.
6247  Assert(spacedim == 3, ExcNotImplemented());
6248 
6249  // Check whether a new level is needed. We have to check for
6250  // this on the highest level only
6251  for (const auto &cell : triangulation.active_cell_iterators_on_level(
6252  triangulation.levels.size() - 1))
6253  if (cell->refine_flag_set())
6254  {
6255  triangulation.levels.push_back(
6256  std::make_unique<
6258  break;
6259  }
6260 
6261 
6262  // first clear user flags for quads and lines; we're going to
6263  // use them to flag which lines and quads need refinement
6264  triangulation.faces->quads.clear_user_data();
6265 
6266  for (typename Triangulation<dim, spacedim>::line_iterator line =
6267  triangulation.begin_line();
6268  line != triangulation.end_line();
6269  ++line)
6270  line->clear_user_flag();
6271  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
6272  triangulation.begin_quad();
6273  quad != triangulation.end_quad();
6274  ++quad)
6275  quad->clear_user_flag();
6276 
6277  // create an array of face refine cases. User indices of faces
6278  // will be set to values corresponding with indices in this
6279  // array.
6280  const RefinementCase<dim - 1> face_refinement_cases[4] = {
6281  RefinementCase<dim - 1>::no_refinement,
6282  RefinementCase<dim - 1>::cut_x,
6283  RefinementCase<dim - 1>::cut_y,
6284  RefinementCase<dim - 1>::cut_xy};
6285 
6286  // check how much space is needed on every level. We need not
6287  // check the highest level since either
6288  // - on the highest level no cells are flagged for refinement
6289  // - there are, but prepare_refinement added another empty
6290  // level which then is the highest level
6291 
6292  // variables to hold the number of newly to be created
6293  // vertices, lines and quads. as these are stored globally,
6294  // declare them outside the loop over al levels. we need lines
6295  // and quads in pairs for refinement of old ones and lines and
6296  // quads, that can be stored as single ones, as they are newly
6297  // created in the inside of an existing cell
6298  unsigned int needed_vertices = 0;
6299  unsigned int needed_lines_single = 0;
6300  unsigned int needed_quads_single = 0;
6301  unsigned int needed_lines_pair = 0;
6302  unsigned int needed_quads_pair = 0;
6303  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
6304  {
6305  // count number of flagged cells on this level and compute
6306  // how many new vertices and new lines will be needed
6307  unsigned int new_cells = 0;
6308 
6309  for (const auto &acell :
6310  triangulation.active_cell_iterators_on_level(level))
6311  if (acell->refine_flag_set())
6312  {
6313  RefinementCase<dim> ref_case = acell->refine_flag_set();
6314 
6315  // now for interior vertices, lines and quads, which
6316  // are needed in any case
6317  if (ref_case == RefinementCase<dim>::cut_x ||
6318  ref_case == RefinementCase<dim>::cut_y ||
6319  ref_case == RefinementCase<dim>::cut_z)
6320  {
6321  ++needed_quads_single;
6322  new_cells += 2;
6323  triangulation.anisotropic_refinement = true;
6324  }
6325  else if (ref_case == RefinementCase<dim>::cut_xy ||
6326  ref_case == RefinementCase<dim>::cut_xz ||
6327  ref_case == RefinementCase<dim>::cut_yz)
6328  {
6329  ++needed_lines_single;
6330  needed_quads_single += 4;
6331  new_cells += 4;
6332  triangulation.anisotropic_refinement = true;
6333  }
6334  else if (ref_case == RefinementCase<dim>::cut_xyz)
6335  {
6336  ++needed_vertices;
6337  needed_lines_single += 6;
6338  needed_quads_single += 12;
6339  new_cells += 8;
6340  }
6341  else
6342  {
6343  // we should never get here
6344  Assert(false, ExcInternalError());
6345  }
6346 
6347  // mark all faces for refinement; checking locally
6348  // if and how the neighbor would like to refine
6349  // these is difficult so we only flag them and after
6350  // visiting all cells, we decide which faces need
6351  // which refinement;
6352  for (const unsigned int face :
6354  {
6356  aface = acell->face(face);
6357  // get the RefineCase this faces has for the
6358  // given RefineCase of the cell
6359  RefinementCase<dim - 1> face_ref_case =
6361  ref_case,
6362  face,
6363  acell->face_orientation(face),
6364  acell->face_flip(face),
6365  acell->face_rotation(face));
6366  // only do something, if this face has to be
6367  // refined
6368  if (face_ref_case)
6369  {
6370  if (face_ref_case ==
6372  {
6373  if (aface->n_active_descendants() < 4)
6374  // we use user_flags to denote needed
6375  // isotropic refinement
6376  aface->set_user_flag();
6377  }
6378  else if (aface->refinement_case() != face_ref_case)
6379  // we use user_indices to denote needed
6380  // anisotropic refinement. note, that we
6381  // can have at most one anisotropic
6382  // refinement case for this face, as
6383  // otherwise prepare_refinement() would
6384  // have changed one of the cells to yield
6385  // isotropic refinement at this
6386  // face. therefore we set the user_index
6387  // uniquely
6388  {
6389  Assert(aface->refinement_case() ==
6391  dim - 1>::isotropic_refinement ||
6392  aface->refinement_case() ==
6394  ExcInternalError());
6395  aface->set_user_index(face_ref_case);
6396  }
6397  }
6398  } // for all faces
6399 
6400  // flag all lines, that have to be refined
6401  for (unsigned int line = 0;
6402  line < GeometryInfo<dim>::lines_per_cell;
6403  ++line)
6405  line) &&
6406  !acell->line(line)->has_children())
6407  acell->line(line)->set_user_flag();
6408 
6409  } // if refine_flag set and for all cells on this level
6410 
6411 
6412  // count number of used cells on the next higher level
6413  const unsigned int used_cells =
6414  std::count(triangulation.levels[level + 1]->cells.used.begin(),
6415  triangulation.levels[level + 1]->cells.used.end(),
6416  true);
6417 
6418 
6419  // reserve space for the used_cells cells already existing
6420  // on the next higher level as well as for the
6421  // 8*flagged_cells that will be created on that level
6422  reserve_space(*triangulation.levels[level + 1],
6423  used_cells + new_cells,
6424  3,
6425  spacedim);
6426  // reserve space for 8*flagged_cells new hexes on the next
6427  // higher level
6428  reserve_space(triangulation.levels[level + 1]->cells, new_cells);
6429  } // for all levels
6430  // now count the quads and lines which were flagged for
6431  // refinement
6432  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
6433  triangulation.begin_quad();