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