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