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