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