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