Reference documentation for deal.II version Git eab6000ff4 2020-05-30 14:50:38 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
19 
20 #include <deal.II/fe/mapping_q1.h>
21 
24 #include <deal.II/grid/manifold.h>
25 #include <deal.II/grid/tria.h>
30 
32 #include <deal.II/lac/vector.h>
33 
34 #include <algorithm>
35 #include <array>
36 #include <cmath>
37 #include <functional>
38 #include <list>
39 #include <map>
40 #include <memory>
41 #include <numeric>
42 
43 
45 
46 
47 namespace internal
48 {
49  namespace TriangulationImplementation
50  {
52  : n_levels(0)
53  , n_lines(0)
54  , n_active_lines(0)
55  // all other fields are
56  // default constructed
57  {}
58 
59 
60 
61  std::size_t
63  {
64  return (MemoryConsumption::memory_consumption(n_levels) +
68  MemoryConsumption::memory_consumption(n_active_lines_level));
69  }
70 
71 
73  : n_quads(0)
74  , n_active_quads(0)
75  // all other fields are
76  // default constructed
77  {}
78 
79 
80 
81  std::size_t
83  {
88  MemoryConsumption::memory_consumption(n_active_quads_level));
89  }
90 
91 
92 
94  : n_hexes(0)
95  , n_active_hexes(0)
96  // all other fields are
97  // default constructed
98  {}
99 
100 
101 
102  std::size_t
104  {
108  MemoryConsumption::memory_consumption(n_active_hexes) +
109  MemoryConsumption::memory_consumption(n_active_hexes_level));
110  }
111  } // namespace TriangulationImplementation
112 } // namespace internal
113 
114 // anonymous namespace for internal helper functions
115 namespace
116 {
117  // return whether the given cell is
118  // patch_level_1, i.e. determine
119  // whether either all or none of
120  // its children are further
121  // refined. this function can only
122  // be called for non-active cells.
123  template <int dim, int spacedim>
124  bool
125  cell_is_patch_level_1(
127  {
128  Assert(cell->is_active() == false, ExcInternalError());
129 
130  unsigned int n_active_children = 0;
131  for (unsigned int i = 0; i < cell->n_children(); ++i)
132  if (cell->child(i)->is_active())
133  ++n_active_children;
134 
135  return (n_active_children == 0) ||
136  (n_active_children == cell->n_children());
137  }
138 
139 
140 
141  // return, whether a given @p cell will be
142  // coarsened, which is the case if all
143  // children are active and have their coarsen
144  // flag set. In case only part of the coarsen
145  // flags are set, remove them.
146  template <int dim, int spacedim>
147  bool
148  cell_will_be_coarsened(
150  {
151  // only cells with children should be
152  // considered for coarsening
153 
154  if (cell->has_children())
155  {
156  unsigned int children_to_coarsen = 0;
157  const unsigned int n_children = cell->n_children();
158 
159  for (unsigned int c = 0; c < n_children; ++c)
160  if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
161  ++children_to_coarsen;
162  if (children_to_coarsen == n_children)
163  return true;
164  else
165  for (unsigned int c = 0; c < n_children; ++c)
166  if (cell->child(c)->is_active())
167  cell->child(c)->clear_coarsen_flag();
168  }
169  // no children, so no coarsening
170  // possible. however, no children also
171  // means that this cell will be in the same
172  // state as if it had children and was
173  // coarsened. So, what should we return -
174  // false or true?
175  // make sure we do not have to do this at
176  // all...
177  Assert(cell->has_children(), ExcInternalError());
178  // ... and then simply return false
179  return false;
180  }
181 
182 
183  // return, whether the face @p face_no of the
184  // given @p cell will be refined after the
185  // current refinement step, considering
186  // refine and coarsen flags and considering
187  // only those refinemnts that will be caused
188  // by the neighboring cell.
189 
190  // this function is used on both active cells
191  // and cells with children. on cells with
192  // children it also of interest to know 'how'
193  // the face will be refined. thus there is an
194  // additional third argument @p
195  // expected_face_ref_case returning just
196  // that. be aware, that this vriable will
197  // only contain useful information if this
198  // function is called for an active cell.
199  //
200  // thus, this is an internal function, users
201  // should call one of the two alternatives
202  // following below.
203  template <int dim, int spacedim>
204  bool
205  face_will_be_refined_by_neighbor_internal(
207  const unsigned int face_no,
208  RefinementCase<dim - 1> &expected_face_ref_case)
209  {
210  // first of all: set the default value for
211  // expected_face_ref_case, which is no
212  // refinement at all
213  expected_face_ref_case = RefinementCase<dim - 1>::no_refinement;
214 
215  const typename Triangulation<dim, spacedim>::cell_iterator neighbor =
216  cell->neighbor(face_no);
217 
218  // If we are at the boundary, there is no
219  // neighbor which could refine the face
220  if (neighbor.state() != IteratorState::valid)
221  return false;
222 
223  if (neighbor->has_children())
224  {
225  // if the neighbor is refined, it may be
226  // coarsened. if so, then it won't refine
227  // the face, no matter what else happens
228  if (cell_will_be_coarsened(neighbor))
229  return false;
230  else
231  // if the neighbor is refined, then it
232  // is also refined at our current
233  // face. He will stay so without
234  // coarsening, so return true in that
235  // case.
236  {
237  expected_face_ref_case = cell->face(face_no)->refinement_case();
238  return true;
239  }
240  }
241 
242  // now, the neighbor is not refined, but
243  // perhaps it will be
244  const RefinementCase<dim> nb_ref_flag = neighbor->refine_flag_set();
245  if (nb_ref_flag != RefinementCase<dim>::no_refinement)
246  {
247  // now we need to know, which of the
248  // neighbors faces points towards us
249  const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
250  // check, whether the cell will be
251  // refined in a way that refines our
252  // face
253  const RefinementCase<dim - 1> face_ref_case =
255  nb_ref_flag,
256  neighbor_neighbor,
257  neighbor->face_orientation(neighbor_neighbor),
258  neighbor->face_flip(neighbor_neighbor),
259  neighbor->face_rotation(neighbor_neighbor));
260  if (face_ref_case != RefinementCase<dim - 1>::no_refinement)
261  {
263  neighbor_face = neighbor->face(neighbor_neighbor);
264  const int this_face_index = cell->face_index(face_no);
265 
266  // there are still two basic
267  // possibilities here: the neighbor
268  // might be coarser or as coarse
269  // as we are
270  if (neighbor_face->index() == this_face_index)
271  // the neighbor is as coarse as
272  // we are and will be refined at
273  // the face of consideration, so
274  // return true
275  {
276  expected_face_ref_case = face_ref_case;
277  return true;
278  }
279  else
280  {
281  // the neighbor is coarser.
282  // this is the most complicated
283  // case. It might be, that the
284  // neighbor's face will be
285  // refined, but that we will
286  // not see this, as we are
287  // refined in a similar way.
288 
289  // so, the neighbor's face must
290  // have children. check, if our
291  // cell's face is one of these
292  // (it could also be a
293  // grand_child)
294  for (unsigned int c = 0; c < neighbor_face->n_children(); ++c)
295  if (neighbor_face->child_index(c) == this_face_index)
296  {
297  // if the flagged refine
298  // case of the face is a
299  // subset or the same as
300  // the current refine case,
301  // then the face, as seen
302  // from our cell, won't be
303  // refined by the neighbor
304  if ((neighbor_face->refinement_case() | face_ref_case) ==
305  neighbor_face->refinement_case())
306  return false;
307  else
308  {
309  // if we are active, we
310  // must be an
311  // anisotropic child
312  // and the coming
313  // face_ref_case is
314  // isotropic. Thus,
315  // from our cell we
316  // will see exactly the
317  // opposite refine case
318  // that the face has
319  // now...
320  Assert(
321  face_ref_case ==
323  ExcInternalError());
324  expected_face_ref_case =
325  ~neighbor_face->refinement_case();
326  return true;
327  }
328  }
329 
330  // so, obviously we were not
331  // one of the children, but a
332  // grandchild. This is only
333  // possible in 3d.
334  Assert(dim == 3, ExcInternalError());
335  // In that case, however, no
336  // matter what the neighbor
337  // does, it won't be finer
338  // after the next refinement
339  // step.
340  return false;
341  }
342  } // if face will be refined
343  } // if neighbor is flagged for refinement
344 
345  // no cases left, so the neighbor will not
346  // refine the face
347  return false;
348  }
349 
350  // version of above function for both active
351  // and non-active cells
352  template <int dim, int spacedim>
353  bool
354  face_will_be_refined_by_neighbor(
356  const unsigned int face_no)
357  {
358  RefinementCase<dim - 1> dummy = RefinementCase<dim - 1>::no_refinement;
359  return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
360  }
361 
362  // version of above function for active cells
363  // only. Additionally returning the refine
364  // case (to come) of the face under
365  // consideration
366  template <int dim, int spacedim>
367  bool
368  face_will_be_refined_by_neighbor(
370  const unsigned int face_no,
371  RefinementCase<dim - 1> &expected_face_ref_case)
372  {
373  return face_will_be_refined_by_neighbor_internal(cell,
374  face_no,
375  expected_face_ref_case);
376  }
377 
378 
379 
380  template <int dim, int spacedim>
381  bool
382  satisfies_level1_at_vertex_rule(
384  {
385  std::vector<unsigned int> min_adjacent_cell_level(
386  triangulation.n_vertices(), triangulation.n_levels());
387  std::vector<unsigned int> max_adjacent_cell_level(
388  triangulation.n_vertices(), 0);
389 
390  for (const auto &cell : triangulation.active_cell_iterators())
391  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
392  {
393  min_adjacent_cell_level[cell->vertex_index(v)] =
395  min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
396  max_adjacent_cell_level[cell->vertex_index(v)] =
398  min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
399  }
400 
401  for (unsigned int k = 0; k < triangulation.n_vertices(); ++k)
402  if (triangulation.vertex_used(k))
403  if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
404  return false;
405  return true;
406  }
407 
408 
409 
416  template <int dim, int spacedim>
417  std::vector<unsigned int>
418  count_cells_bounded_by_line(const Triangulation<dim, spacedim> &triangulation)
419  {
420  if (dim >= 2)
421  {
422  std::vector<unsigned int> line_cell_count(triangulation.n_raw_lines(),
423  0);
425  cell = triangulation.begin(),
426  endc = triangulation.end();
427  for (; cell != endc; ++cell)
428  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
429  ++line_cell_count[cell->line_index(l)];
430  return line_cell_count;
431  }
432  else
433  return std::vector<unsigned int>();
434  }
435 
436 
437 
444  template <int dim, int spacedim>
445  std::vector<unsigned int>
446  count_cells_bounded_by_quad(const Triangulation<dim, spacedim> &triangulation)
447  {
448  if (dim >= 3)
449  {
450  std::vector<unsigned int> quad_cell_count(triangulation.n_raw_quads(),
451  0);
453  cell = triangulation.begin(),
454  endc = triangulation.end();
455  for (; cell != endc; ++cell)
456  for (unsigned int q : GeometryInfo<dim>::face_indices())
457  ++quad_cell_count[cell->quad_index(q)];
458  return quad_cell_count;
459  }
460  else
461  return std::vector<unsigned int>();
462  }
463 
464 
465 
477  void
478  reorder_compatibility(const std::vector<CellData<1>> &, const SubCellData &)
479  {
480  // nothing to do here: the format
481  // hasn't changed for 1d
482  }
483 
484 
485  void reorder_compatibility(std::vector<CellData<2>> &cells,
486  const SubCellData &)
487  {
488  for (auto &cell : cells)
489  std::swap(cell.vertices[2], cell.vertices[3]);
490  }
491 
492 
493  void reorder_compatibility(std::vector<CellData<3>> &cells,
494  SubCellData & subcelldata)
495  {
496  unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
497  for (auto &cell : cells)
498  {
499  for (const unsigned int i : GeometryInfo<3>::vertex_indices())
500  tmp[i] = cell.vertices[i];
501  for (const unsigned int i : GeometryInfo<3>::vertex_indices())
502  cell.vertices[GeometryInfo<3>::ucd_to_deal[i]] = tmp[i];
503  }
504 
505  // now points in boundary quads
506  std::vector<CellData<2>>::iterator boundary_quad =
507  subcelldata.boundary_quads.begin();
508  std::vector<CellData<2>>::iterator end_quad =
509  subcelldata.boundary_quads.end();
510  for (unsigned int quad_no = 0; boundary_quad != end_quad;
511  ++boundary_quad, ++quad_no)
512  std::swap(boundary_quad->vertices[2], boundary_quad->vertices[3]);
513  }
514 
515 
516 
534  template <int dim, int spacedim>
535  unsigned int
536  middle_vertex_index(
537  const typename Triangulation<dim, spacedim>::line_iterator &line)
538  {
539  if (line->has_children())
540  return line->child(0)->vertex_index(1);
542  }
543 
544 
545  template <int dim, int spacedim>
546  unsigned int
547  middle_vertex_index(
548  const typename Triangulation<dim, spacedim>::quad_iterator &quad)
549  {
550  switch (static_cast<unsigned char>(quad->refinement_case()))
551  {
553  return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
554  break;
556  return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
557  break;
559  return quad->child(0)->vertex_index(3);
560  break;
561  default:
562  break;
563  }
565  }
566 
567 
568  template <int dim, int spacedim>
569  unsigned int
570  middle_vertex_index(
571  const typename Triangulation<dim, spacedim>::hex_iterator &hex)
572  {
573  switch (static_cast<unsigned char>(hex->refinement_case()))
574  {
576  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
577  break;
579  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
580  break;
582  return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
583  break;
585  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
586  break;
588  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
589  break;
591  return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
592  break;
594  return hex->child(0)->vertex_index(7);
595  break;
596  default:
597  break;
598  }
600  }
601 
602 
615  template <class TRIANGULATION>
616  inline typename TRIANGULATION::DistortedCellList
617  collect_distorted_coarse_cells(const TRIANGULATION &)
618  {
619  return typename TRIANGULATION::DistortedCellList();
620  }
621 
622 
623 
632  template <int dim>
634  collect_distorted_coarse_cells(const Triangulation<dim, dim> &triangulation)
635  {
636  typename Triangulation<dim, dim>::DistortedCellList distorted_cells;
637  for (typename Triangulation<dim, dim>::cell_iterator cell =
638  triangulation.begin(0);
639  cell != triangulation.end(0);
640  ++cell)
641  {
643  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
644  vertices[i] = cell->vertex(i);
645 
647  GeometryInfo<dim>::alternating_form_at_vertices(vertices, determinants);
648 
649  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
650  if (determinants[i] <= 1e-9 * std::pow(cell->diameter(), 1. * dim))
651  {
652  distorted_cells.distorted_cells.push_back(cell);
653  break;
654  }
655  }
656 
657  return distorted_cells;
658  }
659 
660 
667  template <int dim>
668  bool
669  has_distorted_children(
670  const typename Triangulation<dim, dim>::cell_iterator &cell,
671  std::integral_constant<int, dim>,
672  std::integral_constant<int, dim>)
673  {
674  Assert(cell->has_children(), ExcInternalError());
675 
676  for (unsigned int c = 0; c < cell->n_children(); ++c)
677  {
679  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
680  vertices[i] = cell->child(c)->vertex(i);
681 
683  GeometryInfo<dim>::alternating_form_at_vertices(vertices, determinants);
684 
685  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
686  if (determinants[i] <=
687  1e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
688  return true;
689  }
690 
691  return false;
692  }
693 
694 
702  template <int dim, int spacedim>
703  bool
704  has_distorted_children(
706  std::integral_constant<int, dim>,
707  std::integral_constant<int, spacedim>)
708  {
709  return false;
710  }
711 
712 
713 
718  template <int spacedim>
719  void update_neighbors(Triangulation<1, spacedim> &)
720  {}
721 
722 
723  template <int dim, int spacedim>
724  void
725  update_neighbors(Triangulation<dim, spacedim> &triangulation)
726  {
727  // each face can be neighbored on two sides
728  // by cells. according to the face's
729  // intrinsic normal we define the left
730  // neighbor as the one for which the face
731  // normal points outward, and store that
732  // one first; the second one is then
733  // the right neighbor for which the
734  // face normal points inward. This
735  // information depends on the type of cell
736  // and local number of face for the
737  // 'standard ordering and orientation' of
738  // faces and then on the face_orientation
739  // information for the real mesh. Set up a
740  // table to have fast access to those
741  // offsets (0 for left and 1 for
742  // right). Some of the values are invalid
743  // as they reference too large face
744  // numbers, but we just leave them at a
745  // zero value.
746  //
747  // Note, that in 2d for lines as faces the
748  // normal direction given in the
749  // GeometryInfo class is not consistent. We
750  // thus define here that the normal for a
751  // line points to the right if the line
752  // points upwards.
753  //
754  // There is one more point to
755  // consider, however: if we have
756  // dim<spacedim, then we may have
757  // cases where cells are
758  // inverted. In effect, both
759  // cells think they are the left
760  // neighbor of an edge, for
761  // example, which leads us to
762  // forget neighborship
763  // information (a case that shows
764  // this is
765  // codim_one/hanging_nodes_02). We
766  // store whether a cell is
767  // inverted using the
768  // direction_flag, so if a cell
769  // has a false direction_flag,
770  // then we need to invert our
771  // selection whether we are a
772  // left or right neighbor in all
773  // following computations.
774  //
775  // first index: dimension (minus 2)
776  // second index: local face index
777  // third index: face_orientation (false and true)
778  static const unsigned int left_right_offset[2][6][2] = {
779  // quadrilateral
780  {{0, 1}, // face 0, face_orientation = false and true
781  {1, 0}, // face 1, face_orientation = false and true
782  {1, 0}, // face 2, face_orientation = false and true
783  {0, 1}, // face 3, face_orientation = false and true
784  {0, 0}, // face 4, invalid face
785  {0, 0}}, // face 5, invalid face
786  // hexahedron
787  {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
788 
789  // now create a vector of the two active
790  // neighbors (left and right) for each face
791  // and fill it by looping over all cells. For
792  // cases with anisotropic refinement and more
793  // then one cell neighboring at a given side
794  // of the face we will automatically get the
795  // active one on the highest level as we loop
796  // over cells from lower levels first.
798  std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
799  adjacent_cells(2 * triangulation.n_raw_faces(), dummy);
800 
801  typename Triangulation<dim, spacedim>::cell_iterator cell = triangulation
802  .begin(),
803  endc =
804  triangulation.end();
805  for (; cell != endc; ++cell)
806  for (auto f : GeometryInfo<dim>::face_indices())
807  {
808  const typename Triangulation<dim, spacedim>::face_iterator face =
809  cell->face(f);
810 
811  const unsigned int offset =
812  (cell->direction_flag() ?
813  left_right_offset[dim - 2][f][cell->face_orientation(f)] :
814  1 - left_right_offset[dim - 2][f][cell->face_orientation(f)]);
815 
816  adjacent_cells[2 * face->index() + offset] = cell;
817 
818  // if this cell is not refined, but the
819  // face is, then we'll have to set our
820  // cell as neighbor for the child faces
821  // as well. Fortunately the normal
822  // orientation of children will be just
823  // the same.
824  if (dim == 2)
825  {
826  if (cell->is_active() && face->has_children())
827  {
828  adjacent_cells[2 * face->child(0)->index() + offset] = cell;
829  adjacent_cells[2 * face->child(1)->index() + offset] = cell;
830  }
831  }
832  else // -> dim == 3
833  {
834  // We need the same as in 2d
835  // here. Furthermore, if the face is
836  // refined with cut_x or cut_y then
837  // those children again in the other
838  // direction, and if this cell is
839  // refined isotropically (along the
840  // face) then the neighbor will
841  // (probably) be refined as cut_x or
842  // cut_y along the face. For those
843  // neighboring children cells, their
844  // neighbor will be the current,
845  // inactive cell, as our children are
846  // too fine to be neighbors. Catch that
847  // case by also acting on inactive
848  // cells with isotropic refinement
849  // along the face. If the situation
850  // described is not present, the data
851  // will be overwritten later on when we
852  // visit cells on finer levels, so no
853  // harm will be done.
854  if (face->has_children() &&
855  (cell->is_active() ||
857  cell->refinement_case(), f) ==
859  {
860  for (unsigned int c = 0; c < face->n_children(); ++c)
861  adjacent_cells[2 * face->child(c)->index() + offset] = cell;
862  if (face->child(0)->has_children())
863  {
864  adjacent_cells[2 * face->child(0)->child(0)->index() +
865  offset] = cell;
866  adjacent_cells[2 * face->child(0)->child(1)->index() +
867  offset] = cell;
868  }
869  if (face->child(1)->has_children())
870  {
871  adjacent_cells[2 * face->child(1)->child(0)->index() +
872  offset] = cell;
873  adjacent_cells[2 * face->child(1)->child(1)->index() +
874  offset] = cell;
875  }
876  } // if cell active and face refined
877  } // else -> dim==3
878  } // for all faces of all cells
879 
880  // now loop again over all cells and set the
881  // corresponding neighbor cell. Note, that we
882  // have to use the opposite of the
883  // left_right_offset in this case as we want
884  // the offset of the neighbor, not our own.
885  for (cell = triangulation.begin(); cell != endc; ++cell)
886  for (auto f : GeometryInfo<dim>::face_indices())
887  {
888  const unsigned int offset =
889  (cell->direction_flag() ?
890  left_right_offset[dim - 2][f][cell->face_orientation(f)] :
891  1 - left_right_offset[dim - 2][f][cell->face_orientation(f)]);
892  cell->set_neighbor(
893  f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
894  }
895  }
896 
897 
898  template <int dim, int spacedim>
899  void
900  update_periodic_face_map_recursively(
901  const typename Triangulation<dim, spacedim>::cell_iterator &cell_1,
902  const typename Triangulation<dim, spacedim>::cell_iterator &cell_2,
903  unsigned int n_face_1,
904  unsigned int n_face_2,
905  const std::bitset<3> & orientation,
906  typename std::map<
908  unsigned int>,
909  std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
910  unsigned int>,
911  std::bitset<3>>> &periodic_face_map)
912  {
913  using FaceIterator = typename Triangulation<dim, spacedim>::face_iterator;
914  const FaceIterator face_1 = cell_1->face(n_face_1);
915  const FaceIterator face_2 = cell_2->face(n_face_2);
916 
917  const bool face_orientation = orientation[0];
918  const bool face_flip = orientation[1];
919  const bool face_rotation = orientation[2];
920 
921  Assert((dim != 1) || (face_orientation == true && face_flip == false &&
922  face_rotation == false),
923  ExcMessage("The supplied orientation "
924  "(face_orientation, face_flip, face_rotation) "
925  "is invalid for 1D"));
926 
927  Assert((dim != 2) || (face_orientation == true && face_rotation == false),
928  ExcMessage("The supplied orientation "
929  "(face_orientation, face_flip, face_rotation) "
930  "is invalid for 2D"));
931 
932  Assert(face_1 != face_2, ExcMessage("face_1 and face_2 are equal!"));
933 
934  Assert(face_1->at_boundary() && face_2->at_boundary(),
935  ExcMessage("Periodic faces must be on the boundary"));
936 
937  // Check if the requirement that each edge can only have at most one hanging
938  // node, and as a consequence neighboring cells can differ by at most
939  // one refinement level is enforced. In 1d, there are no hanging nodes and
940  // so neighboring cells can differ by more than one refinement level.
941  Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
942  ExcInternalError());
943 
944  // insert periodic face pair for both cells
945  using CellFace =
946  std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
947  unsigned int>;
948  const CellFace cell_face_1(cell_1, n_face_1);
949  const CellFace cell_face_2(cell_2, n_face_2);
950  const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
951  cell_face_2, orientation);
952 
953  const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
954  periodic_faces(cell_face_1, cell_face_orientation_2);
955 
956  // Only one periodic neighbor is allowed
957  Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
958  periodic_face_map.insert(periodic_faces);
959 
960  if (dim == 1)
961  {
962  if (cell_1->has_children())
963  {
964  if (cell_2->has_children())
965  {
966  update_periodic_face_map_recursively<dim, spacedim>(
967  cell_1->child(n_face_1),
968  cell_2->child(n_face_2),
969  n_face_1,
970  n_face_2,
971  orientation,
972  periodic_face_map);
973  }
974  else // only face_1 has children
975  {
976  update_periodic_face_map_recursively<dim, spacedim>(
977  cell_1->child(n_face_1),
978  cell_2,
979  n_face_1,
980  n_face_2,
981  orientation,
982  periodic_face_map);
983  }
984  }
985  }
986  else // dim == 2 || dim == 3
987  {
988  // A lookup table on how to go through the child cells depending on the
989  // orientation:
990  // see Documentation of GeometryInfo for details
991 
992  static const int lookup_table_2d[2][2] =
993  // flip:
994  {
995  {0, 1}, // false
996  {1, 0} // true
997  };
998 
999  static const int lookup_table_3d[2][2][2][4] =
1000  // orientation flip rotation
1001  {{{
1002  {0, 2, 1, 3}, // false false false
1003  {2, 3, 0, 1} // false false true
1004  },
1005  {
1006  {3, 1, 2, 0}, // false true false
1007  {1, 0, 3, 2} // false true true
1008  }},
1009  {{
1010  {0, 1, 2, 3}, // true false false
1011  {1, 3, 0, 2} // true false true
1012  },
1013  {
1014  {3, 2, 1, 0}, // true true false
1015  {2, 0, 3, 1} // true true true
1016  }}};
1017 
1018  if (cell_1->has_children())
1019  {
1020  if (cell_2->has_children())
1021  {
1022  // In the case that both faces have children, we loop over all
1023  // children and apply update_periodic_face_map_recursively
1024  // recursively:
1025 
1026  Assert(face_1->n_children() ==
1028  face_2->n_children() ==
1030  ExcNotImplemented());
1031 
1032  for (unsigned int i = 0;
1033  i < GeometryInfo<dim>::max_children_per_face;
1034  ++i)
1035  {
1036  // Lookup the index for the second face
1037  unsigned int j = 0;
1038  switch (dim)
1039  {
1040  case 2:
1041  j = lookup_table_2d[face_flip][i];
1042  break;
1043  case 3:
1044  j = lookup_table_3d[face_orientation][face_flip]
1045  [face_rotation][i];
1046  break;
1047  default:
1048  AssertThrow(false, ExcNotImplemented());
1049  }
1050 
1051  // find subcell ids that belong to the subface indices
1052  unsigned int child_cell_1 =
1054  cell_1->refinement_case(),
1055  n_face_1,
1056  i,
1057  cell_1->face_orientation(n_face_1),
1058  cell_1->face_flip(n_face_1),
1059  cell_1->face_rotation(n_face_1),
1060  face_1->refinement_case());
1061  unsigned int child_cell_2 =
1063  cell_2->refinement_case(),
1064  n_face_2,
1065  j,
1066  cell_2->face_orientation(n_face_2),
1067  cell_2->face_flip(n_face_2),
1068  cell_2->face_rotation(n_face_2),
1069  face_2->refinement_case());
1070 
1071  Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
1072  face_1->child(i),
1073  ExcInternalError());
1074  Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
1075  face_2->child(j),
1076  ExcInternalError());
1077 
1078  // precondition: subcell has the same orientation as cell
1079  // (so that the face numbers coincide) recursive call
1080  update_periodic_face_map_recursively<dim, spacedim>(
1081  cell_1->child(child_cell_1),
1082  cell_2->child(child_cell_2),
1083  n_face_1,
1084  n_face_2,
1085  orientation,
1086  periodic_face_map);
1087  }
1088  }
1089  else // only face_1 has children
1090  {
1091  for (unsigned int i = 0;
1092  i < GeometryInfo<dim>::max_children_per_face;
1093  ++i)
1094  {
1095  // find subcell ids that belong to the subface indices
1096  unsigned int child_cell_1 =
1098  cell_1->refinement_case(),
1099  n_face_1,
1100  i,
1101  cell_1->face_orientation(n_face_1),
1102  cell_1->face_flip(n_face_1),
1103  cell_1->face_rotation(n_face_1),
1104  face_1->refinement_case());
1105 
1106  // recursive call
1107  update_periodic_face_map_recursively<dim, spacedim>(
1108  cell_1->child(child_cell_1),
1109  cell_2,
1110  n_face_1,
1111  n_face_2,
1112  orientation,
1113  periodic_face_map);
1114  }
1115  }
1116  }
1117  }
1118  }
1119 
1120 
1121 } // end of anonymous namespace
1122 
1123 
1124 namespace internal
1125 {
1126  namespace TriangulationImplementation
1127  {
1128  // make sure that if in the following we
1129  // write Triangulation<dim,spacedim>
1130  // we mean the *class*
1131  // ::Triangulation, not the
1132  // enclosing namespace
1133  // internal::TriangulationImplementation
1135 
1141  int,
1142  << "Something went wrong when making cell " << arg1
1143  << ". Read the docs and the source code "
1144  << "for more information.");
1150  int,
1151  << "Something went wrong upon construction of cell "
1152  << arg1);
1163  int,
1164  << "Cell " << arg1
1165  << " has negative measure. This typically "
1166  << "indicates some distortion in the cell, or a mistakenly "
1167  << "swapped pair of vertices in the input to "
1168  << "Triangulation::create_triangulation().");
1177  int,
1178  int,
1179  int,
1180  << "Error while creating cell " << arg1
1181  << ": the vertex index " << arg2 << " must be between 0 and "
1182  << arg3 << ".");
1188  int,
1189  int,
1190  << "While trying to assign a boundary indicator to a line: "
1191  << "the line with end vertices " << arg1 << " and " << arg2
1192  << " does not exist.");
1198  int,
1199  int,
1200  int,
1201  int,
1202  << "While trying to assign a boundary indicator to a quad: "
1203  << "the quad with bounding lines " << arg1 << ", " << arg2
1204  << ", " << arg3 << ", " << arg4 << " does not exist.");
1211  int,
1212  int,
1214  << "The input data for creating a triangulation contained "
1215  << "information about a line with indices " << arg1 << " and " << arg2
1216  << " that is described to have boundary indicator "
1217  << static_cast<int>(arg3)
1218  << ". However, this is an internal line not located on the "
1219  << "boundary. You cannot assign a boundary indicator to it." << std::endl
1220  << std::endl
1221  << "If this happened at a place where you call "
1222  << "Triangulation::create_triangulation() yourself, you need "
1223  << "to check the SubCellData object you pass to this function."
1224  << std::endl
1225  << std::endl
1226  << "If this happened in a place where you are reading a mesh "
1227  << "from a file, then you need to investigate why such a line "
1228  << "ended up in the input file. A typical case is a geometry "
1229  << "that consisted of multiple parts and for which the mesh "
1230  << "generator program assumes that the interface between "
1231  << "two parts is a boundary when that isn't supposed to be "
1232  << "the case, or where the mesh generator simply assigns "
1233  << "'geometry indicators' to lines at the perimeter of "
1234  << "a part that are not supposed to be interpreted as "
1235  << "'boundary indicators'.");
1242  int,
1243  int,
1244  int,
1245  int,
1247  << "The input data for creating a triangulation contained "
1248  << "information about a quad with indices " << arg1 << ", " << arg2
1249  << ", " << arg3 << ", and " << arg4
1250  << " that is described to have boundary indicator "
1251  << static_cast<int>(arg5)
1252  << ". However, this is an internal quad not located on the "
1253  << "boundary. You cannot assign a boundary indicator to it." << std::endl
1254  << std::endl
1255  << "If this happened at a place where you call "
1256  << "Triangulation::create_triangulation() yourself, you need "
1257  << "to check the SubCellData object you pass to this function."
1258  << std::endl
1259  << std::endl
1260  << "If this happened in a place where you are reading a mesh "
1261  << "from a file, then you need to investigate why such a quad "
1262  << "ended up in the input file. A typical case is a geometry "
1263  << "that consisted of multiple parts and for which the mesh "
1264  << "generator program assumes that the interface between "
1265  << "two parts is a boundary when that isn't supposed to be "
1266  << "the case, or where the mesh generator simply assigns "
1267  << "'geometry indicators' to quads at the surface of "
1268  << "a part that are not supposed to be interpreted as "
1269  << "'boundary indicators'.");
1276  int,
1277  int,
1278  << "In SubCellData the line info of the line with vertex indices " << arg1
1279  << " and " << arg2 << " appears more than once. "
1280  << "This is not allowed.");
1287  int,
1288  int,
1289  std::string,
1290  << "In SubCellData the line info of the line with vertex indices " << arg1
1291  << " and " << arg2 << " appears multiple times with different (valid) "
1292  << arg3 << ". This is not allowed.");
1299  int,
1300  int,
1301  int,
1302  int,
1303  std::string,
1304  << "In SubCellData the quad info of the quad with line indices " << arg1
1305  << ", " << arg2 << ", " << arg3 << " and " << arg4
1306  << " appears multiple times with different (valid) " << arg5
1307  << ". This is not allowed.");
1308 
1405  {
1417  template <int dim, int spacedim>
1418  static void
1421  const unsigned int level_objects,
1423  {
1424  using line_iterator =
1426 
1427  number_cache.n_levels = 0;
1428  if (level_objects > 0)
1429  // find the last level on which there are used cells
1430  for (unsigned int level = 0; level < level_objects; ++level)
1431  if (triangulation.begin(level) != triangulation.end(level))
1432  number_cache.n_levels = level + 1;
1433 
1434  // no cells at all?
1435  Assert(number_cache.n_levels > 0, ExcInternalError());
1436 
1438  // update the number of lines on the different levels in the
1439  // cache
1440  number_cache.n_lines = 0;
1441  number_cache.n_active_lines = 0;
1442 
1443  // for 1d, lines have levels so take count the objects per
1444  // level and globally
1445  if (dim == 1)
1446  {
1447  number_cache.n_lines_level.resize(number_cache.n_levels);
1448  number_cache.n_active_lines_level.resize(number_cache.n_levels);
1449 
1450  for (unsigned int level = 0; level < number_cache.n_levels; ++level)
1451  {
1452  // count lines on this level
1453  number_cache.n_lines_level[level] = 0;
1454  number_cache.n_active_lines_level[level] = 0;
1455 
1456  line_iterator line = triangulation.begin_line(level),
1457  endc =
1458  (level == number_cache.n_levels - 1 ?
1459  line_iterator(triangulation.end_line()) :
1460  triangulation.begin_line(level + 1));
1461  for (; line != endc; ++line)
1462  {
1463  ++number_cache.n_lines_level[level];
1464  if (line->has_children() == false)
1465  ++number_cache.n_active_lines_level[level];
1466  }
1467 
1468  // update total number of lines
1469  number_cache.n_lines += number_cache.n_lines_level[level];
1470  number_cache.n_active_lines +=
1471  number_cache.n_active_lines_level[level];
1472  }
1473  }
1474  else
1475  {
1476  // for dim>1, there are no levels for lines
1477  number_cache.n_lines_level.clear();
1478  number_cache.n_active_lines_level.clear();
1479 
1480  line_iterator line = triangulation.begin_line(),
1481  endc = triangulation.end_line();
1482  for (; line != endc; ++line)
1483  {
1484  ++number_cache.n_lines;
1485  if (line->has_children() == false)
1486  ++number_cache.n_active_lines;
1487  }
1488  }
1489  }
1490 
1505  template <int dim, int spacedim>
1506  static void
1509  const unsigned int level_objects,
1511  {
1512  // update lines and n_levels in number_cache. since we don't
1513  // access any of these numbers, we can do this in the
1514  // background
1515  Threads::Task<void> update_lines = Threads::new_task(
1516  static_cast<
1517  void (*)(const Triangulation<dim, spacedim> &,
1518  const unsigned int,
1520  &compute_number_cache<dim, spacedim>),
1521  triangulation,
1522  level_objects,
1524  number_cache));
1525 
1526  using quad_iterator =
1528 
1530  // update the number of quads on the different levels in the
1531  // cache
1532  number_cache.n_quads = 0;
1533  number_cache.n_active_quads = 0;
1534 
1535  // for 2d, quads have levels so take count the objects per
1536  // level and globally
1537  if (dim == 2)
1538  {
1539  // count the number of levels; the function we called above
1540  // on a separate Task for lines also does this and puts it into
1541  // number_cache.n_levels, but this datum may not yet be
1542  // available as we call the function on a separate task
1543  unsigned int n_levels = 0;
1544  if (level_objects > 0)
1545  // find the last level on which there are used cells
1546  for (unsigned int level = 0; level < level_objects; ++level)
1547  if (triangulation.begin(level) != triangulation.end(level))
1548  n_levels = level + 1;
1549 
1550  number_cache.n_quads_level.resize(n_levels);
1551  number_cache.n_active_quads_level.resize(n_levels);
1552 
1553  for (unsigned int level = 0; level < n_levels; ++level)
1554  {
1555  // count quads on this level
1556  number_cache.n_quads_level[level] = 0;
1557  number_cache.n_active_quads_level[level] = 0;
1558 
1559  quad_iterator quad = triangulation.begin_quad(level),
1560  endc =
1561  (level == n_levels - 1 ?
1562  quad_iterator(triangulation.end_quad()) :
1563  triangulation.begin_quad(level + 1));
1564  for (; quad != endc; ++quad)
1565  {
1566  ++number_cache.n_quads_level[level];
1567  if (quad->has_children() == false)
1568  ++number_cache.n_active_quads_level[level];
1569  }
1570 
1571  // update total number of quads
1572  number_cache.n_quads += number_cache.n_quads_level[level];
1573  number_cache.n_active_quads +=
1574  number_cache.n_active_quads_level[level];
1575  }
1576  }
1577  else
1578  {
1579  // for dim>2, there are no levels for quads
1580  number_cache.n_quads_level.clear();
1581  number_cache.n_active_quads_level.clear();
1582 
1583  quad_iterator quad = triangulation.begin_quad(),
1584  endc = triangulation.end_quad();
1585  for (; quad != endc; ++quad)
1586  {
1587  ++number_cache.n_quads;
1588  if (quad->has_children() == false)
1589  ++number_cache.n_active_quads;
1590  }
1591  }
1592 
1593  // wait for the background computation for lines
1594  update_lines.join();
1595  }
1596 
1612  template <int dim, int spacedim>
1613  static void
1616  const unsigned int level_objects,
1618  {
1619  // update quads, lines and n_levels in number_cache. since we
1620  // don't access any of these numbers, we can do this in the
1621  // background
1622  Threads::Task<void> update_quads_and_lines = Threads::new_task(
1623  static_cast<
1624  void (*)(const Triangulation<dim, spacedim> &,
1625  const unsigned int,
1627  &compute_number_cache<dim, spacedim>),
1628  triangulation,
1629  level_objects,
1631  number_cache));
1632 
1633  using hex_iterator =
1635 
1637  // update the number of hexes on the different levels in the
1638  // cache
1639  number_cache.n_hexes = 0;
1640  number_cache.n_active_hexes = 0;
1641 
1642  // for 3d, hexes have levels so take count the objects per
1643  // level and globally
1644  if (dim == 3)
1645  {
1646  // count the number of levels; the function we called
1647  // above on a separate Task for quads (recursively, via
1648  // the lines function) also does this and puts it into
1649  // number_cache.n_levels, but this datum may not yet be
1650  // available as we call the function on a separate task
1651  unsigned int n_levels = 0;
1652  if (level_objects > 0)
1653  // find the last level on which there are used cells
1654  for (unsigned int level = 0; level < level_objects; ++level)
1655  if (triangulation.begin(level) != triangulation.end(level))
1656  n_levels = level + 1;
1657 
1658  number_cache.n_hexes_level.resize(n_levels);
1659  number_cache.n_active_hexes_level.resize(n_levels);
1660 
1661  for (unsigned int level = 0; level < n_levels; ++level)
1662  {
1663  // count hexes on this level
1664  number_cache.n_hexes_level[level] = 0;
1665  number_cache.n_active_hexes_level[level] = 0;
1666 
1667  hex_iterator hex = triangulation.begin_hex(level),
1668  endc = (level == n_levels - 1 ?
1669  hex_iterator(triangulation.end_hex()) :
1670  triangulation.begin_hex(level + 1));
1671  for (; hex != endc; ++hex)
1672  {
1673  ++number_cache.n_hexes_level[level];
1674  if (hex->has_children() == false)
1675  ++number_cache.n_active_hexes_level[level];
1676  }
1677 
1678  // update total number of hexes
1679  number_cache.n_hexes += number_cache.n_hexes_level[level];
1680  number_cache.n_active_hexes +=
1681  number_cache.n_active_hexes_level[level];
1682  }
1683  }
1684  else
1685  {
1686  // for dim>3, there are no levels for hexes
1687  number_cache.n_hexes_level.clear();
1688  number_cache.n_active_hexes_level.clear();
1689 
1690  hex_iterator hex = triangulation.begin_hex(),
1691  endc = triangulation.end_hex();
1692  for (; hex != endc; ++hex)
1693  {
1694  ++number_cache.n_hexes;
1695  if (hex->has_children() == false)
1696  ++number_cache.n_active_hexes;
1697  }
1698  }
1699 
1700  // wait for the background computation for quads
1701  update_quads_and_lines.join();
1702  }
1703 
1704 
1712  template <int spacedim>
1713  static void
1714  create_triangulation(const std::vector<Point<spacedim>> &v,
1715  const std::vector<CellData<1>> & cells,
1716  const SubCellData & /*subcelldata*/,
1718  {
1719  AssertThrow(v.size() > 0, ExcMessage("No vertices given"));
1720  AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
1721 
1722  // note: since no boundary
1723  // information can be given in one
1724  // dimension, the @p{subcelldata}
1725  // field is ignored. (only used for
1726  // error checking, which is a good
1727  // idea in any case)
1728  const unsigned int dim = 1;
1729 
1730  // copy vertices
1731  triangulation.vertices = v;
1732  triangulation.vertices_used = std::vector<bool>(v.size(), true);
1733 
1734  // Check that all cells have positive volume. This check is not run in
1735  // the codimension one or two cases since cell_measure is not
1736  // implemented for those.
1737 #ifndef _MSC_VER
1738  // TODO: The following code does not compile with MSVC. Find a way
1739  // around it
1740  if (dim == spacedim)
1741  {
1742  for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
1743  {
1744  // If we should check for distorted cells, then we permit them
1745  // to exist. If a cell has negative measure, then it must be
1746  // distorted (the converse is not necessarily true); hence
1747  // throw an exception if no such cells should exist.
1748  if (!triangulation.check_for_distorted_cells)
1749  {
1750  const double cell_measure =
1751  GridTools::cell_measure<1>(triangulation.vertices,
1752  cells[cell_no].vertices);
1753  AssertThrow(cell_measure > 0,
1754  ExcGridHasInvalidCell(cell_no));
1755  }
1756  }
1757  }
1758 #endif
1759 
1760 
1761  // store the indices of the lines
1762  // which are adjacent to a given
1763  // vertex
1764  std::vector<std::vector<int>> lines_at_vertex(v.size());
1765 
1766  // reserve enough space
1767  triangulation.levels.push_back(
1768  std::make_unique<internal::TriangulationImplementation::TriaLevel>(
1769  dim));
1770  triangulation.levels[0]->reserve_space(cells.size(), dim, spacedim);
1771  triangulation.levels[0]->cells.reserve_space(0, cells.size());
1772 
1773  // make up cells
1775  next_free_line = triangulation.begin_raw_line();
1776  for (unsigned int cell = 0; cell < cells.size(); ++cell)
1777  {
1778  while (next_free_line->used())
1779  ++next_free_line;
1780 
1781  next_free_line->set(
1783  cells[cell].vertices[0], cells[cell].vertices[1]));
1784  next_free_line->set_used_flag();
1785  next_free_line->set_material_id(cells[cell].material_id);
1786  next_free_line->set_manifold_id(cells[cell].manifold_id);
1787  next_free_line->clear_user_data();
1788  next_free_line->set_subdomain_id(0);
1789 
1790  // note that this cell is
1791  // adjacent to these vertices
1792  lines_at_vertex[cells[cell].vertices[0]].push_back(cell);
1793  lines_at_vertex[cells[cell].vertices[1]].push_back(cell);
1794  }
1795 
1796 
1797  // some security tests
1798  {
1799  unsigned int boundary_nodes = 0;
1800  for (const auto &line : lines_at_vertex)
1801  switch (line.size())
1802  {
1803  case 1:
1804  // this vertex has only
1805  // one adjacent line
1806  ++boundary_nodes;
1807  break;
1808  case 2:
1809  break;
1810  default:
1811  AssertThrow(
1812  false,
1813  ExcMessage(
1814  "You have a vertex in your triangulation "
1815  "at which more than two cells come together. "
1816  "(For one dimensional triangulation, cells are "
1817  "line segments.)"
1818  "\n\n"
1819  "This is not currently supported because the "
1820  "Triangulation class makes the assumption that "
1821  "every cell has zero or one neighbors behind "
1822  "each face (here, behind each vertex), but in your "
1823  "situation there would be more than one."
1824  "\n\n"
1825  "Support for this is not currently implemented. "
1826  "If you need to work with triangulations where "
1827  "more than two cells come together at a vertex, "
1828  "duplicate the vertices once per cell (i.e., put "
1829  "multiple vertices at the same physical location, "
1830  "but using different vertex indices for each) "
1831  "and then ensure continuity of the solution by "
1832  "explicitly creating constraints that the degrees "
1833  "of freedom at these vertices have the same "
1834  "value, using the AffineConstraints class."));
1835  }
1836  }
1837 
1838 
1839 
1840  // update neighborship info
1842  triangulation.begin_active_line();
1843  // for all lines
1844  for (; line != triangulation.end(); ++line)
1845  // for each of the two vertices
1846  for (const unsigned int vertex : GeometryInfo<dim>::vertex_indices())
1847  // if first cell adjacent to
1848  // this vertex is the present
1849  // one, then the neighbor is
1850  // the second adjacent cell and
1851  // vice versa
1852  if (lines_at_vertex[line->vertex_index(vertex)][0] == line->index())
1853  if (lines_at_vertex[line->vertex_index(vertex)].size() == 2)
1854  {
1856  neighbor(&triangulation,
1857  0, // level
1858  lines_at_vertex[line->vertex_index(vertex)][1]);
1859  line->set_neighbor(vertex, neighbor);
1860  }
1861  else
1862  // no second adjacent cell
1863  // entered -> cell at
1864  // boundary
1865  line->set_neighbor(vertex, triangulation.end());
1866  else
1867  // present line is not first
1868  // adjacent one -> first
1869  // adjacent one is neighbor
1870  {
1872  neighbor(&triangulation,
1873  0, // level
1874  lines_at_vertex[line->vertex_index(vertex)][0]);
1875  line->set_neighbor(vertex, neighbor);
1876  }
1877 
1878  // finally set the
1879  // vertex_to_boundary_id_map_1d
1880  // and vertex_to_manifold_id_map_1d
1881  // maps
1882  triangulation.vertex_to_boundary_id_map_1d->clear();
1883  triangulation.vertex_to_manifold_id_map_1d->clear();
1884  for (const auto &cell : triangulation.active_cell_iterators())
1885  for (auto f : GeometryInfo<dim>::face_indices())
1886  {
1887  (*triangulation.vertex_to_manifold_id_map_1d)
1888  [cell->face(f)->vertex_index()] = numbers::flat_manifold_id;
1889 
1890  if (cell->at_boundary(f))
1891  (*triangulation.vertex_to_boundary_id_map_1d)
1892  [cell->face(f)->vertex_index()] = f;
1893  }
1894  }
1895 
1896 
1904  template <int spacedim>
1905  static void
1906  create_triangulation(const std::vector<Point<spacedim>> &v,
1907  const std::vector<CellData<2>> & cells,
1908  const SubCellData & subcelldata,
1910  {
1911  AssertThrow(v.size() > 0, ExcMessage("No vertices given"));
1912  AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
1913 
1914  const unsigned int dim = 2;
1915 
1916  // copy vertices
1917  triangulation.vertices = v;
1918  triangulation.vertices_used = std::vector<bool>(v.size(), true);
1919 
1920  // Check that all cells have positive volume. This check is not run in
1921  // the codimension one or two cases since cell_measure is not
1922  // implemented for those.
1923 #ifndef _MSC_VER
1924  // TODO: The following code does not compile with MSVC. Find a way
1925  // around it
1926  if (dim == spacedim)
1927  {
1928  for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
1929  {
1930  // See the note in the 1D function on this if statement.
1931  if (!triangulation.check_for_distorted_cells)
1932  {
1933  const double cell_measure =
1934  GridTools::cell_measure<2>(triangulation.vertices,
1935  cells[cell_no].vertices);
1936  AssertThrow(cell_measure > 0,
1937  ExcGridHasInvalidCell(cell_no));
1938  }
1939  }
1940  }
1941 #endif
1942 
1943  // make up a list of the needed
1944  // lines each line is a pair of
1945  // vertices. The list is kept
1946  // sorted and it is guaranteed that
1947  // each line is inserted only once.
1948  // While the key of such an entry
1949  // is the pair of vertices, the
1950  // thing it points to is an
1951  // iterator pointing to the line
1952  // object itself. In the first run,
1953  // these iterators are all invalid
1954  // ones, but they are filled
1955  // afterwards
1956  std::map<std::pair<int, int>,
1958  needed_lines;
1959  for (unsigned int cell = 0; cell < cells.size(); ++cell)
1960  {
1961  for (const auto vertex : cells[cell].vertices)
1962  AssertThrow(vertex < triangulation.vertices.size(),
1963  ExcInvalidVertexIndex(cell,
1964  vertex,
1965  triangulation.vertices.size()));
1966 
1967  for (const unsigned int line : GeometryInfo<dim>::face_indices())
1968  {
1969  // given a line vertex number (0,1) on a specific line
1970  // we get the cell vertex number (0-4) through the
1971  // line_to_cell_vertices function
1972  std::pair<int, int> line_vertices(
1973  cells[cell].vertices[GeometryInfo<dim>::line_to_cell_vertices(
1974  line, 0)],
1975  cells[cell].vertices[GeometryInfo<dim>::line_to_cell_vertices(
1976  line, 1)]);
1977 
1978  // assert that the line was not already inserted in
1979  // reverse order. This happens in spite of the vertex
1980  // rotation above, if the sense of the cell was
1981  // incorrect.
1982  //
1983  // Here is what usually happened when this exception
1984  // is thrown: consider these two cells and the
1985  // vertices
1986  // 3---4---5
1987  // | | |
1988  // 0---1---2
1989  // If in the input vector the two cells are given with
1990  // vertices <0 1 3 4> and <4 1 5 2>, in the first cell
1991  // the middle line would have direction 1->4, while in
1992  // the second it would be 4->1. This will cause the
1993  // exception.
1994  AssertThrow(needed_lines.find(std::make_pair(
1995  line_vertices.second, line_vertices.first)) ==
1996  needed_lines.end(),
1997  ExcGridHasInvalidCell(cell));
1998 
1999  // insert line, with
2000  // invalid iterator if line
2001  // already exists, then
2002  // nothing bad happens here
2003  needed_lines[line_vertices] = triangulation.end_line();
2004  }
2005  }
2006 
2007 
2008  // check that every vertex has at
2009  // least two adjacent lines
2010  {
2011  std::vector<unsigned short int> vertex_touch_count(v.size(), 0);
2012  typename std::map<
2013  std::pair<int, int>,
2014  typename Triangulation<dim, spacedim>::line_iterator>::iterator i;
2015  for (i = needed_lines.begin(); i != needed_lines.end(); ++i)
2016  {
2017  // touch the vertices of
2018  // this line
2019  ++vertex_touch_count[i->first.first];
2020  ++vertex_touch_count[i->first.second];
2021  }
2022 
2023  // assert minimum touch count
2024  // is at least two. if not so,
2025  // then clean triangulation and
2026  // exit with an exception
2027  AssertThrow(*(std::min_element(vertex_touch_count.begin(),
2028  vertex_touch_count.end())) >= 2,
2029  ExcMessage(
2030  "During creation of a triangulation, a part of the "
2031  "algorithm encountered a vertex that is part of only "
2032  "a single adjacent line. However, in 2d, every vertex "
2033  "needs to be at least part of two lines."));
2034  }
2035 
2036  // reserve enough space
2037  triangulation.levels.push_back(
2038  std::make_unique<internal::TriangulationImplementation::TriaLevel>(
2039  dim));
2040  triangulation.faces =
2041  std::make_unique<internal::TriangulationImplementation::TriaFaces>(
2042  dim);
2043  triangulation.levels[0]->reserve_space(cells.size(), dim, spacedim);
2044  triangulation.faces->lines.reserve_space(0, needed_lines.size());
2045  triangulation.levels[0]->cells.reserve_space(0, cells.size());
2046 
2047  // make up lines
2048  {
2050  triangulation.begin_raw_line();
2051  typename std::map<
2052  std::pair<int, int>,
2053  typename Triangulation<dim, spacedim>::line_iterator>::iterator i;
2054  for (i = needed_lines.begin(); line != triangulation.end_line();
2055  ++line, ++i)
2056  {
2058  i->first.first, i->first.second));
2059  line->set_used_flag();
2060  line->clear_user_flag();
2061  line->clear_user_data();
2062  i->second = line;
2063  }
2064  }
2065 
2066 
2067  // store for each line index
2068  // the adjacent cells
2069  std::map<
2070  int,
2071  std::vector<typename Triangulation<dim, spacedim>::cell_iterator>>
2072  adjacent_cells;
2073 
2074  // finally make up cells
2075  {
2077  triangulation.begin_raw_quad();
2078  for (unsigned int c = 0; c < cells.size(); ++c, ++cell)
2079  {
2082  for (unsigned int line = 0;
2083  line < GeometryInfo<dim>::lines_per_cell;
2084  ++line)
2085  lines[line] = needed_lines[std::make_pair(
2087  line, 0)],
2088  cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
2089  line, 1)])];
2090 
2092  lines[0]->index(),
2093  lines[1]->index(),
2094  lines[2]->index(),
2095  lines[3]->index()));
2096 
2097  cell->set_used_flag();
2098  cell->set_material_id(cells[c].material_id);
2099  cell->set_manifold_id(cells[c].manifold_id);
2100  cell->clear_user_data();
2101  cell->set_subdomain_id(0);
2102 
2103  // note that this cell is
2104  // adjacent to the four
2105  // lines
2106  for (const auto &line : lines)
2107  adjacent_cells[line->index()].push_back(cell);
2108  }
2109  }
2110 
2111 
2112  for (typename Triangulation<dim, spacedim>::line_iterator line =
2113  triangulation.begin_line();
2114  line != triangulation.end_line();
2115  ++line)
2116  {
2117  const unsigned int n_adj_cells =
2118  adjacent_cells[line->index()].size();
2119 
2120  // assert that every line has one or two adjacent cells.
2121  // this has to be the case for 2d triangulations in 2d.
2122  // in higher dimensions, this may happen but is not
2123  // implemented
2124  if (spacedim == 2)
2125  {
2126  AssertThrow((n_adj_cells >= 1) && (n_adj_cells <= 2),
2127  ExcInternalError());
2128  }
2129  else
2130  {
2131  AssertThrow(
2132  (n_adj_cells >= 1) && (n_adj_cells <= 2),
2133  ExcMessage("You have a line in your triangulation at which "
2134  "more than two cells come together."
2135  "\n\n"
2136  "This is not currently supported because the "
2137  "Triangulation class makes the assumption that "
2138  "every cell has zero or one neighbors behind each "
2139  "face (here, behind each line), but in your "
2140  "situation there would be more than one."
2141  "\n\n"
2142  "Support for this is not currently implemented. "
2143  "If you need to work with triangulations where "
2144  "more than two cells come together at a line, "
2145  "duplicate the vertices once per cell (i.e., put "
2146  "multiple vertices at the same physical location, "
2147  "but using different vertex indices for each) "
2148  "and then ensure continuity of the solution by "
2149  "explicitly creating constraints that the degrees "
2150  "of freedom at these lines have the same "
2151  "value, using the AffineConstraints class."));
2152  }
2153 
2154  // if only one cell: line is at boundary -> give it the boundary
2155  // indicator zero by default
2156  line->set_boundary_id_internal(
2157  (n_adj_cells == 1) ? 0 : numbers::internal_face_boundary_id);
2158  line->set_manifold_id(numbers::flat_manifold_id);
2159  }
2160 
2161  // set boundary indicators where given
2162  for (const auto &subcell_line : subcelldata.boundary_lines)
2163  {
2165  std::pair<int, int> line_vertices(
2166  std::make_pair(subcell_line.vertices[0],
2167  subcell_line.vertices[1]));
2168  if (needed_lines.find(line_vertices) != needed_lines.end())
2169  // line found in this direction
2170  line = needed_lines[line_vertices];
2171  else
2172  {
2173  // look whether it exists in reverse direction
2174  std::swap(line_vertices.first, line_vertices.second);
2175  if (needed_lines.find(line_vertices) != needed_lines.end())
2176  line = needed_lines[line_vertices];
2177  else
2178  // line does not exist
2179  AssertThrow(false,
2180  ExcLineInexistant(line_vertices.first,
2181  line_vertices.second));
2182  }
2183 
2184  // assert that we only set boundary info once
2185  AssertThrow(!(line->boundary_id() != 0 &&
2186  line->boundary_id() !=
2188  ExcMultiplySetLineInfoOfLine(line_vertices.first,
2189  line_vertices.second));
2190 
2191  // assert that the manifold id is not yet set or consistent
2192  // with the previous id
2193  AssertThrow(line->manifold_id() == numbers::flat_manifold_id ||
2194  line->manifold_id() == subcell_line.manifold_id,
2195  ExcInconsistentLineInfoOfLine(line_vertices.first,
2196  line_vertices.second,
2197  "manifold ids"));
2198  line->set_manifold_id(subcell_line.manifold_id);
2199 
2200  // assert that only exterior lines are given a boundary
2201  // indicator
2202  if (subcell_line.boundary_id != numbers::internal_face_boundary_id)
2203  {
2204  AssertThrow(
2205  line->boundary_id() != numbers::internal_face_boundary_id,
2206  ExcInteriorLineCantBeBoundary(line->vertex_index(0),
2207  line->vertex_index(1),
2208  subcell_line.boundary_id));
2209  line->set_boundary_id_internal(subcell_line.boundary_id);
2210  }
2211  }
2212 
2213 
2214  // finally update neighborship info
2215  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
2216  triangulation.begin();
2217  cell != triangulation.end();
2218  ++cell)
2219  for (unsigned int side = 0; side < 4; ++side)
2220  if (adjacent_cells[cell->line(side)->index()][0] == cell)
2221  // first adjacent cell is
2222  // this one
2223  {
2224  if (adjacent_cells[cell->line(side)->index()].size() == 2)
2225  // there is another
2226  // adjacent cell
2227  cell->set_neighbor(
2228  side, adjacent_cells[cell->line(side)->index()][1]);
2229  }
2230  // first adjacent cell is not this
2231  // one, -> it must be the neighbor
2232  // we are looking for
2233  else
2234  cell->set_neighbor(side,
2235  adjacent_cells[cell->line(side)->index()][0]);
2236  }
2237 
2238 
2250  {
2251  inline bool
2255  {
2256  // here is room to
2257  // optimize the repeated
2258  // equality test of the
2259  // previous lines; the
2260  // compiler will probably
2261  // take care of most of
2262  // it anyway
2263  if ((q1.face(0) < q2.face(0)) ||
2264  ((q1.face(0) == q2.face(0)) && (q1.face(1) < q2.face(1))) ||
2265  ((q1.face(0) == q2.face(0)) && (q1.face(1) == q2.face(1)) &&
2266  (q1.face(2) < q2.face(2))) ||
2267  ((q1.face(0) == q2.face(0)) && (q1.face(1) == q2.face(1)) &&
2268  (q1.face(2) == q2.face(2)) && (q1.face(3) < q2.face(3))))
2269  return true;
2270  else
2271  return false;
2272  }
2273  };
2274 
2275 
2283  template <int spacedim>
2284  static void
2285  create_triangulation(const std::vector<Point<spacedim>> &v,
2286  const std::vector<CellData<3>> & cells,
2287  const SubCellData & subcelldata,
2289  {
2290  AssertThrow(v.size() > 0, ExcMessage("No vertices given"));
2291  AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
2292 
2293  const unsigned int dim = 3;
2294 
2295  // copy vertices
2296  triangulation.vertices = v;
2297  triangulation.vertices_used = std::vector<bool>(v.size(), true);
2298 
2299  // Check that all cells have positive volume.
2300 #ifndef _MSC_VER
2301  // TODO: The following code does not compile with MSVC. Find a way
2302  // around it
2303  for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2304  {
2305  // See the note in the 1D function on this if statement.
2306  if (!triangulation.check_for_distorted_cells)
2307  {
2308  const double cell_measure =
2309  GridTools::cell_measure<3>(triangulation.vertices,
2310  cells[cell_no].vertices);
2311  AssertThrow(cell_measure > 0, ExcGridHasInvalidCell(cell_no));
2312  }
2313  }
2314 #endif
2315 
2317  // first set up some collections of data
2318  //
2319  // make up a list of the needed
2320  // lines
2321  //
2322  // each line is a pair of
2323  // vertices. The list is kept
2324  // sorted and it is guaranteed that
2325  // each line is inserted only once.
2326  // While the key of such an entry
2327  // is the pair of vertices, the
2328  // thing it points to is an
2329  // iterator pointing to the line
2330  // object itself. In the first run,
2331  // these iterators are all invalid
2332  // ones, but they are filled
2333  // afterwards same applies for the
2334  // quads
2335  typename std::map<std::pair<int, int>,
2337  needed_lines;
2338  for (unsigned int cell = 0; cell < cells.size(); ++cell)
2339  {
2340  // check whether vertex indices
2341  // are valid ones
2342  for (const auto vertex : cells[cell].vertices)
2343  AssertThrow(vertex < triangulation.vertices.size(),
2344  ExcInvalidVertexIndex(cell,
2345  vertex,
2346  triangulation.vertices.size()));
2347 
2348  for (unsigned int line = 0;
2349  line < GeometryInfo<dim>::lines_per_cell;
2350  ++line)
2351  {
2352  // given a line vertex number
2353  // (0,1) on a specific line we
2354  // get the cell vertex number
2355  // (0-7) through the
2356  // line_to_cell_vertices
2357  // function
2358  std::pair<int, int> line_vertices(
2359  cells[cell].vertices[GeometryInfo<dim>::line_to_cell_vertices(
2360  line, 0)],
2361  cells[cell].vertices[GeometryInfo<dim>::line_to_cell_vertices(
2362  line, 1)]);
2363 
2364  // if that line was already inserted
2365  // in reverse order do nothing, else
2366  // insert the line
2367  if ((needed_lines.find(std::make_pair(line_vertices.second,
2368  line_vertices.first)) ==
2369  needed_lines.end()))
2370  {
2371  // insert line, with
2372  // invalid iterator. if line
2373  // already exists, then
2374  // nothing bad happens here
2375  needed_lines[line_vertices] = triangulation.end_line();
2376  }
2377  }
2378  }
2379 
2380 
2382  // now for some sanity-checks:
2383  //
2384  // check that every vertex has at
2385  // least tree adjacent lines
2386  {
2387  std::vector<unsigned short int> vertex_touch_count(v.size(), 0);
2388  typename std::map<
2389  std::pair<int, int>,
2390  typename Triangulation<dim, spacedim>::line_iterator>::iterator i;
2391  for (i = needed_lines.begin(); i != needed_lines.end(); ++i)
2392  {
2393  // touch the vertices of
2394  // this line
2395  ++vertex_touch_count[i->first.first];
2396  ++vertex_touch_count[i->first.second];
2397  }
2398 
2399  // assert minimum touch count
2400  // is at least three. if not so,
2401  // then clean triangulation and
2402  // exit with an exception
2403  AssertThrow(
2404  *(std::min_element(vertex_touch_count.begin(),
2405  vertex_touch_count.end())) >= 3,
2406  ExcMessage(
2407  "During creation of a triangulation, a part of the "
2408  "algorithm encountered a vertex that is part of only "
2409  "one or two adjacent lines. However, in 3d, every vertex "
2410  "needs to be at least part of three lines."));
2411  }
2412 
2413 
2415  // actually set up data structures
2416  // for the lines
2417  // reserve enough space
2418  triangulation.levels.push_back(
2419  std::make_unique<internal::TriangulationImplementation::TriaLevel>(
2420  dim));
2421  triangulation.faces =
2422  std::make_unique<internal::TriangulationImplementation::TriaFaces>(
2423  dim);
2424  triangulation.levels[0]->reserve_space(cells.size(), dim, spacedim);
2425  triangulation.faces->lines.reserve_space(0, needed_lines.size());
2426 
2427  // make up lines
2428  {
2430  triangulation.begin_raw_line();
2431  typename std::map<
2432  std::pair<int, int>,
2433  typename Triangulation<dim, spacedim>::line_iterator>::iterator i;
2434  for (i = needed_lines.begin(); line != triangulation.end_line();
2435  ++line, ++i)
2436  {
2438  i->first.first, i->first.second));
2439  line->set_used_flag();
2440  line->clear_user_flag();
2441  line->clear_user_data();
2442 
2443  // now set the iterator for
2444  // this line
2445  i->second = line;
2446  }
2447  }
2448 
2449 
2451  // make up the quads of this triangulation
2452  //
2453  // same thing: the iterators are
2454  // set to the invalid value at
2455  // first, we only collect the data
2456  // now
2457 
2458  // the bool array stores, whether the lines
2459  // are in the standard orientation or not
2460 
2461  // note that QuadComparator is a
2462  // class declared and defined in
2463  // this file
2464  std::map<internal::TriangulationImplementation::TriaObject<2>,
2465  std::pair<typename Triangulation<dim, spacedim>::quad_iterator,
2466  std::array<bool, GeometryInfo<dim>::lines_per_face>>,
2468  needed_quads;
2469  for (const auto &cell : cells)
2470  {
2471  // the faces are quads which
2472  // consist of four numbers
2473  // denoting the index of the
2474  // four lines bounding the
2475  // quad. we can get this index
2476  // by asking @p{needed_lines}
2477  // for an iterator to this
2478  // line, dereferencing it and
2479  // thus return an iterator into
2480  // the @p{lines} array of the
2481  // triangulation, which is
2482  // already set up. we can then
2483  // ask this iterator for its
2484  // index within the present
2485  // level (the level is zero, of
2486  // course)
2487  //
2488  // to make things easier, we
2489  // don't create the lines
2490  // (pairs of their vertex
2491  // indices) in place, but
2492  // before they are really
2493  // needed.
2494  std::pair<int, int> line_list[GeometryInfo<dim>::lines_per_cell],
2495  inverse_line_list[GeometryInfo<dim>::lines_per_cell];
2496  unsigned int face_line_list[GeometryInfo<dim>::lines_per_face];
2497  std::array<bool, GeometryInfo<dim>::lines_per_face> orientation;
2498 
2499  for (unsigned int line = 0;
2500  line < GeometryInfo<dim>::lines_per_cell;
2501  ++line)
2502  {
2503  line_list[line] = std::pair<int, int>(
2504  cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(line,
2505  0)],
2506  cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(line,
2507  1)]);
2508  inverse_line_list[line] = std::pair<int, int>(
2509  cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(line,
2510  1)],
2511  cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(line,
2512  0)]);
2513  }
2514 
2515  for (const unsigned int face : GeometryInfo<dim>::face_indices())
2516  {
2517  // set up a list of the lines to be
2518  // used for this face. check the
2519  // direction for each line
2520  //
2521  // given a face line number (0-3) on
2522  // a specific face we get the cell
2523  // line number (0-11) through the
2524  // face_to_cell_lines function
2525  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_face;
2526  ++l)
2527  if (needed_lines.find(
2528  inverse_line_list[GeometryInfo<dim>::face_to_cell_lines(
2529  face, l)]) == needed_lines.end())
2530  {
2531  face_line_list[l] =
2532  needed_lines[line_list[GeometryInfo<
2533  dim>::face_to_cell_lines(face, l)]]
2534  ->index();
2535  orientation[l] = true;
2536  }
2537  else
2538  {
2539  face_line_list[l] =
2540  needed_lines[inverse_line_list[GeometryInfo<
2541  dim>::face_to_cell_lines(face, l)]]
2542  ->index();
2543  orientation[l] = false;
2544  }
2545 
2546 
2548  face_line_list[0],
2549  face_line_list[1],
2550  face_line_list[2],
2551  face_line_list[3]);
2552 
2553  // insert quad, with
2554  // invalid iterator
2555  //
2556  // if quad already exists,
2557  // then nothing bad happens
2558  // here, as this will then
2559  // simply become an
2560  // interior face of the
2561  // triangulation. however,
2562  // we will run into major
2563  // trouble if the face was
2564  // already inserted in the
2565  // opposite
2566  // direction. there are
2567  // really only two
2568  // orientations for a face
2569  // to be in, since the edge
2570  // directions are already
2571  // set. thus, vertex 0 is
2572  // the one from which two
2573  // edges originate, and
2574  // vertex 3 is the one to
2575  // which they converge. we
2576  // are then left with
2577  // orientations 0-1-2-3 and
2578  // 2-3-0-1 for the order of
2579  // lines. the
2580  // corresponding quad can
2581  // be easily constructed by
2582  // exchanging lines. we do
2583  // so here, just to check
2584  // that that flipped quad
2585  // isn't already in the
2586  // triangulation. if it is,
2587  // then don't insert the
2588  // new one and instead
2589  // later set the
2590  // face_orientation flag
2592  test_quad_1(quad.face(2),
2593  quad.face(3),
2594  quad.face(0),
2595  quad.face(
2596  1)), // face_orientation=false, face_flip=false,
2597  // face_rotation=false
2598  test_quad_2(quad.face(0),
2599  quad.face(1),
2600  quad.face(3),
2601  quad.face(
2602  2)), // face_orientation=false, face_flip=false,
2603  // face_rotation=true
2604  test_quad_3(quad.face(3),
2605  quad.face(2),
2606  quad.face(1),
2607  quad.face(
2608  0)), // face_orientation=false, face_flip=true,
2609  // face_rotation=false
2610  test_quad_4(quad.face(1),
2611  quad.face(0),
2612  quad.face(2),
2613  quad.face(
2614  3)), // face_orientation=false, face_flip=true,
2615  // face_rotation=true
2616  test_quad_5(quad.face(2),
2617  quad.face(3),
2618  quad.face(1),
2619  quad.face(
2620  0)), // face_orientation=true, face_flip=false,
2621  // face_rotation=true
2622  test_quad_6(quad.face(1),
2623  quad.face(0),
2624  quad.face(3),
2625  quad.face(
2626  2)), // face_orientation=true, face_flip=true,
2627  // face_rotation=false
2628  test_quad_7(quad.face(3),
2629  quad.face(2),
2630  quad.face(0),
2631  quad.face(
2632  1)); // face_orientation=true, face_flip=true,
2633  // face_rotation=true
2634  if (needed_quads.find(test_quad_1) == needed_quads.end() &&
2635  needed_quads.find(test_quad_2) == needed_quads.end() &&
2636  needed_quads.find(test_quad_3) == needed_quads.end() &&
2637  needed_quads.find(test_quad_4) == needed_quads.end() &&
2638  needed_quads.find(test_quad_5) == needed_quads.end() &&
2639  needed_quads.find(test_quad_6) == needed_quads.end() &&
2640  needed_quads.find(test_quad_7) == needed_quads.end())
2641  needed_quads[quad] =
2642  std::make_pair(triangulation.end_quad(), orientation);
2643  }
2644  }
2645 
2646 
2648  // enter the resulting quads into
2649  // the arrays of the Triangulation
2650  //
2651  // first reserve enough space
2652  triangulation.faces->reserve_space(0, needed_quads.size());
2653  triangulation.faces->quads.reserve_space(0, needed_quads.size());
2654 
2655  {
2657  triangulation.begin_raw_quad();
2658  typename std::map<
2660  std::pair<typename Triangulation<dim, spacedim>::quad_iterator,
2661  std::array<bool, GeometryInfo<dim>::lines_per_face>>,
2662  QuadComparator>::iterator q;
2663  for (q = needed_quads.begin(); quad != triangulation.end_quad();
2664  ++quad, ++q)
2665  {
2666  quad->set(q->first);
2667  quad->set_used_flag();
2668  quad->clear_user_flag();
2669  quad->clear_user_data();
2670  // set the line orientation
2671  quad->set_line_orientation(0, q->second.second[0]);
2672  quad->set_line_orientation(1, q->second.second[1]);
2673  quad->set_line_orientation(2, q->second.second[2]);
2674  quad->set_line_orientation(3, q->second.second[3]);
2675 
2676 
2677  // now set the iterator for
2678  // this quad
2679  q->second.first = quad;
2680  }
2681  }
2682 
2684  // finally create the cells
2685  triangulation.levels[0]->cells.reserve_space(cells.size());
2686 
2687  // store for each quad index the
2688  // adjacent cells
2689  std::map<
2690  int,
2691  std::vector<typename Triangulation<dim, spacedim>::cell_iterator>>
2692  adjacent_cells;
2693 
2694  // finally make up cells
2695  {
2697  triangulation.begin_raw_hex();
2698  for (unsigned int c = 0; c < cells.size(); ++c, ++cell)
2699  {
2700  // first find for each of
2701  // the cells the quad
2702  // iterator of the
2703  // respective faces.
2704  //
2705  // to this end, set up the
2706  // lines of this cell and
2707  // find the quads that are
2708  // bounded by these lines;
2709  // these are then the faces
2710  // of the present cell
2711  std::pair<int, int> line_list[GeometryInfo<dim>::lines_per_cell],
2712  inverse_line_list[GeometryInfo<dim>::lines_per_cell];
2713  unsigned int face_line_list[4];
2714  for (unsigned int line = 0;
2715  line < GeometryInfo<dim>::lines_per_cell;
2716  ++line)
2717  {
2718  line_list[line] = std::make_pair(
2720  line, 0)],
2722  line, 1)]);
2723  inverse_line_list[line] = std::pair<int, int>(
2724  cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
2725  line, 1)],
2726  cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
2727  line, 0)]);
2728  }
2729 
2730  // get the iterators
2731  // corresponding to the
2732  // faces. also store
2733  // whether they are
2734  // reversed or not
2736  face_iterator[GeometryInfo<dim>::faces_per_cell];
2737  bool face_orientation[GeometryInfo<dim>::faces_per_cell];
2738  bool face_flip[GeometryInfo<dim>::faces_per_cell];
2739  bool face_rotation[GeometryInfo<dim>::faces_per_cell];
2740  for (const unsigned int face : GeometryInfo<dim>::face_indices())
2741  {
2742  for (unsigned int l = 0;
2743  l < GeometryInfo<dim>::lines_per_face;
2744  ++l)
2745  if (needed_lines.find(inverse_line_list[GeometryInfo<
2746  dim>::face_to_cell_lines(face, l)]) ==
2747  needed_lines.end())
2748  face_line_list[l] =
2749  needed_lines[line_list[GeometryInfo<
2750  dim>::face_to_cell_lines(face, l)]]
2751  ->index();
2752  else
2753  face_line_list[l] =
2754  needed_lines[inverse_line_list[GeometryInfo<
2755  dim>::face_to_cell_lines(face, l)]]
2756  ->index();
2757 
2759  face_line_list[0],
2760  face_line_list[1],
2761  face_line_list[2],
2762  face_line_list[3]);
2763 
2764  if (needed_quads.find(quad) != needed_quads.end())
2765  {
2766  // face is in standard
2767  // orientation (and not
2768  // flipped or rotated). this
2769  // must be true for at least
2770  // one of the two cells
2771  // containing this face
2772  // (i.e. for the cell which
2773  // originally inserted the
2774  // face)
2775  face_iterator[face] = needed_quads[quad].first;
2776  face_orientation[face] = true;
2777  face_flip[face] = false;
2778  face_rotation[face] = false;
2779  }
2780  else
2781  {
2782  // face must be available in
2783  // reverse order
2784  // then. construct all
2785  // possibilities and check
2786  // them one after the other
2788  test_quad_1(
2789  quad.face(2),
2790  quad.face(3),
2791  quad.face(0),
2792  quad.face(1)), // face_orientation=false,
2793  // face_flip=false, face_rotation=false
2794  test_quad_2(
2795  quad.face(0),
2796  quad.face(1),
2797  quad.face(3),
2798  quad.face(2)), // face_orientation=false,
2799  // face_flip=false, face_rotation=true
2800  test_quad_3(
2801  quad.face(3),
2802  quad.face(2),
2803  quad.face(1),
2804  quad.face(0)), // face_orientation=false,
2805  // face_flip=true, face_rotation=false
2806  test_quad_4(quad.face(1),
2807  quad.face(0),
2808  quad.face(2),
2809  quad.face(
2810  3)), // face_orientation=false,
2811  // face_flip=true, face_rotation=true
2812  test_quad_5(
2813  quad.face(2),
2814  quad.face(3),
2815  quad.face(1),
2816  quad.face(0)), // face_orientation=true,
2817  // face_flip=false, face_rotation=true
2818  test_quad_6(
2819  quad.face(1),
2820  quad.face(0),
2821  quad.face(3),
2822  quad.face(2)), // face_orientation=true,
2823  // face_flip=true, face_rotation=false
2824  test_quad_7(quad.face(3),
2825  quad.face(2),
2826  quad.face(0),
2827  quad.face(
2828  1)); // face_orientation=true,
2829  // face_flip=true, face_rotation=true
2830  if (needed_quads.find(test_quad_1) != needed_quads.end())
2831  {
2832  face_iterator[face] = needed_quads[test_quad_1].first;
2833  face_orientation[face] = false;
2834  face_flip[face] = false;
2835  face_rotation[face] = false;
2836  }
2837  else if (needed_quads.find(test_quad_2) !=
2838  needed_quads.end())
2839  {
2840  face_iterator[face] = needed_quads[test_quad_2].first;
2841  face_orientation[face] = false;
2842  face_flip[face] = false;
2843  face_rotation[face] = true;
2844  }
2845  else if (needed_quads.find(test_quad_3) !=
2846  needed_quads.end())
2847  {
2848  face_iterator[face] = needed_quads[test_quad_3].first;
2849  face_orientation[face] = false;
2850  face_flip[face] = true;
2851  face_rotation[face] = false;
2852  }
2853  else if (needed_quads.find(test_quad_4) !=
2854  needed_quads.end())
2855  {
2856  face_iterator[face] = needed_quads[test_quad_4].first;
2857  face_orientation[face] = false;
2858  face_flip[face] = true;
2859  face_rotation[face] = true;
2860  }
2861  else if (needed_quads.find(test_quad_5) !=
2862  needed_quads.end())
2863  {
2864  face_iterator[face] = needed_quads[test_quad_5].first;
2865  face_orientation[face] = true;
2866  face_flip[face] = false;
2867  face_rotation[face] = true;
2868  }
2869  else if (needed_quads.find(test_quad_6) !=
2870  needed_quads.end())
2871  {
2872  face_iterator[face] = needed_quads[test_quad_6].first;
2873  face_orientation[face] = true;
2874  face_flip[face] = true;
2875  face_rotation[face] = false;
2876  }
2877  else if (needed_quads.find(test_quad_7) !=
2878  needed_quads.end())
2879  {
2880  face_iterator[face] = needed_quads[test_quad_7].first;
2881  face_orientation[face] = true;
2882  face_flip[face] = true;
2883  face_rotation[face] = true;
2884  }
2885 
2886  else
2887  // we didn't find the
2888  // face in any direction,
2889  // so something went
2890  // wrong above
2891  Assert(false, ExcInternalError());
2892  }
2893  } // for all faces
2894 
2895  // make the cell out of
2896  // these iterators
2898  face_iterator[0]->index(),
2899  face_iterator[1]->index(),
2900  face_iterator[2]->index(),
2901  face_iterator[3]->index(),
2902  face_iterator[4]->index(),
2903  face_iterator[5]->index()));
2904 
2905  cell->set_used_flag();
2906  cell->set_material_id(cells[c].material_id);
2907  cell->set_manifold_id(cells[c].manifold_id);
2908  cell->clear_user_flag();
2909  cell->clear_user_data();
2910  cell->set_subdomain_id(0);
2911 
2912  // set orientation flag for
2913  // each of the faces
2914  for (const unsigned int quad : GeometryInfo<dim>::face_indices())
2915  {
2916  cell->set_face_orientation(quad, face_orientation[quad]);
2917  cell->set_face_flip(quad, face_flip[quad]);
2918  cell->set_face_rotation(quad, face_rotation[quad]);
2919  }
2920 
2921 
2922  // note that this cell is
2923  // adjacent to the six
2924  // quads
2925  for (const auto &quad : face_iterator)
2926  adjacent_cells[quad->index()].push_back(cell);
2927 
2928 #ifdef DEBUG
2929  // make some checks on the
2930  // lines and their
2931  // ordering
2932 
2933  // first map all cell lines
2934  // to the two face lines
2935  // which should
2936  // coincide. all face lines
2937  // are included with a cell
2938  // line number (0-11)
2939  // key. At the end all keys
2940  // will be included twice
2941  // (for each of the two
2942  // coinciding lines once)
2943  std::multimap<unsigned int, std::pair<unsigned int, unsigned int>>
2944  cell_to_face_lines;
2945  for (const unsigned int face : GeometryInfo<dim>::face_indices())
2946  for (unsigned int line = 0;
2947  line < GeometryInfo<dim>::lines_per_face;
2948  ++line)
2949  cell_to_face_lines.insert(
2950  std::pair<unsigned int,
2951  std::pair<unsigned int, unsigned int>>(
2953  std::pair<unsigned int, unsigned int>(face, line)));
2954  std::multimap<unsigned int,
2955  std::pair<unsigned int, unsigned int>>::
2956  const_iterator map_iter = cell_to_face_lines.begin();
2957 
2958  for (; map_iter != cell_to_face_lines.end(); ++map_iter)
2959  {
2960  const unsigned int cell_line = map_iter->first;
2961  const unsigned int face1 = map_iter->second.first;
2962  const unsigned int line1 = map_iter->second.second;
2963  ++map_iter;
2964  Assert(map_iter != cell_to_face_lines.end(),
2966  Assert(map_iter->first == cell_line,
2968  const unsigned int face2 = map_iter->second.first;
2969  const unsigned int line2 = map_iter->second.second;
2970 
2971  // check that the pair
2972  // of lines really
2973  // coincide. Take care
2974  // about the face
2975  // orientation;
2976  Assert(face_iterator[face1]->line(
2978  line1,
2979  face_orientation[face1],
2980  face_flip[face1],
2981  face_rotation[face1])) ==
2982  face_iterator[face2]->line(
2984  line2,
2985  face_orientation[face2],
2986  face_flip[face2],
2987  face_rotation[face2])),
2989  }
2990 #endif
2991  }
2992  }
2993 
2994 
2996  // find those quads which are at the
2997  // boundary and mark them appropriately
2998  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
2999  triangulation.begin_quad();
3000  quad != triangulation.end_quad();
3001  ++quad)
3002  {
3003  const unsigned int n_adj_cells =
3004  adjacent_cells[quad->index()].size();
3005  // assert that every quad has
3006  // one or two adjacent cells
3007  AssertThrow((n_adj_cells >= 1) && (n_adj_cells <= 2),
3008  ExcInternalError());
3009 
3010  // if only one cell: quad is at boundary -> give it the boundary
3011  // indicator zero by default
3012  quad->set_boundary_id_internal(
3013  (n_adj_cells == 1) ? 0 : numbers::internal_face_boundary_id);
3014 
3015  // Manifold ids are set independently of where they are
3016  quad->set_manifold_id(numbers::flat_manifold_id);
3017  }
3018 
3020  // next find those lines which are at
3021  // the boundary and mark all others as
3022  // interior ones
3023  //
3024  // for this: first mark all lines as interior. use this loop
3025  // to also set all manifold ids of all lines
3026  for (typename Triangulation<dim, spacedim>::line_iterator line =
3027  triangulation.begin_line();
3028  line != triangulation.end_line();
3029  ++line)
3030  {
3031  line->set_boundary_id_internal(numbers::internal_face_boundary_id);
3032  line->set_manifold_id(numbers::flat_manifold_id);
3033  }
3034 
3035  // next reset all lines bounding
3036  // boundary quads as on the
3037  // boundary also. note that since
3038  // we are in 3d, there are cases
3039  // where one or more lines of a
3040  // quad that is not on the
3041  // boundary, are actually boundary
3042  // lines. they will not be marked
3043  // when visiting this
3044  // face. however, since we do not
3045  // support dim-2 dimensional
3046  // boundaries (i.e. internal lines
3047  // constituting boundaries), every
3048  // such line is also part of a face
3049  // that is actually on the
3050  // boundary, so sooner or later we
3051  // get to mark that line for being
3052  // on the boundary
3053  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
3054  triangulation.begin_quad();
3055  quad != triangulation.end_quad();
3056  ++quad)
3057  if (quad->at_boundary())
3058  {
3059  for (unsigned int l = 0; l < 4; ++l)
3060  {
3062  quad->line(l);
3063  line->set_boundary_id_internal(0);
3064  }
3065  }
3066 
3068  // now set boundary indicators
3069  // where given
3070  //
3071  // first do so for lines
3072  for (const auto &subcell_line : subcelldata.boundary_lines)
3073  {
3075  std::pair<int, int> line_vertices(
3076  std::make_pair(subcell_line.vertices[0],
3077  subcell_line.vertices[1]));
3078  if (needed_lines.find(line_vertices) != needed_lines.end())
3079  // line found in this
3080  // direction
3081  line = needed_lines[line_vertices];
3082 
3083  else
3084  {
3085  // look whether it exists in
3086  // reverse direction
3087  std::swap(line_vertices.first, line_vertices.second);
3088  if (needed_lines.find(line_vertices) != needed_lines.end())
3089  line = needed_lines[line_vertices];
3090  else
3091  // line does not exist
3092  AssertThrow(false,
3093  ExcLineInexistant(line_vertices.first,
3094  line_vertices.second));
3095  }
3096  // Only exterior lines can be given a boundary indicator
3097  if (line->at_boundary())
3098  {
3099  // make sure that we don't attempt to reset the boundary
3100  // indicator to a different than the previously set value
3101  AssertThrow(line->boundary_id() == 0 ||
3102  line->boundary_id() == subcell_line.boundary_id,
3103  ExcInconsistentLineInfoOfLine(line_vertices.first,
3104  line_vertices.second,
3105  "boundary ids"));
3106  // If the boundary id provided in subcell_line
3107  // is anything other than the default
3108  // (internal_face_boundary_id), then set it in the new
3109  // triangulation.
3110  if (subcell_line.boundary_id !=
3112  line->set_boundary_id(subcell_line.boundary_id);
3113  }
3114  // Set manifold id if given
3115  AssertThrow(line->manifold_id() == numbers::flat_manifold_id ||
3116  line->manifold_id() == subcell_line.manifold_id,
3117  ExcInconsistentLineInfoOfLine(line_vertices.first,
3118  line_vertices.second,
3119  "manifold ids"));
3120  line->set_manifold_id(subcell_line.manifold_id);
3121  }
3122 
3123 
3124  // now go on with the faces
3125  for (const auto &subcell_quad : subcelldata.boundary_quads)
3126  {
3129 
3130  // first find the lines that
3131  // are made up of the given
3132  // vertices, then build up a
3133  // quad from these lines
3134  // finally use the find
3135  // function of the map template
3136  // to find the quad
3137  for (unsigned int i = 0; i < 4; ++i)
3138  {
3139  std::pair<int, int> line_vertices(
3140  subcell_quad
3142  0)],
3143  subcell_quad
3145  1)]);
3146 
3147  // check whether line
3148  // already exists
3149  if (needed_lines.find(line_vertices) != needed_lines.end())
3150  line[i] = needed_lines[line_vertices];
3151  else
3152  // look whether it exists
3153  // in reverse direction
3154  {
3155  std::swap(line_vertices.first, line_vertices.second);
3156  if (needed_lines.find(line_vertices) != needed_lines.end())
3157  line[i] = needed_lines[line_vertices];
3158  else
3159  // line does
3160  // not exist
3161  AssertThrow(false,
3162  ExcLineInexistant(line_vertices.first,
3163  line_vertices.second));
3164  }
3165  }
3166 
3167 
3168  // Set up 2 quads that are
3169  // built up from the lines for
3170  // reasons of comparison to
3171  // needed_quads. The second
3172  // quad is the reversed version
3173  // of the first quad in order
3174  // find the quad regardless of
3175  // its orientation. This is
3176  // introduced for convenience
3177  // and because boundary quad
3178  // orientation does not carry
3179  // any information.
3181  line[0]->index(),
3182  line[1]->index(),
3183  line[2]->index(),
3184  line[3]->index());
3186  line[2]->index(),
3187  line[3]->index(),
3188  line[0]->index(),
3189  line[1]->index());
3190 
3191  // try to find the quad with
3192  // lines situated as
3193  // constructed above. if it
3194  // could not be found, rotate
3195  // the boundary lines 3 times
3196  // until it is found or it does
3197  // not exist.
3198 
3199  // mapping from counterclock to
3200  // lexicographic ordering of
3201  // quad lines
3202  static const unsigned int lex2cclock[4] = {3, 1, 0, 2};
3203  // copy lines from
3204  // lexicographic to
3205  // counterclock ordering, as
3206  // rotation is much simpler in
3207  // counterclock ordering
3209  line_counterclock[4];
3210  for (unsigned int i = 0; i < 4; ++i)
3211  line_counterclock[lex2cclock[i]] = line[i];
3212  unsigned int n_rotations = 0;
3213  bool not_found_quad_1;
3214  while ((not_found_quad_1 = (needed_quads.find(quad_compare_1) ==
3215  needed_quads.end())) &&
3216  (needed_quads.find(quad_compare_2) == needed_quads.end()) &&
3217  (n_rotations < 4))
3218  {
3219  // use the rotate defined
3220  // in <algorithms>
3221  std::rotate(line_counterclock,
3222  line_counterclock + 1,
3223  line_counterclock + 4);
3224  // update the quads with
3225  // rotated lines (i runs in
3226  // lexicographic ordering)
3227  for (unsigned int i = 0; i < 4; ++i)
3228  {
3229  quad_compare_1.set_face(
3230  i, line_counterclock[lex2cclock[i]]->index());
3231  quad_compare_2.set_face(
3232  (i + 2) % 4, line_counterclock[lex2cclock[i]]->index());
3233  }
3234 
3235  ++n_rotations;
3236  }
3237 
3238  AssertThrow(n_rotations != 4,
3239  ExcQuadInexistant(line[0]->index(),
3240  line[1]->index(),
3241  line[2]->index(),
3242  line[3]->index()));
3243 
3244  if (not_found_quad_1)
3245  quad = needed_quads[quad_compare_2].first;
3246  else
3247  quad = needed_quads[quad_compare_1].first;
3248 
3249  // check whether this face is
3250  // really an exterior one
3251  if (quad->at_boundary())
3252  {
3253  // and make sure that we don't attempt to reset the boundary
3254  // indicator to a different than the previously set value
3255  AssertThrow(quad->boundary_id() == 0 ||
3256  quad->boundary_id() == subcell_quad.boundary_id,
3257  ExcInconsistentQuadInfoOfQuad(line[0]->index(),
3258  line[1]->index(),
3259  line[2]->index(),
3260  line[3]->index(),
3261  "boundary ids"));
3262  // If the boundary id provided in subcell_line
3263  // is anything other than the default
3264  // (internal_face_boundary_id), then set it in the new
3265  // triangulation.
3266  if (subcell_quad.boundary_id !=
3268  quad->set_boundary_id(subcell_quad.boundary_id);
3269  }
3270  // Set manifold id if given
3271  if (quad->manifold_id() != numbers::flat_manifold_id)
3272  AssertThrow(quad->manifold_id() == subcell_quad.manifold_id,
3273  ExcInconsistentQuadInfoOfQuad(line[0]->index(),
3274  line[1]->index(),
3275  line[2]->index(),
3276  line[3]->index(),
3277  "manifold ids"));
3278 
3279  quad->set_manifold_id(subcell_quad.manifold_id);
3280  }
3281 
3282 
3284  // finally update neighborship info
3285  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3286  triangulation.begin();
3287  cell != triangulation.end();
3288  ++cell)
3289  for (unsigned int face = 0; face < 6; ++face)
3290  if (adjacent_cells[cell->quad(face)->index()][0] == cell)
3291  // first adjacent cell is
3292  // this one
3293  {
3294  if (adjacent_cells[cell->quad(face)->index()].size() == 2)
3295  // there is another
3296  // adjacent cell
3297  cell->set_neighbor(
3298  face, adjacent_cells[cell->quad(face)->index()][1]);
3299  }
3300  // first adjacent cell is not this
3301  // one, -> it must be the neighbor
3302  // we are looking for
3303  else
3304  cell->set_neighbor(face,
3305  adjacent_cells[cell->quad(face)->index()][0]);
3306  }
3307 
3308 
3324  template <int spacedim>
3325  static void delete_children(
3328  std::vector<unsigned int> &,
3329  std::vector<unsigned int> &)
3330  {
3331  const unsigned int dim = 1;
3332 
3333  // first we need to reset the
3334  // neighbor pointers of the
3335  // neighbors of this cell's
3336  // children to this cell. This is
3337  // different for one dimension,
3338  // since there neighbors can have a
3339  // refinement level differing from
3340  // that of this cell's children by
3341  // more than one level.
3342 
3343  Assert(!cell->child(0)->has_children() &&
3344  !cell->child(1)->has_children(),
3345  ExcInternalError());
3346 
3347  // first do it for the cells to the
3348  // left
3349  if (cell->neighbor(0).state() == IteratorState::valid)
3350  if (cell->neighbor(0)->has_children())
3351  {
3353  cell->neighbor(0);
3354  Assert(neighbor->level() == cell->level(), ExcInternalError());
3355 
3356  // right child
3357  neighbor = neighbor->child(1);
3358  while (true)
3359  {
3360  Assert(neighbor->neighbor(1) == cell->child(0),
3361  ExcInternalError());
3362  neighbor->set_neighbor(1, cell);
3363 
3364  // move on to further
3365  // children on the
3366  // boundary between this
3367  // cell and its neighbor
3368  if (neighbor->has_children())
3369  neighbor = neighbor->child(1);
3370  else
3371  break;
3372  }
3373  }
3374 
3375  // now do it for the cells to the
3376  // left
3377  if (cell->neighbor(1).state() == IteratorState::valid)
3378  if (cell->neighbor(1)->has_children())
3379  {
3381  cell->neighbor(1);
3382  Assert(neighbor->level() == cell->level(), ExcInternalError());
3383 
3384  // left child
3385  neighbor = neighbor->child(0);
3386  while (true)
3387  {
3388  Assert(neighbor->neighbor(0) == cell->child(1),
3389  ExcInternalError());
3390  neighbor->set_neighbor(0, cell);
3391 
3392  // move on to further
3393  // children on the
3394  // boundary between this
3395  // cell and its neighbor
3396  if (neighbor->has_children())
3397  neighbor = neighbor->child(0);
3398  else
3399  break;
3400  }
3401  }
3402 
3403 
3404  // delete the vertex which will not
3405  // be needed anymore. This vertex
3406  // is the second of the first child
3407  triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
3408 
3409  // invalidate children. clear user
3410  // pointers, to avoid that they may
3411  // appear at unwanted places later
3412  // on...
3413  for (unsigned int child = 0; child < cell->n_children(); ++child)
3414  {
3415  cell->child(child)->clear_user_data();
3416  cell->child(child)->clear_user_flag();
3417  cell->child(child)->clear_used_flag();
3418  }
3419 
3420 
3421  // delete pointer to children
3422  cell->clear_children();
3423  cell->clear_user_flag();
3424  }
3425 
3426 
3427 
3428  template <int spacedim>
3429  static void delete_children(
3432  std::vector<unsigned int> & line_cell_count,
3433  std::vector<unsigned int> &)
3434  {
3435  const unsigned int dim = 2;
3436  const RefinementCase<dim> ref_case = cell->refinement_case();
3437 
3438  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
3439  ExcInternalError());
3440 
3441  // vectors to hold all lines which
3442  // may be deleted
3443  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3444  lines_to_delete(0);
3445 
3446  lines_to_delete.reserve(4 * 2 + 4);
3447 
3448  // now we decrease the counters for
3449  // lines contained in the child
3450  // cells
3451  for (unsigned int c = 0; c < cell->n_children(); ++c)
3452  {
3454  cell->child(c);
3455  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3456  --line_cell_count[child->line_index(l)];
3457  }
3458 
3459 
3460  // delete the vertex which will not
3461  // be needed anymore. This vertex
3462  // is the second of the second line
3463  // of the first child, if the cell
3464  // is refined with cut_xy, else there
3465  // is no inner vertex.
3466  // additionally delete unneeded inner
3467  // lines
3468  if (ref_case == RefinementCase<dim>::cut_xy)
3469  {
3470  triangulation
3471  .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
3472 
3473  lines_to_delete.push_back(cell->child(0)->line(1));
3474  lines_to_delete.push_back(cell->child(0)->line(3));
3475  lines_to_delete.push_back(cell->child(3)->line(0));
3476  lines_to_delete.push_back(cell->child(3)->line(2));
3477  }
3478  else
3479  {
3480  unsigned int inner_face_no =
3481  ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
3482 
3483  // the inner line will not be
3484  // used any more
3485  lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
3486  }
3487 
3488  // invalidate children
3489  for (unsigned int child = 0; child < cell->n_children(); ++child)
3490  {
3491  cell->child(child)->clear_user_data();
3492  cell->child(child)->clear_user_flag();
3493  cell->child(child)->clear_used_flag();
3494  }
3495 
3496 
3497  // delete pointer to children
3498  cell->clear_children();
3499  cell->clear_refinement_case();
3500  cell->clear_user_flag();
3501 
3502  // look at the refinement of outer
3503  // lines. if nobody needs those
3504  // anymore we can add them to the
3505  // list of lines to be deleted.
3506  for (unsigned int line_no = 0;
3507  line_no < GeometryInfo<dim>::lines_per_cell;
3508  ++line_no)
3509  {
3511  cell->line(line_no);
3512 
3513  if (line->has_children())
3514  {
3515  // if one of the cell counters is
3516  // zero, the other has to be as well
3517 
3518  Assert((line_cell_count[line->child_index(0)] == 0 &&
3519  line_cell_count[line->child_index(1)] == 0) ||
3520  (line_cell_count[line->child_index(0)] > 0 &&
3521  line_cell_count[line->child_index(1)] > 0),
3522  ExcInternalError());
3523 
3524  if (line_cell_count[line->child_index(0)] == 0)
3525  {
3526  for (unsigned int c = 0; c < 2; ++c)
3527  Assert(!line->child(c)->has_children(),
3528  ExcInternalError());
3529 
3530  // we may delete the line's
3531  // children and the middle vertex
3532  // as no cell references them
3533  // anymore
3534  triangulation
3535  .vertices_used[line->child(0)->vertex_index(1)] = false;
3536 
3537  lines_to_delete.push_back(line->child(0));
3538  lines_to_delete.push_back(line->child(1));
3539 
3540  line->clear_children();
3541  }
3542  }
3543  }
3544 
3545  // finally, delete unneeded lines
3546 
3547  // clear user pointers, to avoid that
3548  // they may appear at unwanted places
3549  // later on...
3550  // same for user flags, then finally
3551  // delete the lines
3552  typename std::vector<
3554  line = lines_to_delete.begin(),
3555  endline = lines_to_delete.end();
3556  for (; line != endline; ++line)
3557  {
3558  (*line)->clear_user_data();
3559  (*line)->clear_user_flag();
3560  (*line)->clear_used_flag();
3561  }
3562  }
3563 
3564 
3565 
3566  template <int spacedim>
3567  static void delete_children(
3570  std::vector<unsigned int> & line_cell_count,
3571  std::vector<unsigned int> & quad_cell_count)
3572  {
3573  const unsigned int dim = 3;
3574 
3575  Assert(line_cell_count.size() == triangulation.n_raw_lines(),
3576  ExcInternalError());
3577  Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
3578  ExcInternalError());
3579 
3580  // first of all, we store the RefineCase of
3581  // this cell
3582  const RefinementCase<dim> ref_case = cell->refinement_case();
3583  // vectors to hold all lines and quads which
3584  // may be deleted
3585  std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3586  lines_to_delete(0);
3587  std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3588  quads_to_delete(0);
3589 
3590  lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3591  quads_to_delete.reserve(6 * 4 + 12);
3592 
3593  // now we decrease the counters for lines and
3594  // quads contained in the child cells
3595  for (unsigned int c = 0; c < cell->n_children(); ++c)
3596  {
3598  cell->child(c);
3599  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3600  --line_cell_count[child->line_index(l)];
3601  for (auto f : GeometryInfo<dim>::face_indices())
3602  --quad_cell_count[child->quad_index(f)];
3603  }
3604 
3606  // delete interior quads and lines and the
3607  // interior vertex, depending on the
3608  // refinement case of the cell
3609  //
3610  // for append quads and lines: only append
3611  // them to the list of objects to be deleted
3612 
3613  switch (ref_case)
3614  {
3616  quads_to_delete.push_back(cell->child(0)->face(1));
3617  break;
3619  quads_to_delete.push_back(cell->child(0)->face(3));
3620  break;
3622  quads_to_delete.push_back(cell->child(0)->face(5));
3623  break;
3625  quads_to_delete.push_back(cell->child(0)->face(1));
3626  quads_to_delete.push_back(cell->child(0)->face(3));
3627  quads_to_delete.push_back(cell->child(3)->face(0));
3628  quads_to_delete.push_back(cell->child(3)->face(2));
3629 
3630  lines_to_delete.push_back(cell->child(0)->line(11));
3631  break;
3633  quads_to_delete.push_back(cell->child(0)->face(1));
3634  quads_to_delete.push_back(cell->child(0)->face(5));
3635  quads_to_delete.push_back(cell->child(3)->face(0));
3636  quads_to_delete.push_back(cell->child(3)->face(4));
3637 
3638  lines_to_delete.push_back(cell->child(0)->line(5));
3639  break;
3641  quads_to_delete.push_back(cell->child(0)->face(3));
3642  quads_to_delete.push_back(cell->child(0)->face(5));
3643  quads_to_delete.push_back(cell->child(3)->face(2));
3644  quads_to_delete.push_back(cell->child(3)->face(4));
3645 
3646  lines_to_delete.push_back(cell->child(0)->line(7));
3647  break;
3649  quads_to_delete.push_back(cell->child(0)->face(1));
3650  quads_to_delete.push_back(cell->child(2)->face(1));
3651  quads_to_delete.push_back(cell->child(4)->face(1));
3652  quads_to_delete.push_back(cell->child(6)->face(1));
3653 
3654  quads_to_delete.push_back(cell->child(0)->face(3));
3655  quads_to_delete.push_back(cell->child(1)->face(3));
3656  quads_to_delete.push_back(cell->child(4)->face(3));
3657  quads_to_delete.push_back(cell->child(5)->face(3));
3658 
3659  quads_to_delete.push_back(cell->child(0)->face(5));
3660  quads_to_delete.push_back(cell->child(1)->face(5));
3661  quads_to_delete.push_back(cell->child(2)->face(5));
3662  quads_to_delete.push_back(cell->child(3)->face(5));
3663 
3664  lines_to_delete.push_back(cell->child(0)->line(5));
3665  lines_to_delete.push_back(cell->child(0)->line(7));
3666  lines_to_delete.push_back(cell->child(0)->line(11));
3667  lines_to_delete.push_back(cell->child(7)->line(0));
3668  lines_to_delete.push_back(cell->child(7)->line(2));
3669  lines_to_delete.push_back(cell->child(7)->line(8));
3670  // delete the vertex which will not
3671  // be needed anymore. This vertex
3672  // is the vertex at the heart of
3673  // this cell, which is the sixth of
3674  // the first child
3675  triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3676  false;
3677  break;
3678  default:
3679  // only remaining case is
3680  // no_refinement, thus an error
3681  Assert(false, ExcInternalError());
3682  break;
3683  }
3684 
3685 
3686  // invalidate children
3687  for (unsigned int child = 0; child < cell->n_children(); ++child)
3688  {
3689  cell->child(child)->clear_user_data();
3690  cell->child(child)->clear_user_flag();
3691 
3692  for (auto f : GeometryInfo<dim>::face_indices())
3693  {
3694  // set flags denoting deviations from
3695  // standard orientation of faces back
3696  // to initialization values
3697  cell->child(child)->set_face_orientation(f, true);
3698  cell->child(child)->set_face_flip(f, false);
3699  cell->child(child)->set_face_rotation(f, false);
3700  }
3701 
3702  cell->child(child)->clear_used_flag();
3703  }
3704 
3705 
3706  // delete pointer to children
3707  cell->clear_children();
3708  cell->clear_refinement_case();
3709  cell->clear_user_flag();
3710 
3711  // so far we only looked at inner quads,
3712  // lines and vertices. Now we have to
3713  // consider outer ones as well. here, we have
3714  // to check, whether there are other cells
3715  // still needing these objects. otherwise we
3716  // can delete them. first for quads (and
3717  // their inner lines).
3718 
3719  for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
3720  {
3722  cell->face(quad_no);
3723 
3724  Assert(
3725  (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
3726  quad->has_children()) ||
3727  GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
3728  RefinementCase<dim - 1>::no_refinement,
3729  ExcInternalError());
3730 
3731  switch (quad->refinement_case())
3732  {
3733  case RefinementCase<dim - 1>::no_refinement:
3734  // nothing to do as the quad
3735  // is not refined
3736  break;
3737  case RefinementCase<dim - 1>::cut_x:
3738  case RefinementCase<dim - 1>::cut_y:
3739  {
3740  // if one of the cell counters is
3741  // zero, the other has to be as
3742  // well
3743  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3744  quad_cell_count[quad->child_index(1)] == 0) ||
3745  (quad_cell_count[quad->child_index(0)] > 0 &&
3746  quad_cell_count[quad->child_index(1)] > 0),
3747  ExcInternalError());
3748  // it might be, that the quad is
3749  // refined twice anisotropically,
3750  // first check, whether we may
3751  // delete possible grand_children
3752  unsigned int deleted_grandchildren = 0;
3753  unsigned int number_of_child_refinements = 0;
3754 
3755  for (unsigned int c = 0; c < 2; ++c)
3756  if (quad->child(c)->has_children())
3757  {
3758  ++number_of_child_refinements;
3759  // if one of the cell counters is
3760  // zero, the other has to be as
3761  // well
3762  Assert(
3763  (quad_cell_count[quad->child(c)->child_index(0)] ==
3764  0 &&
3765  quad_cell_count[quad->child(c)->child_index(1)] ==
3766  0) ||
3767  (quad_cell_count[quad->child(c)->child_index(0)] >
3768  0 &&
3769  quad_cell_count[quad->child(c)->child_index(1)] >
3770  0),
3771  ExcInternalError());
3772  if (quad_cell_count[quad->child(c)->child_index(0)] ==
3773  0)
3774  {
3775  // Assert, that the two
3776  // anisotropic
3777  // refinements add up to
3778  // isotropic refinement
3779  Assert(quad->refinement_case() +
3780  quad->child(c)->refinement_case() ==
3782  ExcInternalError());
3783  // we may delete the
3784  // quad's children and
3785  // the inner line as no
3786  // cell references them
3787  // anymore
3788  quads_to_delete.push_back(
3789  quad->child(c)->child(0));
3790  quads_to_delete.push_back(
3791  quad->child(c)->child(1));
3792  if (quad->child(c)->refinement_case() ==
3794  lines_to_delete.push_back(
3795  quad->child(c)->child(0)->line(1));
3796  else
3797  lines_to_delete.push_back(
3798  quad->child(c)->child(0)->line(3));
3799  quad->child(c)->clear_children();
3800  quad->child(c)->clear_refinement_case();
3801  ++deleted_grandchildren;
3802  }
3803  }
3804  // if no grandchildren are left, we
3805  // may as well delete the
3806  // refinement of the inner line
3807  // between our children and the
3808  // corresponding vertex
3809  if (number_of_child_refinements > 0 &&
3810  deleted_grandchildren == number_of_child_refinements)
3811  {
3813  middle_line;
3814  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3815  middle_line = quad->child(0)->line(1);
3816  else
3817  middle_line = quad->child(0)->line(3);
3818 
3819  lines_to_delete.push_back(middle_line->child(0));
3820  lines_to_delete.push_back(middle_line->child(1));
3821  triangulation
3822  .vertices_used[middle_vertex_index<dim, spacedim>(
3823  middle_line)] = false;
3824  middle_line->clear_children();
3825  }
3826 
3827  // now consider the direct children
3828  // of the given quad
3829  if (quad_cell_count[quad->child_index(0)] == 0)
3830  {
3831  // we may delete the quad's
3832  // children and the inner line
3833  // as no cell references them
3834  // anymore
3835  quads_to_delete.push_back(quad->child(0));
3836  quads_to_delete.push_back(quad->child(1));
3837  if (quad->refinement_case() == RefinementCase<2>::cut_x)
3838  lines_to_delete.push_back(quad->child(0)->line(1));
3839  else
3840  lines_to_delete.push_back(quad->child(0)->line(3));
3841 
3842  // if the counters just dropped
3843  // to zero, otherwise the
3844  // children would have been
3845  // deleted earlier, then this
3846  // cell's children must have
3847  // contained the anisotropic
3848  // quad children. thus, if
3849  // those have again anisotropic
3850  // children, which are in
3851  // effect isotropic children of
3852  // the original quad, those are
3853  // still needed by a
3854  // neighboring cell and we
3855  // cannot delete them. instead,
3856  // we have to reset this quad's
3857  // refine case to isotropic and
3858  // set the children
3859  // accordingly.
3860  if (quad->child(0)->has_children())
3861  if (quad->refinement_case() ==
3863  {
3864  // now evereything is
3865  // quite complicated. we
3866  // have the children
3867  // numbered according to
3868  //
3869  // *---*---*
3870  // |n+1|m+1|
3871  // *---*---*
3872  // | n | m |
3873  // *---*---*
3874  //
3875  // from the original
3876  // anisotropic
3877  // refinement. we have to
3878  // reorder them as
3879  //
3880  // *---*---*
3881  // | m |m+1|
3882  // *---*---*
3883  // | n |n+1|
3884  // *---*---*
3885  //
3886  // for isotropic refinement.
3887  //
3888  // this is a bit ugly, of
3889  // course: loop over all
3890  // cells on all levels
3891  // and look for faces n+1
3892  // (switch_1) and m
3893  // (switch_2).
3894  const typename Triangulation<dim, spacedim>::
3895  quad_iterator switch_1 =
3896  quad->child(0)->child(1),
3897  switch_2 =
3898  quad->child(1)->child(0);
3899 
3900  Assert(!switch_1->has_children(),
3901  ExcInternalError());
3902  Assert(!switch_2->has_children(),
3903  ExcInternalError());
3904 
3905  const int switch_1_index = switch_1->index();
3906  const int switch_2_index = switch_2->index();
3907  for (unsigned int l = 0;
3908  l < triangulation.levels.size();
3909  ++l)
3910  for (unsigned int h = 0;
3911  h <
3912  triangulation.levels[l]->cells.n_objects();
3913  ++h)
3914  for (const unsigned int q :
3916  {
3917  const int index = triangulation.levels[l]
3918  ->cells.get_object(h)
3919  .face(q);
3920  if (index == switch_1_index)
3921  triangulation.levels[l]
3922  ->cells.get_object(h)
3923  .set_face(q, switch_2_index);
3924  else if (index == switch_2_index)
3925  triangulation.levels[l]
3926  ->cells.get_object(h)
3927  .set_face(q, switch_1_index);
3928  }
3929  // now we have to copy
3930  // all information of the
3931  // two quads
3932  const int switch_1_lines[4] = {
3933  static_cast<signed int>(
3934  switch_1->line_index(0)),
3935  static_cast<signed int>(
3936  switch_1->line_index(1)),
3937  static_cast<signed int>(
3938  switch_1->line_index(2)),
3939  static_cast<signed int>(
3940  switch_1->line_index(3))};
3941  const bool switch_1_line_orientations[4] = {
3942  switch_1->line_orientation(0),
3943  switch_1->line_orientation(1),
3944  switch_1->line_orientation(2),
3945  switch_1->line_orientation(3)};
3946  const types::boundary_id switch_1_boundary_id =
3947  switch_1->boundary_id();
3948  const unsigned int switch_1_user_index =
3949  switch_1->user_index();
3950  const bool switch_1_user_flag =
3951  switch_1->user_flag_set();
3952 
3953  switch_1->set(
3955  TriaObject<2>(switch_2->line_index(0),
3956  switch_2->line_index(1),
3957  switch_2->line_index(2),
3958  switch_2->line_index(3)));
3959  switch_1->set_line_orientation(
3960  0, switch_2->line_orientation(0));
3961  switch_1->set_line_orientation(
3962  1, switch_2->line_orientation(1));
3963  switch_1->set_line_orientation(
3964  2, switch_2->line_orientation(2));
3965  switch_1->set_line_orientation(
3966  3, switch_2->line_orientation(3));
3967  switch_1->set_boundary_id_internal(
3968  switch_2->boundary_id());
3969  switch_1->set_manifold_id(
3970  switch_2->manifold_id());
3971  switch_1->set_user_index(switch_2->user_index());
3972  if (switch_2->user_flag_set())
3973  switch_1->set_user_flag();
3974  else
3975  switch_1->clear_user_flag();
3976 
3977  switch_2->set(
3979  TriaObject<2>(switch_1_lines[0],
3980  switch_1_lines[1],
3981  switch_1_lines[2],
3982  switch_1_lines[3]));
3983  switch_2->set_line_orientation(
3984  0, switch_1_line_orientations[0]);
3985  switch_2->set_line_orientation(
3986  1, switch_1_line_orientations[1]);
3987  switch_2->set_line_orientation(
3988  2, switch_1_line_orientations[2]);
3989  switch_2->set_line_orientation(
3990  3, switch_1_line_orientations[3]);
3991  switch_2->set_boundary_id_internal(
3992  switch_1_boundary_id);
3993  switch_2->set_manifold_id(
3994  switch_1->manifold_id());
3995  switch_2->set_user_index(switch_1_user_index);
3996  if (switch_1_user_flag)
3997  switch_2->set_user_flag();
3998  else
3999  switch_2->clear_user_flag();
4000 
4001  const unsigned int child_0 =
4002  quad->child(0)->child_index(0);
4003  const unsigned int child_2 =
4004  quad->child(1)->child_index(0);
4005  quad->clear_children();
4006  quad->clear_refinement_case();
4007  quad->set_refinement_case(
4009  quad->set_children(0, child_0);
4010  quad->set_children(2, child_2);
4011  std::swap(quad_cell_count[child_0 + 1],
4012  quad_cell_count[child_2]);
4013  }
4014  else
4015  {
4016  // the face was refined
4017  // with cut_y, thus the
4018  // children are already
4019  // in correct order. we
4020  // only have to set them
4021  // correctly, deleting
4022  // the indirection of two
4023  // anisotropic refinement
4024  // and going directly
4025  // from the quad to
4026  // isotropic children
4027  const unsigned int child_0 =
4028  quad->child(0)->child_index(0);
4029  const unsigned int child_2 =
4030  quad->child(1)->child_index(0);
4031  quad->clear_children();
4032  quad->clear_refinement_case();
4033  quad->set_refinement_case(
4035  quad->set_children(0, child_0);
4036  quad->set_children(2, child_2);
4037  }
4038  else
4039  {
4040  quad->clear_children();
4041  quad->clear_refinement_case();
4042  }
4043  }
4044  break;
4045  }
4046  case RefinementCase<dim - 1>::cut_xy:
4047  {
4048  // if one of the cell counters is
4049  // zero, the others have to be as
4050  // well
4051 
4052  Assert((quad_cell_count[quad->child_index(0)] == 0 &&
4053  quad_cell_count[quad->child_index(1)] == 0 &&
4054  quad_cell_count[quad->child_index(2)] == 0 &&
4055  quad_cell_count[quad->child_index(3)] == 0) ||
4056  (quad_cell_count[quad->child_index(0)] > 0 &&
4057  quad_cell_count[quad->child_index(1)] > 0 &&
4058  quad_cell_count[quad->child_index(2)] > 0 &&
4059  quad_cell_count[quad->child_index(3)] > 0),
4060  ExcInternalError());
4061 
4062  if (quad_cell_count[quad->child_index(0)] == 0)
4063  {
4064  // we may delete the quad's
4065  // children, the inner lines
4066  // and the middle vertex as no
4067  // cell references them anymore
4068  lines_to_delete.push_back(quad->child(0)->line(1));
4069  lines_to_delete.push_back(quad->child(3)->line(0));
4070  lines_to_delete.push_back(quad->child(0)->line(3));
4071  lines_to_delete.push_back(quad->child(3)->line(2));
4072 
4073  for (unsigned int child = 0; child < quad->n_children();
4074  ++child)
4075  quads_to_delete.push_back(quad->child(child));
4076 
4077  triangulation
4078  .vertices_used[quad->child(0)->vertex_index(3)] =
4079  false;
4080 
4081  quad->clear_children();
4082  quad->clear_refinement_case();
4083  }
4084  }
4085  break;
4086 
4087  default:
4088  Assert(false, ExcInternalError());
4089  break;
4090  }
4091  }
4092 
4093  // now we repeat a similar procedure
4094  // for the outer lines of this cell.
4095 
4096  // if in debug mode: check that each
4097  // of the lines for which we consider
4098  // deleting the children in fact has
4099  // children (the bits/coarsening_3d
4100  // test tripped over this initially)
4101  for (unsigned int line_no = 0;
4102  line_no < GeometryInfo<dim>::lines_per_cell;
4103  ++line_no)
4104  {
4106  cell->line(line_no);
4107 
4108  Assert(
4109  (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
4110  line->has_children()) ||
4111  GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
4113  ExcInternalError());
4114 
4115  if (line->has_children())
4116  {
4117  // if one of the cell counters is
4118  // zero, the other has to be as well
4119 
4120  Assert((line_cell_count[line->child_index(0)] == 0 &&
4121  line_cell_count[line->child_index(1)] == 0) ||
4122  (line_cell_count[line->child_index(0)] > 0 &&
4123  line_cell_count[line->child_index(1)] > 0),
4124  ExcInternalError());
4125 
4126  if (line_cell_count[line->child_index(0)] == 0)
4127  {
4128  for (unsigned int c = 0; c < 2; ++c)
4129  Assert(!line->child(c)->has_children(),
4130  ExcInternalError());
4131 
4132  // we may delete the line's
4133  // children and the middle vertex
4134  // as no cell references them
4135  // anymore
4136  triangulation
4137  .vertices_used[line->child(0)->vertex_index(1)] = false;
4138 
4139  lines_to_delete.push_back(line->child(0));
4140  lines_to_delete.push_back(line->child(1));
4141 
4142  line->clear_children();
4143  }
4144  }
4145  }
4146 
4147  // finally, delete unneeded quads and lines
4148 
4149  // clear user pointers, to avoid that
4150  // they may appear at unwanted places
4151  // later on...
4152  // same for user flags, then finally
4153  // delete the quads and lines
4154  typename std::vector<
4156  line = lines_to_delete.begin(),
4157  endline = lines_to_delete.end();
4158  for (; line != endline; ++line)
4159  {
4160  (*line)->clear_user_data();
4161  (*line)->clear_user_flag();
4162  (*line)->clear_used_flag();
4163  }
4164 
4165  typename std::vector<
4167  quad = quads_to_delete.begin(),
4168  endquad = quads_to_delete.end();
4169  for (; quad != endquad; ++quad)
4170  {
4171  (*quad)->clear_user_data();
4172  (*quad)->clear_children();
4173  (*quad)->clear_refinement_case();
4174  (*quad)->clear_user_flag();
4175  (*quad)->clear_used_flag();
4176  }
4177  }
4178 
4179 
4197  template <int spacedim>
4198  static void create_children(
4200  unsigned int & next_unused_vertex,
4202  &next_unused_line,
4204  & next_unused_cell,
4206  {
4207  const unsigned int dim = 2;
4208  // clear refinement flag
4209  const RefinementCase<dim> ref_case = cell->refine_flag_set();
4210  cell->clear_refine_flag();
4211 
4212  /* For the refinement process: since we go the levels up from the
4213  lowest, there are (unlike above) only two possibilities: a neighbor
4214  cell is on the same level or one level up (in both cases, it may or
4215  may not be refined later on, but we don't care here).
4216 
4217  First:
4218  Set up an array of the 3x3 vertices, which are distributed on the
4219  cell (the array consists of indices into the @p{vertices} std::vector
4220 
4221  2--7--3
4222  | | |
4223  4--8--5
4224  | | |
4225  0--6--1
4226 
4227  note: in case of cut_x or cut_y not all these vertices are needed for
4228  the new cells
4229 
4230  Second:
4231  Set up an array of the new lines (the array consists of iterator
4232  pointers into the lines arrays)
4233 
4234  .-6-.-7-. The directions are: .->-.->-.
4235  1 9 3 ^ ^ ^
4236  .-10.11-. .->-.->-.
4237  0 8 2 ^ ^ ^
4238  .-4-.-5-. .->-.->-.
4239 
4240  cut_x:
4241  .-4-.-5-.
4242  | | |
4243  0 6 1
4244  | | |
4245  .-2-.-3-.
4246 
4247  cut_y:
4248  .---5---.
4249  1 3
4250  .---6---.
4251  0 2
4252  .---4---.
4253 
4254 
4255  Third:
4256  Set up an array of neighbors:
4257 
4258  6 7
4259  .--.--.
4260  1| | |3
4261  .--.--.
4262  0| | |2
4263  .--.--.
4264  4 5
4265 
4266  We need this array for two reasons: first to get the lines which will
4267  bound the four subcells (if the neighboring cell is refined, these
4268  lines already exist), and second to update neighborship information.
4269  Since if a neighbor is not refined, its neighborship record only
4270  points to the present, unrefined, cell rather than the children we
4271  are presently creating, we only need the neighborship information
4272  if the neighbor cells are refined. In all other cases, we store
4273  the unrefined neighbor address
4274 
4275  We also need for every neighbor (if refined) which number among its
4276  neighbors the present (unrefined) cell has, since that number is to
4277  be replaced and because that also is the number of the subline which
4278  will be the interface between that neighbor and the to be created
4279  cell. We will store this number (between 0 and 3) in the field
4280  @p{neighbors_neighbor}.
4281 
4282  It would be sufficient to use the children of the common line to the
4283  neighbor, if we only wanted to get the new sublines and the new
4284  vertex, but because we need to update the neighborship information of
4285  the two refined subcells of the neighbor, we need to search these
4286  anyway.
4287 
4288  Convention:
4289  The created children are numbered like this:
4290 
4291  .--.--.
4292  |2 . 3|
4293  .--.--.
4294  |0 | 1|
4295  .--.--.
4296  */
4297  // collect the
4298  // indices of the
4299  // eight
4300  // surrounding
4301  // vertices
4302  // 2--7--3
4303  // | | |
4304  // 4--9--5
4305  // | | |
4306  // 0--6--1
4307  int new_vertices[9];
4308  for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
4309  new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4310  for (unsigned int line_no = 0; line_no < 4; ++line_no)
4311  if (cell->line(line_no)->has_children())
4312  new_vertices[4 + line_no] =
4313  cell->line(line_no)->child(0)->vertex_index(1);
4314 
4315  if (ref_case == RefinementCase<dim>::cut_xy)
4316  {
4317  // find the next
4318  // unused vertex and
4319  // allocate it for
4320  // the new vertex we
4321  // need here
4322  while (triangulation.vertices_used[next_unused_vertex] == true)
4323  ++next_unused_vertex;
4324  Assert(
4325  next_unused_vertex < triangulation.vertices.size(),
4326  ExcMessage(
4327  "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."));
4328  triangulation.vertices_used[next_unused_vertex] = true;
4329 
4330  new_vertices[8] = next_unused_vertex;
4331 
4332  // if this quad lives
4333  // in 2d, then we can
4334  // compute the new
4335  // central vertex
4336  // location just from
4337  // the surrounding
4338  // ones. If this is
4339  // not the case, then
4340  // we need to ask a
4341  // boundary object
4342  if (dim == spacedim)
4343  {
4344  // triangulation.vertices[next_unused_vertex] = new_point;
4345  triangulation.vertices[next_unused_vertex] = cell->center(true);
4346 
4347  // if the user_flag is set, i.e. if the cell is at the
4348  // boundary, use a different calculation of the middle vertex
4349  // here. this is of advantage if the boundary is strongly
4350  // curved (whereas the cell is not) and the cell has a high
4351  // aspect ratio.
4352  if (cell->user_flag_set())
4353  {
4354  // first reset the user_flag and then refine
4355  cell->clear_user_flag();
4356  triangulation.vertices[next_unused_vertex] =
4357  cell->center(true, true);
4358  }
4359  }
4360  else
4361  {
4362  // if this quad lives in a higher dimensional space
4363  // then we don't need to worry if it is at the
4364  // boundary of the manifold -- we always have to use
4365  // the boundary object anyway; so ignore whether the
4366  // user flag is set or not
4367  cell->clear_user_flag();
4368 
4369  // determine middle vertex by transfinite interpolation to be
4370  // consistent with what happens to quads in a Triangulation<3,
4371  // 3> when they are refined
4372  triangulation.vertices[next_unused_vertex] =
4373  cell->center(true, true);
4374  }
4375  }
4376 
4377 
4378  // Now the lines:
4379  typename Triangulation<dim, spacedim>::raw_line_iterator new_lines[12];
4380  unsigned int lmin = 8;
4381  unsigned int lmax = 12;
4382  if (ref_case != RefinementCase<dim>::cut_xy)
4383  {
4384  lmin = 6;
4385  lmax = 7;
4386  }
4387 
4388  for (unsigned int l = lmin; l < lmax; ++l)
4389  {
4390  while (next_unused_line->used() == true)
4391  ++next_unused_line;
4392  new_lines[l] = next_unused_line;
4393  ++next_unused_line;
4394 
4395  Assert(
4396  new_lines[l]->used() == false,
4397  ExcMessage(
4398  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4399  }
4400 
4401  if (ref_case == RefinementCase<dim>::cut_xy)
4402  {
4403  // .-6-.-7-.
4404  // 1 9 3
4405  // .-10.11-.
4406  // 0 8 2
4407  // .-4-.-5-.
4408 
4409  // lines 0-7 already exist, create only the four interior
4410  // lines 8-11
4411  unsigned int l = 0;
4412  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
4413  for (unsigned int c = 0; c < 2; ++c, ++l)
4414  new_lines[l] = cell->line(face_no)->child(c);
4415  Assert(l == 8, ExcInternalError());
4416 
4417  new_lines[8]->set(
4419  new_vertices[6], new_vertices[8]));
4420  new_lines[9]->set(
4422  new_vertices[8], new_vertices[7]));
4423  new_lines[10]->set(
4425  new_vertices[4], new_vertices[8]));
4426  new_lines[11]->set(
4428  new_vertices[8], new_vertices[5]));
4429  }
4430  else if (ref_case == RefinementCase<dim>::cut_x)
4431  {
4432  // .-4-.-5-.
4433  // | | |
4434  // 0 6 1
4435  // | | |
4436  // .-2-.-3-.
4437  new_lines[0] = cell->line(0);
4438  new_lines[1] = cell->line(1);
4439  new_lines[2] = cell->line(2)->child(0);
4440  new_lines[3] = cell->line(2)->child(1);
4441  new_lines[4] = cell->line(3)->child(0);
4442  new_lines[5] = cell->line(3)->child(1);
4443  new_lines[6]->set(
4445  new_vertices[6], new_vertices[7]));
4446  }
4447  else
4448  {
4450  // .---5---.
4451  // 1 3
4452  // .---6---.
4453  // 0 2
4454  // .---4---.
4455  new_lines[0] = cell->line(0)->child(0);
4456  new_lines[1] = cell->line(0)->child(1);
4457  new_lines[2] = cell->line(1)->child(0);
4458  new_lines[3] = cell->line(1)->child(1);
4459  new_lines[4] = cell->line(2);
4460  new_lines[5] = cell->line(3);
4461  new_lines[6]->set(
4463  new_vertices[4], new_vertices[5]));
4464  }
4465 
4466  for (unsigned int l = lmin; l < lmax; ++l)
4467  {
4468  new_lines[l]->set_used_flag();
4469  new_lines[l]->clear_user_flag();
4470  new_lines[l]->clear_user_data();
4471  new_lines[l]->clear_children();
4472  // interior line
4473  new_lines[l]->set_boundary_id_internal(
4475  new_lines[l]->set_manifold_id(cell->manifold_id());
4476  }
4477 
4478  // Now add the four (two)
4479  // new cells!
4482  while (next_unused_cell->used() == true)
4483  ++next_unused_cell;
4484 
4485  const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
4486  for (unsigned int i = 0; i < n_children; ++i)
4487  {
4488  Assert(
4489  next_unused_cell->used() == false,
4490  ExcMessage(
4491  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4492  subcells[i] = next_unused_cell;
4493  ++next_unused_cell;
4494  if (i % 2 == 1 && i < n_children - 1)
4495  while (next_unused_cell->used() == true)
4496  ++next_unused_cell;
4497  }
4498 
4499  if (ref_case == RefinementCase<dim>::cut_xy)
4500  {
4501  // children:
4502  // .--.--.
4503  // |2 . 3|
4504  // .--.--.
4505  // |0 | 1|
4506  // .--.--.
4507  // lines:
4508  // .-6-.-7-.
4509  // 1 9 3
4510  // .-10.11-.
4511  // 0 8 2
4512  // .-4-.-5-.
4513  subcells[0]->set(
4515  new_lines[0]->index(),
4516  new_lines[8]->index(),
4517  new_lines[4]->index(),
4518  new_lines[10]->index()));
4519  subcells[1]->set(
4521  new_lines[8]->index(),
4522  new_lines[2]->index(),
4523  new_lines[5]->index(),
4524  new_lines[11]->index()));
4525  subcells[2]->set(
4527  new_lines[1]->index(),
4528  new_lines[9]->index(),
4529  new_lines[10]->index(),
4530  new_lines[6]->index()));
4531  subcells[3]->set(
4533  new_lines[9]->index(),
4534  new_lines[3]->index(),
4535  new_lines[11]->index(),
4536  new_lines[7]->index()));
4537  }
4538  else if (ref_case == RefinementCase<dim>::cut_x)
4539  {
4540  // children:
4541  // .--.--.
4542  // | . |
4543  // .0 . 1.
4544  // | | |
4545  // .--.--.
4546  // lines:
4547  // .-4-.-5-.
4548  // | | |
4549  // 0 6 1
4550  // | | |
4551  // .-2-.-3-.
4552  subcells[0]->set(
4554  new_lines[0]->index(),
4555  new_lines[6]->index(),
4556  new_lines[2]->index(),
4557  new_lines[4]->index()));
4558  subcells[1]->set(
4560  new_lines[6]->index(),
4561  new_lines[1]->index(),
4562  new_lines[3]->index(),
4563  new_lines[5]->index()));
4564  }
4565  else
4566  {
4568  // children:
4569  // .-----.
4570  // | 1 |
4571  // .-----.
4572  // | 0 |
4573  // .-----.
4574  // lines:
4575  // .---5---.
4576  // 1 3
4577  // .---6---.
4578  // 0 2
4579  // .---4---.
4580  subcells[0]->set(
4582  new_lines[0]->index(),
4583  new_lines[2]->index(),
4584  new_lines[4]->index(),
4585  new_lines[6]->index()));
4586  subcells[1]->set(
4588  new_lines[1]->index(),
4589  new_lines[3]->index(),
4590  new_lines[6]->index(),
4591  new_lines[5]->index()));
4592  }
4593 
4594  types::subdomain_id subdomainid = cell->subdomain_id();
4595 
4596  for (unsigned int i = 0; i < n_children; ++i)
4597  {
4598  subcells[i]->set_used_flag();
4599  subcells[i]->clear_refine_flag();
4600  subcells[i]->clear_user_flag();
4601  subcells[i]->clear_user_data();
4602  subcells[i]->clear_children();
4603  // inherit material
4604  // properties
4605  subcells[i]->set_material_id(cell->material_id());
4606  subcells[i]->set_manifold_id(cell->manifold_id());
4607  subcells[i]->set_subdomain_id(subdomainid);
4608 
4609  if (i % 2 == 0)
4610  subcells[i]->set_parent(cell->index());
4611  }
4612 
4613 
4614 
4615  // set child index for
4616  // even children children
4617  // i=0,2 (0)
4618  for (unsigned int i = 0; i < n_children / 2; ++i)
4619  cell->set_children(2 * i, subcells[2 * i]->index());
4620  // set the refine case
4621  cell->set_refinement_case(ref_case);
4622 
4623  // note that the
4624  // refinement flag was
4625  // already cleared at the
4626  // beginning of this function
4627 
4628  if (dim < spacedim)
4629  for (unsigned int c = 0; c < n_children; ++c)
4630  cell->child(c)->set_direction_flag(cell->direction_flag());
4631  }
4632 
4633 
4634 
4639  template <int spacedim>
4642  const bool /*check_for_distorted_cells*/)
4643  {
4644  const unsigned int dim = 1;
4645 
4646  // check whether a new level is needed we have to check for
4647  // this on the highest level only (on this, all used cells are
4648  // also active, so we only have to check for this)
4649  {
4651  cell = triangulation.begin_active(triangulation.levels.size() - 1),
4652  endc = triangulation.end();
4653  for (; cell != endc; ++cell)
4654  if (cell->used())
4655  if (cell->refine_flag_set())
4656  {
4657  triangulation.levels.push_back(
4658  std::make_unique<
4660  break;
4661  }
4662  }
4663 
4664 
4665  // check how much space is needed on every level we need not
4666  // check the highest level since either - on the highest level
4667  // no cells are flagged for refinement - there are, but
4668  // prepare_refinement added another empty level
4669  unsigned int needed_vertices = 0;
4670  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4671  {
4672  // count number of flagged
4673  // cells on this level
4674  unsigned int flagged_cells = 0;
4676  acell = triangulation.begin_active(level),
4677  aendc = triangulation.begin_active(level + 1);
4678  for (; acell != aendc; ++acell)
4679  if (acell->refine_flag_set())
4680  ++flagged_cells;
4681 
4682  // count number of used cells
4683  // on the next higher level
4684  const unsigned int used_cells =
4685  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4686  triangulation.levels[level + 1]->cells.used.end(),
4687  true);
4688 
4689  // reserve space for the used_cells cells already existing
4690  // on the next higher level as well as for the
4691  // 2*flagged_cells that will be created on that level
4692  triangulation.levels[level + 1]->reserve_space(
4693  used_cells +
4695  1,
4696  spacedim);
4697  // reserve space for 2*flagged_cells new lines on the next
4698  // higher level
4699  triangulation.levels[level + 1]->cells.reserve_space(
4700  GeometryInfo<1>::max_children_per_cell * flagged_cells, 0);
4701 
4702  needed_vertices += flagged_cells;
4703  }
4704 
4705  // add to needed vertices how many
4706  // vertices are already in use
4707  needed_vertices += std::count(triangulation.vertices_used.begin(),
4708  triangulation.vertices_used.end(),
4709  true);
4710  // if we need more vertices: create them, if not: leave the
4711  // array as is, since shrinking is not really possible because
4712  // some of the vertices at the end may be in use
4713  if (needed_vertices > triangulation.vertices.size())
4714  {
4715  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4716  triangulation.vertices_used.resize(needed_vertices, false);
4717  }
4718 
4719 
4720  // Do REFINEMENT on every level; exclude highest level as
4721  // above
4722 
4723  // index of next unused vertex
4724  unsigned int next_unused_vertex = 0;
4725 
4726  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4727  {
4729  cell = triangulation.begin_active(level),
4730  endc = triangulation.begin_active(level + 1);
4731 
4733  next_unused_cell = triangulation.begin_raw(level + 1);
4734 
4735  for (; (cell != endc) && (cell->level() == level); ++cell)
4736  if (cell->refine_flag_set())
4737  {
4738  // clear refinement flag
4739  cell->clear_refine_flag();
4740 
4741  // search for next unused
4742  // vertex
4743  while (triangulation.vertices_used[next_unused_vertex] ==
4744  true)
4745  ++next_unused_vertex;
4746  Assert(
4747  next_unused_vertex < triangulation.vertices.size(),
4748  ExcMessage(
4749  "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."));
4750 
4751  // Now we always ask the cell itself where to put
4752  // the new point. The cell in turn will query the
4753  // manifold object internally.
4754  triangulation.vertices[next_unused_vertex] =
4755  cell->center(true);
4756 
4757  triangulation.vertices_used[next_unused_vertex] = true;
4758 
4759  // search for next two unused cell (++ takes care of
4760  // the end of the vector)
4762  first_child,
4763  second_child;
4764  while (next_unused_cell->used() == true)
4765  ++next_unused_cell;
4766  first_child = next_unused_cell;
4767  first_child->set_used_flag();
4768  first_child->clear_user_data();
4769  ++next_unused_cell;
4770  Assert(
4771  next_unused_cell->used() == false,
4772  ExcMessage(
4773  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4774  second_child = next_unused_cell;
4775  second_child->set_used_flag();
4776  second_child->clear_user_data();
4777 
4778  types::subdomain_id subdomainid = cell->subdomain_id();
4779 
4780  // insert first child
4781  cell->set_children(0, first_child->index());
4782  first_child->clear_children();
4783  first_child->set(
4785  cell->vertex_index(0), next_unused_vertex));
4786  first_child->set_material_id(cell->material_id());
4787  first_child->set_manifold_id(cell->manifold_id());
4788  first_child->set_subdomain_id(subdomainid);
4789  first_child->set_direction_flag(cell->direction_flag());
4790 
4791  first_child->set_parent(cell->index());
4792 
4793  // Set manifold id of the right face. Only do this
4794  // on the first child.
4795  first_child->face(1)->set_manifold_id(cell->manifold_id());
4796 
4797  // reset neighborship info (refer to
4798  // internal::TriangulationImplementation::TriaLevel<0> for
4799  // details)
4800  first_child->set_neighbor(1, second_child);
4801  if (cell->neighbor(0).state() != IteratorState::valid)
4802  first_child->set_neighbor(0, cell->neighbor(0));
4803  else if (cell->neighbor(0)->is_active())
4804  {
4805  // since the neighbors level is always <=level,
4806  // if the cell is active, then there are no
4807  // cells to the left which may want to know
4808  // about this new child cell.
4809  Assert(cell->neighbor(0)->level() <= cell->level(),
4810  ExcInternalError());
4811  first_child->set_neighbor(0, cell->neighbor(0));
4812  }
4813  else
4814  // left neighbor is refined
4815  {
4816  // set neighbor to cell on same level
4817  const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4818  first_child->set_neighbor(0,
4819  cell->neighbor(0)->child(nbnb));
4820 
4821  // reset neighbor info of all right descendant
4822  // of the left neighbor of cell
4824  left_neighbor = cell->neighbor(0);
4825  while (left_neighbor->has_children())
4826  {
4827  left_neighbor = left_neighbor->child(nbnb);
4828  left_neighbor->set_neighbor(nbnb, first_child);
4829  }
4830  }
4831 
4832  // insert second child
4833  second_child->clear_children();
4834  second_child->set(
4836  next_unused_vertex, cell->vertex_index(1)));
4837  second_child->set_neighbor(0, first_child);
4838  second_child->set_material_id(cell->material_id());
4839  second_child->set_manifold_id(cell->manifold_id());
4840  second_child->set_subdomain_id(subdomainid);
4841  second_child->set_direction_flag(cell->direction_flag());
4842 
4843  if (cell->neighbor(1).state() != IteratorState::valid)
4844  second_child->set_neighbor(1, cell->neighbor(1));
4845  else if (cell->neighbor(1)->is_active())
4846  {
4847  Assert(cell->neighbor(1)->level() <= cell->level(),
4848  ExcInternalError());
4849  second_child->set_neighbor(1, cell->neighbor(1));
4850  }
4851  else
4852  // right neighbor is refined same as above
4853  {
4854  const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4855  second_child->set_neighbor(
4856  1, cell->neighbor(1)->child(nbnb));
4857 
4859  right_neighbor = cell->neighbor(1);
4860  while (right_neighbor->has_children())
4861  {
4862  right_neighbor = right_neighbor->child(nbnb);
4863  right_neighbor->set_neighbor(nbnb, second_child);
4864  }
4865  }
4866  // inform all listeners that cell refinement is done
4867  triangulation.signals.post_refinement_on_cell(cell);
4868  }
4869  }
4870 
4871  // in 1d, we can not have distorted children unless the parent
4872  // was already distorted (that is because we don't use
4873  // boundary information for 1d triangulations). so return an
4874  // empty list
4876  }
4877 
4878 
4883  template <int spacedim>
4886  const bool check_for_distorted_cells)
4887  {
4888  const unsigned int dim = 2;
4889 
4890  // check whether a new level is needed we have to check for
4891  // this on the highest level only (on this, all used cells are
4892  // also active, so we only have to check for this)
4893  {
4895  cell = triangulation.begin_active(triangulation.levels.size() - 1),
4896  endc = triangulation.end();
4897  for (; cell != endc; ++cell)
4898  if (cell->used())
4899  if (cell->refine_flag_set())
4900  {
4901  triangulation.levels.push_back(
4902  std::make_unique<
4904  break;
4905  }
4906  }
4907 
4908  // TODO[WB]: we clear user flags and pointers of lines; we're going
4909  // to use them to flag which lines need refinement
4910  for (typename Triangulation<dim, spacedim>::line_iterator line =
4911  triangulation.begin_line();
4912  line != triangulation.end_line();
4913  ++line)
4914  {
4915  line->clear_user_flag();
4916  line->clear_user_data();
4917  }
4918  // running over all cells and lines count the number
4919  // n_single_lines of lines which can be stored as single
4920  // lines, e.g. inner lines
4921  unsigned int n_single_lines = 0;
4922 
4923  // New lines to be created: number lines which are stored in
4924  // pairs (the children of lines must be stored in pairs)
4925  unsigned int n_lines_in_pairs = 0;
4926 
4927  // check how much space is needed on every level we need not
4928  // check the highest level since either - on the highest level
4929  // no cells are flagged for refinement - there are, but
4930  // prepare_refinement added another empty level
4931  unsigned int needed_vertices = 0;
4932  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4933  {
4934  // count number of flagged cells on this level and compute
4935  // how many new vertices and new lines will be needed
4936  unsigned int needed_cells = 0;
4937 
4939  cell = triangulation.begin_active(level),
4940  endc = triangulation.begin_active(level + 1);
4941  for (; cell != endc; ++cell)
4942  if (cell->refine_flag_set())
4943  {
4944  if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
4945  {
4946  needed_cells += 4;
4947 
4948  // new vertex at center of cell is needed in any
4949  // case
4950  ++needed_vertices;
4951 
4952  // the four inner lines can be stored as singles
4953  n_single_lines += 4;
4954  }
4955  else // cut_x || cut_y
4956  {
4957  // set the flag showing that anisotropic
4958  // refinement is used for at least one cell
4959  triangulation.anisotropic_refinement = true;
4960 
4961  needed_cells += 2;
4962  // no vertex at center
4963 
4964  // the inner line can be stored as single
4965  n_single_lines += 1;
4966  }
4967 
4968  // mark all faces (lines) for refinement; checking
4969  // locally whether the neighbor would also like to
4970  // refine them is rather difficult for lines so we
4971  // only flag them and after visiting all cells, we
4972  // decide which lines need refinement;
4973  for (const unsigned int line_no :
4975  {
4977  cell->refine_flag_set(), line_no) ==
4979  {
4981  line = cell->line(line_no);
4982  if (line->has_children() == false)
4983  line->set_user_flag();
4984  }
4985  }
4986  }
4987 
4988 
4989  // count number of used cells on the next higher level
4990  const unsigned int used_cells =
4991  std::count(triangulation.levels[level + 1]->cells.used.begin(),
4992  triangulation.levels[level + 1]->cells.used.end(),
4993  true);
4994 
4995 
4996  // reserve space for the used_cells cells already existing
4997  // on the next higher level as well as for the
4998  // needed_cells that will be created on that level
4999  triangulation.levels[level + 1]->reserve_space(
5000  used_cells + needed_cells, 2, spacedim);
5001 
5002  // reserve space for needed_cells new quads on the next
5003  // higher level
5004  triangulation.levels[level + 1]->cells.reserve_space(needed_cells,
5005  0);
5006  }
5007 
5008  // now count the lines which were flagged for refinement
5009  for (typename Triangulation<dim, spacedim>::line_iterator line =
5010  triangulation.begin_line();
5011  line != triangulation.end_line();
5012  ++line)
5013  if (line->user_flag_set())
5014  {
5015  Assert(line->has_children() == false, ExcInternalError());
5016  n_lines_in_pairs += 2;
5017  needed_vertices += 1;
5018  }
5019  // reserve space for n_lines_in_pairs new lines. note, that
5020  // we can't reserve space for the single lines here as well,
5021  // as all the space reserved for lines in pairs would be
5022  // counted as unused and we would end up with too little space
5023  // to store all lines. memory reservation for n_single_lines
5024  // can only be done AFTER we refined the lines of the current
5025  // cells
5026  triangulation.faces->lines.reserve_space(n_lines_in_pairs, 0);
5027 
5028  // add to needed vertices how many vertices are already in use
5029  needed_vertices += std::count(triangulation.vertices_used.begin(),
5030  triangulation.vertices_used.end(),
5031  true);
5032  // if we need more vertices: create them, if not: leave the
5033  // array as is, since shrinking is not really possible because
5034  // some of the vertices at the end may be in use
5035  if (needed_vertices > triangulation.vertices.size())
5036  {
5037  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5038  triangulation.vertices_used.resize(needed_vertices, false);
5039  }
5040 
5041 
5042  // Do REFINEMENT on every level; exclude highest level as
5043  // above
5044 
5045  // index of next unused vertex
5046  unsigned int next_unused_vertex = 0;
5047 
5048  // first the refinement of lines. children are stored
5049  // pairwise
5050  {
5051  // only active objects can be refined further
5053  line = triangulation.begin_active_line(),
5054  endl = triangulation.end_line();
5056  next_unused_line = triangulation.begin_raw_line();
5057 
5058  for (; line != endl; ++line)
5059  if (line->user_flag_set())
5060  {
5061  // this line needs to be refined
5062 
5063  // find the next unused vertex and set it
5064  // appropriately
5065  while (triangulation.vertices_used[next_unused_vertex] == true)
5066  ++next_unused_vertex;
5067  Assert(
5068  next_unused_vertex < triangulation.vertices.size(),
5069  ExcMessage(
5070  "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."));
5071  triangulation.vertices_used[next_unused_vertex] = true;
5072 
5073  triangulation.vertices[next_unused_vertex] = line->center(true);
5074 
5075  // now that we created the right point, make up the
5076  // two child lines. To this end, find a pair of
5077  // unused lines
5078  bool pair_found = false;
5079  (void)pair_found;
5080  for (; next_unused_line != endl; ++next_unused_line)
5081  if (!next_unused_line->used() &&
5082  !(++next_unused_line)->used())
5083  {
5084  // go back to the first of the two unused
5085  // lines
5086  --next_unused_line;
5087  pair_found = true;
5088  break;
5089  }
5090  Assert(pair_found, ExcInternalError());
5091 
5092  // there are now two consecutive unused lines, such
5093  // that the children of a line will be consecutive.
5094  // then set the child pointer of the present line
5095  line->set_children(0, next_unused_line->index());
5096 
5097  // set the two new lines
5099  children[2] = {next_unused_line, ++next_unused_line};
5100  // some tests; if any of the iterators should be
5101  // invalid, then already dereferencing will fail
5102  Assert(
5103  children[0]->used() == false,
5104  ExcMessage(
5105  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5106  Assert(
5107  children[1]->used() == false,
5108  ExcMessage(
5109  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5110 
5111  children[0]->set(
5113  line->vertex_index(0), next_unused_vertex));
5114  children[1]->set(
5116  next_unused_vertex, line->vertex_index(1)));
5117 
5118  children[0]->set_used_flag();
5119  children[1]->set_used_flag();
5120  children[0]->clear_children();
5121  children[1]->clear_children();
5122  children[0]->clear_user_data();
5123  children[1]->clear_user_data();
5124  children[0]->clear_user_flag();
5125  children[1]->clear_user_flag();
5126 
5127 
5128  children[0]->set_boundary_id_internal(line->boundary_id());
5129  children[1]->set_boundary_id_internal(line->boundary_id());
5130 
5131  children[0]->set_manifold_id(line->manifold_id());
5132  children[1]->set_manifold_id(line->manifold_id());
5133 
5134  // finally clear flag indicating the need for
5135  // refinement
5136  line->clear_user_flag();
5137  }
5138  }
5139 
5140 
5141  // Now set up the new cells
5142 
5143  // reserve space for inner lines (can be stored as single
5144  // lines)
5145  triangulation.faces->lines.reserve_space(0, n_single_lines);
5146 
5148  cells_with_distorted_children;
5149 
5150  // reset next_unused_line, as now also single empty places in
5151  // the vector can be used
5153  next_unused_line = triangulation.begin_raw_line();
5154 
5155  for (int level = 0;
5156  level < static_cast<int>(triangulation.levels.size()) - 1;
5157  ++level)
5158  {
5159  // Remember: as we don't operate on the finest level,
5160  // begin_*(level+1) is allowed
5162  cell = triangulation.begin_active(level),
5163  endc = triangulation.begin_active(level + 1);
5164 
5166  next_unused_cell = triangulation.begin_raw(level + 1);
5167 
5168  for (; cell != endc; ++cell)
5169  if (cell->refine_flag_set())
5170  {
5171  // set the user flag to indicate, that at least one
5172  // line is at the boundary
5173 
5174  // TODO[Tobias Leicht] find a better place to set
5175  // this flag, so that we do not need so much time to
5176  // check each cell here
5177  if (cell->at_boundary())
5178  cell->set_user_flag();
5179 
5180  // actually set up the children and update neighbor
5181  // information
5182  create_children(triangulation,
5183  next_unused_vertex,
5184  next_unused_line,
5185  next_unused_cell,
5186  cell);
5187 
5188  if ((check_for_distorted_cells == true) &&
5189  has_distorted_children(
5190  cell,
5191  std::integral_constant<int, dim>(),
5192  std::integral_constant<int, spacedim>()))
5193  cells_with_distorted_children.distorted_cells.push_back(
5194  cell);
5195  // inform all listeners that cell refinement is done
5196  triangulation.signals.post_refinement_on_cell(cell);
5197  }
5198  }
5199 
5200  return cells_with_distorted_children;
5201  }
5202 
5203 
5208  template <int spacedim>
5211  const bool check_for_distorted_cells)
5212  {
5213  const unsigned int dim = 3;
5214 
5215  // this function probably also works for spacedim>3 but it
5216  // isn't tested. it will probably be necessary to pull new
5217  // vertices onto the manifold just as we do for the other
5218  // functions above.
5219  Assert(spacedim == 3, ExcNotImplemented());
5220 
5221  // check whether a new level is needed we have to check for
5222  // this on the highest level only (on this, all used cells are
5223  // also active, so we only have to check for this)
5224  {
5226  cell = triangulation.begin_active(triangulation.levels.size() - 1),
5227  endc = triangulation.end();
5228  for (; cell != endc; ++cell)
5229  if (cell->used())
5230  if (cell->refine_flag_set())
5231  {
5232  triangulation.levels.push_back(
5233  std::make_unique<
5235  break;
5236  }
5237  }
5238 
5239 
5240  // first clear user flags for quads and lines; we're going to
5241  // use them to flag which lines and quads need refinement
5242  triangulation.faces->quads.clear_user_data();
5243 
5244  for (typename Triangulation<dim, spacedim>::line_iterator line =
5245  triangulation.begin_line();
5246  line != triangulation.end_line();
5247  ++line)
5248  line->clear_user_flag();
5249  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
5250  triangulation.begin_quad();
5251  quad != triangulation.end_quad();
5252  ++quad)
5253  quad->clear_user_flag();
5254 
5255  // create an array of face refine cases. User indices of faces
5256  // will be set to values corresponding with indices in this
5257  // array.
5258  const RefinementCase<dim - 1> face_refinement_cases[4] = {
5259  RefinementCase<dim - 1>::no_refinement,
5260  RefinementCase<dim - 1>::cut_x,
5261  RefinementCase<dim - 1>::cut_y,
5262  RefinementCase<dim - 1>::cut_xy};
5263 
5264  // check how much space is needed on every level we need not
5265  // check the highest level since either
5266  // - on the highest level no cells are flagged for refinement
5267  // - there are, but prepare_refinement added another empty
5268  // level which then is the highest level
5269 
5270  // variables to hold the number of newly to be created
5271  // vertices, lines and quads. as these are stored globally,
5272  // declare them outside the loop over al levels. we need lines
5273  // and quads in pairs for refinement of old ones and lines and
5274  // quads, that can be stored as single ones, as they are newly
5275  // created in the inside of an existing cell
5276  unsigned int needed_vertices = 0;
5277  unsigned int needed_lines_single = 0;
5278  unsigned int needed_quads_single = 0;
5279  unsigned int needed_lines_pair = 0;
5280  unsigned int needed_quads_pair = 0;
5281  for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5282  {
5283  // count number of flagged cells on this level and compute
5284  // how many new vertices and new lines will be needed
5285  unsigned int new_cells = 0;
5286 
5288  acell = triangulation.begin_active(level),
5289  aendc = triangulation.begin_active(level + 1);
5290  for (; acell != aendc; ++acell)
5291  if (acell->refine_flag_set())
5292  {
5293  RefinementCase<dim> ref_case = acell->refine_flag_set();
5294 
5295  // now for interior vertices, lines and quads, which
5296  // are needed in any case
5297  if (ref_case == RefinementCase<dim>::cut_x ||
5298  ref_case == RefinementCase<dim>::cut_y ||
5299  ref_case == RefinementCase<dim>::cut_z)
5300  {
5301  ++needed_quads_single;
5302  new_cells += 2;
5303  triangulation.anisotropic_refinement = true;
5304  }
5305  else if (ref_case == RefinementCase<dim>::cut_xy ||
5306  ref_case == RefinementCase<dim>::cut_xz ||
5307  ref_case == RefinementCase<dim>::cut_yz)
5308  {
5309  ++needed_lines_single;
5310  needed_quads_single += 4;
5311  new_cells += 4;
5312  triangulation.anisotropic_refinement = true;
5313  }
5314  else if (ref_case == RefinementCase<dim>::cut_xyz)
5315  {
5316  ++needed_vertices;
5317  needed_lines_single += 6;
5318  needed_quads_single += 12;
5319  new_cells += 8;
5320  }
5321  else
5322  {
5323  // we should never get here
5324  Assert(false, ExcInternalError());
5325  }
5326 
5327  // mark all faces for refinement; checking locally
5328  // if and how the neighbor would like to refine
5329  // these is difficult so we only flag them and after
5330  // visiting all cells, we decide which faces need
5331  // which refinement;
5332  for (const unsigned int face :
5334  {
5336  aface = acell->face(face);
5337  // get the RefineCase this faces has for the
5338  // given RefineCase of the cell
5339  RefinementCase<dim - 1> face_ref_case =
5341  ref_case,
5342  face,
5343  acell->face_orientation(face),
5344  acell->face_flip(face),
5345  acell->face_rotation(face));
5346  // only do something, if this face has to be
5347  // refined
5348  if (face_ref_case)
5349  {
5350  if (face_ref_case ==
5352  {
5353  if (aface->number_of_children() < 4)
5354  // we use user_flags to denote needed
5355  // isotropic refinement
5356  aface->set_user_flag();
5357  }
5358  else if (aface->refinement_case() != face_ref_case)
5359  // we use user_indices to denote needed
5360  // anisotropic refinement. note, that we
5361  // can have at most one anisotropic
5362  // refinement case for this face, as
5363  // otherwise prepare_refinement() would
5364  // have changed one of the cells to yield
5365  // isotropic refinement at this
5366  // face. therefore we set the user_index
5367  // uniquely
5368  {
5369  Assert(aface->refinement_case() ==
5371  dim - 1>::isotropic_refinement ||
5372  aface->refinement_case() ==
5373  RefinementCase<dim - 1>::no_refinement,
5374  ExcInternalError());
5375  aface->set_user_index(face_ref_case);
5376  }
5377  }
5378  } // for all faces
5379 
5380  // flag all lines, that have to be refined
5381  for (unsigned int line = 0;
5382  line < GeometryInfo<dim>::lines_per_cell;
5383  ++line)
5385  line) &&
5386  !acell->line(line)->has_children())
5387  acell->line(line)->set_user_flag();
5388 
5389  } // if refine_flag set and for all cells on this level
5390 
5391 
5392  // count number of used cells on the next higher level
5393  const unsigned int used_cells =
5394  std::count(triangulation.levels[level + 1]->cells.used.begin(),
5395  triangulation.levels[level + 1]->cells.used.end(),
5396  true);
5397 
5398 
5399  // reserve space for the used_cells cells already existing
5400  // on the next higher level as well as for the
5401  // 8*flagged_cells that will be created on that level
5402  triangulation.levels[level + 1]->reserve_space(
5403  used_cells + new_cells, 3, spacedim);
5404  // reserve space for 8*flagged_cells new hexes on the next
5405  // higher level
5406  triangulation.levels[level + 1]->cells.reserve_space(new_cells);
5407  } // for all levels
5408  // now count the quads and lines which were flagged for
5409  // refinement
5410  for (typename Triangulation<dim, spacedim>::quad_iterator quad =
5411  triangulation.begin_quad();
5412  quad != triangulation.end_quad();
5413  ++quad)
5414  {
5415  if (quad->user_flag_set())
5416  {
5417  // isotropic refinement: 1 interior vertex, 4 quads
5418  // and 4 interior lines. we store the interior lines
5419  // in pairs in case the face is already or will be
5420  // refined anisotropically
5421  needed_quads_pair += 4;
5422  needed_lines_pair += 4;
5423  needed_vertices += 1;
5424  }
5425  if (quad->user_index())
5426  {
5427  // anisotropic refinement: 1 interior
5428  // line and two quads
5429  needed_quads_pair += 2;
5430  needed_lines_single += 1;
5431  // there is a kind of complicated situation here which
5432  // requires our attention. if the quad is refined
5433  // isotropcally, two of the interior lines will get a
5434  // new mother line - the interior line of our
5435  // anisotropically refined quad. if those two lines
5436  // are not consecutive, we cannot do so and have to
5437  // replace them by two lines that are consecutive. we
5438  // try to avoid that situation, but it may happen
5439  // nevertheless through repeated refinement and
5440  // coarsening. thus we have to check here, as we will
5441  // need some additional space to store those new lines
5442  // in case we need them...
5443  if (quad->has_children())
5444  {
5445  Assert(quad->refinement_case() ==
5446  RefinementCase<dim - 1>::isotropic_refinement,
5447  ExcInternalError());
5448  if ((face_refinement_cases[quad->user_index()] ==
5449  RefinementCase<dim - 1>::cut_x &&
5450  (quad->child(0)->line_index(1) + 1 !=
5451  quad->child(2)->line_index(1))) ||
5452  (face_refinement_cases[quad->user_index()] ==
5453  RefinementCase<dim - 1>::cut_y &&
5454  (quad->child(0)->line_index(3) + 1 !=
5455  quad->child(1)->line_index(3))))
5456  needed_lines_pair += 2;
5457  }
5458  }
5459  }
5460 
5461  for (typename Triangulation<dim, spacedim>::line_iterator line =
5462  triangulation.begin_line();
5463  line != triangulation.end_line();
5464  ++line)
5465  if (line->user_flag_set())
5466  {
5467  needed_lines_pair += 2;
5468  needed_vertices += 1;
5469  }
5470 
5471  // reserve space for needed_lines new lines stored in pairs
5472  triangulation.faces->lines.reserve_space(needed_lines_pair,
5473  needed_lines_single);
5474  // reserve space for needed_quads new quads stored in pairs
5475  triangulation.faces->reserve_space(needed_quads_pair,
5476  needed_quads_single);
5477  triangulation.faces->quads.reserve_space(needed_quads_pair,
5478  needed_quads_single);
5479 
5480 
5481  // add to needed vertices how many vertices are already in use
5482  needed_vertices += std::count(triangulation.vertices_used.begin(),
5483  triangulation.vertices_used.end(),
5484  true);
5485  // if we need more vertices: create them, if not: leave the
5486  // array as is, since shrinking is not really possible because
5487  // some of the vertices at the end may be in use
5488  if (needed_vertices > triangulation.vertices.size())
5489  {
5490  triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5491  triangulation.vertices_used.resize(needed_vertices, false);
5492  }
5493 
5494 
5496  // Before we start with the actual refinement, we do some
5497  // sanity checks if in debug mode. especially, we try to catch
5498  // the notorious problem with lines being twice refined,
5499  // i.e. there are cells adjacent at one line ("around the
5500  // edge", but not at a face), with two cells differing by more
5501  // than one refinement level
5502  //
5503  // this check is very simple to implement here, since we have
5504  // all lines flagged if they shall be refined
5505 #ifdef DEBUG
5506  for (const auto &cell : triangulation.active_cell_iterators())
5507  if (!cell->refine_flag_set())
5508  for (unsigned int line = 0;
5509  line < GeometryInfo<dim>::lines_per_cell;
5510  ++line)
5511  if (cell->line(line)->has_children())
5512  for (unsigned int c = 0; c < 2; ++c)
5513  Assert(cell->line(line)->child(c)->user_flag_set() == false,
5514  ExcInternalError());
5515 #endif
5516 
5518  // Do refinement on every level
5519  //
5520  // To make life a bit easier, we first refine those lines and
5521  // quads that were flagged for refinement and then compose the
5522  // newly to be created cells.
5523  //
5524  // index of next unused vertex
5525  unsigned int next_unused_vertex = 0;
5526 
5527  // first for lines
5528  {
5529  // only active objects can be refined further
5531  line = triangulation.begin_active_line(),
5532  endl = triangulation.end_line();
5534  next_unused_line = triangulation.begin_raw_line();
5535 
5536  for (; line != endl; ++line)
5537  if (line->user_flag_set())
5538  {
5539  // this line needs to be refined
5540 
5541  // find the next unused vertex and set it
5542  // appropriately
5543  while (triangulation.vertices_used[next_unused_vertex] == true)
5544  ++next_unused_vertex;
5545  Assert(
5546  next_unused_vertex < triangulation.vertices.size(),
5547  ExcMessage(
5548  "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."));
5549  triangulation.vertices_used[next_unused_vertex] = true;
5550 
5551  triangulation.vertices[next_unused_vertex] = line->center(true);
5552 
5553  // now that we created the right point, make up the
5554  // two child lines (++ takes care of the end of the
5555  // vector)
5556  next_unused_line =
5557  triangulation.faces->lines.template next_free_pair_object<1>(
5558  triangulation);
5559  Assert(next_unused_line.state() == IteratorState::valid,
5560  ExcInternalError());
5561 
5562  // now we found two consecutive unused lines, such
5563  // that the children of a line will be consecutive.
5564  // then set the child pointer of the present line
5565  line->set_children(0, next_unused_line->index());
5566 
5567  // set the two new lines
5569  children[2] = {next_unused_line, ++next_unused_line};
5570 
5571  // some tests; if any of the iterators should be
5572  // invalid, then already dereferencing will fail
5573  Assert(
5574  children[0]->used() == false,
5575  ExcMessage(
5576  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5577  Assert(
5578  children[1]->used() == false,
5579  ExcMessage(
5580  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5581 
5582  children[0]->set(
5584  line->vertex_index(0), next_unused_vertex));
5585  children[1]->set(
5587  next_unused_vertex, line->vertex_index(1)));
5588 
5589  children[0]->set_used_flag();
5590  children[1]->set_used_flag();
5591  children[0]->clear_children();
5592  children[1]->clear_children();
5593  children[0]->clear_user_data();
5594  children[1]->clear_user_data();
5595  children[0]->clear_user_flag();
5596  children[1]->clear_user_flag();
5597 
5598  children[0]->set_boundary_id_internal(line->boundary_id());
5599  children[1]->set_boundary_id_internal(line->boundary_id());
5600 
5601  children[0]->set_manifold_id(line->manifold_id());
5602  children[1]->set_manifold_id(line->manifold_id());
5603 
5604  // finally clear flag
5605  // indicating the need
5606  // for refinement
5607  line->clear_user_flag();
5608  }
5609  }
5610 
5611 
5613  // now refine marked quads
5615 
5616  // here we encounter several cases:
5617 
5618  // a) the quad is unrefined and shall be refined isotropically
5619 
5620  // b) the quad is unrefined and shall be refined
5621  // anisotropically
5622 
5623  // c) the quad is unrefined and shall be refined both
5624  // anisotropically and isotropically (this is reduced to case
5625  // b) and then case b) for the children again)
5626 
5627  // d) the quad is refined anisotropically and shall be refined
5628  // isotropically (this is reduced to case b) for the
5629  // anisotropic children)
5630 
5631  // e) the quad is refined isotropically and shall be refined
5632  // anisotropically (this is transformed to case c), however we
5633  // might have to renumber/rename children...)
5634 
5635  // we need a loop in cases c) and d), as the anisotropic
5636  // children migt have a lower index than the mother quad
5637  for (unsigned int loop = 0; loop < 2; ++loop)
5638  {
5639  // usually, only active objects can be refined
5640  // further. however, in cases d) and e) that is not true,
5641  // so we have to use 'normal' iterators here
5643  quad = triangulation.begin_quad(),
5644  endq = triangulation.end_quad();
5646  next_unused_line = triangulation.begin_raw_line();
5648  next_unused_quad = triangulation.begin_raw_quad();
5649 
5650  for (; quad != endq; ++quad)
5651  {
5652  if (quad->user_index())
5653  {
5654  RefinementCase<dim - 1> aniso_quad_ref_case =
5655  face_refinement_cases[quad->user_index()];
5656  // there is one unlikely event here, where we
5657  // already have refind the face: if the face was
5658  // refined anisotropically and we want to refine
5659  // it isotropically, both children are flagged for
5660  // anisotropic refinement. however, if those
5661  // children were already flagged for anisotropic
5662  // refinement, they might already be processed and
5663  // refined.
5664  if (aniso_quad_ref_case == quad->refinement_case())
5665  continue;
5666 
5667  Assert(quad->refinement_case() ==
5668  RefinementCase<dim - 1>::cut_xy ||
5669  quad->refinement_case() ==
5670  RefinementCase<dim - 1>::no_refinement,
5671  ExcInternalError());
5672 
5673  // this quad needs to be refined anisotropically
5674  Assert(quad->user_index() ==
5675  RefinementCase<dim - 1>::cut_x ||
5676  quad->user_index() ==
5677  RefinementCase<dim - 1>::cut_y,
5678  ExcInternalError());
5679 
5680  // make the new line interior to the quad
5682  new_line;
5683 
5684  new_line =
5685  triangulation.faces->lines
5686  .template next_free_single_object<1>(triangulation);
5687  Assert(
5688  new_line->used() == false,
5689  ExcMessage(
5690  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5691 
5692  // first collect the
5693  // indices of the vertices:
5694  // *--1--*
5695  // | | |
5696  // | | | cut_x
5697  // | | |
5698  // *--0--*
5699  //
5700  // *-----*
5701  // | |
5702  // 0-----1 cut_y
5703  // | |
5704  // *-----*
5705  unsigned int vertex_indices[2];
5706  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5707  {
5708  vertex_indices[0] =
5709  quad->line(2)->child(0)->vertex_index(1);
5710  vertex_indices[1] =
5711  quad->line(3)->child(0)->vertex_index(1);
5712  }
5713  else
5714  {
5715  vertex_indices[0] =
5716  quad->line(0)->child(0)->vertex_index(1);
5717  vertex_indices[1] =
5718  quad->line(1)->child(0)->vertex_index(1);
5719  }
5720 
5721  new_line->set(
5723  vertex_indices[0], vertex_indices[1]));
5724  new_line->set_used_flag();
5725  new_line->clear_user_flag();
5726  new_line->clear_user_data();
5727  new_line->clear_children();
5728  new_line->set_boundary_id_internal(quad->boundary_id());
5729  new_line->set_manifold_id(quad->manifold_id());
5730 
5731  // child 0 and 1 of a line are switched if the
5732  // line orientation is false. set up a miniature
5733  // table, indicating which child to take for line
5734  // orientations false and true. first index: child
5735  // index in standard orientation, second index:
5736  // line orientation
5737  const unsigned int index[2][2] = {
5738  {1, 0}, // child 0, line_orientation=false and true
5739  {0, 1}}; // child 1, line_orientation=false and true
5740 
5741  // find some space (consecutive) for the two newly
5742  // to be created quads.
5744  new_quads[2];
5745 
5746  next_unused_quad =
5747  triangulation.faces->quads
5748  .template next_free_pair_object<2>(triangulation);
5749  new_quads[0] = next_unused_quad;
5750  Assert(
5751  new_quads[0]->used() == false,
5752  ExcMessage(
5753  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5754 
5755  ++next_unused_quad;
5756  new_quads[1] = next_unused_quad;
5757  Assert(
5758  new_quads[1]->used() == false,
5759  ExcMessage(
5760  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5761 
5762 
5763  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5764  {
5765  new_quads[0]->set(
5767  quad->line_index(0),
5768  new_line->index(),
5769  quad->line(2)
5770  ->child(index[0][quad->line_orientation(2)])
5771  ->index(),
5772  quad->line(3)
5773  ->child(index[0][quad->line_orientation(3)])
5774  ->index()));
5775  new_quads[1]->set(
5777  new_line->index(),
5778  quad->line_index(1),
5779  quad->line(2)
5780  ->child(index[1][quad->line_orientation(2)])
5781  ->index(),
5782  quad->line(3)
5783  ->child(index[1][quad->line_orientation(3)])
5784  ->index()));
5785  }
5786  else
5787  {
5788  new_quads[0]->set(
5790  quad->line(0)
5791  ->child(index[0][quad->line_orientation(0)])
5792  ->index(),
5793  quad->line(1)
5794  ->child(index[0][quad->line_orientation(1)])
5795  ->index(),
5796  quad->line_index(2),
5797  new_line->index()));
5798  new_quads[1]->set(
5800  quad->line(0)
5801  ->child(index[1][quad->line_orientation(0)])
5802  ->index(),
5803  quad->line(1)
5804  ->child(index[1][quad->line_orientation(1)])
5805  ->index(),
5806  new_line->index(),
5807  quad->line_index(3)));
5808  }
5809 
5810  for (const auto &new_quad : new_quads)
5811  {
5812  new_quad->set_used_flag();
5813  new_quad->clear_user_flag();
5814  new_quad->clear_user_data();
5815  new_quad->clear_children();
5816  new_quad->set_boundary_id_internal(quad->boundary_id());
5817  new_quad->set_manifold_id(quad->manifold_id());
5818  // set all line orientations to true, change
5819  // this after the loop, as we have to consider
5820  // different lines for each child
5821  for (unsigned int j = 0;
5822  j < GeometryInfo<dim>::lines_per_face;
5823  ++j)
5824  new_quad->set_line_orientation(j, true);
5825  }
5826  // now set the line orientation of children of
5827  // outer lines correctly, the lines in the
5828  // interior of the refined quad are automatically
5829  // oriented conforming to the standard
5830  new_quads[0]->set_line_orientation(
5831  0, quad->line_orientation(0));
5832  new_quads[0]->set_line_orientation(
5833  2, quad->line_orientation(2));
5834  new_quads[1]->set_line_orientation(
5835  1, quad->line_orientation(1));
5836  new_quads[1]->set_line_orientation(
5837  3, quad->line_orientation(3));
5838  if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5839  {
5840  new_quads[0]->set_line_orientation(
5841  3, quad->line_orientation(3));
5842  new_quads[1]->set_line_orientation(
5843  2, quad->line_orientation(2));
5844  }
5845  else
5846  {
5847  new_quads[0]->set_line_orientation(
5848  1, quad->line_orientation(1));
5849  new_quads[1]->set_line_orientation(
5850  0, quad->line_orientation(0));
5851  }
5852 
5853  // test, whether this face is refined
5854  // isotropically already. if so, set the correct
5855  // children pointers.
5856  if (quad->refinement_case() ==
5857  RefinementCase<dim - 1>::cut_xy)
5858  {
5859  // we will put a new refinemnt level of
5860  // anisotropic refinement between the
5861  // unrefined and isotropically refined quad
5862  // ending up with the same fine quads but
5863  // introducing anisotropically refined ones as
5864  // children of the unrefined quad and mother
5865  // cells of the original fine ones.
5866 
5867  // this process includes the creation of a new
5868  // middle line which we will assign as the
5869  // mother line of two of the existing inner
5870  // lines. If those inner lines are not
5871  // consecutive in memory, we won't find them
5872  // later on, so we have to create new ones
5873  // instead and replace all occurrences of the
5874  // old ones with those new ones. As this is
5875  // kind of ugly, we hope we don't have to do
5876  // it often...
5878  old_child[2];
5879  if (aniso_quad_ref_case ==
5881  {
5882  old_child[0] = quad->child(0)->line(1);
5883  old_child[1] = quad->child(2)->line(1);
5884  }
5885  else
5886  {
5887  Assert(aniso_quad_ref_case ==
5889  ExcInternalError());
5890 
5891  old_child[0] = quad->child(0)->line(3);
5892  old_child[1] = quad->child(1)->line(3);
5893  }
5894 
5895  if (old_child[0]->index() + 1 != old_child[1]->index())
5896  {
5897  // this is exactly the ugly case we taked
5898  // about. so, no coimplaining, lets get
5899  // two new lines and copy all info
5900  typename Triangulation<dim,
5901  spacedim>::raw_line_iterator
5902  new_child[2];
5903 
5904  new_child[0] = new_child[1] =
5905  triangulation.faces->lines
5906  .template next_free_pair_object<1>(
5907  triangulation);
5908  ++new_child[1];
5909 
5910  new_child[0]->set_used_flag();
5911  new_child[1]->set_used_flag();
5912 
5913  const int old_index_0 = old_child[0]->index(),
5914  old_index_1 = old_child[1]->index(),
5915  new_index_0 = new_child[0]->index(),
5916  new_index_1 = new_child[1]->index();
5917 
5918  // loop over all quads and replace the old
5919  // lines
5920  for (unsigned int q = 0;
5921  q < triangulation.faces->quads.n_objects();
5922  ++q)
5923  for (unsigned int l = 0;
5924  l < GeometryInfo<dim>::lines_per_face;
5925  ++l)
5926  {
5927  const int this_index =
5928  triangulation.faces->quads.get_object(q)
5929  .face(l);
5930  if (this_index == old_index_0)
5931  triangulation.faces->quads.get_object(q)
5932  .set_face(l, new_index_0);
5933  else if (this_index == old_index_1)
5934  triangulation.faces->quads.get_object(q)
5935  .set_face(l, new_index_1);
5936  }
5937  // now we have to copy all information of
5938  // the two lines
5939  for (unsigned int i = 0; i < 2; ++i)
5940  {
5941  Assert(!old_child[i]->has_children(),
5942  ExcInternalError());
5943 
5944  new_child[i]->set(
5946  TriaObject<1>(old_child[i]->vertex_index(0),
5947  old_child[i]->vertex_index(
5948  1)));
5949  new_child[i]->set_boundary_id_internal(
5950  old_child[i]->boundary_id());
5951  new_child[i]->set_manifold_id(
5952  old_child[i]->manifold_id());
5953  new_child[i]->set_user_index(
5954  old_child[i]->user_index());
5955  if (old_child[i]->user_flag_set())
5956  new_child[i]->set_user_flag();
5957  else
5958  new_child[i]->clear_user_flag();
5959 
5960  new_child[i]->clear_children();
5961 
5962  old_child[i]->clear_user_flag();
5963  old_child[i]->clear_user_index();
5964  old_child[i]->clear_used_flag();
5965  }
5966  }
5967  // now that we cared about the lines, go on
5968  // with the quads themselves, where we might
5969  // encounter similar situations...
5970  if (aniso_quad_ref_case ==
5972  {
5973  new_line->set_children(
5974  0, quad->child(0)->line_index(1));
5975  Assert(new_line->child(1) ==
5976  quad->child(2)->line(1),
5977  ExcInternalError());
5978  // now evereything is quite
5979  // complicated. we have the children
5980  // numbered according to
5981  //
5982  // *---*---*
5983  // |n+2|n+3|
5984  // *---*---*
5985  // | n |n+1|
5986  // *---*---*
5987  //
5988  // from the original isotropic
5989  // refinement. we have to reorder them as
5990  //
5991  // *---*---*
5992  // |n+1|n+3|
5993  // *---*---*
5994  // | n |n+2|
5995  // *---*---*
5996  //
5997  // such that n and n+1 are consecutive
5998  // children of m and n+2 and n+3 are
5999  // consecutive children of m+1, where m
6000  // and m+1 are given as in
6001  //
6002  // *---*---*
6003  // | | |
6004  // | m |m+1|
6005  // | | |
6006  // *---*---*
6007  //
6008  // this is a bit ugly, of course: loop
6009  // over all cells on all levels and look
6010  // for faces n+1 (switch_1) and n+2
6011  // (switch_2).
6012  const typename Triangulation<dim, spacedim>::
6013  quad_iterator switch_1 = quad->child(1),
6014  switch_2 = quad->child(2);
6015  const int switch_1_index = switch_1->index();
6016  const int switch_2_index = switch_2->index();
6017  for (unsigned int l = 0;
6018  l < triangulation.levels.size();
6019  ++l)
6020  for (unsigned int h = 0;
6021  h <
6022  triangulation.levels[l]->cells.n_objects();
6023  ++h)
6024  for (const unsigned int q :
6026  {
6027  const int face_index =
6028  triangulation.levels[l]
6029  ->cells.get_object(h)
6030  .face(q);
6031  if (face_index == switch_1_index)
6032  triangulation.levels[l]
6033  ->cells.get_object(h)
6034  .set_face(q, switch_2_index);
6035  else if (face_index == switch_2_index)
6036  triangulation.levels[l]
6037  ->cells.get_object(h)
6038  .set_face(q, switch_1_index);
6039  }
6040  // now we have to copy all information of
6041  // the two quads
6042  const unsigned int switch_1_lines[4] = {
6043  switch_1->line_index(0),
6044  switch_1->line_index(1),
6045  switch_1->line_index(2),
6046  switch_1->line_index(3)};
6047  const bool switch_1_line_orientations[4] = {
6048  switch_1->line_orientation(0),
6049  switch_1->line_orientation(1),
6050  switch_1->line_orientation(2),
6051  switch_1->line_orientation(3)};
6052  const types::boundary_id switch_1_boundary_id =
6053  switch_1->boundary_id();
6054  const unsigned int switch_1_user_index =
6055  switch_1->user_index();
6056  const bool switch_1_user_flag =
6057  switch_1->user_flag_set();
6058  const RefinementCase<dim - 1>
6059  switch_1_refinement_case =
6060  switch_1->refinement_case();
6061  const int switch_1_first_child_pair =
6062  (switch_1_refinement_case ?
6063  switch_1->child_index(0) :
6064  -1);
6065  const int switch_1_second_child_pair =
6066  (switch_1_refinement_case ==
6067  RefinementCase<dim - 1>::cut_xy ?
6068  switch_1->child_index(2) :
6069  -1);
6070 
6071  switch_1->set(
6073  2>(switch_2->line_index(0),
6074  switch_2->line_index(1),
6075  switch_2->line_index(2),
6076  switch_2->line_index(3)));
6077  switch_1->set_line_orientation(
6078  0, switch_2->line_orientation(0));
6079  switch_1->set_line_orientation(
6080  1, switch_2->line_orientation(1));
6081  switch_1->set_line_orientation(
6082  2, switch_2->line_orientation(2));
6083  switch_1->set_line_orientation(
6084  3, switch_2->line_orientation(3));
6085  switch_1->set_boundary_id_internal(
6086  switch_2->boundary_id());
6087  switch_1->set_manifold_id(switch_2->manifold_id());
6088  switch_1->set_user_index(switch_2->user_index());
6089  if (switch_2->user_flag_set())
6090  switch_1->set_user_flag();
6091  else
6092  switch_1->clear_user_flag();
6093  switch_1->clear_refinement_case();
6094  switch_1->set_refinement_case(
6095  switch_2->refinement_case());
6096  switch_1->clear_children();
6097  if (switch_2->refinement_case())
6098  switch_1->set_children(0,
6099  switch_2->child_index(0));
6100  if (switch_2->refinement_case() ==
6101  RefinementCase<dim - 1>::cut_xy)
6102  switch_1->set_children(2,
6103  switch_2->child_index(2));
6104 
6105  switch_2->set(
6107  2>(switch_1_lines[0],
6108  switch_1_lines[1],
6109  switch_1_lines[2],
6110  switch_1_lines[3]));
6111  switch_2->set_line_orientation(
6112  0, switch_1_line_orientations[0]);
6113  switch_2->set_line_orientation(
6114  1, switch_1_line_orientations[1]);
6115  switch_2->set_line_orientation(
6116  2, switch_1_line_orientations[2]);
6117  switch_2->set_line_orientation(
6118  3, switch_1_line_orientations[3]);
6119  switch_2->set_boundary_id_internal(
6120  switch_1_boundary_id);
6121  switch_2->set_manifold_id(switch_1->manifold_id());
6122  switch_2->set_user_index(switch_1_user_index);
6123  if (switch_1_user_flag)
6124  switch_2->set_user_flag();
6125  else
6126  switch_2->clear_user_flag();
6127  switch_2->clear_refinement_case();
6128  switch_2->set_refinement_case(
6129  switch_1_refinement_case);
6130  switch_2->clear_children();
6131  switch_2->set_children(0,
6132  switch_1_first_child_pair);
6133  switch_2->set_children(2,
6134  switch_1_second_child_pair);
6135 
6136  new_quads[0]->set_refinement_case(
6138  new_quads[0]->set_children(0, quad->child_index(0));
6139  new_quads[1]->set_refinement_case(
6141  new_quads[1]->set_children(0, quad->child_index(2));
6142  }
6143  else
6144  {
6145  new_quads[0]->set_refinement_case(
6147  new_quads[0]->set_children(0, quad->child_index(0));
6148  new_quads[1]->set_refinement_case(
6150  new_quads[1]->set_children(0, quad->child_index(2));
6151  new_line->set_children(
6152  0, quad->child(0)->line_index(3));
6153  Assert(new_line->child(1) ==
6154  quad->child(1)->line(3),
6155  ExcInternalError());
6156  }
6157  quad->clear_children();
6158  }
6159 
6160  // note these quads as children to the present one
6161  quad->set_children(0, new_quads[0]->index());
6162 
6163  quad->set_refinement_case(aniso_quad_ref_case);
6164 
6165  // finally clear flag indicating the need for
6166  // refinement
6167  quad->clear_user_data();
6168  } // if (anisotropic refinement)
6169 
6170  if (quad->user_flag_set())
6171  {
6172  // this quad needs to be refined isotropically
6173 
6174  // first of all: we only get here in the first run
6175  // of the loop
6176  Assert(loop == 0, ExcInternalError());
6177 
6178  // find the next unused vertex. we'll need this in
6179  // any case
6180  while (triangulation.vertices_used[next_unused_vertex] ==
6181  true)
6182  ++next_unused_vertex;
6183  Assert(
6184  next_unused_vertex < triangulation.vertices.size(),
6185  ExcMessage(
6186  "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."));
6187 
6188  // now: if the quad is refined anisotropically
6189  // already, set the anisotropic refinement flag
6190  // for both children. Additionally, we have to
6191  // refine the inner line, as it is an outer line
6192  // of the two (anisotropic) children
6193  const RefinementCase<dim - 1> quad_ref_case =
6194  quad->refinement_case();
6195 
6196  if (quad_ref_case == RefinementCase<dim - 1>::cut_x ||
6197  quad_ref_case == RefinementCase<dim - 1>::cut_y)
6198  {
6199  // set the 'opposite' refine case for children
6200  quad->child(0)->set_user_index(
6201  RefinementCase<dim - 1>::cut_xy - quad_ref_case);
6202  quad->child(1)->set_user_index(
6203  RefinementCase<dim - 1>::cut_xy - quad_ref_case);
6204  // refine the inner line
6206  middle_line;
6207  if (quad_ref_case == RefinementCase<dim - 1>::cut_x)
6208  middle_line = quad->child(0)->line(1);
6209  else
6210  middle_line = quad->child(0)->line(3);
6211 
6212  // if the face has been refined
6213  // anisotropically in the last refinement step
6214  // it might be, that it is flagged already and
6215  // that the middle line is thus refined
6216  // already. if not create children.
6217  if (!middle_line->has_children())
6218  {
6219  // set the middle vertex
6220  // appropriately. double refinement of
6221  // quads can only happen in the interior
6222  // of the domain, so we need not care
6223  // about boundary quads here
6224  triangulation.vertices[next_unused_vertex] =
6225  middle_line->center(true);
6226  triangulation.vertices_used[next_unused_vertex] =
6227  true;
6228 
6229  // now search a slot for the two
6230  // child lines
6231  next_unused_line =
6232  triangulation.faces->lines
6233  .template next_free_pair_object<1>(
6234  triangulation);
6235 
6236  // set the child pointer of the present
6237  // line
6238  middle_line->set_children(
6239  0, next_unused_line->index());
6240 
6241  // set the two new lines
6242  const typename Triangulation<dim, spacedim>::
6243  raw_line_iterator children[2] = {
6244  next_unused_line, ++next_unused_line};
6245 
6246  // some tests; if any of the iterators
6247  // should be invalid, then already
6248  // dereferencing will fail
6249  Assert(
6250  children[0]->used() == false,
6251  ExcMessage(
6252  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6253  Assert(
6254  children[1]->used() == false,
6255  ExcMessage(
6256  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6257 
6258  children[0]->set(
6260  1>(middle_line->vertex_index(0),
6261  next_unused_vertex));
6262  children[1]->set(
6264  1>(next_unused_vertex,
6265  middle_line->vertex_index(1)));
6266 
6267  children[0]->set_used_flag();
6268  children[1]->set_used_flag();
6269  children[0]->clear_children();
6270  children[1]->clear_children();
6271  children[0]->clear_user_data();
6272  children[1]->clear_user_data();
6273  children[0]->clear_user_flag();
6274  children[1]->clear_user_flag();
6275 
6276  children[0]->set_boundary_id_internal(
6277  middle_line->boundary_id());
6278  children[1]->set_boundary_id_internal(
6279  middle_line->boundary_id());
6280 
6281  children[0]->set_manifold_id(
6282  middle_line->manifold_id());
6283  children[1]->set_manifold_id(
6284  middle_line->manifold_id());
6285  }
6286  // now remove the flag from the quad and go to
6287  // the next quad, the actual refinement of the
6288  // quad takes place later on in this pass of
6289  // the loop or in the next one
6290  quad->clear_user_flag();
6291  continue;
6292  } // if (several refinement cases)
6293 
6294  // if we got here, we have an unrefined quad and
6295  // have to do the usual work like in an purely
6296  // isotropic refinement
6297  Assert(quad_ref_case ==
6299  ExcInternalError());
6300 
6301  // set the middle vertex appropriately: it might be that
6302  // the quad itself is not at the boundary, but that one of
6303  // its lines actually is. in this case, the newly created
6304  // vertices at the centers of the lines are not
6305  // necessarily the mean values of the adjacent vertices,
6306  // so do not compute the new vertex as the mean value of
6307  // the 4 vertices of the face, but rather as a weighted
6308  // mean value of the 8 vertices which we already have (the
6309  // four old ones, and the four ones inserted as middle
6310  // points for the four lines). summing up some more points
6311  // is generally cheaper than first asking whether one of
6312  // the lines is at the boundary
6313  //
6314  // note that the exact weights are chosen such as to
6315  // minimize the distortion of the four new quads from the
6316  // optimal shape. their description uses the formulas
6317  // underlying the TransfiniteInterpolationManifold
6318  // implementation
6319  triangulation.vertices[next_unused_vertex] =
6320  quad->center(true, true);
6321  triangulation.vertices_used[next_unused_vertex] = true;
6322 
6323  // now that we created the right point, make up
6324  // the four lines interior to the quad (++ takes
6325  // care of the end of the vector)
6327  new_lines[4];
6328 
6329  for (unsigned int i = 0; i < 4; ++i)
6330  {
6331  if (i % 2 == 0)
6332  // search a free pair of lines for 0. and
6333  // 2. line, so that two of them end up
6334  // together, which is necessary if later on
6335  // we want to refine the quad
6336  // anisotropically and the two lines end up
6337  // as children of new line
6338  next_unused_line =
6339  triangulation.faces->lines
6340  .template next_free_pair_object<1>(triangulation);
6341 
6342  new_lines[i] = next_unused_line;
6343  ++next_unused_line;
6344 
6345  Assert(
6346  new_lines[i]->used() == false,
6347  ExcMessage(
6348  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6349  }
6350 
6351  // set the data of the four lines. first collect
6352  // the indices of the five vertices:
6353  //
6354  // *--3--*
6355  // | | |
6356  // 0--4--1
6357  // | | |
6358  // *--2--*
6359  //
6360  // the lines are numbered as follows:
6361  //
6362  // *--*--*
6363  // | 1 |
6364  // *2-*-3*
6365  // | 0 |
6366  // *--*--*
6367 
6368  const unsigned int vertex_indices[5] = {
6369  quad->line(0)->child(0)->vertex_index(1),
6370  quad->line(1)->child(0)->vertex_index(1),
6371  quad->line(2)->child(0)->vertex_index(1),
6372  quad->line(3)->child(0)->vertex_index(1),
6373  next_unused_vertex};
6374 
6375  new_lines[0]->set(
6377  vertex_indices[2], vertex_indices[4]));
6378  new_lines[1]->set(
6380  vertex_indices[4], vertex_indices[3]));
6381  new_lines[2]->set(
6383  vertex_indices[0], vertex_indices[4]));
6384  new_lines[3]->set(
6386  vertex_indices[4], vertex_indices[1]));
6387 
6388  for (const auto &new_line : new_lines)
6389  {
6390  new_line->set_used_flag();
6391  new_line->clear_user_flag();
6392  new_line->clear_user_data();
6393  new_line->clear_children();
6394  new_line->set_boundary_id_internal(quad->boundary_id());
6395  new_line->set_manifold_id(quad->manifold_id());
6396  }
6397 
6398  // now for the quads. again, first collect some
6399  // data about the indices of the lines, with the
6400  // following numbering:
6401  //
6402  // .-6-.-7-.
6403  // 1 9 3
6404  // .-10.11-.
6405  // 0 8 2
6406  // .-4-.-5-.
6407 
6408  // child 0 and 1 of a line are switched if the
6409  // line orientation is false. set up a miniature
6410  // table, indicating which child to take for line
6411  // orientations false and true. first index: child
6412  // index in standard orientation, second index:
6413  // line orientation
6414  const unsigned int index[2][2] = {
6415  {1, 0}, // child 0, line_orientation=false and true
6416  {0, 1}}; // child 1, line_orientation=false and true
6417 
6418  const int line_indices[12] = {
6419  quad->line(0)
6420  ->child(index[0][quad->line_orientation(0)])
6421  ->index(),
6422  quad->line(0)
6423  ->child(index[1][quad->line_orientation(0)])
6424  ->index(),
6425  quad->line(1)
6426  ->child(index[0][quad->line_orientation(1)])
6427  ->index(),
6428  quad->line(1)
6429  ->child(index[1][quad->line_orientation(1)])
6430  ->index(),
6431  quad->line(2)
6432  ->child(index[0][quad->line_orientation(2)])
6433  ->index(),
6434  quad->line(2)
6435  ->child(index[1][quad->line_orientation(2)])
6436  ->index(),
6437  quad->line(3)
6438  ->child(index[0][quad->line_orientation(3)])
6439  ->index(),
6440  quad->line(3)
6441  ->child(index[1][quad->line_orientation(3)])
6442  ->index(),
6443  new_lines[0]->index(),
6444  new_lines[1]->index(),
6445  new_lines[2]->index(),
6446  new_lines[3]->index()};
6447 
6448  // find some space (consecutive)
6449  // for the first two newly to be
6450  // created quads.
6452  new_quads[4];
6453 
6454  next_unused_quad =
6455  triangulation.faces->quads
6456  .template next_free_pair_object<2>(triangulation);
6457 
6458  new_quads[0] = next_unused_quad;
6459  Assert(
6460  new_quads[0]->used() == false,
6461  ExcMessage(
6462  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6463 
6464  ++next_unused_quad;
6465  new_quads[1] = next_unused_quad;
6466  Assert(
6467  new_quads[1]->used() == false,
6468  ExcMessage(
6469  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6470 
6471  next_unused_quad =
6472  triangulation.faces->quads
6473  .template next_free_pair_object<2>(triangulation);
6474  new_quads[2] = next_unused_quad;
6475  Assert(
6476  new_quads[2]->used() == false,
6477  ExcMessage(
6478  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6479 
6480  ++next_unused_quad;
6481  new_quads[3] = next_unused_quad;
6482  Assert(
6483  new_quads[3]->used() == false,
6484  ExcMessage(
6485  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6486 
6487  // note these quads as children to the present one
6488  quad->set_children(0, new_quads[0]->index());
6489  quad->set_children(2, new_quads[2]->index());
6490  new_quads[0]->set(
6492  line_indices[0],
6493  line_indices[8],
6494  line_indices[4],
6495  line_indices[10]));
6496 
6497  quad->set_refinement_case(RefinementCase<2>::cut_xy);
6498 
6499  new_quads[0]->set(
6501  line_indices[0],
6502  line_indices[8],
6503  line_indices[4],
6504  line_indices[10]));
6505  new_quads[1]->set(
6507  line_indices[8],
6508  line_indices[2],
6509  line_indices[5],
6510  line_indices[11]));
6511  new_quads[2]->set(
6513  line_indices[1],
6514  line_indices[9],
6515  line_indices[10],
6516  line_indices[6]));
6517  new_quads[3]->set(
6519  line_indices[9],
6520  line_indices[3],
6521  line_indices[11],
6522  line_indices[7]));
6523  for (const auto &new_quad : new_quads)
6524  {
6525  new_quad->set_used_flag();
6526  new_quad->clear_user_flag();
6527  new_quad->clear_user_data();
6528  new_quad->clear_children();
6529  new_quad->set_boundary_id_internal(quad->boundary_id());
6530  new_quad->set_manifold_id(quad->manifold_id());
6531  // set all line orientations to true, change
6532  // this after the loop, as we have to consider
6533  // different lines for each child
6534  for (unsigned int j = 0;
6535  j < GeometryInfo<dim>::lines_per_face;
6536  ++j)
6537  new_quad->set_line_orientation(j, true);
6538  }
6539  // now set the line orientation of children of
6540  // outer lines correctly, the lines in the
6541  // interior of the refined quad are automatically
6542  // oriented conforming to the standard
6543  new_quads[0]->set_line_orientation(
6544  0, quad->line_orientation(0));
6545  new_quads[0]->set_line_orientation(
6546  2, quad->line_orientation(2));
6547  new_quads[1]->set_line_orientation(
6548  1, quad->line_orientation(1));
6549  new_quads[1]->set_line_orientation(
6550  2, quad->line_orientation(2));
6551  new_quads[2]->set_line_orientation(
6552  0, quad->line_orientation(0));
6553  new_quads[2]->set_line_orientation(
6554  3, quad->line_orientation(3));
6555  new_quads[3]->set_line_orientation(
6556  1, quad->line_orientation(1));
6557  new_quads[3]->set_line_orientation(
6558  3, quad->line_orientation(3));
6559 
6560  // finally clear flag indicating the need for
6561  // refinement
6562  quad->clear_user_flag();
6563  } // if (isotropic refinement)
6564  } // for all quads
6565  } // looped two times over all quads, all quads refined now
6566 
6568  // Now, finally, set up the new
6569  // cells
6571 
6573  cells_with_distorted_children;
6574 
6575  for (unsigned int level = 0; level != triangulation.levels.size() - 1;
6576  ++level)
6577  {
6578  // only active objects can be refined further; remember
6579  // that we won't operate on the finest level, so
6580  // triangulation.begin_*(level+1) is allowed
6582  hex = triangulation.begin_active_hex(level),
6583  endh = triangulation.begin_active_hex(level + 1);
6585  next_unused_hex = triangulation.begin_raw_hex(level + 1);
6586 
6587  for (; hex != endh; ++hex)
6588  if (hex->refine_flag_set())
6589  {
6590  // this hex needs to be refined
6591 
6592  // clear flag indicating the need for refinement. do
6593  // it here already, since we can't do it anymore
6594  // once the cell has children
6595  const RefinementCase<dim> ref_case = hex->refine_flag_set();
6596  hex->clear_refine_flag();
6597  hex->set_refinement_case(ref_case);
6598 
6599  // depending on the refine case we might have to
6600  // create additional vertices, lines and quads
6601  // interior of the hex before the actual children
6602  // can be set up.
6603 
6604  // in a first step: reserve the needed space for
6605  // lines, quads and hexes and initialize them
6606  // correctly
6607 
6608  unsigned int n_new_lines = 0;
6609  unsigned int n_new_quads = 0;
6610  unsigned int n_new_hexes = 0;
6611  switch (ref_case)
6612  {
6616  n_new_lines = 0;
6617  n_new_quads = 1;
6618  n_new_hexes = 2;
6619  break;
6623  n_new_lines = 1;
6624  n_new_quads = 4;
6625  n_new_hexes = 4;
6626  break;
6628  n_new_lines = 6;
6629  n_new_quads = 12;
6630  n_new_hexes = 8;
6631  break;
6632  default:
6633  Assert(false, ExcInternalError());
6634  break;
6635  }
6636 
6637  // find some space for the newly to be created
6638  // interior lines and initialize them.
6639  std::vector<
6641  new_lines(n_new_lines);
6642  for (unsigned int i = 0; i < n_new_lines; ++i)
6643  {
6644  new_lines[i] =
6645  triangulation.faces->lines
6646  .template next_free_single_object<1>(triangulation);
6647 
6648  Assert(
6649  new_lines[i]->used() == false,
6650  ExcMessage(
6651  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6652  new_lines[i]->set_used_flag();
6653  new_lines[i]->clear_user_flag();
6654  new_lines[i]->clear_user_data();
6655  new_lines[i]->clear_children();
6656  // interior line
6657  new_lines[i]->set_boundary_id_internal(
6659  // they inherit geometry description of the hex they
6660  // belong to
6661  new_lines[i]->set_manifold_id(hex->manifold_id());
6662  }
6663 
6664  // find some space for the newly to be created
6665  // interior quads and initialize them.
6666  std::vector<
6668  new_quads(n_new_quads);
6669  for (unsigned int i = 0; i < n_new_quads; ++i)
6670  {
6671  new_quads[i] =
6672  triangulation.faces->quads
6673  .template next_free_single_object<2>(triangulation);
6674 
6675  Assert(
6676  new_quads[i]->used() == fa