Reference documentation for deal.II version GIT 725517848b 2022-05-23 21:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
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 
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 
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  1u);
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  1u);
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
1599 
1603  virtual void
1607  std::vector<unsigned int> & line_cell_count,
1608  std::vector<unsigned int> &quad_cell_count) = 0;
1609 
1615  const bool check_for_distorted_cells) = 0;
1616 
1620  virtual void
1623 
1627  virtual void
1630 
1634  virtual bool
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 
1841  //---------------------------------
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 
1933  //---------------------------------
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 
2040  //---------------------------------
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) ==
2250  RefinementCase<dim - 1>::isotropic_refinement))
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.
2317  {
2318  const double cell_measure = GridTools::cell_measure<spacedim>(
2319  vertices,
2320  ArrayView<const unsigned int>(cells[cell_no].vertices));
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
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!?
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  vertices);
2527 
2528  // SubCellData: quad
2529  if (dim == 3)
2530  process_subcelldata(connectivity.entity_to_entities(2, 0),
2531  tria.faces->quads,
2532  subcelldata.boundary_quads,
2533  vertices);
2534  }
2535 
2536 
2537  template <int structdim, int spacedim, typename T>
2538  static void
2540  const CRS<T> & crs,
2541  TriaObjects & obj,
2542  const std::vector<CellData<structdim>> &boundary_objects_in,
2543  const std::vector<Point<spacedim>> & vertex_locations)
2544  {
2545  AssertDimension(obj.structdim, structdim);
2546 
2547  if (boundary_objects_in.size() == 0)
2548  return; // empty subcelldata -> nothing to do
2549 
2550  // pre-sort subcelldata
2551  auto boundary_objects = boundary_objects_in;
2552 
2553  // ... sort vertices
2554  for (auto &boundary_object : boundary_objects)
2555  std::sort(boundary_object.vertices.begin(),
2556  boundary_object.vertices.end());
2557 
2558  // ... sort cells
2559  std::sort(boundary_objects.begin(),
2560  boundary_objects.end(),
2561  [](const auto &a, const auto &b) {
2562  return a.vertices < b.vertices;
2563  });
2564 
2565  unsigned int counter = 0;
2566 
2567  std::vector<unsigned int> key;
2569 
2570  for (unsigned int o = 0; o < obj.n_objects(); ++o)
2571  {
2572  auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
2573  auto &manifold_id = obj.manifold_id[o];
2574 
2575  // assert that object has not been visited yet and its value
2576  // has not been modified yet
2577  AssertThrow(boundary_id == 0 ||
2579  ExcNotImplemented());
2581  ExcNotImplemented());
2582 
2583  // create key
2584  key.assign(crs.col.data() + crs.ptr[o],
2585  crs.col.data() + crs.ptr[o + 1]);
2586  std::sort(key.begin(), key.end());
2587 
2588  // is subcelldata provided? -> binary search
2589  const auto subcell_object =
2590  std::lower_bound(boundary_objects.begin(),
2591  boundary_objects.end(),
2592  key,
2593  [&](const auto &cell, const auto &key) {
2594  return cell.vertices < key;
2595  });
2596 
2597  // no subcelldata provided for this object
2598  if (subcell_object == boundary_objects.end() ||
2599  subcell_object->vertices != key)
2600  continue;
2601 
2602  counter++;
2603 
2604  // set manifold id
2605  manifold_id = subcell_object->manifold_id;
2606 
2607  // set boundary id
2608  if (subcell_object->boundary_id !=
2610  {
2611  (void)vertex_locations;
2612  AssertThrow(
2614  ExcMessage(
2615  "The input arguments for creating a triangulation "
2616  "specified a boundary id for an internal face. This "
2617  "is not allowed."
2618  "\n\n"
2619  "The object in question has vertex indices " +
2620  [subcell_object]() {
2621  std::string s;
2622  for (const auto v : subcell_object->vertices)
2623  s += std::to_string(v) + ',';
2624  return s;
2625  }() +
2626  " which are located at positions " +
2627  [vertex_locations, subcell_object]() {
2628  std::ostringstream s;
2629  for (const auto v : subcell_object->vertices)
2630  s << '(' << vertex_locations[v] << ')';
2631  return s.str();
2632  }() +
2633  "."));
2634  boundary_id = subcell_object->boundary_id;
2635  }
2636  }
2637 
2638  // make sure that all subcelldata entries have been processed
2639  // TODO: this is not guaranteed, why?
2640  // AssertDimension(counter, boundary_objects_in.size());
2641  }
2642 
2643 
2644 
2645  static void
2647  const unsigned structdim,
2648  const unsigned int size)
2649  {
2650  const unsigned int dim = faces.dim;
2651 
2652  const unsigned int max_faces_per_cell = 2 * structdim;
2653 
2654  if (dim == 3 && structdim == 2)
2655  {
2656  // quad entity types
2657  faces.quad_reference_cell.assign(size,
2659 
2660  // quad line orientations
2661  faces.quads_line_orientations.assign(size * max_faces_per_cell, -1);
2662  }
2663  }
2664 
2665 
2666 
2667  static void
2669  const unsigned int spacedim,
2670  const unsigned int size,
2671  const bool orientation_needed)
2672  {
2673  const unsigned int dim = level.dim;
2674 
2675  const unsigned int max_faces_per_cell = 2 * dim;
2676 
2677  level.active_cell_indices.assign(size, -1);
2678  level.subdomain_ids.assign(size, 0);
2679  level.level_subdomain_ids.assign(size, 0);
2680 
2681  level.refine_flags.assign(size, 0u);
2682  level.coarsen_flags.assign(size, false);
2683 
2684  level.parents.assign((size + 1) / 2, -1);
2685 
2686  if (dim < spacedim)
2687  level.direction_flags.assign(size, true);
2688 
2689  level.neighbors.assign(size * max_faces_per_cell, {-1, -1});
2690 
2691  level.reference_cell.assign(size, ::ReferenceCells::Invalid);
2692 
2693  if (orientation_needed)
2694  level.face_orientations.assign(size * max_faces_per_cell, -1);
2695 
2696  level.global_active_cell_indices.assign(size,
2698  level.global_level_cell_indices.assign(size,
2700  }
2701 
2702 
2703 
2704  static void
2705  reserve_space_(TriaObjects &obj, const unsigned int size)
2706  {
2707  const unsigned int structdim = obj.structdim;
2708 
2709  const unsigned int max_children_per_cell = 1 << structdim;
2710  const unsigned int max_faces_per_cell = 2 * structdim;
2711 
2712  obj.used.assign(size, true);
2713  obj.boundary_or_material_id.assign(
2714  size,
2716  BoundaryOrMaterialId());
2717  obj.manifold_id.assign(size, -1);
2718  obj.user_flags.assign(size, false);
2719  obj.user_data.resize(size);
2720 
2721  if (structdim > 1) // TODO: why?
2722  obj.refinement_cases.assign(size, 0);
2723 
2724  obj.children.assign(max_children_per_cell / 2 * size, -1);
2725 
2726  obj.cells.assign(max_faces_per_cell * size, -1);
2727 
2728  if (structdim <= 2)
2729  {
2730  obj.next_free_single = size - 1;
2731  obj.next_free_pair = 0;
2732  obj.reverse_order_next_free_single = true;
2733  }
2734  else
2735  {
2736  obj.next_free_single = obj.next_free_pair = 0;
2737  }
2738  }
2739 
2740 
2756  template <int spacedim>
2757  static void
2760  std::vector<unsigned int> &,
2761  std::vector<unsigned int> &)
2762  {
2763  const unsigned int dim = 1;
2764 
2765  // first we need to reset the
2766  // neighbor pointers of the
2767  // neighbors of this cell's
2768  // children to this cell. This is
2769  // different for one dimension,
2770  // since there neighbors can have a
2771  // refinement level differing from
2772  // that of this cell's children by
2773  // more than one level.
2774 
2775  Assert(!cell->child(0)->has_children() &&
2776  !cell->child(1)->has_children(),
2777  ExcInternalError());
2778 
2779  // first do it for the cells to the
2780  // left
2781  if (cell->neighbor(0).state() == IteratorState::valid)
2782  if (cell->neighbor(0)->has_children())
2783  {
2785  cell->neighbor(0);
2786  Assert(neighbor->level() == cell->level(), ExcInternalError());
2787 
2788  // right child
2789  neighbor = neighbor->child(1);
2790  while (true)
2791  {
2792  Assert(neighbor->neighbor(1) == cell->child(0),
2793  ExcInternalError());
2794  neighbor->set_neighbor(1, cell);
2795 
2796  // move on to further
2797  // children on the
2798  // boundary between this
2799  // cell and its neighbor
2800  if (neighbor->has_children())
2801  neighbor = neighbor->child(1);
2802  else
2803  break;
2804  }
2805  }
2806 
2807  // now do it for the cells to the
2808  // left
2809  if (cell->neighbor(1).state() == IteratorState::valid)
2810  if (cell->neighbor(1)->has_children())
2811  {
2813  cell->neighbor(1);
2814  Assert(neighbor->level() == cell->level(), ExcInternalError());
2815 
2816  // left child
2817  neighbor = neighbor->child(0);
2818  while (true)
2819  {
2820  Assert(neighbor->neighbor(0) == cell->child(1),
2821  ExcInternalError());
2822  neighbor->set_neighbor(0, cell);
2823 
2824  // move on to further
2825  // children on the
2826  // boundary between this
2827  // cell and its neighbor
2828  if (neighbor->has_children())
2829  neighbor = neighbor->child(0);
2830  else
2831  break;
2832  }
2833  }
2834 
2835 
2836  // delete the vertex which will not
2837  // be needed anymore. This vertex
2838  // is the second of the first child
2839  triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
2840 
2841  // invalidate children. clear user
2842  // pointers, to avoid that they may
2843  // appear at unwanted places later
2844  // on...
2845  for (unsigned int child = 0; child < cell->n_children(); ++child)
2846  {
2847  cell->child(child)->clear_user_data();
2848  cell->child(child)->clear_user_flag();
2849  cell->child(child)->clear_used_flag();
2850  }
2851 
2852 
2853  // delete pointer to children
2854  cell->clear_children();
2855  cell->clear_user_flag();
2856  }
2857 
2858 
2859 
2860  template <int spacedim>
2861  static void
2864  std::vector<unsigned int> &line_cell_count,
2865  std::vector<unsigned int> &)
2866  {
2867  const unsigned int dim = 2;
2868  const RefinementCase<dim> ref_case = cell->refinement_case();
2869 
2870  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2871  ExcInternalError());
2872 
2873  // vectors to hold all lines which
2874  // may be deleted
2875  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2876  lines_to_delete(0);
2877 
2878  lines_to_delete.reserve(4 * 2 + 4);
2879 
2880  // now we decrease the counters for
2881  // lines contained in the child
2882  // cells
2883  for (unsigned int c = 0; c < cell->n_children(); ++c)
2884  {
2886  cell->child(c);
2887  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2888  --line_cell_count[child->line_index(l)];
2889  }
2890 
2891 
2892  // delete the vertex which will not
2893  // be needed anymore. This vertex
2894  // is the second of the second line
2895  // of the first child, if the cell
2896  // is refined with cut_xy, else there
2897  // is no inner vertex.
2898  // additionally delete unneeded inner
2899  // lines
2900  if (ref_case == RefinementCase<dim>::cut_xy)
2901  {
2903  .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
2904 
2905  lines_to_delete.push_back(cell->child(0)->line(1));
2906  lines_to_delete.push_back(cell->child(0)->line(3));
2907  lines_to_delete.push_back(cell->child(3)->line(0));
2908  lines_to_delete.push_back(cell->child(3)->line(2));
2909  }
2910  else
2911  {
2912  unsigned int inner_face_no =
2913  ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
2914 
2915  // the inner line will not be
2916  // used any more
2917  lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2918  }
2919 
2920  // invalidate children
2921  for (unsigned int child = 0; child < cell->n_children(); ++child)
2922  {
2923  cell->child(child)->clear_user_data();
2924  cell->child(child)->clear_user_flag();
2925  cell->child(child)->clear_used_flag();
2926  }
2927 
2928 
2929  // delete pointer to children
2930  cell->clear_children();
2931  cell->clear_refinement_case();
2932  cell->clear_user_flag();
2933 
2934  // look at the refinement of outer
2935  // lines. if nobody needs those
2936  // anymore we can add them to the
2937  // list of lines to be deleted.
2938  for (unsigned int line_no = 0;
2939  line_no < GeometryInfo<dim>::lines_per_cell;
2940  ++line_no)
2941  {
2943  cell->line(line_no);
2944 
2945  if (line->has_children())
2946  {
2947  // if one of the cell counters is
2948  // zero, the other has to be as well
2949 
2950  Assert((line_cell_count[line->child_index(0)] == 0 &&
2951  line_cell_count[line->child_index(1)] == 0) ||
2952  (line_cell_count[line->child_index(0)] > 0 &&
2953  line_cell_count[line->child_index(1)] > 0),
2954  ExcInternalError());
2955 
2956  if (line_cell_count[line->child_index(0)] == 0)
2957  {
2958  for (unsigned int c = 0; c < 2; ++c)
2959  Assert(!line->child(c)->has_children(),
2960  ExcInternalError());
2961 
2962  // we may delete the line's
2963  // children and the middle vertex
2964  // as no cell references them
2965  // anymore
2967  .vertices_used[line->child(0)->vertex_index(1)] = false;
2968 
2969  lines_to_delete.push_back(line->child(0));
2970  lines_to_delete.push_back(line->child(1));
2971 
2972  line->clear_children();
2973  }
2974  }
2975  }
2976 
2977  // finally, delete unneeded lines
2978 
2979  // clear user pointers, to avoid that
2980  // they may appear at unwanted places
2981  // later on...
2982  // same for user flags, then finally
2983  // delete the lines
2984  typename std::vector<
2986  line = lines_to_delete.begin(),
2987  endline = lines_to_delete.end();
2988  for (; line != endline; ++line)
2989  {
2990  (*line)->clear_user_data();
2991  (*line)->clear_user_flag();
2992  (*line)->clear_used_flag();
2993  }
2994  }
2995 
2996 
2997 
2998  template <int spacedim>
2999  static void
3002  std::vector<unsigned int> &line_cell_count,
3003  std::vector<unsigned int> &quad_cell_count)
3004  {
3005  const unsigned int dim = 3;
3006 
3007  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
3008  ExcInternalError());
3009  Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
3010  ExcInternalError());
3011 
3012  // first of all, we store the RefineCase of
3013  // this cell
3014  const RefinementCase<dim> ref_case = cell->refinement_case();
3015  // vectors to hold all lines and quads which
3016  // may be deleted
3017  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3018  lines_to_delete(0);
3019  std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3020  quads_to_delete(0);
3021 
3022  lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3023  quads_to_delete.reserve(6 * 4 + 12);
3024 
3025  // now we decrease the counters for lines and
3026  // quads contained in the child cells
3027  for (unsigned int c = 0; c < cell->n_children(); ++c)
3028  {
3030  cell->child(c);
3031  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3032  --line_cell_count[child->line_index(l)];
3033  for (auto f : GeometryInfo<dim>::face_indices())
3034  --quad_cell_count[child->quad_index(f)];
3035  }
3036 
3037  //-------------------------------------
3038  // delete interior quads and lines and the
3039  // interior vertex, depending on the
3040  // refinement case of the cell
3041  //
3042  // for append quads and lines: only append
3043  // them to the list of objects to be deleted
3044 
3045  switch (ref_case)
3046  {
3048  quads_to_delete.push_back(cell->child(0)->face(1));
3049  break;
3051  quads_to_delete.push_back(cell->child(0)->face(3));
3052  break;
3054  quads_to_delete.push_back(cell->child(0)->face(5));
3055  break;
3057  quads_to_delete.push_back(cell->child(0)->face(1));
3058  quads_to_delete.push_back(cell->child(0)->face(3));
3059  quads_to_delete.push_back(cell->child(3)->face(0));
3060  quads_to_delete.push_back(cell->child(3)->face(2));
3061 
3062  lines_to_delete.push_back(cell->child(0)->line(11));
3063  break;
3065  quads_to_delete.push_back(cell->child(0)->face(1));
3066  quads_to_delete.push_back(cell->child(0)->face(5));
3067  quads_to_delete.push_back(cell->child(3)->face(0));
3068  quads_to_delete.push_back(cell->child(3)->face(4));
3069 
3070  lines_to_delete.push_back(cell->child(0)->line(5));
3071  break;
3073  quads_to_delete.push_back(cell->child(0)->face(3));
3074  quads_to_delete.push_back(cell->child(0)->face(5));
3075  quads_to_delete.push_back(cell->child(3)->face(2));
3076  quads_to_delete.push_back(cell->child(3)->face(4));
3077 
3078  lines_to_delete.push_back(cell->child(0)->line(7));
3079  break;
3081  quads_to_delete.push_back(cell->child(0)->face(1));
3082  quads_to_delete.push_back(cell->child(2)->face(1));
3083  quads_to_delete.push_back(cell->child(4)->face(1));
3084  quads_to_delete.push_back(cell->child(6)->face(1));
3085 
3086  quads_to_delete.push_back(cell->child(0)->face(3));
3087  quads_to_delete.push_back(cell->child(1)->face(3));
3088  quads_to_delete.push_back(cell->child(4)->face(3));
3089  quads_to_delete.push_back(cell->child(5)->face(3));
3090 
3091  quads_to_delete.push_back(cell->child(0)->face(5));
3092  quads_to_delete.push_back(cell->child(1)->face(5));
3093  quads_to_delete.push_back(cell->child(2)->face(5));
3094  quads_to_delete.push_back(cell->child(3)->face(5));
3095 
3096  lines_to_delete.push_back(cell->child(0)->line(5));
3097  lines_to_delete.push_back(cell->child(0)->line(7));
3098  lines_to_delete.push_back(cell->child(0)->line(11));
3099  lines_to_delete.push_back(cell->child(7)->line(0));
3100  lines_to_delete.push_back(cell->child(7)->line(2));
3101  lines_to_delete.push_back(cell->child(7)->line(8));
3102  // delete the vertex which will not
3103  // be needed anymore. This vertex
3104  // is the vertex at the heart of
3105  // this cell, which is the sixth of
3106  // the first child
3107  triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3108  false;
3109  break;
3110  default:
3111  // only remaining case is
3112  // no_refinement, thus an error
3113  Assert(false, ExcInternalError());
3114  break;
3115  }
3116 
3117 
3118  // invalidate children
3119  for (unsigned int child = 0; child < cell->n_children(); ++child)
3120  {
3121  cell->child(child)->clear_user_data();
3122  cell->child(child)->clear_user_flag();
3123 
3124  for (auto f : GeometryInfo<dim>::face_indices())
3125  {
3126  // set flags denoting deviations from
3127  // standard orientation of faces back
3128  // to initialization values
3129  cell->child(child)->set_face_orientation(f, true);
3130  cell->child(child)->set_face_flip(f, false);
3131  cell->child(child)->set_face_rotation(f, false);
3132  }
3133 
3134  cell->child(child)->clear_used_flag();
3135  }
3136 
3137 
3138  // delete pointer to children
3139  cell->clear_children();
3140  cell->clear_refinement_case();
3141  cell->clear_user_flag();
3142 
3143  // so far we only looked at inner quads,
3144  // lines and vertices. Now we have to
3145  // consider outer ones as well. here, we have
3146  // to check, whether there are other cells
3147  // still needing these objects. otherwise we
3148  // can delete them. first for quads (and
3149  // their inner lines).
3150 
3151  for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
3152  {
3154  cell->face(quad_no);
3155 
3156  Assert(
3157  (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
3158  quad->has_children()) ||
3159  GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
3161  ExcInternalError());
3162 
3163  switch (quad->refinement_case())
3164  {
3165  case RefinementCase<dim - 1>::no_refinement:
3166  // nothing to do as the quad
3167  // is not refined
3168  break;
3169  case RefinementCase<dim - 1>::cut_x:
3170  case RefinementCase<dim - 1>::cut_y:
3171  {
3172  // if one of the cell counters is
3173  // zero, the other has to be as
3174  // well
3175  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3176  quad_cell_count[quad->child_index(1)] == 0) ||
3177  (quad_cell_count[quad->child_index(0)] > 0 &&
3178  quad_cell_count[quad->child_index(1)] > 0),
3179  ExcInternalError());
3180  // it might be, that the quad is
3181  // refined twice anisotropically,
3182  // first check, whether we may
3183  // delete possible grand_children
3184  unsigned int deleted_grandchildren = 0;
3185  unsigned int number_of_child_refinements = 0;
3186 
3187  for (unsigned int c = 0; c < 2; ++c)
3188  if (quad->child(c)->has_children())
3189  {
3190  ++number_of_child_refinements;
3191  // if one of the cell counters is
3192  // zero, the other has to be as
3193  // well
3194  Assert(
3195  (quad_cell_count[quad->child(c)->child_index(0)] ==
3196  0 &&
3197  quad_cell_count[quad->child(c)->child_index(1)] ==
3198  0) ||
3199  (quad_cell_count[quad->child(c)->child_index(0)] >
3200  0 &&
3201  quad_cell_count[quad->child(c)->child_index(1)] >
3202  0),
3203  ExcInternalError());
3204  if (quad_cell_count[quad->child(c)->child_index(0)] ==
3205  0)
3206  {
3207  // Assert, that the two
3208  // anisotropic
3209  // refinements add up to
3210  // isotropic refinement
3211  Assert(quad->refinement_case() +
3212  quad->child(c)->refinement_case() ==
3214  ExcInternalError());
3215  // we may delete the
3216  // quad's children and
3217  // the inner line as no
3218  // cell references them
3219  // anymore
3220  quads_to_delete.push_back(
3221  quad->child(c)->child(0));
3222  quads_to_delete.push_back(
3223  quad->child(c)->child(1));
3224  if (quad->child(c)->refinement_case() ==
3226  lines_to_delete.push_back(
3227  quad->child(c)->child(0)->line(1));
3228  else
3229  lines_to_delete.push_back(
3230  quad->child(c)->child(0)->line(3));
3231  quad->child(c)->clear_children();
3232  quad->child(c)->clear_refinement_case();
3233  ++deleted_grandchildren;
3234  }
3235  }
3236  // if no grandchildren are left, we
3237  // may as well delete the
3238  // refinement of the inner line
3239  // between our children and the
3240  // corresponding vertex
3241  if (number_of_child_refinements > 0 &&
3242  deleted_grandchildren == number_of_child_refinements)
3243  {
3245  middle_line;
3246  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3247  middle_line = quad->child(0)->line(1);
3248  else
3249  middle_line = quad->child(0)->line(3);
3250 
3251  lines_to_delete.push_back(middle_line->child(0));
3252  lines_to_delete.push_back(middle_line->child(1));
3254  .vertices_used[middle_vertex_index<dim, spacedim>(
3255  middle_line)] = false;
3256  middle_line->clear_children();
3257  }
3258 
3259  // now consider the direct children
3260  // of the given quad
3261  if (quad_cell_count[quad->child_index(0)] == 0)
3262  {
3263  // we may delete the quad's
3264  // children and the inner line
3265  // as no cell references them
3266  // anymore
3267  quads_to_delete.push_back(quad->child(0));
3268  quads_to_delete.push_back(quad->child(1));
3269  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3270  lines_to_delete.push_back(quad->child(0)->line(1));
3271  else
3272  lines_to_delete.push_back(quad->child(0)->line(3));
3273 
3274  // if the counters just dropped
3275  // to zero, otherwise the
3276  // children would have been
3277  // deleted earlier, then this
3278  // cell's children must have
3279  // contained the anisotropic
3280  // quad children. thus, if
3281  // those have again anisotropic
3282  // children, which are in
3283  // effect isotropic children of
3284  // the original quad, those are
3285  // still needed by a
3286  // neighboring cell and we
3287  // cannot delete them. instead,
3288  // we have to reset this quad's
3289  // refine case to isotropic and
3290  // set the children
3291  // accordingly.
3292  if (quad->child(0)->has_children())
3293  if (quad->refinement_case() ==
3295  {
3296  // now evereything is
3297  // quite complicated. we
3298  // have the children
3299  // numbered according to
3300  //
3301  // *---*---*
3302  // |n+1|m+1|
3303  // *---*---*
3304  // | n | m |
3305  // *---*---*
3306  //
3307  // from the original
3308  // anisotropic
3309  // refinement. we have to
3310  // reorder them as
3311  //
3312  // *---*---*
3313  // | m |m+1|
3314  // *---*---*
3315  // | n |n+1|
3316  // *---*---*
3317  //
3318  // for isotropic refinement.
3319  //
3320  // this is a bit ugly, of
3321  // course: loop over all
3322  // cells on all levels
3323  // and look for faces n+1
3324  // (switch_1) and m
3325  // (switch_2).
3326  const typename Triangulation<dim, spacedim>::
3327  quad_iterator switch_1 =
3328  quad->child(0)->child(1),
3329  switch_2 =
3330  quad->child(1)->child(0);
3331 
3332  Assert(!switch_1->has_children(),
3333  ExcInternalError());
3334  Assert(!switch_2->has_children(),
3335  ExcInternalError());
3336 
3337  const int switch_1_index = switch_1->index();
3338  const int switch_2_index = switch_2->index();
3339  for (unsigned int l = 0;
3340  l < triangulation.levels.size();
3341  ++l)
3342  for (unsigned int h = 0;
3343  h <
3344  triangulation.levels[l]->cells.n_objects();
3345  ++h)
3346  for (const unsigned int q :
3348  {
3349  const int index =
3350  triangulation.levels[l]
3351  ->cells.get_bounding_object_indices(
3352  h)[q];
3353  if (index == switch_1_index)
3354  triangulation.levels[l]
3355  ->cells.get_bounding_object_indices(
3356  h)[q] = switch_2_index;
3357  else if (index == switch_2_index)
3358  triangulation.levels[l]
3359  ->cells.get_bounding_object_indices(
3360  h)[q] = switch_1_index;
3361  }
3362  // now we have to copy
3363  // all information of the
3364  // two quads
3365  const int switch_1_lines[4] = {
3366  static_cast<signed int>(
3367  switch_1->line_index(0)),
3368  static_cast<signed int>(
3369  switch_1->line_index(1)),
3370  static_cast<signed int>(
3371  switch_1->line_index(2)),
3372  static_cast<signed int>(
3373  switch_1->line_index(3))};
3374  const bool switch_1_line_orientations[4] = {
3375  switch_1->line_orientation(0),
3376  switch_1->line_orientation(1),
3377  switch_1->line_orientation(2),
3378  switch_1->line_orientation(3)};
3379  const types::boundary_id switch_1_boundary_id =
3380  switch_1->boundary_id();
3381  const unsigned int switch_1_user_index =
3382  switch_1->user_index();
3383  const bool switch_1_user_flag =
3384  switch_1->user_flag_set();
3385 
3386  switch_1->set_bounding_object_indices(
3387  {switch_2->line_index(0),
3388  switch_2->line_index(1),
3389  switch_2->line_index(2),
3390  switch_2->line_index(3)});
3391  switch_1->set_line_orientation(
3392  0, switch_2->line_orientation(0));
3393  switch_1->set_line_orientation(
3394  1, switch_2->line_orientation(1));
3395  switch_1->set_line_orientation(
3396  2, switch_2->line_orientation(2));
3397  switch_1->set_line_orientation(
3398  3, switch_2->line_orientation(3));
3399  switch_1->set_boundary_id_internal(
3400  switch_2->boundary_id());
3401  switch_1->set_manifold_id(
3402  switch_2->manifold_id());
3403  switch_1->set_user_index(switch_2->user_index());
3404  if (switch_2->user_flag_set())
3405  switch_1->set_user_flag();
3406  else
3407  switch_1->clear_user_flag();
3408 
3409  switch_2->set_bounding_object_indices(
3410  {switch_1_lines[0],
3411  switch_1_lines[1],
3412  switch_1_lines[2],
3413  switch_1_lines[3]});
3414  switch_2->set_line_orientation(
3415  0, switch_1_line_orientations[0]);
3416  switch_2->set_line_orientation(
3417  1, switch_1_line_orientations[1]);
3418  switch_2->set_line_orientation(
3419  2, switch_1_line_orientations[2]);
3420  switch_2->set_line_orientation(
3421  3, switch_1_line_orientations[3]);
3422  switch_2->set_boundary_id_internal(
3423  switch_1_boundary_id);
3424  switch_2->set_manifold_id(
3425  switch_1->manifold_id());
3426  switch_2->set_user_index(switch_1_user_index);
3427  if (switch_1_user_flag)
3428  switch_2->set_user_flag();
3429  else
3430  switch_2->clear_user_flag();
3431 
3432  const unsigned int child_0 =
3433  quad->child(0)->child_index(0);
3434  const unsigned int child_2 =
3435  quad->child(1)->child_index(0);
3436  quad->clear_children();
3437  quad->clear_refinement_case();
3438  quad->set_refinement_case(
3440  quad->set_children(0, child_0);
3441  quad->set_children(2, child_2);
3442  std::swap(quad_cell_count[child_0 + 1],
3443  quad_cell_count[child_2]);
3444  }
3445  else
3446  {
3447  // the face was refined
3448  // with cut_y, thus the
3449  // children are already
3450  // in correct order. we
3451  // only have to set them
3452  // correctly, deleting
3453  // the indirection of two
3454  // anisotropic refinement
3455  // and going directly
3456  // from the quad to
3457  // isotropic children
3458  const unsigned int child_0 =
3459  quad->child(0)->child_index(0);
3460  const unsigned int child_2 =
3461  quad->child(1)->child_index(0);
3462  quad->clear_children();
3463  quad->clear_refinement_case();
3464  quad->set_refinement_case(
3466  quad->set_children(0, child_0);
3467  quad->set_children(2, child_2);
3468  }
3469  else
3470  {
3471  quad->clear_children();
3472  quad->clear_refinement_case();
3473  }
3474  }
3475  break;
3476  }
3477  case RefinementCase<dim - 1>::cut_xy:
3478  {
3479  // if one of the cell counters is
3480  // zero, the others have to be as
3481  // well
3482 
3483  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3484  quad_cell_count[quad->child_index(1)] == 0 &&
3485  quad_cell_count[quad->child_index(2)] == 0 &&
3486  quad_cell_count[quad->child_index(3)] == 0) ||
3487  (quad_cell_count[quad->child_index(0)] > 0 &&
3488  quad_cell_count[quad->child_index(1)] > 0 &&
3489  quad_cell_count[quad->child_index(2)] > 0 &&
3490  quad_cell_count[quad->child_index(3)] > 0),
3491  ExcInternalError());
3492 
3493  if (quad_cell_count[quad->child_index(0)] == 0)
3494  {
3495  // we may delete the quad's
3496  // children, the inner lines
3497  // and the middle vertex as no
3498  // cell references them anymore
3499  lines_to_delete.push_back(quad->child(0)->line(1));
3500  lines_to_delete.push_back(quad->child(3)->line(0));
3501  lines_to_delete.push_back(quad->child(0)->line(3));
3502  lines_to_delete.push_back(quad->child(3)->line(2));
3503 
3504  for (unsigned int child = 0; child < quad->n_children();
3505  ++child)
3506  quads_to_delete.push_back(quad->child(child));
3507 
3509  .vertices_used[quad->child(0)->vertex_index(3)] =
3510  false;
3511 
3512  quad->clear_children();
3513  quad->clear_refinement_case();
3514  }
3515  }
3516  break;
3517 
3518  default:
3519  Assert(false, ExcInternalError());
3520  break;
3521  }
3522  }
3523 
3524  // now we repeat a similar procedure
3525  // for the outer lines of this cell.
3526 
3527  // if in debug mode: check that each
3528  // of the lines for which we consider
3529  // deleting the children in fact has
3530  // children (the bits/coarsening_3d
3531  // test tripped over this initially)
3532  for (unsigned int line_no = 0;
3533  line_no < GeometryInfo<dim>::lines_per_cell;
3534  ++line_no)
3535  {
3537  cell->line(line_no);
3538 
3539  Assert(
3540  (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
3541  line->has_children()) ||
3542  GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
3544  ExcInternalError());
3545 
3546  if (line->has_children())
3547  {
3548  // if one of the cell counters is
3549  // zero, the other has to be as well
3550 
3551  Assert((line_cell_count[line->child_index(0)] == 0 &&
3552  line_cell_count[line->child_index(1)] == 0) ||
3553  (line_cell_count[line->child_index(0)] > 0 &&
3554  line_cell_count[line->child_index(1)] > 0),
3555  ExcInternalError());
3556 
3557  if (line_cell_count[line->child_index(0)] == 0)
3558  {
3559  for (unsigned int c = 0; c < 2; ++c)
3560  Assert(!line->child(c)->has_children(),
3561  ExcInternalError());
3562 
3563  // we may delete the line's
3564  // children and the middle vertex
3565  // as no cell references them
3566  // anymore
3568  .vertices_used[line->child(0)->vertex_index(1)] = false;
3569 
3570  lines_to_delete.push_back(line->child(0));
3571  lines_to_delete.push_back(line->child(1));
3572 
3573  line->clear_children();
3574  }
3575  }
3576  }
3577 
3578  // finally, delete unneeded quads and lines
3579 
3580  // clear user pointers, to avoid that
3581  // they may appear at unwanted places
3582  // later on...
3583  // same for user flags, then finally
3584  // delete the quads and lines
3585  typename std::vector<
3587  line = lines_to_delete.begin(),
3588  endline = lines_to_delete.end();
3589  for (; line != endline; ++line)
3590  {
3591  (*line)->clear_user_data();
3592  (*line)->clear_user_flag();
3593  (*line)->clear_used_flag();
3594  }
3595 
3596  typename std::vector<
3598  quad = quads_to_delete.begin(),
3599  endquad = quads_to_delete.end();
3600  for (; quad != endquad; ++quad)
3601  {
3602  (*quad)->clear_user_data();
3603  (*quad)->clear_children();
3604  (*quad)->clear_refinement_case();
3605  (*quad)->clear_user_flag();
3606  (*quad)->clear_used_flag();
3607  }
3608  }
3609 
3610 
3628  template <int spacedim>
3629  static void
3632  unsigned int & next_unused_vertex,
3634  &next_unused_line,
3636  &next_unused_cell,
3637  const typename Triangulation<2, spacedim>::cell_iterator &cell)
3638  {
3639  const unsigned int dim = 2;
3640  // clear refinement flag
3641  const RefinementCase<dim> ref_case = cell->refine_flag_set();
3642  cell->clear_refine_flag();
3643 
3644  /* For the refinement process: since we go the levels up from the
3645  lowest, there are (unlike above) only two possibilities: a neighbor
3646  cell is on the same level or one level up (in both cases, it may or
3647  may not be refined later on, but we don't care here).
3648 
3649  First:
3650  Set up an array of the 3x3 vertices, which are distributed on the
3651  cell (the array consists of indices into the @p{vertices} std::vector
3652 
3653  2--7--3
3654  | | |
3655  4--8--5
3656  | | |
3657  0--6--1
3658 
3659  note: in case of cut_x or cut_y not all these vertices are needed for
3660  the new cells
3661 
3662  Second:
3663  Set up an array of the new lines (the array consists of iterator
3664  pointers into the lines arrays)
3665 
3666  .-6-.-7-. The directions are: .->-.->-.
3667  1 9 3 ^ ^ ^
3668  .-10.11-. .->-.->-.
3669  0 8 2 ^ ^ ^
3670  .-4-.-5-. .->-.->-.
3671 
3672  cut_x:
3673  .-4-.-5-.
3674  | | |
3675  0 6 1
3676  | | |
3677  .-2-.-3-.
3678 
3679  cut_y:
3680  .---5---.
3681  1 3
3682  .---6---.
3683  0 2
3684  .---4---.
3685 
3686 
3687  Third:
3688  Set up an array of neighbors:
3689 
3690  6 7
3691  .--.--.
3692  1| | |3
3693  .--.--.
3694  0| | |2
3695  .--.--.
3696  4 5
3697 
3698  We need this array for two reasons: first to get the lines which will
3699  bound the four subcells (if the neighboring cell is refined, these
3700  lines already exist), and second to update neighborship information.
3701  Since if a neighbor is not refined, its neighborship record only
3702  points to the present, unrefined, cell rather than the children we
3703  are presently creating, we only need the neighborship information
3704  if the neighbor cells are refined. In all other cases, we store
3705  the unrefined neighbor address
3706 
3707  We also need for every neighbor (if refined) which number among its
3708  neighbors the present (unrefined) cell has, since that number is to
3709  be replaced and because that also is the number of the subline which
3710  will be the interface between that neighbor and the to be created
3711  cell. We will store this number (between 0 and 3) in the field
3712  @p{neighbors_neighbor}.
3713 
3714  It would be sufficient to use the children of the common line to the
3715  neighbor, if we only wanted to get the new sublines and the new
3716  vertex, but because we need to update the neighborship information of
3717  the two refined subcells of the neighbor, we need to search these
3718  anyway.
3719 
3720  Convention:
3721  The created children are numbered like this:
3722 
3723  .--.--.
3724  |2 . 3|
3725  .--.--.
3726  |0 | 1|
3727  .--.--.
3728  */
3729  // collect the indices of the eight surrounding vertices
3730  // 2--7--3
3731  // | | |
3732  // 4--8--5
3733  // | | |
3734  // 0--6--1
3735  int new_vertices[9];
3736  for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3737  new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3738  for (unsigned int line_no = 0; line_no < 4; ++line_no)
3739  if (cell->line(line_no)->has_children())
3740  new_vertices[4 + line_no] =
3741  cell->line(line_no)->child(0)->vertex_index(1);
3742 
3743  if (ref_case == RefinementCase<dim>::cut_xy)
3744  {
3745  // find the next
3746  // unused vertex and
3747  // allocate it for
3748  // the new vertex we
3749  // need here
3750  while (triangulation.vertices_used[next_unused_vertex] == true)
3751  ++next_unused_vertex;
3752  Assert(next_unused_vertex < triangulation.vertices.size(),
3753  ExcMessage(
3754  "Internal error: During refinement, the triangulation "
3755  "wants to access an element of the 'vertices' array "
3756  "but it turns out that the array is not large enough."));
3757  triangulation.vertices_used[next_unused_vertex] = true;
3758 
3759  new_vertices[8] = next_unused_vertex;
3760 
3761  // determine middle vertex by transfinite interpolation to be
3762  // consistent with what happens to quads in a
3763  // Triangulation<3,3> when they are refined
3764  triangulation.vertices[next_unused_vertex] =
3765  cell->center(true, true);
3766  }
3767 
3768 
3769  // Now the lines:
3770  typename Triangulation<dim, spacedim>::raw_line_iterator new_lines[12];
3771  unsigned int lmin = 8;
3772  unsigned int lmax = 12;
3773  if (ref_case != RefinementCase<dim>::cut_xy)
3774  {
3775  lmin = 6;
3776  lmax = 7;
3777  }
3778 
3779  for (unsigned int l = lmin; l < lmax; ++l)
3780  {
3781  while (next_unused_line->used() == true)
3782  ++next_unused_line;
3783  new_lines[l] = next_unused_line;
3784  ++next_unused_line;
3785 
3786  AssertIsNotUsed(new_lines[l]);
3787  }
3788 
3789  if (ref_case == RefinementCase<dim>::cut_xy)
3790  {
3791  // .-6-.-7-.
3792  // 1 9 3
3793  // .-10.11-.
3794  // 0 8 2
3795  // .-4-.-5-.
3796 
3797  // lines 0-7 already exist, create only the four interior
3798  // lines 8-11
3799  unsigned int l = 0;
3800  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
3801  for (unsigned int c = 0; c < 2; ++c, ++l)
3802  new_lines[l] = cell->line(face_no)->child(c);
3803  Assert(l == 8, ExcInternalError());
3804 
3805  new_lines[8]->set_bounding_object_indices(
3806  {new_vertices[6], new_vertices[8]});
3807  new_lines[9]->set_bounding_object_indices(
3808  {new_vertices[8], new_vertices[7]});
3809  new_lines[10]->set_bounding_object_indices(
3810  {new_vertices[4], new_vertices[8]});
3811  new_lines[11]->set_bounding_object_indices(
3812  {new_vertices[8], new_vertices[5]});
3813  }
3814  else if (ref_case == RefinementCase<dim>::cut_x)
3815  {
3816  // .-4-.-5-.
3817  // | | |
3818  // 0 6 1
3819  // | | |
3820  // .-2-.-3-.
3821  new_lines[0] = cell->line(0);
3822  new_lines[1] = cell->line(1);
3823  new_lines[2] = cell->line(2)->child(0);
3824  new_lines[3] = cell->line(2)->child(1);
3825  new_lines[4] = cell->line(3)->child(0);
3826  new_lines[5] = cell->line(3)->child(1);
3827  new_lines[6]->set_bounding_object_indices(
3828  {new_vertices[6], new_vertices[7]});
3829  }
3830  else
3831  {
3833  // .---5---.
3834  // 1 3
3835  // .---6---.
3836  // 0 2
3837  // .---4---.
3838  new_lines[0] = cell->line(0)->child(0);
3839  new_lines[1] = cell->line(0)->child(1);
3840  new_lines[2] = cell->line(1)->child(0);
3841  new_lines[3] = cell->line(1)->child(1);
3842  new_lines[4] = cell->line(2);
3843  new_lines[5] = cell->line(3);
3844  new_lines[6]->set_bounding_object_indices(
3845  {new_vertices[4], new_vertices[5]});
3846  }
3847 
3848  for (unsigned int l = lmin; l < lmax; ++l)
3849  {
3850  new_lines[l]->set_used_flag();
3851  new_lines[l]->clear_user_flag();
3852  new_lines[l]->clear_user_data();
3853  new_lines[l]->clear_children();
3854  // interior line
3855  new_lines[l]->set_boundary_id_internal(
3857  new_lines[l]->set_manifold_id(cell->manifold_id());
3858  }
3859 
3860  // Now add the four (two)
3861  // new cells!
3864  while (next_unused_cell->used() == true)
3865  ++next_unused_cell;
3866 
3867  const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
3868  for (unsigned int i = 0; i < n_children; ++i)
3869  {
3870  AssertIsNotUsed(next_unused_cell);
3871  subcells[i] = next_unused_cell;
3872  ++next_unused_cell;
3873  if (i % 2 == 1 && i < n_children - 1)
3874  while (next_unused_cell->used() == true)
3875  ++next_unused_cell;
3876  }
3877 
3878  if (ref_case == RefinementCase<dim>::cut_xy)
3879  {
3880  // children:
3881  // .--.--.
3882  // |2 . 3|
3883  // .--.--.
3884  // |0 | 1|
3885  // .--.--.
3886  // lines:
3887  // .-6-.-7-.
3888  // 1 9 3
3889  // .-10.11-.
3890  // 0 8 2
3891  // .-4-.-5-.
3892  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3893  new_lines[8]->index(),
3894  new_lines[4]->index(),
3895  new_lines[10]->index()});
3896  subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3897  new_lines[2]->index(),
3898  new_lines[5]->index(),
3899  new_lines[11]->index()});
3900  subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3901  new_lines[9]->index(),
3902  new_lines[10]->index(),
3903  new_lines[6]->index()});
3904  subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3905  new_lines[3]->index(),
3906  new_lines[11]->index(),
3907  new_lines[7]->index()});
3908  }
3909  else if (ref_case == RefinementCase<dim>::cut_x)
3910  {
3911  // children:
3912  // .--.--.
3913  // | . |
3914  // .0 . 1.
3915  // | | |
3916  // .--.--.
3917  // lines:
3918  // .-4-.-5-.
3919  // | | |
3920  // 0 6 1
3921  // | | |
3922  // .-2-.-3-.
3923  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3924  new_lines[6]->index(),
3925  new_lines[2]->index(),
3926  new_lines[4]->index()});
3927  subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3928  new_lines[1]->index(),
3929  new_lines[3]->index(),
3930  new_lines[5]->index()});
3931  }
3932  else
3933  {
3935  // children:
3936  // .-----.
3937  // | 1 |
3938  // .-----.
3939  // | 0 |
3940  // .-----.
3941  // lines:
3942  // .---5---.
3943  // 1 3
3944  // .---6---.
3945  // 0 2
3946  // .---4---.
3947  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3948  new_lines[2]->index(),
3949  new_lines[4]->index(),
3950  new_lines[6]->index()});
3951  subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
3952  new_lines[3]->index(),
3953  new_lines[6]->index(),
3954  new_lines[5]->index()});
3955  }
3956 
3957  types::subdomain_id subdomainid = cell->subdomain_id();
3958 
3959  for (unsigned int i = 0; i < n_children; ++i)
3960  {
3961  subcells[i]->set_used_flag();
3962  subcells[i]->clear_refine_flag();
3963  subcells[i]->clear_user_flag();
3964  subcells[i]->clear_user_data();
3965  subcells[i]->clear_children();
3966  // inherit material properties
3967  subcells[i]->set_material_id(cell->material_id());
3968  subcells[i]->set_manifold_id(cell->manifold_id());
3969  subcells[i]->set_subdomain_id(subdomainid);
3970 
3971  if (i % 2 == 0)
3972  subcells[i]->set_parent(cell->index());
3973  }
3974 
3975 
3976 
3977  // set child index for even children i=0,2 (0)
3978  for (unsigned int i = 0; i < n_children / 2; ++i)
3979  cell->set_children(2 * i, subcells[2 * i]->index());
3980  // set the refine case
3981  cell->set_refinement_case(ref_case);
3982 
3983  // note that the
3984  // refinement flag was
3985  // already cleared at the
3986  // beginning of this function
3987 
3988  if (dim < spacedim)
3989  for (unsigned int c = 0; c < n_children; ++c)
3990  cell->child(c)->set_direction_flag(cell->direction_flag());
3991  }
3992 
3993 
3994 
3995  template <int dim, int spacedim>
3998  const bool check_for_distorted_cells)
3999  {
4000  AssertDimension(dim, 2);
4001 
4002  // Check whether a new level is needed. We have to check for
4003  // this on the highest level only
4004  for (const auto &cell : triangulation.active_cell_iterators_on_level(
4005  triangulation.levels.size() - 1))
4006  if (cell->refine_flag_set())
4007  {
4008  triangulation.levels.push_back(
4009  std::make_unique<
4011  break;
4012  }
4013 
4014  for (typename Triangulation<dim, spacedim>::line_iterator line =
4015  triangulation.begin_line();
4016  line != triangulation.end_line();
4017  ++line)
4018  {
4019  line->clear_user_flag();
4020  line->clear_user_data();
4021  }
4022 
4023  unsigned int n_single_lines = 0;
4024  unsigned int n_lines_in_pairs = 0;
4025  unsigned int needed_vertices = 0;
4026 
4027  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4028  {
4029  // count number of flagged cells on this level and compute
4030  // how many new vertices and new lines will be needed
4031  unsigned int needed_cells = 0;
4032 
4033  for (const auto &cell :
4034  triangulation.active_cell_iterators_on_level(level))
4035  if (cell->refine_flag_set())
4036  {
4037  if (cell->reference_cell() ==
4039  {
4040  needed_cells += 4;
4041  needed_vertices += 0;
4042  n_single_lines += 3;
4043  }
4044  else if (cell->reference_cell() ==
4046  {
4047  needed_cells += 4;
4048  needed_vertices += 1;
4049  n_single_lines += 4;
4050  }
4051  else
4052  {
4053  AssertThrow(false, ExcNotImplemented());
4054  }
4055 
4056  for (const auto line_no : cell->face_indices())
4057  {
4058  auto line = cell->line(line_no);
4059  if (line->has_children() == false)
4060  line->set_user_flag();
4061  }
4062  }
4063 
4064 
4065  const unsigned int used_cells =
4066  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4067  triangulation.levels[level + 1]->cells.used.end(),
4068  true);
4069 
4070 
4071  reserve_space(*triangulation.levels[level + 1],
4072  used_cells + needed_cells,
4073  2,
4074  spacedim);
4075 
4076  reserve_space(triangulation.levels[level + 1]->cells,
4077  needed_cells,
4078  0);
4079  }
4080 
4081  for (auto line = triangulation.begin_line();
4082  line != triangulation.end_line();
4083  ++line)
4084  if (line->user_flag_set())
4085  {
4086  Assert(line->has_children() == false, ExcInternalError());
4087  n_lines_in_pairs += 2;
4088  needed_vertices += 1;
4089  }
4090 
4091  reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4092 
4093  needed_vertices += std::count(triangulation.vertices_used.begin(),
4094  triangulation.vertices_used.end(),
4095  true);
4096 
4097  if (needed_vertices > triangulation.vertices.size())
4098  {
4099  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4100  triangulation.vertices_used.resize(needed_vertices, false);
4101  }
4102 
4103  unsigned int next_unused_vertex = 0;
4104 
4105  {
4107  line = triangulation.begin_active_line(),
4108  endl = triangulation.end_line();
4110  next_unused_line = triangulation.begin_raw_line();
4111 
4112  for (; line != endl; ++line)
4113  if (line->user_flag_set())
4114  {
4115  // this line needs to be refined
4116 
4117  // find the next unused vertex and set it
4118  // appropriately
4119  while (triangulation.vertices_used[next_unused_vertex] == true)
4120  ++next_unused_vertex;
4121  Assert(
4122  next_unused_vertex < triangulation.vertices.size(),
4123  ExcMessage(
4124  "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."));
4125  triangulation.vertices_used[next_unused_vertex] = true;
4126 
4127  triangulation.vertices[next_unused_vertex] = line->center(true);
4128 
4129  bool pair_found = false;
4130  (void)pair_found;
4131  for (; next_unused_line != endl; ++next_unused_line)
4132  if (!next_unused_line->used() &&
4133  !(++next_unused_line)->used())
4134  {
4135  --next_unused_line;
4136  pair_found = true;
4137  break;
4138  }
4139  Assert(pair_found, ExcInternalError());
4140 
4141  line->set_children(0, next_unused_line->index());
4142 
4144  children[2] = {next_unused_line, ++next_unused_line};
4145 
4146  AssertIsNotUsed(children[0]);
4147  AssertIsNotUsed(children[1]);
4148 
4149  children[0]->set_bounding_object_indices(
4150  {line->vertex_index(0), next_unused_vertex});
4151  children[1]->set_bounding_object_indices(
4152  {next_unused_vertex, line->vertex_index(1)});
4153 
4154  children[0]->set_used_flag();
4155  children[1]->set_used_flag();
4156  children[0]->clear_children();
4157  children[1]->clear_children();
4158  children[0]->clear_user_data();
4159  children[1]->clear_user_data();
4160  children[0]->clear_user_flag();
4161  children[1]->clear_user_flag();
4162 
4163 
4164  children[0]->set_boundary_id_internal(line->boundary_id());
4165  children[1]->set_boundary_id_internal(line->boundary_id());
4166 
4167  children[0]->set_manifold_id(line->manifold_id());
4168  children[1]->set_manifold_id(line->manifold_id());
4169 
4170  line->clear_user_flag();
4171  }
4172  }
4173 
4174  reserve_space(triangulation.faces->lines, 0, n_single_lines);
4175 
4177  cells_with_distorted_children;
4178 
4180  next_unused_line = triangulation.begin_raw_line();
4181 
4182  const auto create_children = [](auto & triangulation,
4183  unsigned int &next_unused_vertex,
4184  auto & next_unused_line,
4185  auto & next_unused_cell,
4186  const auto & cell) {
4187  const auto ref_case = cell->refine_flag_set();
4188  cell->clear_refine_flag();
4189 
4190  unsigned int n_new_vertices = 0;
4191 
4192  if (cell->reference_cell() == ::ReferenceCells::Triangle)
4193  n_new_vertices = 6;
4194  else if (cell->reference_cell() ==
4196  n_new_vertices = 9;
4197  else
4198  AssertThrow(false, ExcNotImplemented());
4199 
4200  std::vector<int> new_vertices(n_new_vertices);
4201  for (unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
4202  ++vertex_no)
4203  new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4204  for (unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
4205  if (cell->line(line_no)->has_children())
4206  new_vertices[cell->n_vertices() + line_no] =
4207  cell->line(line_no)->child(0)->vertex_index(1);
4208 
4209  if (cell->reference_cell() == ::ReferenceCells::Quadrilateral)
4210  {
4211  while (triangulation.vertices_used[next_unused_vertex] == true)
4212  ++next_unused_vertex;
4213  Assert(
4214  next_unused_vertex < triangulation.vertices.size(),
4215  ExcMessage(
4216  "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."));
4217  triangulation.vertices_used[next_unused_vertex] = true;
4218 
4219  new_vertices[8] = next_unused_vertex;
4220 
4221  triangulation.vertices[next_unused_vertex] =
4222  cell->center(true, true);
4223  }
4224 
4225  std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
4226  12>
4227  new_lines;
4228  unsigned int lmin = 0;
4229  unsigned int lmax = 0;
4230 
4231  if (cell->reference_cell() == ::ReferenceCells::Triangle)
4232  {
4233  lmin = 6;
4234  lmax = 9;
4235  }
4236  else if (cell->reference_cell() ==
4238  {
4239  lmin = 8;
4240  lmax = 12;
4241  }
4242  else
4243  {
4244  AssertThrow(false, ExcNotImplemented());
4245  }
4246 
4247  for (unsigned int l = lmin; l < lmax; ++l)
4248  {
4249  while (next_unused_line->used() == true)
4250  ++next_unused_line;
4251  new_lines[l] = next_unused_line;
4252  ++next_unused_line;
4253 
4254  AssertIsNotUsed(new_lines[l]);
4255  }
4256 
4257  if (true)
4258  {
4259  if (cell->reference_cell() == ::ReferenceCells::Triangle)
4260  {
4261  // add lines in the right order [TODO: clean up]
4262  const auto ref = [&](const unsigned int face_no,
4263  const unsigned int vertex_no) {
4264  if (cell->line(face_no)->child(0)->vertex_index(0) ==
4265  static_cast<unsigned int>(new_vertices[vertex_no]) ||
4266  cell->line(face_no)->child(0)->vertex_index(1) ==
4267  static_cast<unsigned int>(new_vertices[vertex_no]))
4268  {
4269  new_lines[2 * face_no + 0] =
4270  cell->line(face_no)->child(0);
4271  new_lines[2 * face_no + 1] =
4272  cell->line(face_no)->child(1);
4273  }
4274  else
4275  {
4276  new_lines[2 * face_no + 0] =
4277  cell->line(face_no)->child(1);
4278  new_lines[2 * face_no + 1] =
4279  cell->line(face_no)->child(0);
4280  }
4281  };
4282 
4283  ref(0, 0);
4284  ref(1, 1);
4285  ref(2, 2);
4286 
4287  new_lines[6]->set_bounding_object_indices(
4288  {new_vertices[3], new_vertices[4]});
4289  new_lines[7]->set_bounding_object_indices(
4290  {new_vertices[4], new_vertices[5]});
4291  new_lines[8]->set_bounding_object_indices(
4292  {new_vertices[5], new_vertices[3]});
4293  }
4294  else if (cell->reference_cell() ==
4296  {
4297  unsigned int l = 0;
4298  for (const unsigned int face_no : cell->face_indices())
4299  for (unsigned int c = 0; c < 2; ++c, ++l)
4300  new_lines[l] = cell->line(face_no)->child(c);
4301 
4302  new_lines[8]->set_bounding_object_indices(
4303  {new_vertices[6], new_vertices[8]});
4304  new_lines[9]->set_bounding_object_indices(
4305  {new_vertices[8], new_vertices[7]});
4306  new_lines[10]->set_bounding_object_indices(
4307  {new_vertices[4], new_vertices[8]});
4308  new_lines[11]->set_bounding_object_indices(
4309  {new_vertices[8], new_vertices[5]});
4310  }
4311  else
4312  {
4313  AssertThrow(false, ExcNotImplemented());
4314  }
4315  }
4316 
4317 
4318  for (unsigned int l = lmin; l < lmax; ++l)
4319  {
4320  new_lines[l]->set_used_flag();
4321  new_lines[l]->clear_user_flag();
4322  new_lines[l]->clear_user_data();
4323  new_lines[l]->clear_children();
4324  // interior line
4325  new_lines[l]->set_boundary_id_internal(
4327  new_lines[l]->set_manifold_id(cell->manifold_id());
4328  }
4329 
4332  while (next_unused_cell->used() == true)
4333  ++next_unused_cell;
4334 
4335  unsigned int n_children = 0;
4336 
4337  if (cell->reference_cell() == ::ReferenceCells::Triangle)
4338  n_children = 4;
4339  else if (cell->reference_cell() ==
4341  n_children = 4;
4342  else
4343  AssertThrow(false, ExcNotImplemented());
4344 
4345  for (unsigned int i = 0; i < n_children; ++i)
4346  {
4347  AssertIsNotUsed(next_unused_cell);
4348  subcells[i] = next_unused_cell;
4349  ++next_unused_cell;
4350  if (i % 2 == 1 && i < n_children - 1)
4351  while (next_unused_cell->used() == true)
4352  ++next_unused_cell;
4353  }
4354 
4355  if ((dim == 2) &&
4356  (cell->reference_cell() == ::ReferenceCells::Triangle))
4357  {
4358  subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4359  new_lines[8]->index(),
4360  new_lines[5]->index()});
4361  subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4362  new_lines[2]->index(),
4363  new_lines[6]->index()});
4364  subcells[2]->set_bounding_object_indices({new_lines[7]->index(),
4365  new_lines[3]->index(),
4366  new_lines[4]->index()});
4367  subcells[3]->set_bounding_object_indices({new_lines[6]->index(),
4368  new_lines[7]->index(),
4369  new_lines[8]->index()});
4370 
4371  // subcell 0
4372 
4373  const auto ref = [&](const unsigned int line_no,
4374  const unsigned int vertex_no,
4375  const unsigned int subcell_no,
4376  const unsigned int subcell_line_no) {
4377  if (new_lines[line_no]->vertex_index(1) !=
4378  static_cast<unsigned int>(new_vertices[vertex_no]))
4379  triangulation.levels[subcells[subcell_no]->level()]
4380  ->face_orientations[subcells[subcell_no]->index() *
4382  subcell_line_no] = 0;
4383  };
4384 
4385  ref(0, 3, 0, 0);
4386  ref(8, 5, 0, 1);
4387  ref(5, 0, 0, 2);
4388 
4389  ref(1, 1, 1, 0);
4390  ref(2, 4, 1, 1);
4391  ref(6, 3, 1, 2);
4392 
4393  ref(7, 4, 2, 0);
4394  ref(3, 2, 2, 1);
4395  ref(4, 5, 2, 2);
4396 
4397  ref(6, 4, 3, 0);
4398  ref(7, 5, 3, 1);
4399  ref(8, 3, 3, 2);
4400 
4401  // triangulation.levels[subcells[1]->level()]->face_orientations[subcells[1]->index()
4402  // * GeometryInfo<2>::faces_per_cell + 2] = 0;
4403  // triangulation.levels[subcells[2]->level()]->face_orientations[subcells[2]->index()
4404  // * GeometryInfo<2>::faces_per_cell + 0] = 0;
4405  }
4406  else if ((dim == 2) && (cell->reference_cell() ==
4408  {
4409  subcells[0]->set_bounding_object_indices(
4410  {new_lines[0]->index(),
4411  new_lines[8]->index(),
4412  new_lines[4]->index(),
4413  new_lines[10]->index()});
4414  subcells[1]->set_bounding_object_indices(
4415  {new_lines[8]->index(),
4416  new_lines[2]->index(),
4417  new_lines[5]->index(),
4418  new_lines[11]->index()});
4419  subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
4420  new_lines[9]->index(),
4421  new_lines[10]->index(),
4422  new_lines[6]->index()});
4423  subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
4424  new_lines[3]->index(),
4425  new_lines[11]->index(),
4426  new_lines[7]->index()});
4427  }
4428  else
4429  {
4430  AssertThrow(false, ExcNotImplemented());
4431  }
4432 
4433  types::subdomain_id subdomainid = cell->subdomain_id();
4434 
4435  for (unsigned int i = 0; i < n_children; ++i)
4436  {
4437  subcells[i]->set_used_flag();
4438  subcells[i]->clear_refine_flag();
4439  subcells[i]->clear_user_flag();
4440  subcells[i]->clear_user_data();
4441  subcells[i]->clear_children();
4442  // inherit material
4443  // properties
4444  subcells[i]->set_material_id(cell->material_id());
4445  subcells[i]->set_manifold_id(cell->manifold_id());
4446  subcells[i]->set_subdomain_id(subdomainid);
4447 
4448  // TODO: here we assume that all children have the same reference
4449  // cell type as the parent! This is justified for 2D.
4450  triangulation.levels[subcells[i]->level()]
4451  ->reference_cell[subcells[i]->index()] = cell->reference_cell();
4452 
4453  if (i % 2 == 0)
4454  subcells[i]->set_parent(cell->index());
4455  }
4456 
4457  for (unsigned int i = 0; i < n_children / 2; ++i)
4458  cell->set_children(2 * i, subcells[2 * i]->index());
4459 
4460  cell->set_refinement_case(ref_case);
4461 
4462  if (dim < spacedim)
4463  for (unsigned int c = 0; c < n_children; ++c)
4464  cell->child(c)->set_direction_flag(cell->direction_flag());
4465  };
4466 
4467  for (int level = 0;
4468  level < static_cast<int>(triangulation.levels.size()) - 1;
4469  ++level)
4470  {
4472  next_unused_cell = triangulation.begin_raw(level + 1);
4473 
4474  for (const auto &cell :
4475  triangulation.active_cell_iterators_on_level(level))
4476  if (cell->refine_flag_set())
4477  {
4479  next_unused_vertex,
4480  next_unused_line,
4481  next_unused_cell,
4482  cell);
4483 
4484  if (cell->reference_cell() ==
4486  check_for_distorted_cells &&
4487  has_distorted_children<dim, spacedim>(cell))
4488  cells_with_distorted_children.distorted_cells.push_back(
4489  cell);
4490 
4491  triangulation.signals.post_refinement_on_cell(cell);
4492  }
4493  }
4494 
4495  return cells_with_distorted_children;
4496  }
4497 
4498 
4499 
4504  template <int spacedim>
4507  const bool /*check_for_distorted_cells*/)
4508  {
4509  const unsigned int dim = 1;
4510 
4511  // Check whether a new level is needed. We have to check for
4512  // this on the highest level only
4513  for (const auto &cell : triangulation.active_cell_iterators_on_level(
4514  triangulation.levels.size() - 1))
4515  if (cell->refine_flag_set())
4516  {
4517  triangulation.levels.push_back(
4518  std::make_unique<
4520  break;
4521  }
4522 
4523 
4524  // check how much space is needed on every level. We need not
4525  // check the highest level since either - on the highest level
4526  // no cells are flagged for refinement - there are, but
4527  // prepare_refinement added another empty level
4528  unsigned int needed_vertices = 0;
4529  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4530  {
4531  // count number of flagged
4532  // cells on this level
4533  unsigned int flagged_cells = 0;
4534 
4535  for (const auto &acell :
4536  triangulation.active_cell_iterators_on_level(level))
4537  if (acell->refine_flag_set())
4538  ++flagged_cells;
4539 
4540  // count number of used cells
4541  // on the next higher level
4542  const unsigned int used_cells =
4543  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4544  triangulation.levels[level + 1]->cells.used.end(),
4545  true);
4546 
4547  // reserve space for the used_cells cells already existing
4548  // on the next higher level as well as for the
4549  // 2*flagged_cells that will be created on that level
4550  reserve_space(*triangulation.levels[level + 1],
4552  flagged_cells,
4553  1,
4554  spacedim);
4555  // reserve space for 2*flagged_cells new lines on the next
4556  // higher level
4557  reserve_space(triangulation.levels[level + 1]->cells,
4559  flagged_cells,
4560  0);
4561 
4562  needed_vertices += flagged_cells;
4563  }
4564 
4565  // add to needed vertices how many
4566  // vertices are already in use
4567  needed_vertices += std::count(triangulation.vertices_used.begin(),
4568  triangulation.vertices_used.end(),
4569  true);
4570  // if we need more vertices: create them, if not: leave the
4571  // array as is, since shrinking is not really possible because
4572  // some of the vertices at the end may be in use
4573  if (needed_vertices > triangulation.vertices.size())
4574  {
4575  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4576  triangulation.vertices_used.resize(needed_vertices, false);
4577  }
4578 
4579 
4580  // Do REFINEMENT on every level; exclude highest level as
4581  // above
4582 
4583  // index of next unused vertex
4584  unsigned int next_unused_vertex = 0;
4585 
4586  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4587  {
4589  next_unused_cell = triangulation.begin_raw(level + 1);
4590 
4591  for (const auto &cell :
4592  triangulation.active_cell_iterators_on_level(level))
4593  if (cell->refine_flag_set())
4594  {
4595  // clear refinement flag
4596  cell->clear_refine_flag();
4597 
4598  // search for next unused
4599  // vertex
4600  while (triangulation.vertices_used[next_unused_vertex] ==
4601  true)
4602  ++next_unused_vertex;
4603  Assert(
4604  next_unused_vertex < triangulation.vertices.size(),
4605  ExcMessage(
4606  "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."));
4607 
4608  // Now we always ask the cell itself where to put
4609  // the new point. The cell in turn will query the
4610  // manifold object internally.
4611  triangulation.vertices[next_unused_vertex] =
4612  cell->center(true);
4613 
4614  triangulation.vertices_used[next_unused_vertex] = true;
4615 
4616  // search for next two unused cell (++ takes care of
4617  // the end of the vector)
4619  first_child,
4620  second_child;
4621  while (next_unused_cell->used() == true)
4622  ++next_unused_cell;
4623  first_child = next_unused_cell;
4624  first_child->set_used_flag();
4625  first_child->clear_user_data();
4626  ++next_unused_cell;
4627  AssertIsNotUsed(next_unused_cell);
4628  second_child = next_unused_cell;
4629  second_child->set_used_flag();
4630  second_child->clear_user_data();
4631 
4632  types::subdomain_id subdomainid = cell->subdomain_id();
4633 
4634  // insert first child
4635  cell->set_children(0, first_child->index());
4636  first_child->clear_children();
4637  first_child->set_bounding_object_indices(
4638  {cell->vertex_index(0), next_unused_vertex});
4639  first_child->set_material_id(cell->material_id());
4640  first_child->set_manifold_id(cell->manifold_id());
4641  first_child->set_subdomain_id(subdomainid);
4642  first_child->set_direction_flag(cell->direction_flag());
4643 
4644  first_child->set_parent(cell->index());
4645 
4646  // Set manifold id of the right face. Only do this
4647  // on the first child.
4648  first_child->face(1)->set_manifold_id(cell->manifold_id());
4649 
4650  // reset neighborship info (refer to
4651  // internal::TriangulationImplementation::TriaLevel<0> for
4652  // details)
4653  first_child->set_neighbor(1, second_child);
4654  if (cell->neighbor(0).state() != IteratorState::valid)
4655  first_child->set_neighbor(0, cell->neighbor(0));
4656  else if (cell->neighbor(0)->is_active())
4657  {
4658  // since the neighbors level is always <=level,
4659  // if the cell is active, then there are no
4660  // cells to the left which may want to know
4661  // about this new child cell.
4662  Assert(cell->neighbor(0)->level() <= cell->level(),
4663  ExcInternalError());
4664  first_child->set_neighbor(0, cell->neighbor(0));
4665  }
4666  else
4667  // left neighbor is refined
4668  {
4669  // set neighbor to cell on same level
4670  const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4671  first_child->set_neighbor(0,
4672  cell->neighbor(0)->child(nbnb));
4673 
4674  // reset neighbor info of all right descendant
4675  // of the left neighbor of cell
4677  left_neighbor = cell->neighbor(0);
4678  while (left_neighbor->has_children())
4679  {
4680  left_neighbor = left_neighbor->child(nbnb);
4681  left_neighbor->set_neighbor(nbnb, first_child);
4682  }
4683  }
4684 
4685  // insert second child
4686  second_child->clear_children();
4687  second_child->set_bounding_object_indices(
4688  {next_unused_vertex, cell->vertex_index(1)});
4689  second_child->set_neighbor(0, first_child);
4690  second_child->set_material_id(cell->material_id());
4691  second_child->set_manifold_id(cell->manifold_id());
4692  second_child->set_subdomain_id(subdomainid);
4693  second_child->set_direction_flag(cell->direction_flag());
4694 
4695  if (cell->neighbor(1).state() != IteratorState::valid)
4696  second_child->set_neighbor(1, cell->neighbor(1));
4697  else if (cell->neighbor(1)->is_active())
4698  {
4699  Assert(cell->neighbor(1)->level() <= cell->level(),
4700  ExcInternalError());
4701  second_child->set_neighbor(1, cell->neighbor(1));
4702  }
4703  else
4704  // right neighbor is refined same as above
4705  {
4706  const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4707  second_child->set_neighbor(
4708  1, cell->neighbor(1)->child(nbnb));
4709 
4711  right_neighbor = cell->neighbor(1);
4712  while (right_neighbor->has_children())
4713  {
4714  right_neighbor = right_neighbor->child(nbnb);
4715  right_neighbor->set_neighbor(nbnb, second_child);
4716  }
4717  }
4718  // inform all listeners that cell refinement is done
4719  triangulation.signals.post_refinement_on_cell(cell);
4720  }
4721  }
4722 
4723  // in 1d, we can not have distorted children unless the parent
4724  // was already distorted (that is because we don't use
4725  // boundary information for 1d triangulations). so return an
4726  // empty list
4728  }
4729 
4730 
4735  template <int spacedim>
4738  const bool check_for_distorted_cells)
4739  {
4740  const unsigned int dim = 2;
4741 
4742 
4743  // First check whether we can get away with isotropic refinement, or
4744  // whether we need to run through the full anisotropic algorithm
4745  {
4746  bool do_isotropic_refinement = true;
4747  for (const auto &cell : triangulation.active_cell_iterators())
4748  if (cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
4749  cell->refine_flag_set() == RefinementCase<dim>::cut_y)
4750  {
4751  do_isotropic_refinement = false;
4752  break;
4753  }
4754 
4755  if (do_isotropic_refinement)
4757  check_for_distorted_cells);
4758  }
4759 
4760  // Check whether a new level is needed. We have to check for
4761  // this on the highest level only
4762  for (const auto &cell : triangulation.active_cell_iterators_on_level(
4763  triangulation.levels.size() - 1))
4764  if (cell->refine_flag_set())
4765  {
4766  triangulation.levels.push_back(
4767  std::make_unique<
4769  break;
4770  }
4771 
4772  // TODO[WB]: we clear user flags and pointers of lines; we're going
4773  // to use them to flag which lines need refinement
4774  for (typename Triangulation<dim, spacedim>::line_iterator line =
4775  triangulation.begin_line();
4776  line != triangulation.end_line();
4777  ++line)
4778  {
4779  line->clear_user_flag();
4780  line->clear_user_data();
4781  }
4782  // running over all cells and lines count the number
4783  // n_single_lines of lines which can be stored as single
4784  // lines, e.g. inner lines
4785  unsigned int n_single_lines = 0;
4786 
4787  // New lines to be created: number lines which are stored in
4788  // pairs (the children of lines must be stored in pairs)
4789  unsigned int n_lines_in_pairs = 0;
4790 
4791  // check how much space is needed on every level. We need not
4792  // check the highest level since either - on the highest level
4793  // no cells are flagged for refinement - there are, but
4794  // prepare_refinement added another empty level
4795  unsigned int needed_vertices = 0;
4796  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4797  {
4798  // count number of flagged cells on this level and compute
4799  // how many new vertices and new lines will be needed
4800  unsigned int needed_cells = 0;
4801 
4802  for (const auto &cell :
4803  triangulation.active_cell_iterators_on_level(level))
4804  if (cell->refine_flag_set())
4805  {
4806  if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
4807  {
4808  needed_cells += 4;
4809 
4810  // new vertex at center of cell is needed in any
4811  // case
4812  ++needed_vertices;
4813 
4814  // the four inner lines can be stored as singles
4815  n_single_lines += 4;
4816  }
4817  else // cut_x || cut_y
4818  {
4819  // set the flag showing that anisotropic
4820  // refinement is used for at least one cell
4821  triangulation.anisotropic_refinement = true;
4822 
4823  needed_cells += 2;
4824  // no vertex at center
4825 
4826  // the inner line can be stored as single
4827  n_single_lines += 1;
4828  }
4829 
4830  // mark all faces (lines) for refinement; checking
4831  // locally whether the neighbor would also like to
4832  // refine them is rather difficult for lines so we
4833  // only flag them and after visiting all cells, we
4834  // decide which lines need refinement;
4835  for (const unsigned int line_no :
4837  {
4839  cell->refine_flag_set(), line_no) ==
4841  {
4843  line = cell->line(line_no);
4844  if (line->has_children() == false)
4845  line->set_user_flag();
4846  }
4847  }
4848  }
4849 
4850 
4851  // count number of used cells on the next higher level
4852  const unsigned int used_cells =
4853  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4854  triangulation.levels[level + 1]->cells.used.end(),
4855  true);
4856 
4857 
4858  // reserve space for the used_cells cells already existing
4859  // on the next higher level as well as for the
4860  // needed_cells that will be created on that level
4861  reserve_space(*triangulation.levels[level + 1],
4862  used_cells + needed_cells,
4863  2,
4864  spacedim);
4865 
4866  // reserve space for needed_cells new quads on the next
4867  // higher level
4868  reserve_space(triangulation.levels[level + 1]->cells,
4869  needed_cells,
4870  0);
4871  }
4872 
4873  // now count the lines which were flagged for refinement
4874  for (typename Triangulation<dim, spacedim>::line_iterator line =
4875  triangulation.begin_line();
4876  line != triangulation.end_line();
4877  ++line)
4878  if (line->user_flag_set())
4879  {
4880  Assert(line->has_children() == false, ExcInternalError());
4881  n_lines_in_pairs += 2;
4882  needed_vertices += 1;
4883  }
4884  // reserve space for n_lines_in_pairs new lines. note, that
4885  // we can't reserve space for the single lines here as well,
4886  // as all the space reserved for lines in pairs would be
4887  // counted as unused and we would end up with too little space
4888  // to store all lines. memory reservation for n_single_lines
4889  // can only be done AFTER we refined the lines of the current
4890  // cells
4891  reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4892 
4893  // add to needed vertices how many vertices are already in use
4894  needed_vertices += std::count(triangulation.vertices_used.begin(),
4895  triangulation.vertices_used.end(),
4896  true);
4897  // if we need more vertices: create them, if not: leave the
4898  // array as is, since shrinking is not really possible because
4899  // some of the vertices at the end may be in use
4900  if (needed_vertices > triangulation.vertices.size())
4901  {
4902  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4903  triangulation.vertices_used.resize(needed_vertices, false);
4904  }
4905 
4906 
4907  // Do REFINEMENT on every level; exclude highest level as
4908  // above
4909 
4910  // index of next unused vertex
4911  unsigned int next_unused_vertex = 0;
4912 
4913  // first the refinement of lines. children are stored
4914  // pairwise
4915  {
4916  // only active objects can be refined further
4918  line = triangulation.begin_active_line(),
4919  endl = triangulation.end_line();
4921  next_unused_line = triangulation.begin_raw_line();
4922 
4923  for (; line != endl; ++line)
4924  if (line->user_flag_set())
4925  {
4926  // this line needs to be refined
4927 
4928  // find the next unused vertex and set it
4929  // appropriately
4930  while (triangulation.vertices_used[next_unused_vertex] == true)
4931  ++next_unused_vertex;
4932  Assert(
4933  next_unused_vertex < triangulation.vertices.size(),
4934  ExcMessage(
4935  "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."));
4936  triangulation.vertices_used[next_unused_vertex] = true;
4937 
4938  triangulation.vertices[next_unused_vertex] = line->center(true);
4939 
4940  // now that we created the right point, make up the
4941  // two child lines. To this end, find a pair of
4942  // unused lines
4943  bool pair_found = false;
4944  (void)pair_found;
4945  for (; next_unused_line != endl; ++next_unused_line)
4946  if (!next_unused_line->used() &&
4947  !(++next_unused_line)->used())
4948  {
4949  // go back to the first of the two unused
4950  // lines
4951  --next_unused_line;
4952  pair_found = true;
4953  break;
4954  }
4955  Assert(pair_found, ExcInternalError());
4956 
4957  // there are now two consecutive unused lines, such
4958  // that the children of a line will be consecutive.
4959  // then set the child pointer of the present line
4960  line->set_children(0, next_unused_line->index());
4961 
4962  // set the two new lines
4964  children[2] = {next_unused_line, ++next_unused_line};
4965  // some tests; if any of the iterators should be
4966  // invalid, then already dereferencing will fail
4967  AssertIsNotUsed(children[0]);
4968  AssertIsNotUsed(children[1]);
4969 
4970  children[0]->set_bounding_object_indices(
4971  {line->vertex_index(0), next_unused_vertex});
4972  children[1]->set_bounding_object_indices(
4973  {next_unused_vertex, line->vertex_index(1)});
4974 
4975  children[0]->set_used_flag();
4976  children[1]->set_used_flag();
4977  children[0]->clear_children();
4978  children[1]->clear_children();
4979  children[0]->clear_user_data();
4980  children[1]->clear_user_data();
4981  children[0]->clear_user_flag();
4982  children[1]->clear_user_flag();
4983 
4984 
4985  children[0]->set_boundary_id_internal(line->boundary_id());
4986  children[1]->set_boundary_id_internal(line->boundary_id());
4987 
4988  children[0]->set_manifold_id(line->manifold_id());
4989  children[1]->set_manifold_id(line->manifold_id());
4990 
4991  // finally clear flag indicating the need for
4992  // refinement
4993  line->clear_user_flag();
4994  }
4995  }
4996 
4997 
4998  // Now set up the new cells
4999 
5000  // reserve space for inner lines (can be stored as single
5001  // lines)
5002  reserve_space(triangulation.faces->lines, 0, n_single_lines);
5003 
5005  cells_with_distorted_children;
5006 
5007  // reset next_unused_line, as now also single empty places in
5008  // the vector can be used
5010  next_unused_line = triangulation.begin_raw_line();
5011 
5012  for (int level = 0;
5013  level < static_cast<int>(triangulation.levels.size()) - 1;
5014  ++level)
5015  {
5017  next_unused_cell = triangulation.begin_raw(level + 1);
5018 
5019  for (const auto &cell :
5020  triangulation.active_cell_iterators_on_level(level))
5021  if (cell->refine_flag_set())
5022  {
5023  // actually set up the children and update neighbor
5024  // information
5026  next_unused_vertex,
5027  next_unused_line,
5028  next_unused_cell,
5029  cell);
5030 
5031  if (check_for_distorted_cells &&
5032  has_distorted_children<dim, spacedim>(cell))
5033  cells_with_distorted_children.distorted_cells.push_back(
5034  cell);
5035  // inform all listeners that cell refinement is done
5036  triangulation.signals.post_refinement_on_cell(cell);
5037  }
5038  }
5039 
5040  return cells_with_distorted_children;
5041  }
5042 
5043 
5044  template <int spacedim>
5047  const bool check_for_distorted_cells)
5048  {
5049  static const int dim = 3;
5050  static const unsigned int X = numbers::invalid_unsigned_int;
5051 
5052  Assert(spacedim == 3, ExcNotImplemented());
5053 
5054  Assert(triangulation.vertices.size() ==
5055  triangulation.vertices_used.size(),
5056  ExcInternalError());
5057 
5058  // Check whether a new level is needed. We have to check for
5059  // this on the highest level only
5060  for (const auto &cell : triangulation.active_cell_iterators_on_level(
5061  triangulation.levels.size() - 1))
5062  if (cell->refine_flag_set())
5063  {
5064  triangulation.levels.push_back(
5065  std::make_unique<
5067  break;
5068  }
5069 
5070  // first clear user flags for quads and lines; we're going to
5071  // use them to flag which lines and quads need refinement
5072  triangulation.faces->quads.clear_user_data();
5073 
5074  for (typename Triangulation<dim, spacedim>::line_iterator line =
5075  triangulation.begin_line();
5076  line != triangulation.end_line();
5077  ++line)
5078  line->clear_user_flag();
5079 
5080  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
5081  triangulation.begin_quad();
5082  quad != triangulation.end_quad();
5083  ++quad)
5084  quad->clear_user_flag();
5085 
5086  // check how much space is needed on every level. We need not
5087  // check the highest level since either
5088  // - on the highest level no cells are flagged for refinement
5089  // - there are, but prepare_refinement added another empty
5090  // level which then is the highest level
5091 
5092  // variables to hold the number of newly to be created
5093  // vertices, lines and quads. as these are stored globally,
5094  // declare them outside the loop over al levels. we need lines
5095  // and quads in pairs for refinement of old ones and lines and
5096  // quads, that can be stored as single ones, as they are newly
5097  // created in the inside of an existing cell
5098  unsigned int needed_vertices = 0;
5099  unsigned int needed_lines_single = 0;
5100  unsigned int needed_quads_single = 0;
5101  unsigned int needed_lines_pair = 0;
5102  unsigned int needed_quads_pair = 0;
5103  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5104  {
5105  unsigned int new_cells = 0;
5106 
5107  for (const auto &cell :
5108  triangulation.active_cell_iterators_on_level(level))
5109  if (cell->refine_flag_set())
5110  {
5111  // Only support isotropic refinement
5112  Assert(cell->refine_flag_set() ==
5114  ExcInternalError());
5115 
5116  // Now count up how many new cells, faces, edges, and vertices
5117  // we will need to allocate to do this refinement.
5118  new_cells += cell->reference_cell().n_isotropic_children();
5119 
5120  if (cell->reference_cell() == ReferenceCells::Hexahedron)
5121  {
5122  ++needed_vertices;
5123  needed_lines_single += 6;
5124  needed_quads_single += 12;
5125  }
5126  else if (cell->reference_cell() ==
5128  {
5129  needed_lines_single += 1;
5130  needed_quads_single += 8;
5131  }
5132  else
5133  {
5134  Assert(false, ExcInternalError());
5135  }
5136 
5137  // Also check whether we have to refine any of the faces and
5138  // edges that bound this cell. They may of course already be
5139  // refined, so we only *mark* them for refinement by setting
5140  // the user flags
5141  for (const auto face : cell->face_indices())
5142  if (cell->face(face)->n_children() == 0)
5143  cell->face(face)->set_user_flag();
5144  else
5145  Assert(cell->face(face)->n_children() ==
5146  cell->reference_cell()
5147  .face_reference_cell(face)
5148  .n_isotropic_children(),
5149  ExcInternalError());
5150 
5151  for (const auto line : cell->line_indices())
5152  if (cell->line(line)->has_children() == false)
5153  cell->line(line)->set_user_flag();
5154  else
5155  Assert(cell->line(line)->n_children() == 2,
5156  ExcInternalError());
5157  }
5158 
5159  const unsigned int used_cells =
5160  std::count(triangulation.levels[level + 1]->cells.used.begin(),
5161  triangulation.levels[level + 1]->cells.used.end(),
5162  true);
5163 
5164  reserve_space(*triangulation.levels[level + 1],
5165  used_cells + new_cells,
5166  3,
5167  spacedim);
5168 
5169  reserve_space(triangulation.levels[level + 1]->cells, new_cells);
5170  }
5171 
5172  // now count the quads and lines which were flagged for
5173  // refinement
5174  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
5175  triangulation.begin_quad();
5176  quad != triangulation.end_quad();
5177  ++quad)
5178  {
5179  if (quad->user_flag_set() == false)
5180  continue;
5181 
5182  if (quad->reference_cell() == ReferenceCells::Quadrilateral)
5183  {
5184  needed_quads_pair += 4;
5185  needed_lines_pair += 4;
5186  needed_vertices += 1;
5187  }
5188  else if (quad->reference_cell() == ReferenceCells::Triangle)
5189  {
5190  needed_quads_pair += 4;
5191  needed_lines_single += 3;
5192  }
5193  else
5194  {
5195  Assert(false, ExcInternalError());
5196  }
5197  }
5198 
5199  for (typename Triangulation<dim, spacedim>::line_iterator line =
5200  triangulation.begin_line();
5201  line != triangulation.end_line();
5202  ++line)
5203  {
5204  if (line->user_flag_set() == false)
5205  continue;
5206 
5207  needed_lines_pair += 2;
5208  needed_vertices += 1;
5209  }
5210 
5211  reserve_space(triangulation.faces->lines,
5212  needed_lines_pair,
5213  needed_lines_single);
5215  needed_quads_pair,
5216  needed_quads_single);
5217  reserve_space(triangulation.faces->quads,
5218  needed_quads_pair,
5219  needed_quads_single);
5220 
5221 
5222  // add to needed vertices how many vertices are already in use
5223  needed_vertices += std::count(triangulation.vertices_used.begin(),
5224  triangulation.vertices_used.end(),
5225  true);
5226 
5227  if (needed_vertices > triangulation.vertices.size())
5228  {
5229  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5230  triangulation.vertices_used.resize(needed_vertices, false);
5231  }
5232 
5233  //-----------------------------------------
5234  // Before we start with the actual refinement, we do some
5235  // sanity checks if in debug mode. especially, we try to catch
5236  // the notorious problem with lines being twice refined,
5237  // i.e. there are cells adjacent at one line ("around the
5238  // edge", but not at a face), with two cells differing by more
5239  // than one refinement level
5240  //
5241  // this check is very simple to implement here, since we have
5242  // all lines flagged if they shall be refined
5243 #ifdef DEBUG
5244  for (const auto &cell : triangulation.active_cell_iterators())
5245  if (!cell->refine_flag_set())
5246  for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
5247  if (cell->line(line_n)->has_children())
5248  for (unsigned int c = 0; c < 2; ++c)
5249  Assert(cell->line(line_n)->child(c)->user_flag_set() == false,
5250  ExcInternalError());
5251 #endif
5252 
5253  unsigned int current_vertex = 0;
5254 
5255  // helper function - find the next available vertex number and mark it
5256  // as used.
5257  auto get_next_unused_vertex = [](const unsigned int current_vertex,
5258  std::vector<bool> &vertices_used) {
5259  unsigned int next_vertex = current_vertex;
5260  while (next_vertex < vertices_used.size() &&
5261  vertices_used[next_vertex] == true)
5262  ++next_vertex;
5263  Assert(next_vertex < vertices_used.size(), ExcInternalError());
5264  vertices_used[next_vertex] = true;
5265 
5266  return next_vertex;
5267  };
5268 
5269  // LINES
5270  {
5272  line = triangulation.begin_active_line(),
5273  endl = triangulation.end_line();
5275  next_unused_line = triangulation.begin_raw_line();
5276 
5277  for (; line != endl; ++line)
5278  {
5279  if (line->user_flag_set() == false)
5280  continue;
5281 
5282  current_vertex =
5283  get_next_unused_vertex(current_vertex,
5284  triangulation.vertices_used);
5285  triangulation.vertices[current_vertex] = line->center(true);
5286 
5287  next_unused_line =
5288  triangulation.faces->lines.template next_free_pair_object<1>(
5289  triangulation);
5290  Assert(next_unused_line.state() == IteratorState::valid,
5291  ExcInternalError());
5292 
5293  // now we found two consecutive unused lines, such
5294  // that the children of a line will be consecutive.
5295  // then set the child pointer of the present line
5296  line->set_children(0, next_unused_line->index());
5297 
5299  children[2] = {next_unused_line, ++next_unused_line};
5300 
5301  AssertIsNotUsed(children[0]);
5302  AssertIsNotUsed(children[1]);
5303 
5304  children[0]->set_bounding_object_indices(
5305  {line->vertex_index(0), current_vertex});
5306  children[1]->set_bounding_object_indices(
5307  {current_vertex, line->vertex_index(1)});
5308 
5309  children[0]->set_used_flag();
5310  children[1]->set_used_flag();
5311  children[0]->clear_children();
5312  children[1]->clear_children();
5313  children[0]->clear_user_data();
5314  children[1]->clear_user_data();
5315  children[0]->clear_user_flag();
5316  children[1]->clear_user_flag();
5317 
5318  children[0]->set_boundary_id_internal(line->boundary_id());
5319  children[1]->set_boundary_id_internal(line->boundary_id());
5320 
5321  children[0]->set_manifold_id(line->manifold_id());
5322  children[1]->set_manifold_id(line->manifold_id());
5323 
5324  line->clear_user_flag();
5325  }
5326  }
5327 
5328  // QUADS
5329  {
5331  quad = triangulation.begin_quad(),
5332  endq = triangulation.end_quad();
5334  next_unused_line = triangulation.begin_raw_line();
5336  next_unused_quad = triangulation.begin_raw_quad();
5337 
5338  for (; quad != endq; ++quad)
5339  {
5340  if (quad->user_flag_set() == false)
5341  continue;
5342 
5343  const auto reference_face_type = quad->reference_cell();
5344 
5345  // 1) create new vertex (at the center of the face)
5346  if (reference_face_type == ReferenceCells::Quadrilateral)
5347  {
5348  current_vertex =
5349  get_next_unused_vertex(current_vertex,
5350  triangulation.vertices_used);
5351  triangulation.vertices[current_vertex] =
5352  quad->center(true, true);
5353  }
5354 
5355  // 2) create new lines (property is set later)
5356  boost::container::small_vector<
5359  new_lines(quad->n_lines());
5360  {
5361  for (unsigned int i = 0; i < new_lines.size(); ++i)
5362  {
5363  if (reference_face_type == ReferenceCells::Quadrilateral)
5364  {
5365  if (i % 2 == 0)
5366  next_unused_line =
5367  triangulation.faces->lines
5368  .template next_free_pair_object<1>(triangulation);
5369  }
5370  else if (reference_face_type == ReferenceCells::Triangle)
5371  {
5372  next_unused_line =
5373  triangulation.faces->lines
5374  .template next_free_single_object<1>(triangulation);
5375  }
5376  else
5377  {
5378  Assert(false, ExcNotImplemented());
5379  }
5380 
5381  new_lines[i] = next_unused_line;
5382  ++next_unused_line;
5383  AssertIsNotUsed(new_lines[i]);
5384  }
5385  }
5386 
5387  // 3) create new quads (properties are set below). Both triangles
5388  // and quads are divided in four.
5389  std::array<
5391  4>
5392  new_quads;
5393  {
5394  next_unused_quad =
5395  triangulation.faces->quads.template next_free_pair_object<2>(
5396  triangulation);
5397 
5398  new_quads[0] = next_unused_quad;
5399  AssertIsNotUsed(new_quads[0]);
5400 
5401  ++next_unused_quad;
5402  new_quads[1] = next_unused_quad;
5403  AssertIsNotUsed(new_quads[1]);
5404 
5405  next_unused_quad =
5406  triangulation.faces->quads.template next_free_pair_object<2>(
5407  triangulation);
5408  new_quads[2] = next_unused_quad;
5409  AssertIsNotUsed(new_quads[2]);
5410 
5411  ++next_unused_quad;
5412  new_quads[3] = next_unused_quad;
5413  AssertIsNotUsed(new_quads[3]);
5414 
5415  quad->set_children(0, new_quads[0]->index());
5416  quad->set_children(2, new_quads[2]->index());
5417  quad->set_refinement_case(RefinementCase<2>::cut_xy);
5418  }
5419 
5420  // Maximum of 9 vertices per refined quad (9 for Quadrilateral, 6
5421  // for Triangle)
5422  std::array<unsigned int, 9> vertex_indices = {};
5423  {
5424  unsigned int k = 0;
5425  for (const auto i : quad->vertex_indices())
5426  vertex_indices[k++] = quad->vertex_index(i);
5427 
5428  for (const auto i : quad->line_indices())
5429  vertex_indices[k++] =
5430  quad->line(i)->child(0)->vertex_index(1);
5431 
5432  vertex_indices[k++] = current_vertex;
5433  }
5434 
5435  boost::container::small_vector<
5437  12>
5438  lines(reference_face_type == ReferenceCells::Quadrilateral ?
5439  12 :
5440  9);
5441  {
5442  unsigned int k = 0;
5443 
5444  for (unsigned int l = 0; l < quad->n_lines(); ++l)
5445  for (unsigned int c = 0; c < 2; ++c)
5446  {
5447  static constexpr std::array<std::array<unsigned int, 2>,
5448  2>
5449  index = {// child 0, line_orientation=false and true
5450  {{{1, 0}},
5451  // child 1, line_orientation=false and true
5452  {{0, 1}}}};
5453 
5454  lines[k++] = quad->line(l)->child(
5455  index[c][quad->line_orientation(l)]);
5456  }
5457 
5458  for (unsigned int l = 0; l < new_lines.size(); ++l)
5459  lines[k++] = new_lines[l];
5460  }
5461 
5462  boost::container::small_vector<int, 12> line_indices(
5463  lines.size());
5464  for (unsigned int i = 0; i < line_indices.size(); ++i)
5465  line_indices[i] = lines[i]->index();
5466 
5467  static constexpr std::array<std::array<unsigned int, 2>, 12>
5468  line_vertices_quad{{{{0, 4}},
5469  {{4, 2}},
5470  {{1, 5}},
5471  {{5, 3}},
5472  {{0, 6}},
5473  {{6, 1}},
5474  {{2, 7}},
5475  {{7, 3}},
5476  {{6, 8}},
5477  {{8, 7}},
5478  {{4, 8}},
5479  {{8, 5}}}};
5480 
5481  static constexpr std::array<std::array<unsigned int, 4>, 4>
5482  quad_lines_quad{{{{0, 8, 4, 10}},
5483  {{8, 2, 5, 11}},
5484  {{1, 9, 10, 6}},
5485  {{9, 3, 11, 7}}}};
5486 
5487  static constexpr std::
5488  array<std::array<std::array<unsigned int, 2>, 4>, 4>
5489  quad_line_vertices_quad{
5490  {{{{{0, 4}}, {{6, 8}}, {{0, 6}}, {{4, 8}}}},
5491  {{{{6, 8}}, {{1, 5}}, {{6, 1}}, {{8, 5}}}},
5492  {{{{4, 2}}, {{8, 7}}, {{4, 8}}, {{2, 7}}}},
5493  {{{{8, 7}}, {{5, 3}}, {{8, 5}}, {{7, 3}}}}}};
5494 
5495  static constexpr std::array<std::array<unsigned int, 2>, 12>
5496  line_vertices_tri{{{{0, 3}},
5497  {{3, 1}},
5498  {{1, 4}},
5499  {{4, 2}},
5500  {{2, 5}},
5501  {{5, 0}},
5502  {{3, 4}},
5503  {{4, 5}},
5504  {{3, 5}},
5505  {{X, X}},
5506  {{X, X}},
5507  {{X, X}}}};
5508 
5509  static constexpr std::array<std::array<unsigned int, 4>, 4>
5510  quad_lines_tri{{{{0, 8, 5, X}},
5511  {{1, 2, 6, X}},
5512  {{7, 3, 4, X}},
5513  {{6, 7, 8, X}}}};
5514 
5515  static constexpr std::
5516  array<std::array<std::array<unsigned int, 2>, 4>, 4>
5517  quad_line_vertices_tri{
5518  {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
5519  {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
5520  {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
5521  {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
5522 
5523  const auto &line_vertices =
5524  (reference_face_type == ReferenceCells::Quadrilateral) ?
5525  line_vertices_quad :
5526  line_vertices_tri;
5527  const auto &quad_lines =
5528  (reference_face_type == ReferenceCells::Quadrilateral) ?
5529  quad_lines_quad :
5530  quad_lines_tri;
5531  const auto &quad_line_vertices =
5532  (reference_face_type == ReferenceCells::Quadrilateral) ?
5533  quad_line_vertices_quad :
5534  quad_line_vertices_tri;
5535 
5536  // 4) set properties of lines
5537  for (unsigned int i = 0, j = lines.size() - new_lines.size();
5538  i < new_lines.size();
5539  ++i, ++j)
5540  {
5541  auto &new_line = new_lines[i];
5542  new_line->set_bounding_object_indices(
5543  {vertex_indices[line_vertices[j][0]],
5544  vertex_indices[line_vertices[j][1]]});
5545  new_line->set_used_flag();
5546  new_line->clear_user_flag();
5547  new_line->clear_user_data();
5548  new_line->clear_children();
5549  new_line->set_boundary_id_internal(quad->boundary_id());
5550  new_line->set_manifold_id(quad->manifold_id());
5551  }
5552 
5553  // 5) set properties of quads
5554  for (unsigned int i = 0; i < new_quads.size(); ++i)
5555  {
5556  auto &new_quad = new_quads[i];
5557 
5558  // TODO: we assume here that all children have the same type
5559  // as the parent
5560  triangulation.faces->quad_reference_cell[new_quad->index()] =
5561  reference_face_type;
5562 
5563  if (new_quad->n_lines() == 3)
5564  new_quad->set_bounding_object_indices(
5565  {line_indices[quad_lines[i][0]],
5566  line_indices[quad_lines[i][1]],
5567  line_indices[quad_lines[i][2]]});
5568  else if (new_quad->n_lines() == 4)
5569  new_quad->set_bounding_object_indices(
5570  {line_indices[quad_lines[i][0]],
5571  line_indices[quad_lines[i][1]],
5572  line_indices[quad_lines[i][2]],
5573  line_indices[quad_lines[i][3]]});
5574  else
5575  Assert(false, ExcNotImplemented());
5576 
5577  new_quad->set_used_flag();
5578  new_quad->clear_user_flag();
5579  new_quad->clear_user_data();
5580  new_quad->clear_children();
5581  new_quad->set_boundary_id_internal(quad->boundary_id());
5582  new_quad->set_manifold_id(quad->manifold_id());
5583 
5584 #ifdef DEBUG
5585  std::set<unsigned int> s;
5586 #endif
5587 
5588  // ... and fix orientation of faces (lines) of quad
5589  for (const auto f : new_quad->line_indices())
5590  {
5591  std::array<unsigned int, 2> vertices_0, vertices_1;
5592 
5593  for (unsigned int v = 0; v < 2; ++v)
5594  vertices_0[v] =
5595  lines[quad_lines[i][f]]->vertex_index(v);
5596 
5597  for (unsigned int v = 0; v < 2; ++v)
5598  vertices_1[v] =
5599  vertex_indices[quad_line_vertices[i][f][v]];
5600 
5601  const auto orientation =
5603  vertices_1);
5604 
5605 #ifdef DEBUG
5606  for (const auto i : vertices_0)
5607  s.insert(i);
5608  for (const auto i : vertices_1)
5609  s.insert(i);
5610 #endif
5611 
5612  new_quad->set_line_orientation(f, orientation);
5613  }
5614 #ifdef DEBUG
5616  s.size(),
5617  (reference_face_type == ReferenceCells::Quadrilateral ? 4 :
5618  3));
5619 #endif
5620  }
5621 
5622  quad->clear_user_flag();
5623  }
5624  }
5625 
5627  cells_with_distorted_children;
5628 
5629  for (unsigned int level = 0; level != triangulation.levels.size() - 1;
5630  ++level)
5631  {
5633  hex = triangulation.begin_active_hex(level),
5634  endh = triangulation.begin_active_hex(level + 1);
5636  next_unused_hex = triangulation.begin_raw_hex(level + 1);
5637 
5638  for (; hex != endh; ++hex)
5639  {
5640  if (hex->refine_flag_set() ==
5642  continue;
5643 
5644  const auto &reference_cell_type = hex->reference_cell();
5645 
5646  const RefinementCase<dim> ref_case = hex->refine_flag_set();
5647  hex->clear_refine_flag();
5648  hex->set_refinement_case(ref_case);
5649 
5650  unsigned int n_new_lines = 0;
5651  unsigned int n_new_quads = 0;
5652  unsigned int n_new_hexes = 0;
5653 
5654  if (reference_cell_type == ReferenceCells::Hexahedron)
5655  {
5656  n_new_lines = 6;
5657  n_new_quads = 12;
5658  n_new_hexes = 8;
5659  }
5660  else if (reference_cell_type == ReferenceCells::Tetrahedron)
5661  {
5662  n_new_lines = 1;
5663  n_new_quads = 8;
5664  n_new_hexes = 8;
5665  }
5666  else
5667  Assert(false, ExcNotImplemented());
5668 
5669  // Hexes add a single new internal vertex
5670  if (reference_cell_type == ReferenceCells::Hexahedron)
5671  {
5672  current_vertex =
5673  get_next_unused_vertex(current_vertex,
5674  triangulation.vertices_used);
5675  triangulation.vertices[current_vertex] =
5676  hex->center(true, true);
5677  }
5678 
5679  boost::container::small_vector<
5681  6>
5682  new_lines(n_new_lines);
5683  for (unsigned int i = 0; i < n_new_lines; ++i)
5684  {
5685  new_lines[i] =
5686  triangulation.faces->lines
5687  .template next_free_single_object<1>(triangulation);
5688 
5689  AssertIsNotUsed(new_lines[i]);
5690  new_lines[i]->set_used_flag();
5691  new_lines[i]->clear_user_flag();
5692  new_lines[i]->clear_user_data();
5693  new_lines[i]->clear_children();
5694  new_lines[i]->set_boundary_id_internal(
5696  new_lines[i]->set_manifold_id(hex->manifold_id());
5697  }
5698 
5699  boost::container::small_vector<
5701  12>
5702  new_quads(n_new_quads);
5703  for (unsigned int i = 0; i < n_new_quads; ++i)
5704  {
5705  new_quads[i] =
5706  triangulation.faces->quads
5707  .template next_free_single_object<2>(triangulation);
5708 
5709  auto &new_quad = new_quads[i];
5710 
5711  // TODO: faces of children have the same type as the faces
5712  // of the parent
5713  triangulation.faces
5714  ->quad_reference_cell[new_quad->index()] =
5715  (reference_cell_type == ReferenceCells::Hexahedron) ?
5718 
5719  AssertIsNotUsed(new_quad);
5720  new_quad->set_used_flag();
5721  new_quad->clear_user_flag();
5722  new_quad->clear_user_data();
5723  new_quad->clear_children();
5724  new_quad->set_boundary_id_internal(
5726  new_quad->set_manifold_id(hex->manifold_id());
5727  for (const auto j : new_quads[i]->line_indices())
5728  new_quad->set_line_orientation(j, true);
5729  }
5730 
5731  // we always get 8 children per refined cell
5732  std::array<
5734  8>
5735  new_hexes;
5736  {
5737  for (unsigned int i = 0; i < n_new_hexes; ++i)
5738  {
5739  if (i % 2 == 0)
5740  next_unused_hex =
5741  triangulation.levels[level + 1]->cells.next_free_hex(
5742  triangulation, level + 1);
5743  else
5744  ++next_unused_hex;
5745 
5746  new_hexes[i] = next_unused_hex;
5747 
5748  auto &new_hex = new_hexes[i];
5749 
5750  // TODO: children have the same type as the parent
5751  triangulation.levels[new_hex->level()]
5752  ->reference_cell[new_hex->index()] =
5753  reference_cell_type;
5754 
5755  AssertIsNotUsed(new_hex);
5756  new_hex->set_used_flag();
5757  new_hex->clear_user_flag();
5758  new_hex->clear_user_data();
5759  new_hex->clear_children();
5760  new_hex->set_material_id(hex->material_id());
5761  new_hex->set_manifold_id(hex->manifold_id());
5762  new_hex->set_subdomain_id(hex->subdomain_id());
5763 
5764  if (i % 2)
5765  new_hex->set_parent(hex->index());
5766  // set the face_orientation flag to true for all
5767  // faces initially, as this is the default value
5768  // which is true for all faces interior to the
5769  // hex. later on go the other way round and
5770  // reset faces that are at the boundary of the
5771  // mother cube
5772  //
5773  // the same is true for the face_flip and
5774  // face_rotation flags. however, the latter two
5775  // are set to false by default as this is the
5776  // standard value
5777  for (const auto f : new_hex->face_indices())
5778  {
5779  new_hex->set_face_orientation(f, true);
5780  new_hex->set_face_flip(f, false);
5781  new_hex->set_face_rotation(f, false);
5782  }
5783  }
5784  for (unsigned int i = 0; i < n_new_hexes / 2; ++i)
5785  hex->set_children(2 * i, new_hexes[2 * i]->index());
5786  }
5787 
5788  {
5789  // load vertex indices
5790  std::array<unsigned int, 27> vertex_indices = {};
5791 
5792  {
5793  unsigned int k = 0;
5794 
5795  for (const unsigned int i : hex->vertex_indices())
5796  vertex_indices[k++] = hex->vertex_index(i);
5797 
5798  for (const unsigned int i : hex->line_indices())
5799  vertex_indices[k++] =
5800  hex->line(i)->child(0)->vertex_index(1);
5801 
5802  if (reference_cell_type == ReferenceCells::Hexahedron)
5803  {
5804  for (const unsigned int i : hex->face_indices())
5805  vertex_indices[k++] =
5806  middle_vertex_index<dim, spacedim>(hex->face(i));
5807 
5808  vertex_indices[k++] = current_vertex;
5809  }
5810  }
5811 
5812  // set up new lines
5813  {
5814  static constexpr std::array<std::array<unsigned int, 2>, 6>
5815  new_line_vertices_hex = {{{{22, 26}},
5816  {{26, 23}},
5817  {{20, 26}},
5818  {{26, 21}},
5819  {{24, 26}},
5820  {{26, 25}}}};
5821 
5822  static constexpr std::array<std::array<unsigned int, 2>, 6>
5823  new_line_vertices_tet = {{{{6, 8}},
5824  {{X, X}},
5825  {{X, X}},
5826  {{X, X}},
5827  {{X, X}},
5828  {{X, X}}}};
5829 
5830  const auto &new_line_vertices =
5831  (reference_cell_type == ReferenceCells::Hexahedron) ?
5832  new_line_vertices_hex :
5833  new_line_vertices_tet;
5834 
5835  for (unsigned int i = 0; i < new_lines.size(); ++i)
5836  new_lines[i]->set_bounding_object_indices(
5837  {vertex_indices[new_line_vertices[i][0]],
5838  vertex_indices[new_line_vertices[i][1]]});
5839  }
5840 
5841  // set up new quads
5842  {
5843  boost::container::small_vector<
5845  30>
5846  relevant_lines(0);
5847 
5848  if (reference_cell_type == ReferenceCells::Hexahedron)
5849  {
5850  relevant_lines.resize(30);
5851  for (unsigned int f = 0, k = 0; f < 6; ++f)
5852  for (unsigned int c = 0; c < 4; ++c, ++k)
5853  {
5854  static constexpr std::
5855  array<std::array<unsigned int, 2>, 4>
5856  temp = {
5857  {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}};
5858 
5859  relevant_lines[k] =
5860  hex->face(f)
5861  ->isotropic_child(
5863  standard_to_real_face_vertex(
5864  temp[c][0],
5865  hex->face_orientation(f),
5866  hex->face_flip(f),
5867  hex->face_rotation(f)))
5868  ->line(GeometryInfo<dim>::
5869  standard_to_real_face_line(
5870  temp[c][1],
5871  hex->face_orientation(f),
5872  hex->face_flip(f),
5873  hex->face_rotation(f)));
5874  }
5875 
5876  for (unsigned int i = 0, k = 24; i < 6; ++i, ++k)
5877  relevant_lines[k] = new_lines[i];
5878  }
5879  else if (reference_cell_type == ReferenceCells::Tetrahedron)
5880  {
5881  relevant_lines.resize(13);
5882 
5883  unsigned int k = 0;
5884  for (unsigned int f = 0; f < 4; ++f)
5885  for (unsigned int l = 0; l < 3; ++l, ++k)
5886  {
5887  // TODO: add comment
5888  static const std::
5889  array<std::array<unsigned int, 3>, 6>
5890  table = {{{{1, 0, 2}}, // 0
5891  {{0, 1, 2}},
5892  {{0, 2, 1}}, // 2
5893  {{1, 2, 0}},
5894  {{2, 1, 0}}, // 4
5895  {{2, 0, 1}}}};
5896 
5897  relevant_lines[k] =
5898  hex->face(f)
5899  ->child(3 /*center triangle*/)
5900  ->line(
5901  table[triangulation.levels[hex->level()]
5902  ->face_orientations
5903  [hex->index() *
5904  GeometryInfo<
5905  dim>::faces_per_cell +
5906  f]][l]);
5907  }
5908 
5909  relevant_lines[k++] = new_lines[0];
5910 
5911  AssertDimension(k, 13);
5912  }
5913  else
5914  Assert(false, ExcNotImplemented());
5915 
5916  boost::container::small_vector<unsigned int, 30>
5917  relevant_line_indices(relevant_lines.size());
5918  for (unsigned int i = 0; i < relevant_line_indices.size();
5919  ++i)
5920  relevant_line_indices[i] = relevant_lines[i]->index();
5921 
5922  static constexpr std::array<std::array<unsigned int, 4>, 12>
5923  new_quad_lines_hex = {{{{10, 28, 16, 24}},
5924  {{28, 14, 17, 25}},
5925  {{11, 29, 24, 20}},
5926  {{29, 15, 25, 21}},
5927  {{18, 26, 0, 28}},
5928  {{26, 22, 1, 29}},
5929  {{19, 27, 28, 4}},
5930  {{27, 23, 29, 5}},
5931  {{2, 24, 8, 26}},
5932  {{24, 6, 9, 27}},
5933  {{3, 25, 26, 12}},
5934  {{25, 7, 27, 13}}}};
5935 
5936  static constexpr std::array<std::array<unsigned int, 4>, 12>
5937  new_quad_lines_tet = {{{{2, 3, 8, X}},
5938  {{0, 9, 5, X}},
5939  {{1, 6, 11, X}},
5940  {{4, 10, 7, X}},
5941  {{2, 12, 5, X}},
5942  {{1, 9, 12, X}},
5943  {{4, 8, 12, X}},
5944  {{6, 12, 10, X}},
5945  {{X, X, X, X}},
5946  {{X, X, X, X}},
5947  {{X, X, X, X}},
5948  {{X, X, X, X}}}};
5949 
5950  static constexpr std::
5951  array<std::array<std::array<unsigned int, 2>, 4>, 12>
5952  table_hex = {
5953  {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
5954  {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
5955  {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
5956  {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
5957  {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
5958  {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
5959  {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
5960  {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
5961  {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
5962  {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
5963  {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
5964  {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
5965 
5966  static constexpr std::
5967  array<std::array<std::array<unsigned int, 2>, 4>, 12>
5968  table_tet = {
5969  {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
5970  {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
5971  {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
5972  {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
5973  {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
5974  {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
5975  {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
5976  {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
5977  {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5978  {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5979  {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5980  {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
5981 
5982  const auto &new_quad_lines =
5983  (reference_cell_type == ReferenceCells::Hexahedron) ?
5984  new_quad_lines_hex :
5985  new_quad_lines_tet;
5986 
5987  const auto &table =
5988  (reference_cell_type == ReferenceCells::Hexahedron) ?
5989  table_hex :
5990  table_tet;
5991 
5992  for (unsigned int q = 0; q < new_quads.size(); ++q)
5993  {
5994  for (unsigned int l = 0; l < 3; ++l)
5995  {
5996  std::array<unsigned int, 2> vertices_0, vertices_1;
5997 
5998  for (unsigned int v = 0; v < 2; ++v)
5999  vertices_0[v] =
6000  relevant_lines[new_quad_lines[q][l]]
6001  ->vertex_index(v);
6002 
6003  for (unsigned int v = 0; v < 2; ++v)
6004  vertices_1[v] = vertex_indices[table[q][l][v]];
6005  }
6006  }
6007 
6008  for (unsigned int q = 0; q < new_quads.size(); ++q)
6009  {
6010  auto &new_quad = new_quads[q];
6011 
6012  if (new_quad->n_lines() == 3)
6013  new_quad->set_bounding_object_indices(
6014  {relevant_line_indices[new_quad_lines[q][0]],
6015  relevant_line_indices[new_quad_lines[q][1]],
6016  relevant_line_indices[new_quad_lines[q][2]]});
6017  else if (new_quad->n_lines() == 4)
6018  new_quad->set_bounding_object_indices(
6019  {relevant_line_indices[new_quad_lines[q][0]],
6020  relevant_line_indices[new_quad_lines[q][1]],
6021  relevant_line_indices[new_quad_lines[q][2]],
6022  relevant_line_indices[new_quad_lines[q][3]]});
6023  else
6024  Assert(false, ExcNotImplemented());
6025 
6026  for (const auto l : new_quad->line_indices())
6027  {
6028  std::array<unsigned int, 2> vertices_0, vertices_1;
6029 
6030  for (unsigned int v = 0; v < 2; ++v)
6031  vertices_0[v] =
6032  relevant_lines[new_quad_lines[q][l]]
6033  ->vertex_index(v);
6034 
6035  for (unsigned int v = 0; v < 2; ++v)
6036  vertices_1[v] = vertex_indices[table[q][l][v]];
6037 
6038  const auto orientation =
6040  vertices_0, vertices_1);
6041 
6042  new_quad->set_line_orientation(l, orientation);
6043  }
6044  }
6045  }
6046 
6047  // set up new hex
6048  {
6049  std::array<int, 36> quad_indices;
6050 
6051  if (reference_cell_type == ReferenceCells::Hexahedron)
6052  {
6053  for (unsigned int i = 0; i < new_quads.size(); ++i)
6054  quad_indices[i] = new_quads[i]->index();
6055 
6056  for (unsigned int f = 0, k = new_quads.size(); f < 6;
6057  ++f)
6058  for (unsigned int c = 0; c < 4; ++c, ++k)
6059  quad_indices[k] =
6060  hex->face(f)->isotropic_child_index(
6062  c,
6063  hex->face_orientation(f),
6064  hex->face_flip(f),
6065  hex->face_rotation(f)));
6066  }
6067  else if (reference_cell_type == ReferenceCells::Tetrahedron)
6068  {
6069  for (unsigned int i = 0; i < new_quads.size(); ++i)
6070  quad_indices[i] = new_quads[i]->index();
6071 
6072  for (unsigned int f = 0, k = new_quads.size(); f < 4;
6073  ++f)
6074  for (unsigned int c = 0; c < 4; ++c, ++k)
6075  {
6076  quad_indices[k] = hex->face(f)->child_index(
6077  (c == 3) ?
6078  3 :
6079  reference_cell_type
6080  .standard_to_real_face_vertex(
6081  c,
6082  f,
6083  triangulation.levels[hex->level()]
6084  ->face_orientations
6085  [hex->index() *
6087  f]));
6088  }
6089  }
6090  else
6091  {
6092  Assert(false, ExcNotImplemented());
6093  }
6094 
6095  static constexpr std::array<std::array<unsigned int, 6>, 8>
6096  cell_quads_hex = {{
6097  {{12, 0, 20, 4, 28, 8}}, // bottom children
6098  {{0, 16, 22, 6, 29, 9}}, //
6099  {{13, 1, 4, 24, 30, 10}}, //
6100  {{1, 17, 6, 26, 31, 11}}, //
6101  {{14, 2, 21, 5, 8, 32}}, // top children
6102  {{2, 18, 23, 7, 9, 33}}, //
6103  {{15, 3, 5, 25, 10, 34}}, //
6104  {{3, 19, 7, 27, 11, 35}} //
6105  }};
6106 
6107  static constexpr std::array<std::array<unsigned int, 6>, 8>
6108  cell_quads_tet{{{{8, 13, 16, 0, X, X}},
6109  {{9, 12, 1, 21, X, X}},
6110  {{10, 2, 17, 20, X, X}},
6111  {{3, 14, 18, 22, X, X}},
6112  {{11, 1, 4, 5, X, X}},
6113  {{15, 0, 4, 6, X, X}},
6114  {{19, 7, 6, 3, X, X}},
6115  {{23, 5, 2, 7, X, X}}}};
6116 
6117  static constexpr std::
6118  array<std::array<std::array<unsigned int, 4>, 6>, 8>
6119  cell_face_vertices_hex{{{{{{0, 8, 16, 20}},
6120  {{10, 24, 22, 26}},
6121  {{0, 16, 10, 22}},
6122  {{8, 20, 24, 26}},
6123  {{0, 10, 8, 24}},
6124  {{16, 22, 20, 26}}}},
6125  {{{{10, 24, 22, 26}},
6126  {{1, 9, 17, 21}},
6127  {{10, 22, 1, 17}},
6128  {{24, 26, 9, 21}},
6129  {{10, 1, 24, 9}},
6130  {{22, 17, 26, 21}}}},
6131  {{{{8, 2, 20, 18}},
6132  {{24, 11, 26, 23}},
6133  {{8, 20, 24, 26}},
6134  {{2, 18, 11, 23}},
6135  {{8, 24, 2, 11}},
6136  {{20, 26, 18, 23}}}},
6137  {{{{24, 11, 26, 23}},
6138  {{9, 3, 21, 19}},
6139  {{24, 26, 9, 21}},
6140  {{11, 23, 3, 19}},
6141  {{24, 9, 11, 3}},
6142  {{26, 21, 23, 19}}}},
6143  {{{{16, 20, 4, 12}},
6144  {{22, 26, 14, 25}},
6145  {{16, 4, 22, 14}},
6146  {{20, 12, 26, 25}},
6147  {{16, 22, 20, 26}},
6148  {{4, 14, 12, 25}}}},
6149  {{{{22, 26, 14, 25}},
6150  {{17, 21, 5, 13}},
6151  {{22, 14, 17, 5}},
6152  {{26, 25, 21, 13}},
6153  {{22, 17, 26, 21}},
6154  {{14, 5, 25, 13}}}},
6155  {{{{20, 18, 12, 6}},
6156  {{26, 23, 25, 15}},
6157  {{20, 12, 26, 25}},
6158  {{18, 6, 23, 15}},
6159  {{20, 26, 18, 23}},
6160  {{12, 25, 6, 15}}}},
6161  {{{{26, 23, 25, 15}},
6162  {{21, 19, 13, 7}},
6163  {{26, 25, 21, 13}},
6164  {{23, 15, 19, 7}},
6165  {{26, 21, 23, 19}},
6166  {{25, 13, 15, 7}}}}}};
6167 
6168  static constexpr std::
6169  array<std::array<std::array<unsigned int, 4>, 6>, 8>
6170  cell_face_vertices_tet{{{{{{0, 4, 6, X}},
6171  {{4, 0, 7, X}},
6172  {{0, 6, 7, X}},
6173  {{6, 4, 7, X}},
6174  {{X, X, X, X}},
6175  {{X, X, X, X}}}},
6176  {{{{4, 1, 5, X}},
6177  {{1, 4, 8, X}},
6178  {{4, 5, 8, X}},
6179  {{5, 1, 8, X}},
6180  {{X, X, X, X}},
6181  {{X, X, X, X}}}},
6182  {{{{6, 5, 2, X}},
6183  {{5, 6, 9, X}},
6184  {{6, 2, 9, X}},
6185  {{2, 5, 9, X}},
6186  {{X, X, X, X}},
6187  {{X, X, X, X}}}},
6188  {{{{7, 8, 9, X}},
6189  {{8, 7, 3, X}},
6190  {{7, 9, 3, X}},
6191  {{9, 8, 3, X}},
6192  {{X, X, X, X}},
6193  {{X, X, X, X}}}},
6194  {{{{4, 5, 6, X}},
6195  {{5, 4, 8, X}},
6196  {{4, 6, 8, X}},
6197  {{6, 5, 8, X}},
6198  {{X, X, X, X}},
6199  {{X, X, X, X}}}},
6200  {{{{4, 7, 8, X}},
6201  {{7, 4, 6, X}},
6202  {{4, 8, 6, X}},
6203  {{8, 7, 6, X}},
6204  {{X, X, X, X}},
6205  {{X, X, X, X}}}},
6206  {{{{6, 9, 7, X}},
6207  {{9, 6, 8, X}},
6208  {{6, 7, 8, X}},
6209  {{7, 9, 8, X}},
6210  {{X, X, X, X}},
6211  {{X, X, X, X}}}},
6212  {{{{5, 8, 9, X}},
6213  {{8, 5, 6, X}},
6214  {{5, 9, 6, X}},
6215  {{9, 8, 6, X}},
6216  {{X, X, X, X}},
6217  {{X, X, X, X}}}}}};
6218 
6219  const auto &cell_quads =
6220  (reference_cell_type == ReferenceCells::Hexahedron) ?
6221  cell_quads_hex :
6222  cell_quads_tet;
6223 
6224  const auto &cell_face_vertices =
6225  (reference_cell_type == ReferenceCells::Hexahedron) ?
6226  cell_face_vertices_hex :
6227  cell_face_vertices_tet;
6228 
6229  for (unsigned int c = 0;
6230  c < GeometryInfo<dim>::max_children_per_cell;
6231  ++c)
6232  {
6233  auto &new_hex = new_hexes[c];
6234 
6235  if (new_hex->n_faces() == 4)
6236  new_hex->set_bounding_object_indices(
6237  {quad_indices[cell_quads[c][0]],
6238  quad_indices[cell_quads[c][1]],
6239  quad_indices[cell_quads[c][2]],
6240  quad_indices[cell_quads[c][3]]});
6241  else if (new_hex->n_faces() == 6)
6242  new_hex->set_bounding_object_indices(
6243  {quad_indices[cell_quads[c][0]],
6244  quad_indices[cell_quads[c][1]],
6245  quad_indices[cell_quads[c][2]],
6246  quad_indices[cell_quads[c][3]],
6247  quad_indices[cell_quads[c][4]],
6248  quad_indices[cell_quads[c][5]]});
6249  else
6250  Assert(false, ExcNotImplemented());
6251 
6252  for (const auto f : new_hex->face_indices())
6253  {
6254  std::array<unsigned int, 4> vertices_0, vertices_1;
6255 
6256  const auto &face = new_hex->face(f);
6257 
6258  for (const auto i : face->vertex_indices())
6259  vertices_0[i] = face->vertex_index(i);
6260 
6261  for (const auto i : face->vertex_indices())
6262  vertices_1[i] =
6263  vertex_indices[cell_face_vertices[c][f][i]];
6264 
6265  const auto orientation =
6266  face->reference_cell().compute_orientation(
6267  vertices_1, vertices_0);
6268 
6269  new_hex->set_face_orientation(
6270  f, Utilities::get_bit(orientation, 0));
6271  new_hex->set_face_flip(
6272  f, Utilities::get_bit(orientation, 2));
6273  new_hex->set_face_rotation(
6274  f, Utilities::get_bit(orientation, 1));
6275  }
6276  }
6277  }
6278  }
6279 
6280  if (check_for_distorted_cells &&
6281  has_distorted_children<dim, spacedim>(hex))
6282  cells_with_distorted_children.distorted_cells.push_back(hex);
6283 
6284  triangulation.signals.post_refinement_on_cell(hex);
6285  }
6286  }
6287 
6288  triangulation.faces->quads.clear_user_data();
6289 
6290  return cells_with_distorted_children;
6291  }
6292 
6297  template <int spacedim>
6300  const bool check_for_distorted_cells)
6301  {
6302  const unsigned int dim = 3;
6303 
6304  {
6305  bool flag_isotropic_mesh = true;
6307  cell = triangulation.begin(),
6308  endc = triangulation.end();
6309  for (; cell != endc; ++cell)
6310  if (cell->used())
6311  if (triangulation.get_anisotropic_refinement_flag() ||
6312  cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
6313  cell->refine_flag_set() == RefinementCase<dim>::cut_y ||
6314  cell->refine_flag_set() == RefinementCase<dim>::cut_z ||
6315  cell->refine_flag_set() == RefinementCase<dim>::cut_xy ||
6316  cell->refine_flag_set() == RefinementCase<dim>::cut_xz ||
6317  cell->refine_flag_set() == RefinementCase<dim>::cut_yz)
6318  {
6319  flag_isotropic_mesh = false;
6320  break;
6321  }
6322 
6323  if (flag_isotropic_mesh)
6324  return execute_refinement_isotropic(triangulation,
6325  check_for_distorted_cells);
6326  }
6327 
6328  // this function probably also works for spacedim>3 but it
6329  // isn't tested. it will probably be necessary to pull new
6330  // vertices onto the manifold just as we do for the other
6331  // functions above.
6332  Assert(spacedim == 3, ExcNotImplemented());
6333 
6334  // Check whether a new level is needed. We have to check for
6335  // this on the highest level only
6336  for (const auto &cell : triangulation.active_cell_iterators_on_level(
6337  triangulation.levels.size() - 1))
6338  if (cell->refine_flag_set())
6339  {
6340  triangulation.levels.push_back(
6341  std::make_unique<
6343  break;
6344  }
6345 
6346 
6347  // first clear user flags for quads and lines; we're going to
6348  // use them to flag which lines and quads need refinement
6349  triangulation.faces->quads.clear_user_data();
6350 
6351  for (typename Triangulation<dim, spacedim>::line_iterator line =
6352  triangulation.begin_line();
6353  line != triangulation.end_line();
6354  ++line)
6355  line->clear_user_flag();
6356  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
6357  triangulation.begin_quad();
6358  quad != triangulation.end_quad();
6359  ++quad)
6360  quad->clear_user_flag();
6361 
6362  // create an array of face refine cases. User indices of faces
6363  // will be set to values corresponding with indices in this
6364  // array.
6365  const RefinementCase<dim - 1> face_refinement_cases[4] = {
6366  RefinementCase<dim - 1>::no_refinement,
6367  RefinementCase<dim - 1>::cut_x,
6368  RefinementCase<dim - 1>::cut_y,
6369  RefinementCase<dim - 1>::cut_xy};
6370 
6371  // check how much space is needed on every level. We need not
6372  // check the highest level since either
6373  // - on the highest level no cells are flagged for refinement
6374  // - there are, but prepare_refinement added another empty
6375  // level which then is the highest level
6376 
6377  // variables to hold the number of newly to be created
6378  // vertices, lines and quads. as these are stored globally,
6379  // declare them outside the loop over al levels. we need lines
6380  // and quads in pairs for refinement of old ones and lines and
6381  // quads, that can be stored as single ones, as they are newly
6382  // created in the inside of an existing cell
6383  unsigned int needed_vertices = 0;
6384  unsigned int needed_lines_single = 0;
6385  unsigned int needed_quads_single = 0;
6386  unsigned int needed_lines_pair = 0;
6387  unsigned int needed_quads_pair = 0;
6388  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
6389  {
6390  // count number of flagged cells on this level and compute
6391  // how many new vertices and new lines will be needed
6392  unsigned int new_cells = 0;
6393 
6394  for (const auto &acell :
6395  triangulation.active_cell_iterators_on_level(level))
6396  if (acell->refine_flag_set())
6397  {
6398  RefinementCase<dim> ref_case = acell->refine_flag_set();
6399 
6400  // now for interior vertices, lines and quads, which
6401  // are needed in any case
6402  if (ref_case == RefinementCase<dim>::cut_x ||
6403  ref_case == RefinementCase<dim>::cut_y ||
6404  ref_case == RefinementCase<dim>::cut_z)
6405  {
6406  ++needed_quads_single;
6407  new_cells += 2;
6408  triangulation.anisotropic_refinement = true;
6409  }
6410  else if (ref_case == RefinementCase<dim>::cut_xy ||
6411  ref_case == RefinementCase<dim>::cut_xz ||
6412  ref_case == RefinementCase<dim>::cut_yz)
6413  {
6414  ++needed_lines_single;
6415  needed_quads_single += 4;
6416  new_cells += 4;
6417  triangulation.anisotropic_refinement = true;
6418  }
6419  else if (ref_case == RefinementCase<dim>::cut_xyz)
6420  {
6421  ++needed_vertices;
6422  needed_lines_single += 6;
6423  needed_quads_single += 12;
6424  new_cells += 8;
6425  }
6426  else
6427  {
6428  // we should never get here
6429  Assert(false, ExcInternalError());
6430  }
6431 
6432  // mark all faces for refinement; checking locally
6433  // if and how the neighbor would like to refine
6434  // these is difficult so we only flag them and after
6435  // visiting all cells, we decide which faces need
6436  // which refinement;
6437  for (const unsigned int face :
6439  {
6441  aface = acell->face(face);
6442  // get the RefineCase this faces has for the
6443  // given RefineCase of the cell
6444  RefinementCase<dim - 1> face_ref_case =
6446  ref_case,
6447  face,
6448  acell->face_orientation(face),
6449  acell->face_flip(face),
6450  acell->face_rotation(face));
6451  // only do something, if this face has to be
6452  // refined
6453  if (face_ref_case)
6454  {
6455  if (face_ref_case ==
6457  {
6458  if (aface->n_active_descendants() < 4)
6459  // we use user_flags to denote needed
6460  // isotropic refinement
6461  aface->set_user_flag();
6462  }
6463  else if (aface->refinement_case() != face_ref_case)
6464  // we use user_indices to denote needed
6465  // anisotropic refinement. note, that we
6466  // can have at most one anisotropic
6467  // refinement case for this face, as
6468  // otherwise prepare_refinement() would
6469  // have changed one of the cells to yield
6470  // isotropic refinement at this
6471  // face. therefore we set the user_index
6472  // uniquely
6473  {
6474  Assert(aface->refinement_case() ==
6476  dim - 1>::isotropic_refinement ||
6477  aface->refinement_case() ==
6479  ExcInternalError());
6480  aface->set_user_index(face_ref_case);
6481  }
6482  }
6483  } // for all faces
6484