Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
19
21
26#include <deal.II/grid/tria.h>
31
33#include <deal.II/lac/vector.h>
34
35#include <algorithm>
36#include <array>
37#include <cmath>
38#include <functional>
39#include <list>
40#include <map>
41#include <memory>
42#include <numeric>
43
44
46
47
48namespace internal
49{
50 namespace TriangulationImplementation
51 {
53 : n_levels(0)
54 , n_lines(0)
55 , n_active_lines(0)
56 // all other fields are
57 // default constructed
58 {}
59
60
61
62 std::size_t
64 {
69 MemoryConsumption::memory_consumption(n_active_lines_level));
70 }
71
72
74 : n_quads(0)
75 , n_active_quads(0)
76 // all other fields are
77 // default constructed
78 {}
79
80
81
82 std::size_t
84 {
89 MemoryConsumption::memory_consumption(n_active_quads_level));
90 }
91
92
93
95 : n_hexes(0)
96 , n_active_hexes(0)
97 // all other fields are
98 // default constructed
99 {}
100
101
102
103 std::size_t
105 {
110 MemoryConsumption::memory_consumption(n_active_hexes_level));
111 }
112 } // namespace TriangulationImplementation
113} // namespace internal
114
115// anonymous namespace for internal helper functions
116namespace
117{
118 // return whether the given cell is
119 // patch_level_1, i.e. determine
120 // whether either all or none of
121 // its children are further
122 // refined. this function can only
123 // be called for non-active cells.
124 template <int dim, int spacedim>
125 bool
126 cell_is_patch_level_1(
128 {
129 Assert(cell->is_active() == false, ExcInternalError());
130
131 unsigned int n_active_children = 0;
132 for (unsigned int i = 0; i < cell->n_children(); ++i)
133 if (cell->child(i)->is_active())
134 ++n_active_children;
135
136 return (n_active_children == 0) ||
137 (n_active_children == cell->n_children());
138 }
139
140
141
142 // return, whether a given @p cell will be
143 // coarsened, which is the case if all
144 // children are active and have their coarsen
145 // flag set. In case only part of the coarsen
146 // flags are set, remove them.
147 template <int dim, int spacedim>
148 bool
149 cell_will_be_coarsened(
151 {
152 // only cells with children should be
153 // considered for coarsening
154
155 if (cell->has_children())
156 {
157 unsigned int children_to_coarsen = 0;
158 const unsigned int n_children = cell->n_children();
159
160 for (unsigned int c = 0; c < n_children; ++c)
161 if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
162 ++children_to_coarsen;
163 if (children_to_coarsen == n_children)
164 return true;
165 else
166 for (unsigned int c = 0; c < n_children; ++c)
167 if (cell->child(c)->is_active())
168 cell->child(c)->clear_coarsen_flag();
169 }
170 // no children, so no coarsening
171 // possible. however, no children also
172 // means that this cell will be in the same
173 // state as if it had children and was
174 // coarsened. So, what should we return -
175 // false or true?
176 // make sure we do not have to do this at
177 // all...
178 Assert(cell->has_children(), ExcInternalError());
179 // ... and then simply return false
180 return false;
181 }
182
183
184 // return, whether the face @p face_no of the
185 // given @p cell will be refined after the
186 // current refinement step, considering
187 // refine and coarsen flags and considering
188 // only those refinemnts that will be caused
189 // by the neighboring cell.
190
191 // this function is used on both active cells
192 // and cells with children. on cells with
193 // children it also of interest to know 'how'
194 // the face will be refined. thus there is an
195 // additional third argument @p
196 // expected_face_ref_case returning just
197 // that. be aware, that this variable will
198 // only contain useful information if this
199 // function is called for an active cell.
200 //
201 // thus, this is an internal function, users
202 // should call one of the two alternatives
203 // following below.
204 template <int dim, int spacedim>
205 bool
206 face_will_be_refined_by_neighbor_internal(
208 const unsigned int face_no,
209 RefinementCase<dim - 1> &expected_face_ref_case)
210 {
211 // first of all: set the default value for
212 // expected_face_ref_case, which is no
213 // refinement at all
214 expected_face_ref_case = RefinementCase<dim - 1>::no_refinement;
215
216 const typename Triangulation<dim, spacedim>::cell_iterator neighbor =
217 cell->neighbor(face_no);
218
219 // If we are at the boundary, there is no
220 // neighbor which could refine the face
221 if (neighbor.state() != IteratorState::valid)
222 return false;
223
224 if (neighbor->has_children())
225 {
226 // if the neighbor is refined, it may be
227 // coarsened. if so, then it won't refine
228 // the face, no matter what else happens
229 if (cell_will_be_coarsened(neighbor))
230 return false;
231 else
232 // if the neighbor is refined, then it
233 // is also refined at our current
234 // face. It will stay so without
235 // coarsening, so return true in that
236 // case.
237 {
238 expected_face_ref_case = cell->face(face_no)->refinement_case();
239 return true;
240 }
241 }
242
243 // now, the neighbor is not refined, but
244 // perhaps it will be
245 const RefinementCase<dim> nb_ref_flag = neighbor->refine_flag_set();
246 if (nb_ref_flag != RefinementCase<dim>::no_refinement)
247 {
248 // now we need to know, which of the
249 // neighbors faces points towards us
250 const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
251 // check, whether the cell will be
252 // refined in a way that refines our
253 // face
254 const RefinementCase<dim - 1> face_ref_case =
256 nb_ref_flag,
257 neighbor_neighbor,
258 neighbor->face_orientation(neighbor_neighbor),
259 neighbor->face_flip(neighbor_neighbor),
260 neighbor->face_rotation(neighbor_neighbor));
261 if (face_ref_case != RefinementCase<dim - 1>::no_refinement)
262 {
264 neighbor_face = neighbor->face(neighbor_neighbor);
265 const int this_face_index = cell->face_index(face_no);
266
267 // there are still two basic
268 // possibilities here: the neighbor
269 // might be coarser or as coarse
270 // as we are
271 if (neighbor_face->index() == this_face_index)
272 // the neighbor is as coarse as
273 // we are and will be refined at
274 // the face of consideration, so
275 // return true
276 {
277 expected_face_ref_case = face_ref_case;
278 return true;
279 }
280 else
281 {
282 // the neighbor is coarser.
283 // this is the most complicated
284 // case. It might be, that the
285 // neighbor's face will be
286 // refined, but that we will
287 // not see this, as we are
288 // refined in a similar way.
289
290 // so, the neighbor's face must
291 // have children. check, if our
292 // cell's face is one of these
293 // (it could also be a
294 // grand_child)
295 for (unsigned int c = 0; c < neighbor_face->n_children(); ++c)
296 if (neighbor_face->child_index(c) == this_face_index)
297 {
298 // if the flagged refine
299 // case of the face is a
300 // subset or the same as
301 // the current refine case,
302 // then the face, as seen
303 // from our cell, won't be
304 // refined by the neighbor
305 if ((neighbor_face->refinement_case() | face_ref_case) ==
306 neighbor_face->refinement_case())
307 return false;
308 else
309 {
310 // if we are active, we
311 // must be an
312 // anisotropic child
313 // and the coming
314 // face_ref_case is
315 // isotropic. Thus,
316 // from our cell we
317 // will see exactly the
318 // opposite refine case
319 // that the face has
320 // now...
321 Assert(
322 face_ref_case ==
325 expected_face_ref_case =
326 ~neighbor_face->refinement_case();
327 return true;
328 }
329 }
330
331 // so, obviously we were not
332 // one of the children, but a
333 // grandchild. This is only
334 // possible in 3d.
335 Assert(dim == 3, ExcInternalError());
336 // In that case, however, no
337 // matter what the neighbor
338 // does, it won't be finer
339 // after the next refinement
340 // step.
341 return false;
342 }
343 } // if face will be refined
344 } // if neighbor is flagged for refinement
345
346 // no cases left, so the neighbor will not
347 // refine the face
348 return false;
349 }
350
351 // version of above function for both active
352 // and non-active cells
353 template <int dim, int spacedim>
354 bool
355 face_will_be_refined_by_neighbor(
357 const unsigned int face_no)
358 {
359 RefinementCase<dim - 1> dummy = RefinementCase<dim - 1>::no_refinement;
360 return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
361 }
362
363 // version of above function for active cells
364 // only. Additionally returning the refine
365 // case (to come) of the face under
366 // consideration
367 template <int dim, int spacedim>
368 bool
369 face_will_be_refined_by_neighbor(
371 const unsigned int face_no,
372 RefinementCase<dim - 1> &expected_face_ref_case)
373 {
374 return face_will_be_refined_by_neighbor_internal(cell,
375 face_no,
376 expected_face_ref_case);
377 }
378
379
380
381 template <int dim, int spacedim>
382 bool
383 satisfies_level1_at_vertex_rule(
385 {
386 std::vector<unsigned int> min_adjacent_cell_level(
387 triangulation.n_vertices(), triangulation.n_levels());
388 std::vector<unsigned int> max_adjacent_cell_level(
389 triangulation.n_vertices(), 0);
390
391 for (const auto &cell : triangulation.active_cell_iterators())
392 for (const unsigned int v : cell->vertex_indices())
393 {
394 min_adjacent_cell_level[cell->vertex_index(v)] =
395 std::min<unsigned int>(
396 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
397 max_adjacent_cell_level[cell->vertex_index(v)] =
398 std::max<unsigned int>(
399 min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
400 }
401
402 for (unsigned int k = 0; k < triangulation.n_vertices(); ++k)
403 if (triangulation.vertex_used(k))
404 if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
405 return false;
406 return true;
407 }
408
409
410
417 template <int dim, int spacedim>
418 std::vector<unsigned int>
419 count_cells_bounded_by_line(const Triangulation<dim, spacedim> &triangulation)
420 {
421 if (dim >= 2)
422 {
423 std::vector<unsigned int> line_cell_count(triangulation.n_raw_lines(),
424 0);
425 for (const auto &cell : triangulation.cell_iterators())
426 for (unsigned int l = 0; l < cell->n_lines(); ++l)
427 ++line_cell_count[cell->line_index(l)];
428 return line_cell_count;
429 }
430 else
431 return std::vector<unsigned int>();
432 }
433
434
435
442 template <int dim, int spacedim>
443 std::vector<unsigned int>
444 count_cells_bounded_by_quad(const Triangulation<dim, spacedim> &triangulation)
445 {
446 if (dim >= 3)
447 {
448 std::vector<unsigned int> quad_cell_count(triangulation.n_raw_quads(),
449 0);
450 for (const auto &cell : triangulation.cell_iterators())
451 for (unsigned int q : cell->face_indices())
452 ++quad_cell_count[cell->quad_index(q)];
453 return quad_cell_count;
454 }
455 else
456 return {};
457 }
458
459
460
472 void
473 reorder_compatibility(const std::vector<CellData<1>> &, const SubCellData &)
474 {
475 // nothing to do here: the format
476 // hasn't changed for 1d
477 }
478
479
480 void reorder_compatibility(std::vector<CellData<2>> &cells,
481 const SubCellData &)
482 {
483 for (auto &cell : cells)
484 if (cell.vertices.size() == GeometryInfo<2>::vertices_per_cell)
485 std::swap(cell.vertices[2], cell.vertices[3]);
486 }
487
488
489 void reorder_compatibility(std::vector<CellData<3>> &cells,
490 SubCellData & subcelldata)
491 {
492 unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
493 for (auto &cell : cells)
494 if (cell.vertices.size() == GeometryInfo<3>::vertices_per_cell)
495 {
496 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
497 tmp[i] = cell.vertices[i];
498 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
499 cell.vertices[GeometryInfo<3>::ucd_to_deal[i]] = tmp[i];
500 }
501
502 // now points in boundary quads
503 for (auto &boundary_quad : subcelldata.boundary_quads)
504 if (boundary_quad.vertices.size() == GeometryInfo<2>::vertices_per_cell)
505 std::swap(boundary_quad.vertices[2], boundary_quad.vertices[3]);
506 }
507
508
509
527 template <int dim, int spacedim>
528 unsigned int
529 middle_vertex_index(
531 {
532 if (line->has_children())
533 return line->child(0)->vertex_index(1);
535 }
536
537
538 template <int dim, int spacedim>
539 unsigned int
540 middle_vertex_index(
542 {
543 switch (static_cast<unsigned char>(quad->refinement_case()))
544 {
546 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
547 break;
549 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
550 break;
552 return quad->child(0)->vertex_index(3);
553 break;
554 default:
555 break;
556 }
558 }
559
560
561 template <int dim, int spacedim>
562 unsigned int
563 middle_vertex_index(
565 {
566 switch (static_cast<unsigned char>(hex->refinement_case()))
567 {
569 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
570 break;
572 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
573 break;
575 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
576 break;
578 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
579 break;
581 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
582 break;
584 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
585 break;
587 return hex->child(0)->vertex_index(7);
588 break;
589 default:
590 break;
591 }
593 }
594
595
608 template <class TRIANGULATION>
609 inline typename TRIANGULATION::DistortedCellList
610 collect_distorted_coarse_cells(const TRIANGULATION &)
611 {
612 return typename TRIANGULATION::DistortedCellList();
613 }
614
615
616
625 template <int dim>
627 collect_distorted_coarse_cells(const Triangulation<dim, dim> &triangulation)
628 {
629 typename Triangulation<dim, dim>::DistortedCellList distorted_cells;
630 for (const auto &cell : triangulation.cell_iterators_on_level(0))
631 {
633 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
634 vertices[i] = cell->vertex(i);
635
638
639 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
640 if (determinants[i] <= 1e-9 * std::pow(cell->diameter(), 1. * dim))
641 {
642 distorted_cells.distorted_cells.push_back(cell);
643 break;
644 }
645 }
646
647 return distorted_cells;
648 }
649
650
657 template <int dim>
658 bool
659 has_distorted_children(
660 const typename Triangulation<dim, dim>::cell_iterator &cell)
661 {
662 Assert(cell->has_children(), ExcInternalError());
663
664 for (unsigned int c = 0; c < cell->n_children(); ++c)
665 {
667 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
668 vertices[i] = cell->child(c)->vertex(i);
669
672
673 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
674 if (determinants[i] <=
675 1e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
676 return true;
677 }
678
679 return false;
680 }
681
682
690 template <int dim, int spacedim>
691 bool
692 has_distorted_children(
694 {
695 return false;
696 }
697
698
699 template <int dim, int spacedim>
700 void
701 update_periodic_face_map_recursively(
702 const typename Triangulation<dim, spacedim>::cell_iterator &cell_1,
703 const typename Triangulation<dim, spacedim>::cell_iterator &cell_2,
704 unsigned int n_face_1,
705 unsigned int n_face_2,
706 const std::bitset<3> & orientation,
707 typename std::map<
709 unsigned int>,
710 std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
711 unsigned int>,
712 std::bitset<3>>> &periodic_face_map)
713 {
714 using FaceIterator = typename Triangulation<dim, spacedim>::face_iterator;
715 const FaceIterator face_1 = cell_1->face(n_face_1);
716 const FaceIterator face_2 = cell_2->face(n_face_2);
717
718 const bool face_orientation = orientation[0];
719 const bool face_flip = orientation[1];
720 const bool face_rotation = orientation[2];
721
722 Assert((dim != 1) || (face_orientation == true && face_flip == false &&
723 face_rotation == false),
724 ExcMessage("The supplied orientation "
725 "(face_orientation, face_flip, face_rotation) "
726 "is invalid for 1D"));
727
728 Assert((dim != 2) || (face_orientation == true && face_rotation == false),
729 ExcMessage("The supplied orientation "
730 "(face_orientation, face_flip, face_rotation) "
731 "is invalid for 2D"));
732
733 Assert(face_1 != face_2, ExcMessage("face_1 and face_2 are equal!"));
734
735 Assert(face_1->at_boundary() && face_2->at_boundary(),
736 ExcMessage("Periodic faces must be on the boundary"));
737
738 // Check if the requirement that each edge can only have at most one hanging
739 // node, and as a consequence neighboring cells can differ by at most
740 // one refinement level is enforced. In 1d, there are no hanging nodes and
741 // so neighboring cells can differ by more than one refinement level.
742 Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
744
745 // insert periodic face pair for both cells
746 using CellFace =
747 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
748 unsigned int>;
749 const CellFace cell_face_1(cell_1, n_face_1);
750 const CellFace cell_face_2(cell_2, n_face_2);
751 const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
752 cell_face_2, orientation);
753
754 const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
755 periodic_faces(cell_face_1, cell_face_orientation_2);
756
757 // Only one periodic neighbor is allowed
758 Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
759 periodic_face_map.insert(periodic_faces);
760
761 if (dim == 1)
762 {
763 if (cell_1->has_children())
764 {
765 if (cell_2->has_children())
766 {
767 update_periodic_face_map_recursively<dim, spacedim>(
768 cell_1->child(n_face_1),
769 cell_2->child(n_face_2),
770 n_face_1,
771 n_face_2,
772 orientation,
773 periodic_face_map);
774 }
775 else // only face_1 has children
776 {
777 update_periodic_face_map_recursively<dim, spacedim>(
778 cell_1->child(n_face_1),
779 cell_2,
780 n_face_1,
781 n_face_2,
782 orientation,
783 periodic_face_map);
784 }
785 }
786 }
787 else // dim == 2 || dim == 3
788 {
789 // A lookup table on how to go through the child cells depending on the
790 // orientation:
791 // see Documentation of GeometryInfo for details
792
793 static const int lookup_table_2d[2][2] =
794 // flip:
795 {
796 {0, 1}, // false
797 {1, 0} // true
798 };
799
800 static const int lookup_table_3d[2][2][2][4] =
801 // orientation flip rotation
802 {{{
803 {0, 2, 1, 3}, // false false false
804 {2, 3, 0, 1} // false false true
805 },
806 {
807 {3, 1, 2, 0}, // false true false
808 {1, 0, 3, 2} // false true true
809 }},
810 {{
811 {0, 1, 2, 3}, // true false false
812 {1, 3, 0, 2} // true false true
813 },
814 {
815 {3, 2, 1, 0}, // true true false
816 {2, 0, 3, 1} // true true true
817 }}};
818
819 if (cell_1->has_children())
820 {
821 if (cell_2->has_children())
822 {
823 // In the case that both faces have children, we loop over all
824 // children and apply update_periodic_face_map_recursively
825 // recursively:
826
827 Assert(face_1->n_children() ==
829 face_2->n_children() ==
832
833 for (unsigned int i = 0;
834 i < GeometryInfo<dim>::max_children_per_face;
835 ++i)
836 {
837 // Lookup the index for the second face
838 unsigned int j = 0;
839 switch (dim)
840 {
841 case 2:
842 j = lookup_table_2d[face_flip][i];
843 break;
844 case 3:
845 j = lookup_table_3d[face_orientation][face_flip]
846 [face_rotation][i];
847 break;
848 default:
850 }
851
852 // find subcell ids that belong to the subface indices
853 unsigned int child_cell_1 =
855 cell_1->refinement_case(),
856 n_face_1,
857 i,
858 cell_1->face_orientation(n_face_1),
859 cell_1->face_flip(n_face_1),
860 cell_1->face_rotation(n_face_1),
861 face_1->refinement_case());
862 unsigned int child_cell_2 =
864 cell_2->refinement_case(),
865 n_face_2,
866 j,
867 cell_2->face_orientation(n_face_2),
868 cell_2->face_flip(n_face_2),
869 cell_2->face_rotation(n_face_2),
870 face_2->refinement_case());
871
872 Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
873 face_1->child(i),
875 Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
876 face_2->child(j),
878
879 // precondition: subcell has the same orientation as cell
880 // (so that the face numbers coincide) recursive call
881 update_periodic_face_map_recursively<dim, spacedim>(
882 cell_1->child(child_cell_1),
883 cell_2->child(child_cell_2),
884 n_face_1,
885 n_face_2,
886 orientation,
887 periodic_face_map);
888 }
889 }
890 else // only face_1 has children
891 {
892 for (unsigned int i = 0;
893 i < GeometryInfo<dim>::max_children_per_face;
894 ++i)
895 {
896 // find subcell ids that belong to the subface indices
897 unsigned int child_cell_1 =
899 cell_1->refinement_case(),
900 n_face_1,
901 i,
902 cell_1->face_orientation(n_face_1),
903 cell_1->face_flip(n_face_1),
904 cell_1->face_rotation(n_face_1),
905 face_1->refinement_case());
906
907 // recursive call
908 update_periodic_face_map_recursively<dim, spacedim>(
909 cell_1->child(child_cell_1),
910 cell_2,
911 n_face_1,
912 n_face_2,
913 orientation,
914 periodic_face_map);
915 }
916 }
917 }
918 }
919 }
920
921
922} // end of anonymous namespace
923
924
925namespace internal
926{
927 namespace TriangulationImplementation
928 {
929 // make sure that if in the following we
930 // write Triangulation<dim,spacedim>
931 // we mean the *class*
932 // ::Triangulation, not the
933 // enclosing namespace
934 // internal::TriangulationImplementation
936
942 int,
943 << "Something went wrong when making cell " << arg1
944 << ". Read the docs and the source code "
945 << "for more information.");
951 int,
952 << "Something went wrong upon construction of cell "
953 << arg1);
964 int,
965 << "Cell " << arg1
966 << " has negative measure. This typically "
967 << "indicates some distortion in the cell, or a mistakenly "
968 << "swapped pair of vertices in the input to "
969 << "Triangulation::create_triangulation().");
978 int,
979 int,
980 int,
981 << "Error while creating cell " << arg1
982 << ": the vertex index " << arg2 << " must be between 0 and "
983 << arg3 << ".");
989 int,
990 int,
991 << "While trying to assign a boundary indicator to a line: "
992 << "the line with end vertices " << arg1 << " and " << arg2
993 << " does not exist.");
999 int,
1000 int,
1001 int,
1002 int,
1003 << "While trying to assign a boundary indicator to a quad: "
1004 << "the quad with bounding lines " << arg1 << ", " << arg2
1005 << ", " << arg3 << ", " << arg4 << " does not exist.");
1012 int,
1013 int,
1015 << "The input data for creating a triangulation contained "
1016 << "information about a line with indices " << arg1 << " and " << arg2
1017 << " that is described to have boundary indicator "
1018 << static_cast<int>(arg3)
1019 << ". However, this is an internal line not located on the "
1020 << "boundary. You cannot assign a boundary indicator to it." << std::endl
1021 << std::endl
1022 << "If this happened at a place where you call "
1023 << "Triangulation::create_triangulation() yourself, you need "
1024 << "to check the SubCellData object you pass to this function."
1025 << std::endl
1026 << std::endl
1027 << "If this happened in a place where you are reading a mesh "
1028 << "from a file, then you need to investigate why such a line "
1029 << "ended up in the input file. A typical case is a geometry "
1030 << "that consisted of multiple parts and for which the mesh "
1031 << "generator program assumes that the interface between "
1032 << "two parts is a boundary when that isn't supposed to be "
1033 << "the case, or where the mesh generator simply assigns "
1034 << "'geometry indicators' to lines at the perimeter of "
1035 << "a part that are not supposed to be interpreted as "
1036 << "'boundary indicators'.");
1043 int,
1044 int,
1045 int,
1046 int,
1048 << "The input data for creating a triangulation contained "
1049 << "information about a quad with indices " << arg1 << ", " << arg2
1050 << ", " << arg3 << ", and " << arg4
1051 << " that is described to have boundary indicator "
1052 << static_cast<int>(arg5)
1053 << ". However, this is an internal quad not located on the "
1054 << "boundary. You cannot assign a boundary indicator to it." << std::endl
1055 << std::endl
1056 << "If this happened at a place where you call "
1057 << "Triangulation::create_triangulation() yourself, you need "
1058 << "to check the SubCellData object you pass to this function."
1059 << std::endl
1060 << std::endl
1061 << "If this happened in a place where you are reading a mesh "
1062 << "from a file, then you need to investigate why such a quad "
1063 << "ended up in the input file. A typical case is a geometry "
1064 << "that consisted of multiple parts and for which the mesh "
1065 << "generator program assumes that the interface between "
1066 << "two parts is a boundary when that isn't supposed to be "
1067 << "the case, or where the mesh generator simply assigns "
1068 << "'geometry indicators' to quads at the surface of "
1069 << "a part that are not supposed to be interpreted as "
1070 << "'boundary indicators'.");
1077 int,
1078 int,
1079 << "In SubCellData the line info of the line with vertex indices " << arg1
1080 << " and " << arg2 << " appears more than once. "
1081 << "This is not allowed.");
1088 int,
1089 int,
1090 std::string,
1091 << "In SubCellData the line info of the line with vertex indices " << arg1
1092 << " and " << arg2 << " appears multiple times with different (valid) "
1093 << arg3 << ". This is not allowed.");
1100 int,
1101 int,
1102 int,
1103 int,
1104 std::string,
1105 << "In SubCellData the quad info of the quad with line indices " << arg1
1106 << ", " << arg2 << ", " << arg3 << " and " << arg4
1107 << " appears multiple times with different (valid) " << arg5
1108 << ". This is not allowed.");
1109
1110 /*
1111 * Reserve space for TriaFaces. Details:
1112 *
1113 * Reserve space for line_orientations.
1114 *
1115 * @note Used only for dim=3.
1116 */
1117 void
1119 const unsigned int new_quads_in_pairs,
1120 const unsigned int new_quads_single)
1121 {
1122 AssertDimension(tria_faces.dim, 3);
1123
1124 Assert(new_quads_in_pairs % 2 == 0, ExcInternalError());
1125
1126 unsigned int next_free_single = 0;
1127 unsigned int next_free_pair = 0;
1128
1129 // count the number of objects, of unused single objects and of
1130 // unused pairs of objects
1131 unsigned int n_quads = 0;
1132 unsigned int n_unused_pairs = 0;
1133 unsigned int n_unused_singles = 0;
1134 for (unsigned int i = 0; i < tria_faces.quads.used.size(); ++i)
1135 {
1136 if (tria_faces.quads.used[i])
1137 ++n_quads;
1138 else if (i + 1 < tria_faces.quads.used.size())
1139 {
1140 if (tria_faces.quads.used[i + 1])
1141 {
1142 ++n_unused_singles;
1143 if (next_free_single == 0)
1144 next_free_single = i;
1145 }
1146 else
1147 {
1148 ++n_unused_pairs;
1149 if (next_free_pair == 0)
1150 next_free_pair = i;
1151 ++i;
1152 }
1153 }
1154 else
1155 ++n_unused_singles;
1156 }
1157 Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
1158 tria_faces.quads.used.size(),
1160
1161 // how many single quads are needed in addition to n_unused_quads?
1162 const int additional_single_quads = new_quads_single - n_unused_singles;
1163
1164 unsigned int new_size =
1165 tria_faces.quads.used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
1166 if (additional_single_quads > 0)
1167 new_size += additional_single_quads;
1168
1169 // see above...
1170 if (new_size > tria_faces.quads.n_objects())
1171 {
1172 // reserve the field of the derived class
1173 tria_faces.quads_line_orientations.reserve(
1175 tria_faces.quads_line_orientations.insert(
1176 tria_faces.quads_line_orientations.end(),
1178 tria_faces.quads_line_orientations.size(),
1179 true);
1180
1181 tria_faces.quad_reference_cell.reserve(new_size);
1182 tria_faces.quad_reference_cell.insert(
1183 tria_faces.quad_reference_cell.end(),
1184 new_size - tria_faces.quad_reference_cell.size(),
1186 }
1187 }
1188
1189
1190
1204 void
1206 const unsigned int total_cells,
1207 const unsigned int dimension,
1208 const unsigned int space_dimension)
1209 {
1210 // we need space for total_cells cells. Maybe we have more already
1211 // with those cells which are unused, so only allocate new space if
1212 // needed.
1213 //
1214 // note that all arrays should have equal sizes (checked by
1215 // @p{monitor_memory}
1216 if (total_cells > tria_level.refine_flags.size())
1217 {
1218 tria_level.refine_flags.reserve(total_cells);
1219 tria_level.refine_flags.insert(tria_level.refine_flags.end(),
1220 total_cells -
1221 tria_level.refine_flags.size(),
1222 /*RefinementCase::no_refinement=*/0);
1223
1224 tria_level.coarsen_flags.reserve(total_cells);
1225 tria_level.coarsen_flags.insert(tria_level.coarsen_flags.end(),
1226 total_cells -
1227 tria_level.coarsen_flags.size(),
1228 false);
1229
1230 tria_level.active_cell_indices.reserve(total_cells);
1231 tria_level.active_cell_indices.insert(
1232 tria_level.active_cell_indices.end(),
1233 total_cells - tria_level.active_cell_indices.size(),
1235
1236 tria_level.subdomain_ids.reserve(total_cells);
1237 tria_level.subdomain_ids.insert(tria_level.subdomain_ids.end(),
1238 total_cells -
1239 tria_level.subdomain_ids.size(),
1240 0);
1241
1242 tria_level.level_subdomain_ids.reserve(total_cells);
1243 tria_level.level_subdomain_ids.insert(
1244 tria_level.level_subdomain_ids.end(),
1245 total_cells - tria_level.level_subdomain_ids.size(),
1246 0);
1247
1248 tria_level.global_active_cell_indices.reserve(total_cells);
1249 tria_level.global_active_cell_indices.insert(
1250 tria_level.global_active_cell_indices.end(),
1251 total_cells - tria_level.global_active_cell_indices.size(),
1253
1254 tria_level.global_level_cell_indices.reserve(total_cells);
1255 tria_level.global_level_cell_indices.insert(
1256 tria_level.global_level_cell_indices.end(),
1257 total_cells - tria_level.global_level_cell_indices.size(),
1259
1260 if (dimension < space_dimension)
1261 {
1262 tria_level.direction_flags.reserve(total_cells);
1263 tria_level.direction_flags.insert(
1264 tria_level.direction_flags.end(),
1265 total_cells - tria_level.direction_flags.size(),
1266 true);
1267 }
1268 else
1269 tria_level.direction_flags.clear();
1270
1271 tria_level.parents.reserve((total_cells + 1) / 2);
1272 tria_level.parents.insert(tria_level.parents.end(),
1273 (total_cells + 1) / 2 -
1274 tria_level.parents.size(),
1275 -1);
1276
1277 tria_level.neighbors.reserve(total_cells * (2 * dimension));
1278 tria_level.neighbors.insert(tria_level.neighbors.end(),
1279 total_cells * (2 * dimension) -
1280 tria_level.neighbors.size(),
1281 std::make_pair(-1, -1));
1282
1283 if (tria_level.dim == 2 || tria_level.dim == 3)
1284 {
1285 const unsigned int max_faces_per_cell = 2 * dimension;
1286 tria_level.face_orientations.reserve(total_cells *
1287 max_faces_per_cell);
1288 tria_level.face_orientations.insert(
1289 tria_level.face_orientations.end(),
1290 total_cells * max_faces_per_cell -
1291 tria_level.face_orientations.size(),
1292 true);
1293
1294 tria_level.reference_cell.reserve(total_cells);
1295 tria_level.reference_cell.insert(
1296 tria_level.reference_cell.end(),
1297 total_cells - tria_level.reference_cell.size(),
1298 tria_level.dim == 2 ? ::ReferenceCells::Quadrilateral :
1300 }
1301 }
1302 }
1303
1304
1305
1310 int,
1311 int,
1312 << "The containers have sizes " << arg1 << " and " << arg2
1313 << ", which is not as expected.");
1314
1320 void
1321 monitor_memory(const TriaLevel & tria_level,
1322 const unsigned int true_dimension)
1323 {
1324 (void)tria_level;
1325 (void)true_dimension;
1326 Assert(2 * true_dimension * tria_level.refine_flags.size() ==
1327 tria_level.neighbors.size(),
1328 ExcMemoryInexact(tria_level.refine_flags.size(),
1329 tria_level.neighbors.size()));
1330 Assert(2 * true_dimension * tria_level.coarsen_flags.size() ==
1331 tria_level.neighbors.size(),
1332 ExcMemoryInexact(tria_level.coarsen_flags.size(),
1333 tria_level.neighbors.size()));
1334 }
1335
1336
1337
1350 void
1352 const unsigned int new_objects_in_pairs,
1353 const unsigned int new_objects_single = 0)
1354 {
1355 if (tria_objects.structdim <= 2)
1356 {
1357 Assert(new_objects_in_pairs % 2 == 0, ExcInternalError());
1358
1359 tria_objects.next_free_single = 0;
1360 tria_objects.next_free_pair = 0;
1361 tria_objects.reverse_order_next_free_single = false;
1362
1363 // count the number of objects, of unused single objects and of
1364 // unused pairs of objects
1365 unsigned int n_objects = 0;
1366 unsigned int n_unused_pairs = 0;
1367 unsigned int n_unused_singles = 0;
1368 for (unsigned int i = 0; i < tria_objects.used.size(); ++i)
1369 {
1370 if (tria_objects.used[i])
1371 ++n_objects;
1372 else if (i + 1 < tria_objects.used.size())
1373 {
1374 if (tria_objects.used[i + 1])
1375 {
1376 ++n_unused_singles;
1377 if (tria_objects.next_free_single == 0)
1378 tria_objects.next_free_single = i;
1379 }
1380 else
1381 {
1382 ++n_unused_pairs;
1383 if (tria_objects.next_free_pair == 0)
1384 tria_objects.next_free_pair = i;
1385 ++i;
1386 }
1387 }
1388 else
1389 ++n_unused_singles;
1390 }
1391 Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
1392 tria_objects.used.size(),
1394
1395 // how many single objects are needed in addition to
1396 // n_unused_objects?
1397 const int additional_single_objects =
1398 new_objects_single - n_unused_singles;
1399
1400 unsigned int new_size = tria_objects.used.size() +
1401 new_objects_in_pairs - 2 * n_unused_pairs;
1402 if (additional_single_objects > 0)
1403 new_size += additional_single_objects;
1404
1405 // only allocate space if necessary
1406 if (new_size > tria_objects.n_objects())
1407 {
1408 const unsigned int max_faces_per_cell =
1409 2 * tria_objects.structdim;
1410 const unsigned int max_children_per_cell =
1411 1 << tria_objects.structdim;
1412
1413 tria_objects.cells.reserve(new_size * max_faces_per_cell);
1414 tria_objects.cells.insert(tria_objects.cells.end(),
1415 (new_size - tria_objects.n_objects()) *
1416 max_faces_per_cell,
1417 -1);
1418
1419 tria_objects.used.reserve(new_size);
1420 tria_objects.used.insert(tria_objects.used.end(),
1421 new_size - tria_objects.used.size(),
1422 false);
1423
1424 tria_objects.user_flags.reserve(new_size);
1425 tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1426 new_size -
1427 tria_objects.user_flags.size(),
1428 false);
1429
1430 const unsigned int factor = max_children_per_cell / 2;
1431 tria_objects.children.reserve(factor * new_size);
1432 tria_objects.children.insert(tria_objects.children.end(),
1433 factor * new_size -
1434 tria_objects.children.size(),
1435 -1);
1436
1437 if (tria_objects.structdim > 1)
1438 {
1439 tria_objects.refinement_cases.reserve(new_size);
1440 tria_objects.refinement_cases.insert(
1441 tria_objects.refinement_cases.end(),
1442 new_size - tria_objects.refinement_cases.size(),
1443 /*RefinementCase::no_refinement=*/0);
1444 }
1445
1446 // first reserve, then resize. Otherwise the std library can
1447 // decide to allocate more entries.
1448 tria_objects.boundary_or_material_id.reserve(new_size);
1449 tria_objects.boundary_or_material_id.resize(new_size);
1450
1451 tria_objects.user_data.reserve(new_size);
1452 tria_objects.user_data.resize(new_size);
1453
1454 tria_objects.manifold_id.reserve(new_size);
1455 tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1456 new_size -
1457 tria_objects.manifold_id.size(),
1459 }
1460
1461 if (n_unused_singles == 0)
1462 {
1463 tria_objects.next_free_single = new_size - 1;
1464 tria_objects.reverse_order_next_free_single = true;
1465 }
1466 }
1467 else
1468 {
1469 const unsigned int new_hexes = new_objects_in_pairs;
1470
1471 const unsigned int new_size =
1472 new_hexes + std::count(tria_objects.used.begin(),
1473 tria_objects.used.end(),
1474 true);
1475
1476 // see above...
1477 if (new_size > tria_objects.n_objects())
1478 {
1479 const unsigned int max_faces_per_cell =
1480 2 * tria_objects.structdim;
1481
1482 tria_objects.cells.reserve(new_size * max_faces_per_cell);
1483 tria_objects.cells.insert(tria_objects.cells.end(),
1484 (new_size - tria_objects.n_objects()) *
1485 max_faces_per_cell,
1486 -1);
1487
1488 tria_objects.used.reserve(new_size);
1489 tria_objects.used.insert(tria_objects.used.end(),
1490 new_size - tria_objects.used.size(),
1491 false);
1492
1493 tria_objects.user_flags.reserve(new_size);
1494 tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1495 new_size -
1496 tria_objects.user_flags.size(),
1497 false);
1498
1499 tria_objects.children.reserve(4 * new_size);
1500 tria_objects.children.insert(tria_objects.children.end(),
1501 4 * new_size -
1502 tria_objects.children.size(),
1503 -1);
1504
1505 // for the following fields, we know exactly how many elements
1506 // we need, so first reserve then resize (resize itself, at least
1507 // with some compiler libraries, appears to round up the size it
1508 // actually reserves)
1509 tria_objects.boundary_or_material_id.reserve(new_size);
1510 tria_objects.boundary_or_material_id.resize(new_size);
1511
1512 tria_objects.manifold_id.reserve(new_size);
1513 tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1514 new_size -
1515 tria_objects.manifold_id.size(),
1517
1518 tria_objects.user_data.reserve(new_size);
1519 tria_objects.user_data.resize(new_size);
1520
1521 tria_objects.refinement_cases.reserve(new_size);
1522 tria_objects.refinement_cases.insert(
1523 tria_objects.refinement_cases.end(),
1524 new_size - tria_objects.refinement_cases.size(),
1525 /*RefinementCase::no_refinement=*/0);
1526 }
1527 tria_objects.next_free_single = tria_objects.next_free_pair = 0;
1528 }
1529 }
1530
1531
1532
1538 void
1539 monitor_memory(const TriaObjects &tria_object, const unsigned int)
1540 {
1541 Assert(tria_object.n_objects() == tria_object.used.size(),
1542 ExcMemoryInexact(tria_object.n_objects(),
1543 tria_object.used.size()));
1544 Assert(tria_object.n_objects() == tria_object.user_flags.size(),
1545 ExcMemoryInexact(tria_object.n_objects(),
1546 tria_object.user_flags.size()));
1547 Assert(tria_object.n_objects() ==
1548 tria_object.boundary_or_material_id.size(),
1549 ExcMemoryInexact(tria_object.n_objects(),
1550 tria_object.boundary_or_material_id.size()));
1551 Assert(tria_object.n_objects() == tria_object.manifold_id.size(),
1552 ExcMemoryInexact(tria_object.n_objects(),
1553 tria_object.manifold_id.size()));
1554 Assert(tria_object.n_objects() == tria_object.user_data.size(),
1555 ExcMemoryInexact(tria_object.n_objects(),
1556 tria_object.user_data.size()));
1557
1558 if (tria_object.structdim == 1)
1559 {
1560 Assert(1 * tria_object.n_objects() == tria_object.children.size(),
1561 ExcMemoryInexact(tria_object.n_objects(),
1562 tria_object.children.size()));
1563 }
1564 else if (tria_object.structdim == 2)
1565 {
1566 Assert(2 * tria_object.n_objects() == tria_object.children.size(),
1567 ExcMemoryInexact(tria_object.n_objects(),
1568 tria_object.children.size()));
1569 }
1570 else if (tria_object.structdim == 3)
1571 {
1572 Assert(4 * tria_object.n_objects() == tria_object.children.size(),
1573 ExcMemoryInexact(tria_object.n_objects(),
1574 tria_object.children.size()));
1575 }
1576 }
1577
1578
1579
1584 template <int dim, int spacedim>
1586 {
1587 public:
1591 virtual ~Policy() = default;
1592
1596 virtual void
1598
1602 virtual void
1606 std::vector<unsigned int> & line_cell_count,
1607 std::vector<unsigned int> &quad_cell_count) = 0;
1608
1614 const bool check_for_distorted_cells) = 0;
1615
1619 virtual void
1622
1626 virtual void
1629
1633 virtual bool
1635 const typename Triangulation<dim, spacedim>::cell_iterator &cell) = 0;
1636
1643 virtual std::unique_ptr<Policy<dim, spacedim>>
1644 clone() = 0;
1645 };
1646
1647
1648
1654 template <int dim, int spacedim, typename T>
1655 class PolicyWrapper : public Policy<dim, spacedim>
1656 {
1657 public:
1658 void
1660 {
1661 T::update_neighbors(tria);
1662 }
1663
1664 void
1668 std::vector<unsigned int> & line_cell_count,
1669 std::vector<unsigned int> &quad_cell_count) override
1670 {
1671 T::delete_children(tria, cell, line_cell_count, quad_cell_count);
1672 }
1673
1676 const bool check_for_distorted_cells) override
1677 {
1678 return T::execute_refinement(triangulation, check_for_distorted_cells);
1679 }
1680
1681 void
1684 {
1685 T::prevent_distorted_boundary_cells(triangulation);
1686 }
1687
1688 void
1691 {
1692 T::prepare_refinement_dim_dependent(triangulation);
1693 }
1694
1695 bool
1698 override
1699 {
1700 return T::template coarsening_allowed<dim, spacedim>(cell);
1701 }
1702
1703 std::unique_ptr<Policy<dim, spacedim>>
1704 clone() override
1705 {
1706 return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
1707 }
1708 };
1709
1710
1711
1808 {
1820 template <int dim, int spacedim>
1821 static void
1824 const unsigned int level_objects,
1826 {
1827 using line_iterator =
1829
1830 number_cache.n_levels = 0;
1831 if (level_objects > 0)
1832 // find the last level on which there are used cells
1833 for (unsigned int level = 0; level < level_objects; ++level)
1834 if (triangulation.begin(level) != triangulation.end(level))
1835 number_cache.n_levels = level + 1;
1836
1837 // no cells at all?
1838 Assert(number_cache.n_levels > 0, ExcInternalError());
1839
1841 // update the number of lines on the different levels in the
1842 // cache
1843 number_cache.n_lines = 0;
1844 number_cache.n_active_lines = 0;
1845
1846 // for 1d, lines have levels so take count the objects per
1847 // level and globally
1848 if (dim == 1)
1849 {
1850 number_cache.n_lines_level.resize(number_cache.n_levels);
1851 number_cache.n_active_lines_level.resize(number_cache.n_levels);
1852
1853 for (unsigned int level = 0; level < number_cache.n_levels; ++level)
1854 {
1855 // count lines on this level
1856 number_cache.n_lines_level[level] = 0;
1857 number_cache.n_active_lines_level[level] = 0;
1858
1859 line_iterator line = triangulation.begin_line(level),
1860 endc =
1861 (level == number_cache.n_levels - 1 ?
1862 line_iterator(triangulation.end_line()) :
1863 triangulation.begin_line(level + 1));
1864 for (; line != endc; ++line)
1865 {
1866 ++number_cache.n_lines_level[level];
1867 if (line->has_children() == false)
1868 ++number_cache.n_active_lines_level[level];
1869 }
1870
1871 // update total number of lines
1872 number_cache.n_lines += number_cache.n_lines_level[level];
1873 number_cache.n_active_lines +=
1874 number_cache.n_active_lines_level[level];
1875 }
1876 }
1877 else
1878 {
1879 // for dim>1, there are no levels for lines
1880 number_cache.n_lines_level.clear();
1881 number_cache.n_active_lines_level.clear();
1882
1883 line_iterator line = triangulation.begin_line(),
1884 endc = triangulation.end_line();
1885 for (; line != endc; ++line)
1886 {
1887 ++number_cache.n_lines;
1888 if (line->has_children() == false)
1889 ++number_cache.n_active_lines;
1890 }
1891 }
1892 }
1893
1908 template <int dim, int spacedim>
1909 static void
1912 const unsigned int level_objects,
1914 {
1915 // update lines and n_levels in number_cache. since we don't
1916 // access any of these numbers, we can do this in the
1917 // background
1919 static_cast<
1920 void (*)(const Triangulation<dim, spacedim> &,
1921 const unsigned int,
1923 &compute_number_cache<dim, spacedim>),
1925 level_objects,
1927 number_cache));
1928
1929 using quad_iterator =
1931
1933 // update the number of quads on the different levels in the
1934 // cache
1935 number_cache.n_quads = 0;
1936 number_cache.n_active_quads = 0;
1937
1938 // for 2d, quads have levels so take count the objects per
1939 // level and globally
1940 if (dim == 2)
1941 {
1942 // count the number of levels; the function we called above
1943 // on a separate Task for lines also does this and puts it into
1944 // number_cache.n_levels, but this datum may not yet be
1945 // available as we call the function on a separate task
1946 unsigned int n_levels = 0;
1947 if (level_objects > 0)
1948 // find the last level on which there are used cells
1949 for (unsigned int level = 0; level < level_objects; ++level)
1950 if (triangulation.begin(level) != triangulation.end(level))
1951 n_levels = level + 1;
1952
1953 number_cache.n_quads_level.resize(n_levels);
1954 number_cache.n_active_quads_level.resize(n_levels);
1955
1956 for (unsigned int level = 0; level < n_levels; ++level)
1957 {
1958 // count quads on this level
1959 number_cache.n_quads_level[level] = 0;
1960 number_cache.n_active_quads_level[level] = 0;
1961
1962 quad_iterator quad = triangulation.begin_quad(level),
1963 endc =
1964 (level == n_levels - 1 ?
1965 quad_iterator(triangulation.end_quad()) :
1966 triangulation.begin_quad(level + 1));
1967 for (; quad != endc; ++quad)
1968 {
1969 ++number_cache.n_quads_level[level];
1970 if (quad->has_children() == false)
1971 ++number_cache.n_active_quads_level[level];
1972 }
1973
1974 // update total number of quads
1975 number_cache.n_quads += number_cache.n_quads_level[level];
1976 number_cache.n_active_quads +=
1977 number_cache.n_active_quads_level[level];
1978 }
1979 }
1980 else
1981 {
1982 // for dim>2, there are no levels for quads
1983 number_cache.n_quads_level.clear();
1984 number_cache.n_active_quads_level.clear();
1985
1986 quad_iterator quad = triangulation.begin_quad(),
1987 endc = triangulation.end_quad();
1988 for (; quad != endc; ++quad)
1989 {
1990 ++number_cache.n_quads;
1991 if (quad->has_children() == false)
1992 ++number_cache.n_active_quads;
1993 }
1994 }
1995
1996 // wait for the background computation for lines
1997 update_lines.join();
1998 }
1999
2015 template <int dim, int spacedim>
2016 static void
2019 const unsigned int level_objects,
2021 {
2022 // update quads, lines and n_levels in number_cache. since we
2023 // don't access any of these numbers, we can do this in the
2024 // background
2025 Threads::Task<void> update_quads_and_lines = Threads::new_task(
2026 static_cast<
2027 void (*)(const Triangulation<dim, spacedim> &,
2028 const unsigned int,
2030 &compute_number_cache<dim, spacedim>),
2032 level_objects,
2034 number_cache));
2035
2036 using hex_iterator =
2038
2040 // update the number of hexes on the different levels in the
2041 // cache
2042 number_cache.n_hexes = 0;
2043 number_cache.n_active_hexes = 0;
2044
2045 // for 3d, hexes have levels so take count the objects per
2046 // level and globally
2047 if (dim == 3)
2048 {
2049 // count the number of levels; the function we called
2050 // above on a separate Task for quads (recursively, via
2051 // the lines function) also does this and puts it into
2052 // number_cache.n_levels, but this datum may not yet be
2053 // available as we call the function on a separate task
2054 unsigned int n_levels = 0;
2055 if (level_objects > 0)
2056 // find the last level on which there are used cells
2057 for (unsigned int level = 0; level < level_objects; ++level)
2058 if (triangulation.begin(level) != triangulation.end(level))
2059 n_levels = level + 1;
2060
2061 number_cache.n_hexes_level.resize(n_levels);
2062 number_cache.n_active_hexes_level.resize(n_levels);
2063
2064 for (unsigned int level = 0; level < n_levels; ++level)
2065 {
2066 // count hexes on this level
2067 number_cache.n_hexes_level[level] = 0;
2068 number_cache.n_active_hexes_level[level] = 0;
2069
2070 hex_iterator hex = triangulation.begin_hex(level),
2071 endc = (level == n_levels - 1 ?
2072 hex_iterator(triangulation.end_hex()) :
2073 triangulation.begin_hex(level + 1));
2074 for (; hex != endc; ++hex)
2075 {
2076 ++number_cache.n_hexes_level[level];
2077 if (hex->has_children() == false)
2078 ++number_cache.n_active_hexes_level[level];
2079 }
2080
2081 // update total number of hexes
2082 number_cache.n_hexes += number_cache.n_hexes_level[level];
2083 number_cache.n_active_hexes +=
2084 number_cache.n_active_hexes_level[level];
2085 }
2086 }
2087 else
2088 {
2089 // for dim>3, there are no levels for hexes
2090 number_cache.n_hexes_level.clear();
2091 number_cache.n_active_hexes_level.clear();
2092
2093 hex_iterator hex = triangulation.begin_hex(),
2094 endc = triangulation.end_hex();
2095 for (; hex != endc; ++hex)
2096 {
2097 ++number_cache.n_hexes;
2098 if (hex->has_children() == false)
2099 ++number_cache.n_active_hexes;
2100 }
2101 }
2102
2103 // wait for the background computation for quads
2104 update_quads_and_lines.join();
2105 }
2106
2107
2108
2109 template <int spacedim>
2111 {}
2112
2113
2114 template <int dim, int spacedim>
2115 static void
2117 {
2118 // each face can be neighbored on two sides
2119 // by cells. according to the face's
2120 // intrinsic normal we define the left
2121 // neighbor as the one for which the face
2122 // normal points outward, and store that
2123 // one first; the second one is then
2124 // the right neighbor for which the
2125 // face normal points inward. This
2126 // information depends on the type of cell
2127 // and local number of face for the
2128 // 'standard ordering and orientation' of
2129 // faces and then on the face_orientation
2130 // information for the real mesh. Set up a
2131 // table to have fast access to those
2132 // offsets (0 for left and 1 for
2133 // right). Some of the values are invalid
2134 // as they reference too large face
2135 // numbers, but we just leave them at a
2136 // zero value.
2137 //
2138 // Note, that in 2d for lines as faces the
2139 // normal direction given in the
2140 // GeometryInfo class is not consistent. We
2141 // thus define here that the normal for a
2142 // line points to the right if the line
2143 // points upwards.
2144 //
2145 // There is one more point to
2146 // consider, however: if we have
2147 // dim<spacedim, then we may have
2148 // cases where cells are
2149 // inverted. In effect, both
2150 // cells think they are the left
2151 // neighbor of an edge, for
2152 // example, which leads us to
2153 // forget neighborship
2154 // information (a case that shows
2155 // this is
2156 // codim_one/hanging_nodes_02). We
2157 // store whether a cell is
2158 // inverted using the
2159 // direction_flag, so if a cell
2160 // has a false direction_flag,
2161 // then we need to invert our
2162 // selection whether we are a
2163 // left or right neighbor in all
2164 // following computations.
2165 //
2166 // first index: dimension (minus 2)
2167 // second index: local face index
2168 // third index: face_orientation (false and true)
2169 static const unsigned int left_right_offset[2][6][2] = {
2170 // quadrilateral
2171 {{0, 1}, // face 0, face_orientation = false and true
2172 {1, 0}, // face 1, face_orientation = false and true
2173 {1, 0}, // face 2, face_orientation = false and true
2174 {0, 1}, // face 3, face_orientation = false and true
2175 {0, 0}, // face 4, invalid face
2176 {0, 0}}, // face 5, invalid face
2177 // hexahedron
2178 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
2179
2180 // now create a vector of the two active
2181 // neighbors (left and right) for each face
2182 // and fill it by looping over all cells. For
2183 // cases with anisotropic refinement and more
2184 // then one cell neighboring at a given side
2185 // of the face we will automatically get the
2186 // active one on the highest level as we loop
2187 // over cells from lower levels first.
2189 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
2190 adjacent_cells(2 * triangulation.n_raw_faces(), dummy);
2191
2192 for (const auto &cell : triangulation.cell_iterators())
2193 for (auto f : cell->face_indices())
2194 {
2196 cell->face(f);
2197
2198 const unsigned int offset =
2199 (cell->direction_flag() ?
2200 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2201 1 -
2202 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2203
2204 adjacent_cells[2 * face->index() + offset] = cell;
2205
2206 // if this cell is not refined, but the
2207 // face is, then we'll have to set our
2208 // cell as neighbor for the child faces
2209 // as well. Fortunately the normal
2210 // orientation of children will be just
2211 // the same.
2212 if (dim == 2)
2213 {
2214 if (cell->is_active() && face->has_children())
2215 {
2216 adjacent_cells[2 * face->child(0)->index() + offset] =
2217 cell;
2218 adjacent_cells[2 * face->child(1)->index() + offset] =
2219 cell;
2220 }
2221 }
2222 else // -> dim == 3
2223 {
2224 // We need the same as in 2d
2225 // here. Furthermore, if the face is
2226 // refined with cut_x or cut_y then
2227 // those children again in the other
2228 // direction, and if this cell is
2229 // refined isotropically (along the
2230 // face) then the neighbor will
2231 // (probably) be refined as cut_x or
2232 // cut_y along the face. For those
2233 // neighboring children cells, their
2234 // neighbor will be the current,
2235 // inactive cell, as our children are
2236 // too fine to be neighbors. Catch that
2237 // case by also acting on inactive
2238 // cells with isotropic refinement
2239 // along the face. If the situation
2240 // described is not present, the data
2241 // will be overwritten later on when we
2242 // visit cells on finer levels, so no
2243 // harm will be done.
2244 if (face->has_children() &&
2245 (cell->is_active() ||
2247 cell->refinement_case(), f) ==
2248 RefinementCase<dim - 1>::isotropic_refinement))
2249 {
2250 for (unsigned int c = 0; c < face->n_children(); ++c)
2251 adjacent_cells[2 * face->child(c)->index() + offset] =
2252 cell;
2253 if (face->child(0)->has_children())
2254 {
2255 adjacent_cells[2 * face->child(0)->child(0)->index() +
2256 offset] = cell;
2257 adjacent_cells[2 * face->child(0)->child(1)->index() +
2258 offset] = cell;
2259 }
2260 if (face->child(1)->has_children())
2261 {
2262 adjacent_cells[2 * face->child(1)->child(0)->index() +
2263 offset] = cell;
2264 adjacent_cells[2 * face->child(1)->child(1)->index() +
2265 offset] = cell;
2266 }
2267 } // if cell active and face refined
2268 } // else -> dim==3
2269 } // for all faces of all cells
2270
2271 // now loop again over all cells and set the
2272 // corresponding neighbor cell. Note, that we
2273 // have to use the opposite of the
2274 // left_right_offset in this case as we want
2275 // the offset of the neighbor, not our own.
2276 for (const auto &cell : triangulation.cell_iterators())
2277 for (auto f : cell->face_indices())
2278 {
2279 const unsigned int offset =
2280 (cell->direction_flag() ?
2281 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2282 1 -
2283 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2284 cell->set_neighbor(
2285 f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
2286 }
2287 }
2288
2289
2293 template <int dim, int spacedim>
2294 static void
2296 const std::vector<CellData<dim>> & cells,
2297 const SubCellData & subcelldata,
2299 {
2300 AssertThrow(vertices.size() > 0, ExcMessage("No vertices given"));
2301 AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
2302
2303 // Check that all cells have positive volume.
2304#ifndef _MSC_VER
2305 // TODO: The following code does not compile with MSVC. Find a way
2306 // around it
2307 if (dim == spacedim)
2308 for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2309 {
2310 // If we should check for distorted cells, then we permit them
2311 // to exist. If a cell has negative measure, then it must be
2312 // distorted (the converse is not necessarily true); hence
2313 // throw an exception if no such cells should exist.
2315 {
2316 const double cell_measure = GridTools::cell_measure<spacedim>(
2317 vertices,
2318 ArrayView<const unsigned int>(cells[cell_no].vertices));
2320 }
2321 }
2322#endif
2323
2324 // clear old content
2325 tria.levels.clear();
2326 tria.levels.push_back(
2327 std::make_unique<
2329
2330 if (dim > 1)
2331 tria.faces = std::make_unique<
2333
2334 // copy vertices
2335 tria.vertices = vertices;
2336 tria.vertices_used.assign(vertices.size(), true);
2337
2338 // compute connectivity
2339 const auto connectivity = build_connectivity<unsigned int>(cells);
2340 const unsigned int n_cell = cells.size();
2341
2342 // TriaObjects: lines
2343 if (dim >= 2)
2344 {
2345 auto &lines_0 = tria.faces->lines; // data structure to be filled
2346
2347 // get connectivity between quads and lines
2348 const auto & crs = connectivity.entity_to_entities(1, 0);
2349 const unsigned int n_lines = crs.ptr.size() - 1;
2350
2351 // allocate memory
2352 reserve_space_(lines_0, n_lines);
2353
2354 // loop over lines
2355 for (unsigned int line = 0; line < n_lines; ++line)
2356 for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
2357 ++i, ++j)
2358 lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
2359 crs.col[i]; // set vertex indices
2360 }
2361
2362 // TriaObjects: quads
2363 if (dim == 3)
2364 {
2365 auto &quads_0 = tria.faces->quads; // data structures to be filled
2366 auto &faces = *tria.faces;
2367
2368 // get connectivity between quads and lines
2369 const auto & crs = connectivity.entity_to_entities(2, 1);
2370 const unsigned int n_quads = crs.ptr.size() - 1;
2371
2372 // allocate memory
2373 reserve_space_(quads_0, n_quads);
2374 reserve_space_(faces, 2 /*structdim*/, n_quads);
2375
2376 // loop over all quads -> entity type, line indices/orientations
2377 for (unsigned int q = 0, k = 0; q < n_quads; ++q)
2378 {
2379 // set entity type of quads
2380 faces.quad_reference_cell[q] = connectivity.entity_types(2)[q];
2381
2382 // loop over all its lines
2383 for (unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
2384 ++i, ++j, ++k)
2385 {
2386 // set line index
2387 quads_0.cells[q * GeometryInfo<2>::faces_per_cell + j] =
2388 crs.col[i];
2389
2390 // set line orientations
2391 faces.quads_line_orientations
2393 connectivity.entity_orientations(1)[k];
2394 }
2395 }
2396 }
2397
2398 // TriaObjects/TriaLevel: cell
2399 {
2400 auto &cells_0 = tria.levels[0]->cells; // data structure to be filled
2401 auto &level = *tria.levels[0];
2402
2403 // get connectivity between cells/faces and cells/cells
2404 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2405 const auto &nei = connectivity.entity_to_entities(dim, dim);
2406
2407 // in 2D optional: since in in pure QUAD meshes same line
2408 // orientations can be guaranteed
2409 const bool orientation_needed =
2410 dim == 3 ||
2411 (dim == 2 &&
2412 std::any_of(connectivity.entity_orientations(1).begin(),
2413 connectivity.entity_orientations(1).end(),
2414 [](const auto &i) { return i == 0; }));
2415
2416 // allocate memory
2417 reserve_space_(cells_0, n_cell);
2418 reserve_space_(level, spacedim, n_cell, orientation_needed);
2419
2420 // loop over all cells
2421 for (unsigned int cell = 0; cell < n_cell; ++cell)
2422 {
2423 // set material ids
2424 cells_0.boundary_or_material_id[cell].material_id =
2425 cells[cell].material_id;
2426
2427 // set manifold ids
2428 cells_0.manifold_id[cell] = cells[cell].manifold_id;
2429
2430 // set entity types
2431 level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
2432
2433 // loop over faces
2434 for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2435 ++i, ++j)
2436 {
2437 // set neighbor if not at boundary
2438 if (nei.col[i] != static_cast<unsigned int>(-1))
2439 level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
2440 j] = {0, nei.col[i]};
2441
2442 // set face indices
2443 cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
2444 crs.col[i];
2445
2446 // set face orientation if needed
2447 if (orientation_needed)
2448 level.face_orientations
2450 connectivity.entity_orientations(dim - 1)[i];
2451 }
2452 }
2453 }
2454
2455 // TriaFaces: boundary id of boundary faces
2456 if (dim > 1)
2457 {
2458 auto &bids_face = dim == 3 ?
2459 tria.faces->quads.boundary_or_material_id :
2460 tria.faces->lines.boundary_or_material_id;
2461
2462 // count number of cells a face is belonging to
2463 std::vector<unsigned int> count(bids_face.size(), 0);
2464
2465 // get connectivity between cells/faces
2466 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2467
2468 // count how many cells are adjacent to the same face
2469 for (unsigned int cell = 0; cell < cells.size(); ++cell)
2470 for (unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
2471 count[crs.col[i]]++;
2472
2473 // loop over all faces
2474 for (unsigned int face = 0; face < count.size(); ++face)
2475 {
2476 if (count[face] != 1) // inner face
2477 continue;
2478
2479 // boundary faces ...
2480 bids_face[face].boundary_id = 0;
2481
2482 if (dim != 3)
2483 continue;
2484
2485 // ... and the lines of quads in 3D
2486 const auto &crs = connectivity.entity_to_entities(2, 1);
2487 for (unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
2488 tria.faces->lines.boundary_or_material_id[crs.col[i]]
2489 .boundary_id = 0;
2490 }
2491 }
2492 else // 1D
2493 {
2494 static const unsigned int t_tba = static_cast<unsigned int>(-1);
2495 static const unsigned int t_inner = static_cast<unsigned int>(-2);
2496
2497 std::vector<unsigned int> type(vertices.size(), t_tba);
2498
2499 const auto &crs = connectivity.entity_to_entities(1, 0);
2500
2501 for (unsigned int cell = 0; cell < cells.size(); ++cell)
2502 for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2503 ++i, ++j)
2504 if (type[crs.col[i]] != t_inner)
2505 type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
2506
2507 for (unsigned int face = 0; face < type.size(); ++face)
2508 {
2509 // note: we also treat manifolds here!?
2510 (*tria.vertex_to_manifold_id_map_1d)[face] =
2512 if (type[face] != t_inner && type[face] != t_tba)
2513 (*tria.vertex_to_boundary_id_map_1d)[face] = type[face];
2514 }
2515 }
2516
2517 // SubCellData: line
2518 if (dim >= 2)
2519 process_subcelldata(connectivity.entity_to_entities(1, 0),
2520 tria.faces->lines,
2521 subcelldata.boundary_lines);
2522
2523 // SubCellData: quad
2524 if (dim == 3)
2525 process_subcelldata(connectivity.entity_to_entities(2, 0),
2526 tria.faces->quads,
2527 subcelldata.boundary_quads);
2528 }
2529
2530
2531 template <int structdim, typename T>
2532 static void
2534 const CRS<T> & crs,
2535 TriaObjects & obj,
2536 const std::vector<CellData<structdim>> &boundary_objects_in)
2537 {
2538 AssertDimension(obj.structdim, structdim);
2539
2540 if (boundary_objects_in.size() == 0)
2541 return; // empty subcelldata -> nothing to do
2542
2543 // pre-sort subcelldata
2544 auto boundary_objects = boundary_objects_in;
2545
2546 // ... sort vertices
2547 for (auto &boundary_object : boundary_objects)
2548 std::sort(boundary_object.vertices.begin(),
2549 boundary_object.vertices.end());
2550
2551 // ... sort cells
2552 std::sort(boundary_objects.begin(),
2553 boundary_objects.end(),
2554 [](const auto &a, const auto &b) {
2555 return a.vertices < b.vertices;
2556 });
2557
2558 unsigned int counter = 0;
2559
2560 std::vector<unsigned int> key;
2562
2563 for (unsigned int o = 0; o < obj.n_objects(); ++o)
2564 {
2565 auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
2566 auto &manifold_id = obj.manifold_id[o];
2567
2568 // assert that object has not been visited yet and its value
2569 // has not been modified yet
2570 AssertThrow(boundary_id == 0 ||
2575
2576 // create key
2577 key.assign(crs.col.data() + crs.ptr[o],
2578 crs.col.data() + crs.ptr[o + 1]);
2579 std::sort(key.begin(), key.end());
2580
2581 // is subcelldata provided? -> binary search
2582 const auto subcell_object =
2583 std::lower_bound(boundary_objects.begin(),
2584 boundary_objects.end(),
2585 key,
2586 [&](const auto &cell, const auto &key) {
2587 return cell.vertices < key;
2588 });
2589
2590 // no subcelldata provided for this object
2591 if (subcell_object == boundary_objects.end() ||
2592 subcell_object->vertices != key)
2593 continue;
2594
2595 counter++;
2596
2597 // set manifold id
2598 manifold_id = subcell_object->manifold_id;
2599
2600 // set boundary id
2601 if (subcell_object->boundary_id !=
2603 {
2606 boundary_id = subcell_object->boundary_id;
2607 }
2608 }
2609
2610 // make sure that all subcelldata entries have been processed
2611 // TODO: this is not guaranteed, why?
2612 // AssertDimension(counter, boundary_objects_in.size());
2613 }
2614
2615
2616
2617 static void
2619 const unsigned structdim,
2620 const unsigned int size)
2621 {
2622 const unsigned int dim = faces.dim;
2623
2624 const unsigned int max_faces_per_cell = 2 * structdim;
2625
2626 if (dim == 3 && structdim == 2)
2627 {
2628 // quad entity types
2629 faces.quad_reference_cell.assign(size,
2631
2632 // quad line orientations
2633 faces.quads_line_orientations.assign(size * max_faces_per_cell, -1);
2634 }
2635 }
2636
2637
2638
2639 static void
2641 const unsigned int spacedim,
2642 const unsigned int size,
2643 const bool orientation_needed)
2644 {
2645 const unsigned int dim = level.dim;
2646
2647 const unsigned int max_faces_per_cell = 2 * dim;
2648
2649 level.active_cell_indices.assign(size, -1);
2650 level.subdomain_ids.assign(size, 0);
2651 level.level_subdomain_ids.assign(size, 0);
2652
2653 level.refine_flags.assign(size, false);
2654 level.coarsen_flags.assign(size, false);
2655
2656 level.parents.assign((size + 1) / 2, -1);
2657
2658 if (dim < spacedim)
2659 level.direction_flags.assign(size, true);
2660
2661 level.neighbors.assign(size * max_faces_per_cell, {-1, -1});
2662
2663 level.reference_cell.assign(size, ::ReferenceCells::Invalid);
2664
2665 if (orientation_needed)
2666 level.face_orientations.assign(size * max_faces_per_cell, -1);
2667
2668 level.global_active_cell_indices.assign(size,
2670 level.global_level_cell_indices.assign(size,
2672 }
2673
2674
2675
2676 static void
2677 reserve_space_(TriaObjects &obj, const unsigned int size)
2678 {
2679 const unsigned int structdim = obj.structdim;
2680
2681 const unsigned int max_children_per_cell = 1 << structdim;
2682 const unsigned int max_faces_per_cell = 2 * structdim;
2683
2684 obj.used.assign(size, true);
2685 obj.boundary_or_material_id.assign(
2686 size,
2688 BoundaryOrMaterialId());
2689 obj.manifold_id.assign(size, -1);
2690 obj.user_flags.assign(size, false);
2691 obj.user_data.resize(size);
2692
2693 if (structdim > 1) // TODO: why?
2694 obj.refinement_cases.assign(size, 0);
2695
2696 obj.children.assign(max_children_per_cell / 2 * size, -1);
2697
2698 obj.cells.assign(max_faces_per_cell * size, -1);
2699
2700 if (structdim <= 2)
2701 {
2702 obj.next_free_single = size - 1;
2703 obj.next_free_pair = 0;
2705 }
2706 else
2707 {
2708 obj.next_free_single = obj.next_free_pair = 0;
2709 }
2710 }
2711
2712
2728 template <int spacedim>
2729 static void delete_children(
2732 std::vector<unsigned int> &,
2733 std::vector<unsigned int> &)
2734 {
2735 const unsigned int dim = 1;
2736
2737 // first we need to reset the
2738 // neighbor pointers of the
2739 // neighbors of this cell's
2740 // children to this cell. This is
2741 // different for one dimension,
2742 // since there neighbors can have a
2743 // refinement level differing from
2744 // that of this cell's children by
2745 // more than one level.
2746
2747 Assert(!cell->child(0)->has_children() &&
2748 !cell->child(1)->has_children(),
2750
2751 // first do it for the cells to the
2752 // left
2753 if (cell->neighbor(0).state() == IteratorState::valid)
2754 if (cell->neighbor(0)->has_children())
2755 {
2757 cell->neighbor(0);
2758 Assert(neighbor->level() == cell->level(), ExcInternalError());
2759
2760 // right child
2761 neighbor = neighbor->child(1);
2762 while (true)
2763 {
2764 Assert(neighbor->neighbor(1) == cell->child(0),
2766 neighbor->set_neighbor(1, cell);
2767
2768 // move on to further
2769 // children on the
2770 // boundary between this
2771 // cell and its neighbor
2772 if (neighbor->has_children())
2773 neighbor = neighbor->child(1);
2774 else
2775 break;
2776 }
2777 }
2778
2779 // now do it for the cells to the
2780 // left
2781 if (cell->neighbor(1).state() == IteratorState::valid)
2782 if (cell->neighbor(1)->has_children())
2783 {
2785 cell->neighbor(1);
2786 Assert(neighbor->level() == cell->level(), ExcInternalError());
2787
2788 // left child
2789 neighbor = neighbor->child(0);
2790 while (true)
2791 {
2792 Assert(neighbor->neighbor(0) == cell->child(1),
2794 neighbor->set_neighbor(0, cell);
2795
2796 // move on to further
2797 // children on the
2798 // boundary between this
2799 // cell and its neighbor
2800 if (neighbor->has_children())
2801 neighbor = neighbor->child(0);
2802 else
2803 break;
2804 }
2805 }
2806
2807
2808 // delete the vertex which will not
2809 // be needed anymore. This vertex
2810 // is the second of the first child
2811 triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
2812
2813 // invalidate children. clear user
2814 // pointers, to avoid that they may
2815 // appear at unwanted places later
2816 // on...
2817 for (unsigned int child = 0; child < cell->n_children(); ++child)
2818 {
2819 cell->child(child)->clear_user_data();
2820 cell->child(child)->clear_user_flag();
2821 cell->child(child)->clear_used_flag();
2822 }
2823
2824
2825 // delete pointer to children
2826 cell->clear_children();
2827 cell->clear_user_flag();
2828 }
2829
2830
2831
2832 template <int spacedim>
2833 static void delete_children(
2836 std::vector<unsigned int> & line_cell_count,
2837 std::vector<unsigned int> &)
2838 {
2839 const unsigned int dim = 2;
2840 const RefinementCase<dim> ref_case = cell->refinement_case();
2841
2842 Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2844
2845 // vectors to hold all lines which
2846 // may be deleted
2847 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2848 lines_to_delete(0);
2849
2850 lines_to_delete.reserve(4 * 2 + 4);
2851
2852 // now we decrease the counters for
2853 // lines contained in the child
2854 // cells
2855 for (unsigned int c = 0; c < cell->n_children(); ++c)
2856 {
2858 cell->child(c);
2859 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2860 --line_cell_count[child->line_index(l)];
2861 }
2862
2863
2864 // delete the vertex which will not
2865 // be needed anymore. This vertex
2866 // is the second of the second line
2867 // of the first child, if the cell
2868 // is refined with cut_xy, else there
2869 // is no inner vertex.
2870 // additionally delete unneeded inner
2871 // lines
2872 if (ref_case == RefinementCase<dim>::cut_xy)
2873 {
2875 .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
2876
2877 lines_to_delete.push_back(cell->child(0)->line(1));
2878 lines_to_delete.push_back(cell->child(0)->line(3));
2879 lines_to_delete.push_back(cell->child(3)->line(0));
2880 lines_to_delete.push_back(cell->child(3)->line(2));
2881 }
2882 else
2883 {
2884 unsigned int inner_face_no =
2885 ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
2886
2887 // the inner line will not be
2888 // used any more
2889 lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2890 }
2891
2892 // invalidate children
2893 for (unsigned int child = 0; child < cell->n_children(); ++child)
2894 {
2895 cell->child(child)->clear_user_data();
2896 cell->child(child)->clear_user_flag();
2897 cell->child(child)->clear_used_flag();
2898 }
2899
2900
2901 // delete pointer to children
2902 cell->clear_children();
2903 cell->clear_refinement_case();
2904 cell->clear_user_flag();
2905
2906 // look at the refinement of outer
2907 // lines. if nobody needs those
2908 // anymore we can add them to the
2909 // list of lines to be deleted.
2910 for (unsigned int line_no = 0;
2911 line_no < GeometryInfo<dim>::lines_per_cell;
2912 ++line_no)
2913 {
2915 cell->line(line_no);
2916
2917 if (line->has_children())
2918 {
2919 // if one of the cell counters is
2920 // zero, the other has to be as well
2921
2922 Assert((line_cell_count[line->child_index(0)] == 0 &&
2923 line_cell_count[line->child_index(1)] == 0) ||
2924 (line_cell_count[line->child_index(0)] > 0 &&
2925 line_cell_count[line->child_index(1)] > 0),
2927
2928 if (line_cell_count[line->child_index(0)] == 0)
2929 {
2930 for (unsigned int c = 0; c < 2; ++c)
2931 Assert(!line->child(c)->has_children(),
2933
2934 // we may delete the line's
2935 // children and the middle vertex
2936 // as no cell references them
2937 // anymore
2939 .vertices_used[line->child(0)->vertex_index(1)] = false;
2940
2941 lines_to_delete.push_back(line->child(0));
2942 lines_to_delete.push_back(line->child(1));
2943
2944 line->clear_children();
2945 }
2946 }
2947 }
2948
2949 // finally, delete unneeded lines
2950
2951 // clear user pointers, to avoid that
2952 // they may appear at unwanted places
2953 // later on...
2954 // same for user flags, then finally
2955 // delete the lines
2956 typename std::vector<
2958 line = lines_to_delete.begin(),
2959 endline = lines_to_delete.end();
2960 for (; line != endline; ++line)
2961 {
2962 (*line)->clear_user_data();
2963 (*line)->clear_user_flag();
2964 (*line)->clear_used_flag();
2965 }
2966 }
2967
2968
2969
2970 template <int spacedim>
2971 static void delete_children(
2974 std::vector<unsigned int> & line_cell_count,
2975 std::vector<unsigned int> & quad_cell_count)
2976 {
2977 const unsigned int dim = 3;
2978
2979 Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2981 Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
2983
2984 // first of all, we store the RefineCase of
2985 // this cell
2986 const RefinementCase<dim> ref_case = cell->refinement_case();
2987 // vectors to hold all lines and quads which
2988 // may be deleted
2989 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2990 lines_to_delete(0);
2991 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
2992 quads_to_delete(0);
2993
2994 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
2995 quads_to_delete.reserve(6 * 4 + 12);
2996
2997 // now we decrease the counters for lines and
2998 // quads contained in the child cells
2999 for (unsigned int c = 0; c < cell->n_children(); ++c)
3000 {
3002 cell->child(c);
3003 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3004 --line_cell_count[child->line_index(l)];
3005 for (auto f : GeometryInfo<dim>::face_indices())
3006 --quad_cell_count[child->quad_index(f)];
3007 }
3008
3010 // delete interior quads and lines and the
3011 // interior vertex, depending on the
3012 // refinement case of the cell
3013 //
3014 // for append quads and lines: only append
3015 // them to the list of objects to be deleted
3016
3017 switch (ref_case)
3018 {
3020 quads_to_delete.push_back(cell->child(0)->face(1));
3021 break;
3023 quads_to_delete.push_back(cell->child(0)->face(3));
3024 break;
3026 quads_to_delete.push_back(cell->child(0)->face(5));
3027 break;
3029 quads_to_delete.push_back(cell->child(0)->face(1));
3030 quads_to_delete.push_back(cell->child(0)->face(3));
3031 quads_to_delete.push_back(cell->child(3)->face(0));
3032 quads_to_delete.push_back(cell->child(3)->face(2));
3033
3034 lines_to_delete.push_back(cell->child(0)->line(11));
3035 break;
3037 quads_to_delete.push_back(cell->child(0)->face(1));
3038 quads_to_delete.push_back(cell->child(0)->face(5));
3039 quads_to_delete.push_back(cell->child(3)->face(0));
3040 quads_to_delete.push_back(cell->child(3)->face(4));
3041
3042 lines_to_delete.push_back(cell->child(0)->line(5));
3043 break;
3045 quads_to_delete.push_back(cell->child(0)->face(3));
3046 quads_to_delete.push_back(cell->child(0)->face(5));
3047 quads_to_delete.push_back(cell->child(3)->face(2));
3048 quads_to_delete.push_back(cell->child(3)->face(4));
3049
3050 lines_to_delete.push_back(cell->child(0)->line(7));
3051 break;
3053 quads_to_delete.push_back(cell->child(0)->face(1));
3054 quads_to_delete.push_back(cell->child(2)->face(1));
3055 quads_to_delete.push_back(cell->child(4)->face(1));
3056 quads_to_delete.push_back(cell->child(6)->face(1));
3057
3058 quads_to_delete.push_back(cell->child(0)->face(3));
3059 quads_to_delete.push_back(cell->child(1)->face(3));
3060 quads_to_delete.push_back(cell->child(4)->face(3));
3061 quads_to_delete.push_back(cell->child(5)->face(3));
3062
3063 quads_to_delete.push_back(cell->child(0)->face(5));
3064 quads_to_delete.push_back(cell->child(1)->face(5));
3065 quads_to_delete.push_back(cell->child(2)->face(5));
3066 quads_to_delete.push_back(cell->child(3)->face(5));
3067
3068 lines_to_delete.push_back(cell->child(0)->line(5));
3069 lines_to_delete.push_back(cell->child(0)->line(7));
3070 lines_to_delete.push_back(cell->child(0)->line(11));
3071 lines_to_delete.push_back(cell->child(7)->line(0));
3072 lines_to_delete.push_back(cell->child(7)->line(2));
3073 lines_to_delete.push_back(cell->child(7)->line(8));
3074 // delete the vertex which will not
3075 // be needed anymore. This vertex
3076 // is the vertex at the heart of
3077 // this cell, which is the sixth of
3078 // the first child
3079 triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3080 false;
3081 break;
3082 default:
3083 // only remaining case is
3084 // no_refinement, thus an error
3085 Assert(false, ExcInternalError());
3086 break;
3087 }
3088
3089
3090 // invalidate children
3091 for (unsigned int child = 0; child < cell->n_children(); ++child)
3092 {
3093 cell->child(child)->clear_user_data();
3094 cell->child(child)->clear_user_flag();
3095
3096 for (auto f : GeometryInfo<dim>::face_indices())
3097 {
3098 // set flags denoting deviations from
3099 // standard orientation of faces back
3100 // to initialization values
3101 cell->child(child)->set_face_orientation(f, true);
3102 cell->child(child)->set_face_flip(f, false);
3103 cell->child(child)->set_face_rotation(f, false);
3104 }
3105
3106 cell->child(child)->clear_used_flag();
3107 }
3108
3109
3110 // delete pointer to children
3111 cell->clear_children();
3112 cell->clear_refinement_case();
3113 cell->clear_user_flag();
3114
3115 // so far we only looked at inner quads,
3116 // lines and vertices. Now we have to
3117 // consider outer ones as well. here, we have
3118 // to check, whether there are other cells
3119 // still needing these objects. otherwise we
3120 // can delete them. first for quads (and
3121 // their inner lines).
3122
3123 for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
3124 {
3126 cell->face(quad_no);
3127
3128 Assert(
3129 (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
3130 quad->has_children()) ||
3131 GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
3134
3135 switch (quad->refinement_case())
3136 {
3137 case RefinementCase<dim - 1>::no_refinement:
3138 // nothing to do as the quad
3139 // is not refined
3140 break;
3141 case RefinementCase<dim - 1>::cut_x:
3142 case RefinementCase<dim - 1>::cut_y:
3143 {
3144 // if one of the cell counters is
3145 // zero, the other has to be as
3146 // well
3147 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3148 quad_cell_count[quad->child_index(1)] == 0) ||
3149 (quad_cell_count[quad->child_index(0)] > 0 &&
3150 quad_cell_count[quad->child_index(1)] > 0),
3152 // it might be, that the quad is
3153 // refined twice anisotropically,
3154 // first check, whether we may
3155 // delete possible grand_children
3156 unsigned int deleted_grandchildren = 0;
3157 unsigned int number_of_child_refinements = 0;
3158
3159 for (unsigned int c = 0; c < 2; ++c)
3160 if (quad->child(c)->has_children())
3161 {
3162 ++number_of_child_refinements;
3163 // if one of the cell counters is
3164 // zero, the other has to be as
3165 // well
3166 Assert(
3167 (quad_cell_count[quad->child(c)->child_index(0)] ==
3168 0 &&
3169 quad_cell_count[quad->child(c)->child_index(1)] ==
3170 0) ||
3171 (quad_cell_count[quad->child(c)->child_index(0)] >
3172 0 &&
3173 quad_cell_count[quad->child(c)->child_index(1)] >
3174 0),
3176 if (quad_cell_count[quad->child(c)->child_index(0)] ==
3177 0)
3178 {
3179 // Assert, that the two
3180 // anisotropic
3181 // refinements add up to
3182 // isotropic refinement
3183 Assert(quad->refinement_case() +
3184 quad->child(c)->refinement_case() ==
3187 // we may delete the
3188 // quad's children and
3189 // the inner line as no
3190 // cell references them
3191 // anymore
3192 quads_to_delete.push_back(
3193 quad->child(c)->child(0));
3194 quads_to_delete.push_back(
3195 quad->child(c)->child(1));
3196 if (quad->child(c)->refinement_case() ==
3198 lines_to_delete.push_back(
3199 quad->child(c)->child(0)->line(1));
3200 else
3201 lines_to_delete.push_back(
3202 quad->child(c)->child(0)->line(3));
3203 quad->child(c)->clear_children();
3204 quad->child(c)->clear_refinement_case();
3205 ++deleted_grandchildren;
3206 }
3207 }
3208 // if no grandchildren are left, we
3209 // may as well delete the
3210 // refinement of the inner line
3211 // between our children and the
3212 // corresponding vertex
3213 if (number_of_child_refinements > 0 &&
3214 deleted_grandchildren == number_of_child_refinements)
3215 {
3217 middle_line;
3218 if (quad->refinement_case() == RefinementCase<2>::cut_x)
3219 middle_line = quad->child(0)->line(1);
3220 else
3221 middle_line = quad->child(0)->line(3);
3222
3223 lines_to_delete.push_back(middle_line->child(0));
3224 lines_to_delete.push_back(middle_line->child(1));
3226 .vertices_used[middle_vertex_index<dim, spacedim>(
3227 middle_line)] = false;
3228 middle_line->clear_children();
3229 }
3230
3231 // now consider the direct children
3232 // of the given quad
3233 if (quad_cell_count[quad->child_index(0)] == 0)
3234 {
3235 // we may delete the quad's
3236 // children and the inner line
3237 // as no cell references them
3238 // anymore
3239 quads_to_delete.push_back(quad->child(0));
3240 quads_to_delete.push_back(quad->child(1));
3241 if (quad->refinement_case() == RefinementCase<2>::cut_x)
3242 lines_to_delete.push_back(quad->child(0)->line(1));
3243 else
3244 lines_to_delete.push_back(quad->child(0)->line(3));
3245
3246 // if the counters just dropped
3247 // to zero, otherwise the
3248 // children would have been
3249 // deleted earlier, then this
3250 // cell's children must have
3251 // contained the anisotropic
3252 // quad children. thus, if
3253 // those have again anisotropic
3254 // children, which are in
3255 // effect isotropic children of
3256 // the original quad, those are
3257 // still needed by a
3258 // neighboring cell and we
3259 // cannot delete them. instead,
3260 // we have to reset this quad's
3261 // refine case to isotropic and
3262 // set the children
3263 // accordingly.
3264 if (quad->child(0)->has_children())
3265 if (quad->refinement_case() ==
3267 {
3268 // now evereything is
3269 // quite complicated. we
3270 // have the children
3271 // numbered according to
3272 //
3273 // *---*---*
3274 // |n+1|m+1|
3275 // *---*---*
3276 // | n | m |
3277 // *---*---*
3278 //
3279 // from the original
3280 // anisotropic
3281 // refinement. we have to
3282 // reorder them as
3283 //
3284 // *---*---*
3285 // | m |m+1|
3286 // *---*---*
3287 // | n |n+1|
3288 // *---*---*
3289 //
3290 // for isotropic refinement.
3291 //
3292 // this is a bit ugly, of
3293 // course: loop over all
3294 // cells on all levels
3295 // and look for faces n+1
3296 // (switch_1) and m
3297 // (switch_2).
3298 const typename Triangulation<dim, spacedim>::
3299 quad_iterator switch_1 =
3300 quad->child(0)->child(1),
3301 switch_2 =
3302 quad->child(1)->child(0);
3303
3304 Assert(!switch_1->has_children(),
3306 Assert(!switch_2->has_children(),
3308
3309 const int switch_1_index = switch_1->index();
3310 const int switch_2_index = switch_2->index();
3311 for (unsigned int l = 0;
3312 l < triangulation.levels.size();
3313 ++l)
3314 for (unsigned int h = 0;
3315 h <
3316 triangulation.levels[l]->cells.n_objects();
3317 ++h)
3318 for (const unsigned int q :
3320 {
3321 const int index =
3322 triangulation.levels[l]
3323 ->cells.get_bounding_object_indices(
3324 h)[q];
3325 if (index == switch_1_index)
3326 triangulation.levels[l]
3327 ->cells.get_bounding_object_indices(
3328 h)[q] = switch_2_index;
3329 else if (index == switch_2_index)
3330 triangulation.levels[l]
3331 ->cells.get_bounding_object_indices(
3332 h)[q] = switch_1_index;
3333 }
3334 // now we have to copy
3335 // all information of the
3336 // two quads
3337 const int switch_1_lines[4] = {
3338 static_cast<signed int>(
3339 switch_1->line_index(0)),
3340 static_cast<signed int>(
3341 switch_1->line_index(1)),
3342 static_cast<signed int>(
3343 switch_1->line_index(2)),
3344 static_cast<signed int>(
3345 switch_1->line_index(3))};
3346 const bool switch_1_line_orientations[4] = {
3347 switch_1->line_orientation(0),
3348 switch_1->line_orientation(1),
3349 switch_1->line_orientation(2),
3350 switch_1->line_orientation(3)};
3351 const types::boundary_id switch_1_boundary_id =
3352 switch_1->boundary_id();
3353 const unsigned int switch_1_user_index =
3354 switch_1->user_index();
3355 const bool switch_1_user_flag =
3356 switch_1->user_flag_set();
3357
3358 switch_1->set_bounding_object_indices(
3359 {switch_2->line_index(0),
3360 switch_2->line_index(1),
3361 switch_2->line_index(2),
3362 switch_2->line_index(3)});
3363 switch_1->set_line_orientation(
3364 0, switch_2->line_orientation(0));
3365 switch_1->set_line_orientation(
3366 1, switch_2->line_orientation(1));
3367 switch_1->set_line_orientation(
3368 2, switch_2->line_orientation(2));
3369 switch_1->set_line_orientation(
3370 3, switch_2->line_orientation(3));
3371 switch_1->set_boundary_id_internal(
3372 switch_2->boundary_id());
3373 switch_1->set_manifold_id(
3374 switch_2->manifold_id());
3375 switch_1->set_user_index(switch_2->user_index());
3376 if (switch_2->user_flag_set())
3377 switch_1->set_user_flag();
3378 else
3379 switch_1->clear_user_flag();
3380
3381 switch_2->set_bounding_object_indices(
3382 {switch_1_lines[0],
3383 switch_1_lines[1],
3384 switch_1_lines[2],
3385 switch_1_lines[3]});
3386 switch_2->set_line_orientation(
3387 0, switch_1_line_orientations[0]);
3388 switch_2->set_line_orientation(
3389 1, switch_1_line_orientations[1]);
3390 switch_2->set_line_orientation(
3391 2, switch_1_line_orientations[2]);
3392 switch_2->set_line_orientation(
3393 3, switch_1_line_orientations[3]);
3394 switch_2->set_boundary_id_internal(
3395 switch_1_boundary_id);
3396 switch_2->set_manifold_id(
3397 switch_1->manifold_id());
3398 switch_2->set_user_index(switch_1_user_index);
3399 if (switch_1_user_flag)
3400 switch_2->set_user_flag();
3401 else
3402 switch_2->clear_user_flag();
3403
3404 const unsigned int child_0 =
3405 quad->child(0)->child_index(0);
3406 const unsigned int child_2 =
3407 quad->child(1)->child_index(0);
3408 quad->clear_children();
3409 quad->clear_refinement_case();
3410 quad->set_refinement_case(
3412 quad->set_children(0, child_0);
3413 quad->set_children(2, child_2);
3414 std::swap(quad_cell_count[child_0 + 1],
3415 quad_cell_count[child_2]);
3416 }
3417 else
3418 {
3419 // the face was refined
3420 // with cut_y, thus the
3421 // children are already
3422 // in correct order. we
3423 // only have to set them
3424 // correctly, deleting
3425 // the indirection of two
3426 // anisotropic refinement
3427 // and going directly
3428 // from the quad to
3429 // isotropic children
3430 const unsigned int child_0 =
3431 quad->child(0)->child_index(0);
3432 const unsigned int child_2 =
3433 quad->child(1)->child_index(0);
3434 quad->clear_children();
3435 quad->clear_refinement_case();
3436 quad->set_refinement_case(
3438 quad->set_children(0, child_0);
3439 quad->set_children(2, child_2);
3440 }
3441 else
3442 {
3443 quad->clear_children();
3444 quad->clear_refinement_case();
3445 }
3446 }
3447 break;
3448 }
3449 case RefinementCase<dim - 1>::cut_xy:
3450 {
3451 // if one of the cell counters is
3452 // zero, the others have to be as
3453 // well
3454
3455 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3456 quad_cell_count[quad->child_index(1)] == 0 &&
3457 quad_cell_count[quad->child_index(2)] == 0 &&
3458 quad_cell_count[quad->child_index(3)] == 0) ||
3459 (quad_cell_count[quad->child_index(0)] > 0 &&
3460 quad_cell_count[quad->child_index(1)] > 0 &&
3461 quad_cell_count[quad->child_index(2)] > 0 &&
3462 quad_cell_count[quad->child_index(3)] > 0),
3464
3465 if (quad_cell_count[quad->child_index(0)] == 0)
3466 {
3467 // we may delete the quad's
3468 // children, the inner lines
3469 // and the middle vertex as no
3470 // cell references them anymore
3471 lines_to_delete.push_back(quad->child(0)->line(1));
3472 lines_to_delete.push_back(quad->child(3)->line(0));
3473 lines_to_delete.push_back(quad->child(0)->line(3));
3474 lines_to_delete.push_back(quad->child(3)->line(2));
3475
3476 for (unsigned int child = 0; child < quad->n_children();
3477 ++child)
3478 quads_to_delete.push_back(quad->child(child));
3479
3481 .vertices_used[quad->child(0)->vertex_index(3)] =
3482 false;
3483
3484 quad->clear_children();
3485 quad->clear_refinement_case();
3486 }
3487 }
3488 break;
3489
3490 default:
3491 Assert(false, ExcInternalError());
3492 break;
3493 }
3494 }
3495
3496 // now we repeat a similar procedure
3497 // for the outer lines of this cell.
3498
3499 // if in debug mode: check that each
3500 // of the lines for which we consider
3501 // deleting the children in fact has
3502 // children (the bits/coarsening_3d
3503 // test tripped over this initially)
3504 for (unsigned int line_no = 0;
3505 line_no < GeometryInfo<dim>::lines_per_cell;
3506 ++line_no)
3507 {
3509 cell->line(line_no);
3510
3511 Assert(
3512 (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
3513 line->has_children()) ||
3514 GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
3517
3518 if (line->has_children())
3519 {
3520 // if one of the cell counters is
3521 // zero, the other has to be as well
3522
3523 Assert((line_cell_count[line->child_index(0)] == 0 &&
3524 line_cell_count[line->child_index(1)] == 0) ||
3525 (line_cell_count[line->child_index(0)] > 0 &&
3526 line_cell_count[line->child_index(1)] > 0),
3528
3529 if (line_cell_count[line->child_index(0)] == 0)
3530 {
3531 for (unsigned int c = 0; c < 2; ++c)
3532 Assert(!line->child(c)->has_children(),
3534
3535 // we may delete the line's
3536 // children and the middle vertex
3537 // as no cell references them
3538 // anymore
3540 .vertices_used[line->child(0)->vertex_index(1)] = false;
3541
3542 lines_to_delete.push_back(line->child(0));
3543 lines_to_delete.push_back(line->child(1));
3544
3545 line->clear_children();
3546 }
3547 }
3548 }
3549
3550 // finally, delete unneeded quads and lines
3551
3552 // clear user pointers, to avoid that
3553 // they may appear at unwanted places
3554 // later on...
3555 // same for user flags, then finally
3556 // delete the quads and lines
3557 typename std::vector<
3559 line = lines_to_delete.begin(),
3560 endline = lines_to_delete.end();
3561 for (; line != endline; ++line)
3562 {
3563 (*line)->clear_user_data();
3564 (*line)->clear_user_flag();
3565 (*line)->clear_used_flag();
3566 }
3567
3568 typename std::vector<
3570 quad = quads_to_delete.begin(),
3571 endquad = quads_to_delete.end();
3572 for (; quad != endquad; ++quad)
3573 {
3574 (*quad)->clear_user_data();
3575 (*quad)->clear_children();
3576 (*quad)->clear_refinement_case();
3577 (*quad)->clear_user_flag();
3578 (*quad)->clear_used_flag();
3579 }
3580 }
3581
3582
3600 template <int spacedim>
3601 static void create_children(
3603 unsigned int & next_unused_vertex,
3605 &next_unused_line,
3607 &next_unused_cell,
3608 const typename Triangulation<2, spacedim>::cell_iterator &cell)
3609 {
3610 const unsigned int dim = 2;
3611 // clear refinement flag
3612 const RefinementCase<dim> ref_case = cell->refine_flag_set();
3613 cell->clear_refine_flag();
3614
3615 /* For the refinement process: since we go the levels up from the
3616 lowest, there are (unlike above) only two possibilities: a neighbor
3617 cell is on the same level or one level up (in both cases, it may or
3618 may not be refined later on, but we don't care here).
3619
3620 First:
3621 Set up an array of the 3x3 vertices, which are distributed on the
3622 cell (the array consists of indices into the @p{vertices} std::vector
3623
3624 2--7--3
3625 | | |
3626 4--8--5
3627 | | |
3628 0--6--1
3629
3630 note: in case of cut_x or cut_y not all these vertices are needed for
3631 the new cells
3632
3633 Second:
3634 Set up an array of the new lines (the array consists of iterator
3635 pointers into the lines arrays)
3636
3637 .-6-.-7-. The directions are: .->-.->-.
3638 1 9 3 ^ ^ ^
3639 .-10.11-. .->-.->-.
3640 0 8 2 ^ ^ ^
3641 .-4-.-5-. .->-.->-.
3642
3643 cut_x:
3644 .-4-.-5-.
3645 | | |
3646 0 6 1
3647 | | |
3648 .-2-.-3-.
3649
3650 cut_y:
3651 .---5---.
3652 1 3
3653 .---6---.
3654 0 2
3655 .---4---.
3656
3657
3658 Third:
3659 Set up an array of neighbors:
3660
3661 6 7
3662 .--.--.
3663 1| | |3
3664 .--.--.
3665 0| | |2
3666 .--.--.
3667 4 5
3668
3669 We need this array for two reasons: first to get the lines which will
3670 bound the four subcells (if the neighboring cell is refined, these
3671 lines already exist), and second to update neighborship information.
3672 Since if a neighbor is not refined, its neighborship record only
3673 points to the present, unrefined, cell rather than the children we
3674 are presently creating, we only need the neighborship information
3675 if the neighbor cells are refined. In all other cases, we store
3676 the unrefined neighbor address
3677
3678 We also need for every neighbor (if refined) which number among its
3679 neighbors the present (unrefined) cell has, since that number is to
3680 be replaced and because that also is the number of the subline which
3681 will be the interface between that neighbor and the to be created
3682 cell. We will store this number (between 0 and 3) in the field
3683 @p{neighbors_neighbor}.
3684
3685 It would be sufficient to use the children of the common line to the
3686 neighbor, if we only wanted to get the new sublines and the new
3687 vertex, but because we need to update the neighborship information of
3688 the two refined subcells of the neighbor, we need to search these
3689 anyway.
3690
3691 Convention:
3692 The created children are numbered like this:
3693
3694 .--.--.
3695 |2 . 3|
3696 .--.--.
3697 |0 | 1|
3698 .--.--.
3699 */
3700 // collect the
3701 // indices of the
3702 // eight
3703 // surrounding
3704 // vertices
3705 // 2--7--3
3706 // | | |
3707 // 4--9--5
3708 // | | |
3709 // 0--6--1
3710 int new_vertices[9];
3711 for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3712 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3713 for (unsigned int line_no = 0; line_no < 4; ++line_no)
3714 if (cell->line(line_no)->has_children())
3715 new_vertices[4 + line_no] =
3716 cell->line(line_no)->child(0)->vertex_index(1);
3717
3718 if (ref_case == RefinementCase<dim>::cut_xy)
3719 {
3720 // find the next
3721 // unused vertex and
3722 // allocate it for
3723 // the new vertex we
3724 // need here
3725 while (triangulation.vertices_used[next_unused_vertex] == true)
3726 ++next_unused_vertex;
3727 Assert(
3728 next_unused_vertex < triangulation.vertices.size(),
3729 ExcMessage(
3730 "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."));
3731 triangulation.vertices_used[next_unused_vertex] = true;
3732
3733 new_vertices[8] = next_unused_vertex;
3734
3735 // if this quad lives
3736 // in 2d, then we can
3737 // compute the new
3738 // central vertex
3739 // location just from
3740 // the surrounding
3741 // ones. If this is
3742 // not the case, then
3743 // we need to ask a
3744 // boundary object
3745 if (dim == spacedim)
3746 {
3747 // triangulation.vertices[next_unused_vertex] = new_point;
3748 triangulation.vertices[next_unused_vertex] = cell->center(true);
3749
3750 // if the user_flag is set, i.e. if the cell is at the
3751 // boundary, use a different calculation of the middle vertex
3752 // here. this is of advantage if the boundary is strongly
3753 // curved (whereas the cell is not) and the cell has a high
3754 // aspect ratio.
3755 if (cell->user_flag_set())
3756 {
3757 // first reset the user_flag and then refine
3758 cell->clear_user_flag();
3759 triangulation.vertices[next_unused_vertex] =
3760 cell->center(true, true);
3761 }
3762 }
3763 else
3764 {
3765 // if this quad lives in a higher dimensional space
3766 // then we don't need to worry if it is at the
3767 // boundary of the manifold -- we always have to use
3768 // the boundary object anyway; so ignore whether the
3769 // user flag is set or not
3770 cell->clear_user_flag();
3771
3772 // determine middle vertex by transfinite interpolation to be
3773 // consistent with what happens to quads in a Triangulation<3,
3774 // 3> when they are refined
3775 triangulation.vertices[next_unused_vertex] =
3776 cell->center(true, true);
3777 }
3778 }
3779
3780
3781 // Now the lines:
3783 unsigned int lmin = 8;
3784 unsigned int lmax = 12;
3785 if (ref_case != RefinementCase<dim>::cut_xy)
3786 {
3787 lmin = 6;
3788 lmax = 7;
3789 }
3790
3791 for (unsigned int l = lmin; l < lmax; ++l)
3792 {
3793 while (next_unused_line->used() == true)
3794 ++next_unused_line;
3795 new_lines[l] = next_unused_line;
3796 ++next_unused_line;
3797
3798 Assert(
3799 new_lines[l]->used() == false,
3800 ExcMessage(
3801 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3802 }
3803
3804 if (ref_case == RefinementCase<dim>::cut_xy)
3805 {
3806 // .-6-.-7-.
3807 // 1 9 3
3808 // .-10.11-.
3809 // 0 8 2
3810 // .-4-.-5-.
3811
3812 // lines 0-7 already exist, create only the four interior
3813 // lines 8-11
3814 unsigned int l = 0;
3815 for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
3816 for (unsigned int c = 0; c < 2; ++c, ++l)
3817 new_lines[l] = cell->line(face_no)->child(c);
3818 Assert(l == 8, ExcInternalError());
3819
3820 new_lines[8]->set_bounding_object_indices(
3821 {new_vertices[6], new_vertices[8]});
3822 new_lines[9]->set_bounding_object_indices(
3823 {new_vertices[8], new_vertices[7]});
3824 new_lines[10]->set_bounding_object_indices(
3825 {new_vertices[4], new_vertices[8]});
3826 new_lines[11]->set_bounding_object_indices(
3827 {new_vertices[8], new_vertices[5]});
3828 }
3829 else if (ref_case == RefinementCase<dim>::cut_x)
3830 {
3831 // .-4-.-5-.
3832 // | | |
3833 // 0 6 1
3834 // | | |
3835 // .-2-.-3-.
3836 new_lines[0] = cell->line(0);
3837 new_lines[1] = cell->line(1);
3838 new_lines[2] = cell->line(2)->child(0);
3839 new_lines[3] = cell->line(2)->child(1);
3840 new_lines[4] = cell->line(3)->child(0);
3841 new_lines[5] = cell->line(3)->child(1);
3842 new_lines[6]->set_bounding_object_indices(
3843 {new_vertices[6], new_vertices[7]});
3844 }
3845 else
3846 {
3848 // .---5---.
3849 // 1 3
3850 // .---6---.
3851 // 0 2
3852 // .---4---.
3853 new_lines[0] = cell->line(0)->child(0);
3854 new_lines[1] = cell->line(0)->child(1);
3855 new_lines[2] = cell->line(1)->child(0);
3856 new_lines[3] = cell->line(1)->child(1);
3857 new_lines[4] = cell->line(2);
3858 new_lines[5] = cell->line(3);
3859 new_lines[6]->set_bounding_object_indices(
3860 {new_vertices[4], new_vertices[5]});
3861 }
3862
3863 for (unsigned int l = lmin; l < lmax; ++l)
3864 {
3865 new_lines[l]->set_used_flag();
3866 new_lines[l]->clear_user_flag();
3867 new_lines[l]->clear_user_data();
3868 new_lines[l]->clear_children();
3869 // interior line
3870 new_lines[l]->set_boundary_id_internal(
3872 new_lines[l]->set_manifold_id(cell->manifold_id());
3873 }
3874
3875 // Now add the four (two)
3876 // new cells!
3879 while (next_unused_cell->used() == true)
3880 ++next_unused_cell;
3881
3882 const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
3883 for (unsigned int i = 0; i < n_children; ++i)
3884 {
3885 Assert(
3886 next_unused_cell->used() == false,
3887 ExcMessage(
3888 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
3889 subcells[i] = next_unused_cell;
3890 ++next_unused_cell;
3891 if (i % 2 == 1 && i < n_children - 1)
3892 while (next_unused_cell->used() == true)
3893 ++next_unused_cell;
3894 }
3895
3896 if (ref_case == RefinementCase<dim>::cut_xy)
3897 {
3898 // children:
3899 // .--.--.
3900 // |2 . 3|
3901 // .--.--.
3902 // |0 | 1|
3903 // .--.--.
3904 // lines:
3905 // .-6-.-7-.
3906 // 1 9 3
3907 // .-10.11-.
3908 // 0 8 2
3909 // .-4-.-5-.
3910 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3911 new_lines[8]->index(),
3912 new_lines[4]->index(),
3913 new_lines[10]->index()});
3914 subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3915 new_lines[2]->index(),
3916 new_lines[5]->index(),
3917 new_lines[11]->index()});
3918 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3919 new_lines[9]->index(),
3920 new_lines[10]->index(),
3921 new_lines[6]->index()});
3922 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3923 new_lines[3]->index(),
3924 new_lines[11]->index(),
3925 new_lines[7]->index()});
3926 }
3927 else if (ref_case == RefinementCase<dim>::cut_x)
3928 {
3929 // children:
3930 // .--.--.
3931 // | . |
3932 // .0 . 1.
3933 // | | |
3934 // .--.--.
3935 // lines:
3936 // .-4-.-5-.
3937 // | | |
3938 // 0 6 1
3939 // | | |
3940 // .-2-.-3-.
3941 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3942 new_lines[6]->index(),
3943 new_lines[2]->index(),
3944 new_lines[4]->index()});
3945 subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3946 new_lines[1]->index(),
3947 new_lines[3]->index(),
3948 new_lines[5]->index()});
3949 }
3950 else
3951 {
3953 // children:
3954 // .-----.
3955 // | 1 |
3956 // .-----.
3957 // | 0 |
3958 // .-----.
3959 // lines:
3960 // .---5---.
3961 // 1 3
3962 // .---6---.
3963 // 0 2
3964 // .---4---.
3965 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3966 new_lines[2]->index(),
3967 new_lines[4]->index(),
3968 new_lines[6]->index()});
3969 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
3970 new_lines[3]->index(),
3971 new_lines[6]->index(),
3972 new_lines[5]->index()});
3973 }
3974
3975 types::subdomain_id subdomainid = cell->subdomain_id();
3976
3977 for (unsigned int i = 0; i < n_children; ++i)
3978 {
3979 subcells[i]->set_used_flag();
3980 subcells[i]->clear_refine_flag();
3981 subcells[i]->clear_user_flag();
3982 subcells[i]->clear_user_data();
3983 subcells[i]->clear_children();
3984 // inherit material
3985 // properties
3986 subcells[i]->set_material_id(cell->material_id());
3987 subcells[i]->set_manifold_id(cell->manifold_id());
3988 subcells[i]->set_subdomain_id(subdomainid);
3989
3990 if (i % 2 == 0)
3991 subcells[i]->set_parent(cell->index());
3992 }
3993
3994
3995
3996 // set child index for
3997 // even children children
3998 // i=0,2 (0)
3999 for (unsigned int i = 0; i < n_children / 2; ++i)
4000 cell->set_children(2 * i, subcells[2 * i]->index());
4001 // set the refine case
4002 cell->set_refinement_case(ref_case);
4003
4004 // note that the
4005 // refinement flag was
4006 // already cleared at the
4007 // beginning of this function
4008
4009 if (dim < spacedim)
4010 for (unsigned int c = 0; c < n_children; ++c)
4011 cell->child(c)->set_direction_flag(cell->direction_flag());
4012 }
4013
4014
4015
4016 template <int dim, int spacedim>
4019 const bool check_for_distorted_cells)
4020 {
4021 AssertDimension(dim, 2);
4022
4023 // Check whether a new level is needed. We have to check for
4024 // this on the highest level only
4025 for (const auto &cell : triangulation.active_cell_iterators_on_level(
4026 triangulation.levels.size() - 1))
4027 if (cell->refine_flag_set())
4028 {
4029 triangulation.levels.push_back(
4030 std::make_unique<
4032 break;
4033 }
4034
4036 triangulation.begin_line();
4037 line != triangulation.end_line();
4038 ++line)
4039 {
4040 line->clear_user_flag();
4041 line->clear_user_data();
4042 }
4043
4044 unsigned int n_single_lines = 0;
4045 unsigned int n_lines_in_pairs = 0;
4046 unsigned int needed_vertices = 0;
4047
4048 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4049 {
4050 // count number of flagged cells on this level and compute
4051 // how many new vertices and new lines will be needed
4052 unsigned int needed_cells = 0;
4053
4054 for (const auto &cell :
4055 triangulation.active_cell_iterators_on_level(level))
4056 if (cell->refine_flag_set())
4057 {
4058 if (cell->reference_cell() ==
4060 {
4061 needed_cells += 4;
4062 needed_vertices += 0;
4063 n_single_lines += 3;
4064 }
4065 else if (cell->reference_cell() ==
4067 {
4068 needed_cells += 4;
4069 needed_vertices += 1;
4070 n_single_lines += 4;
4071 }
4072 else
4073 {
4075 }
4076
4077 for (const auto line_no : cell->face_indices())
4078 {
4079 auto line = cell->line(line_no);
4080 if (line->has_children() == false)
4081 line->set_user_flag();
4082 }
4083 }
4084
4085
4086 const unsigned int used_cells =
4087 std::count(triangulation.levels[level + 1]->cells.used.begin(),
4088 triangulation.levels[level + 1]->cells.used.end(),
4089 true);
4090
4091
4092 reserve_space(*triangulation.levels[level + 1],
4093 used_cells + needed_cells,
4094 2,
4095 spacedim);
4096
4097 reserve_space(triangulation.levels[level + 1]->cells,
4098 needed_cells,
4099 0);
4100 }
4101
4102 for (auto line = triangulation.begin_line();
4103 line != triangulation.end_line();
4104 ++line)
4105 if (line->user_flag_set())
4106 {
4107 Assert(line->has_children() == false, ExcInternalError());
4108 n_lines_in_pairs += 2;
4109 needed_vertices += 1;
4110 }
4111
4112 reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4113
4114 needed_vertices += std::count(triangulation.vertices_used.begin(),
4115 triangulation.vertices_used.end(),
4116 true);
4117
4118 if (needed_vertices > triangulation.vertices.size())
4119 {
4120 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4121 triangulation.vertices_used.resize(needed_vertices, false);
4122 }
4123
4124 unsigned int next_unused_vertex = 0;
4125
4126 {
4128 line = triangulation.begin_active_line(),
4129 endl = triangulation.end_line();
4131 next_unused_line = triangulation.begin_raw_line();
4132
4133 for (; line != endl; ++line)
4134 if (line->user_flag_set())
4135 {
4136 // this line needs to be refined
4137
4138 // find the next unused vertex and set it
4139 // appropriately
4140 while (triangulation.vertices_used[next_unused_vertex] == true)
4141 ++next_unused_vertex;
4142 Assert(
4143 next_unused_vertex < triangulation.vertices.size(),
4144 ExcMessage(
4145 "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."));
4146 triangulation.vertices_used[next_unused_vertex] = true;
4147
4148 triangulation.vertices[next_unused_vertex] = line->center(true);
4149
4150 bool pair_found = false;
4151 (void)pair_found;
4152 for (; next_unused_line != endl; ++next_unused_line)
4153 if (!next_unused_line->used() &&
4154 !(++next_unused_line)->used())
4155 {
4156 --next_unused_line;
4157 pair_found = true;
4158 break;
4159 }
4160 Assert(pair_found, ExcInternalError());
4161
4162 line->set_children(0, next_unused_line->index());
4163
4165 children[2] = {next_unused_line, ++next_unused_line};
4166
4167 Assert(
4168 children[0]->used() == false,
4169 ExcMessage(
4170 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4171 Assert(
4172 children[1]->used() == false,
4173 ExcMessage(
4174 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4175
4176 children[0]->set_bounding_object_indices(
4177 {line->vertex_index(0), next_unused_vertex});
4178 children[1]->set_bounding_object_indices(
4179 {next_unused_vertex, line->vertex_index(1)});
4180
4181 children[0]->set_used_flag();
4182 children[1]->set_used_flag();
4183 children[0]->clear_children();
4184 children[1]->clear_children();
4185 children[0]->clear_user_data();
4186 children[1]->clear_user_data();
4187 children[0]->clear_user_flag();
4188 children[1]->clear_user_flag();
4189
4190
4191 children[0]->set_boundary_id_internal(line->boundary_id());
4192 children[1]->set_boundary_id_internal(line->boundary_id());
4193
4194 children[0]->set_manifold_id(line->manifold_id());
4195 children[1]->set_manifold_id(line->manifold_id());
4196
4197 line->clear_user_flag();
4198 }
4199 }
4200
4201 reserve_space(triangulation.faces->lines, 0, n_single_lines);
4202
4204 cells_with_distorted_children;
4205
4207 next_unused_line = triangulation.begin_raw_line();
4208
4209 const auto create_children = [](auto & triangulation,
4210 unsigned int &next_unused_vertex,
4211 auto & next_unused_line,
4212 auto & next_unused_cell,
4213 const auto & cell) {
4214 const auto ref_case = cell->refine_flag_set();
4215 cell->clear_refine_flag();
4216
4217 unsigned int n_new_vertices = 0;
4218
4219 if (cell->reference_cell() == ::ReferenceCells::Triangle)
4220 n_new_vertices = 6;
4221 else if (cell->reference_cell() ==
4223 n_new_vertices = 9;
4224 else
4226
4227 std::vector<int> new_vertices(n_new_vertices);
4228 for (unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
4229 ++vertex_no)
4230 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4231 for (unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
4232 if (cell->line(line_no)->has_children())
4233 new_vertices[cell->n_vertices() + line_no] =
4234 cell->line(line_no)->child(0)->vertex_index(1);
4235
4236 if (cell->reference_cell() == ::ReferenceCells::Quadrilateral)
4237 {
4238 while (triangulation.vertices_used[next_unused_vertex] == true)
4239 ++next_unused_vertex;
4240 Assert(
4241 next_unused_vertex < triangulation.vertices.size(),
4242 ExcMessage(
4243 "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."));
4244 triangulation.vertices_used[next_unused_vertex] = true;
4245
4246 new_vertices[8] = next_unused_vertex;
4247
4248 if (dim == spacedim)
4249 {
4250 triangulation.vertices[next_unused_vertex] =
4251 cell->center(true);
4252
4253 if (cell->user_flag_set())
4254 {
4255 cell->clear_user_flag();
4256 triangulation.vertices[next_unused_vertex] =
4257 cell->center(true, true);
4258 }
4259 }
4260 else
4261 {
4262 cell->clear_user_flag();
4263
4264 triangulation.vertices[next_unused_vertex] =
4265 cell->center(true, true);
4266 }
4267 }
4268
4269 std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
4270 12>
4271 new_lines;
4272 unsigned int lmin = 0;
4273 unsigned int lmax = 0;
4274
4275 if (cell->reference_cell() == ::ReferenceCells::Triangle)
4276 {
4277 lmin = 6;
4278 lmax = 9;
4279 }
4280 else if (cell->reference_cell() ==
4282 {
4283 lmin = 8;
4284 lmax = 12;
4285 }
4286 else
4287 {
4289 }
4290
4291 for (unsigned int l = lmin; l < lmax; ++l)
4292 {
4293 while (next_unused_line->used() == true)
4294 ++next_unused_line;
4295 new_lines[l] = next_unused_line;
4296 ++next_unused_line;
4297
4298 Assert(
4299 new_lines[l]->used() == false,
4300 ExcMessage(
4301 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4302 }
4303
4304 if (true)
4305 {
4306 if (cell->reference_cell() == ::ReferenceCells::Triangle)
4307 {
4308 // add lines in the right order [TODO: clean up]
4309 const auto ref = [&](const unsigned int face_no,
4310 const unsigned int vertex_no) {
4311 if (cell->line(face_no)->child(0)->vertex_index(0) ==
4312 static_cast<unsigned int>(new_vertices[vertex_no]) ||
4313 cell->line(face_no)->child(0)->vertex_index(1) ==
4314 static_cast<unsigned int>(new_vertices[vertex_no]))
4315 {
4316 new_lines[2 * face_no + 0] =
4317 cell->line(face_no)->child(0);
4318 new_lines[2 * face_no + 1] =
4319 cell->line(face_no)->child(1);
4320 }
4321 else
4322 {
4323 new_lines[2 * face_no + 0] =
4324 cell->line(face_no)->child(1);
4325 new_lines[2 * face_no + 1] =
4326 cell->line(face_no)->child(0);
4327 }
4328 };
4329
4330 ref(0, 0);
4331 ref(1, 1);
4332 ref(2, 2);
4333
4334 new_lines[6]->set_bounding_object_indices(
4335 {new_vertices[3], new_vertices[4]});
4336 new_lines[7]->set_bounding_object_indices(
4337 {new_vertices[4], new_vertices[5]});
4338 new_lines[8]->set_bounding_object_indices(
4339 {new_vertices[5], new_vertices[3]});
4340 }
4341 else if (cell->reference_cell() ==
4343 {
4344 unsigned int l = 0;
4345 for (const unsigned int face_no : cell->face_indices())
4346 for (unsigned int c = 0; c < 2; ++c, ++l)
4347 new_lines[l] = cell->line(face_no)->child(c);
4348
4349 new_lines[8]->set_bounding_object_indices(
4350 {new_vertices[6], new_vertices[8]});
4351 new_lines[9]->set_bounding_object_indices(
4352 {new_vertices[8], new_vertices[7]});
4353 new_lines[10]->set_bounding_object_indices(
4354 {new_vertices[4], new_vertices[8]});
4355 new_lines[11]->set_bounding_object_indices(
4356 {new_vertices[8], new_vertices[5]});
4357 }
4358 else
4359 {
4361 }
4362 }
4363
4364
4365 for (unsigned int l = lmin; l < lmax; ++l)
4366 {
4367 new_lines[l]->set_used_flag();
4368 new_lines[l]->clear_user_flag();
4369 new_lines[l]->clear_user_data();
4370 new_lines[l]->clear_children();
4371 // interior line
4372 new_lines[l]->set_boundary_id_internal(
4374 new_lines[l]->set_manifold_id(cell->manifold_id());
4375 }
4376
4379 while (next_unused_cell->used() == true)
4380 ++next_unused_cell;
4381
4382 unsigned int n_children = 0;
4383
4384 if (cell->reference_cell() == ::ReferenceCells::Triangle)
4385 n_children = 4;
4386 else if (cell->reference_cell() ==
4388 n_children = 4;
4389 else
4391
4392 for (unsigned int i = 0; i < n_children; ++i)
4393 {
4394 Assert(
4395 next_unused_cell->used() == false,
4396 ExcMessage(
4397 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4398 subcells[i] = next_unused_cell;
4399 ++next_unused_cell;
4400 if (i % 2 == 1 && i < n_children - 1)
4401 while (next_unused_cell->used() == true)
4402 ++next_unused_cell;
4403 }
4404
4405 if ((dim == 2) &&
4406 (cell->reference_cell() == ::ReferenceCells::Triangle))
4407 {
4408 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4409 new_lines[8]->index(),
4410 new_lines[5]->index()});
4411 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4412 new_lines[2]->index(),
4413 new_lines[6]->index()});
4414 subcells[2]->set_bounding_object_indices({new_lines[7]->index(),
4415 new_lines[3]->index(),
4416 new_lines[4]->index()});
4417 subcells[3]->set_bounding_object_indices({new_lines[6]->index(),
4418 new_lines[7]->index(),
4419 new_lines[8]->index()});
4420
4421 // subcell 0
4422
4423 const auto ref = [&](const unsigned int line_no,
4424 const unsigned int vertex_no,
4425 const unsigned int subcell_no,
4426 const unsigned int subcell_line_no) {
4427 if (new_lines[line_no]->vertex_index(1) !=
4428 static_cast<unsigned int>(new_vertices[vertex_no]))
4429 triangulation.levels[subcells[subcell_no]->level()]
4430 ->face_orientations[subcells[subcell_no]->index() *
4432 subcell_line_no] = 0;
4433 };
4434
4435 ref(0, 3, 0, 0);
4436 ref(8, 5, 0, 1);
4437 ref(5, 0, 0, 2);
4438
4439 ref(1, 1, 1, 0);
4440 ref(2, 4, 1, 1);
4441 ref(6, 3, 1, 2);
4442
4443 ref(7, 4, 2, 0);
4444 ref(3, 2, 2, 1);
4445 ref(4, 5, 2, 2);
4446
4447 ref(6, 4, 3, 0);
4448 ref(7, 5, 3, 1);
4449 ref(8, 3, 3, 2);
4450
4451 // triangulation.levels[subcells[1]->level()]->face_orientations[subcells[1]->index()
4452 // * GeometryInfo<2>::faces_per_cell + 2] = 0;
4453 // triangulation.levels[subcells[2]->level()]->face_orientations[subcells[2]->index()
4454 // * GeometryInfo<2>::faces_per_cell + 0] = 0;
4455 }
4456 else if ((dim == 2) && (cell->reference_cell() ==
4458 {
4459 subcells[0]->set_bounding_object_indices(
4460 {new_lines[0]->index(),
4461 new_lines[8]->index(),
4462 new_lines[4]->index(),
4463 new_lines[10]->index()});
4464 subcells[1]->set_bounding_object_indices(
4465 {new_lines[8]->index(),
4466 new_lines[2]->index(),
4467 new_lines[5]->index(),
4468 new_lines[11]->index()});
4469 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
4470 new_lines[9]->index(),
4471 new_lines[10]->index(),
4472 new_lines[6]->index()});
4473 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
4474 new_lines[3]->index(),
4475 new_lines[11]->index(),
4476 new_lines[7]->index()});
4477 }
4478 else
4479 {
4481 }
4482
4483 types::subdomain_id subdomainid = cell->subdomain_id();
4484
4485 for (unsigned int i = 0; i < n_children; ++i)
4486 {
4487 subcells[i]->set_used_flag();
4488 subcells[i]->clear_refine_flag();
4489 subcells[i]->clear_user_flag();
4490 subcells[i]->clear_user_data();
4491 subcells[i]->clear_children();
4492 // inherit material
4493 // properties
4494 subcells[i]->set_material_id(cell->material_id());
4495 subcells[i]->set_manifold_id(cell->manifold_id());
4496 subcells[i]->set_subdomain_id(subdomainid);
4497
4498 // TODO: here we assume that all children have the same reference
4499 // cell type as the parent! This is justified for 2D.
4500 triangulation.levels[subcells[i]->level()]
4501 ->reference_cell[subcells[i]->index()] = cell->reference_cell();
4502
4503 if (i % 2 == 0)
4504 subcells[i]->set_parent(cell->index());
4505 }
4506
4507 for (unsigned int i = 0; i < n_children / 2; ++i)
4508 cell->set_children(2 * i, subcells[2 * i]->index());
4509
4510 cell->set_refinement_case(ref_case);
4511
4512 if (dim < spacedim)
4513 for (unsigned int c = 0; c < n_children; ++c)
4514 cell->child(c)->set_direction_flag(cell->direction_flag());
4515 };
4516
4517 for (int level = 0;
4518 level < static_cast<int>(triangulation.levels.size()) - 1;
4519 ++level)
4520 {
4522 next_unused_cell = triangulation.begin_raw(level + 1);
4523
4524 for (const auto &cell :
4525 triangulation.active_cell_iterators_on_level(level))
4526 if (cell->refine_flag_set())
4527 {
4528 if (cell->at_boundary())
4529 cell->set_user_flag();
4530
4532 next_unused_vertex,
4533 next_unused_line,
4534 next_unused_cell,
4535 cell);
4536
4537 if (cell->reference_cell() ==
4539 check_for_distorted_cells &&
4540 has_distorted_children<dim, spacedim>(cell))
4541 cells_with_distorted_children.distorted_cells.push_back(
4542 cell);
4543
4544 triangulation.signals.post_refinement_on_cell(cell);
4545 }
4546 }
4547
4548 return cells_with_distorted_children;
4549 }
4550
4551
4552
4557 template <int spacedim>
4560 const bool /*check_for_distorted_cells*/)
4561 {
4562 const unsigned int dim = 1;
4563
4564 // Check whether a new level is needed. We have to check for
4565 // this on the highest level only
4566 for (const auto &cell : triangulation.active_cell_iterators_on_level(
4567 triangulation.levels.size() - 1))
4568 if (cell->refine_flag_set())
4569 {
4570 triangulation.levels.push_back(
4571 std::make_unique<
4573 break;
4574 }
4575
4576
4577 // check how much space is needed on every level we need not
4578 // check the highest level since either - on the highest level
4579 // no cells are flagged for refinement - there are, but
4580 // prepare_refinement added another empty level
4581 unsigned int needed_vertices = 0;
4582 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4583 {
4584 // count number of flagged
4585 // cells on this level
4586 unsigned int flagged_cells = 0;
4587
4588 for (const auto &acell :
4589 triangulation.active_cell_iterators_on_level(level))
4590 if (acell->refine_flag_set())
4591 ++flagged_cells;
4592
4593 // count number of used cells
4594 // on the next higher level
4595 const unsigned int used_cells =
4596 std::count(triangulation.levels[level + 1]->cells.used.begin(),
4597 triangulation.levels[level + 1]->cells.used.end(),
4598 true);
4599
4600 // reserve space for the used_cells cells already existing
4601 // on the next higher level as well as for the
4602 // 2*flagged_cells that will be created on that level
4603 reserve_space(*triangulation.levels[level + 1],
4605 flagged_cells,
4606 1,
4607 spacedim);
4608 // reserve space for 2*flagged_cells new lines on the next
4609 // higher level
4610 reserve_space(triangulation.levels[level + 1]->cells,
4612 flagged_cells,
4613 0);
4614
4615 needed_vertices += flagged_cells;
4616 }
4617
4618 // add to needed vertices how many
4619 // vertices are already in use
4620 needed_vertices += std::count(triangulation.vertices_used.begin(),
4621 triangulation.vertices_used.end(),
4622 true);
4623 // if we need more vertices: create them, if not: leave the
4624 // array as is, since shrinking is not really possible because
4625 // some of the vertices at the end may be in use
4626 if (needed_vertices > triangulation.vertices.size())
4627 {
4628 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4629 triangulation.vertices_used.resize(needed_vertices, false);
4630 }
4631
4632
4633 // Do REFINEMENT on every level; exclude highest level as
4634 // above
4635
4636 // index of next unused vertex
4637 unsigned int next_unused_vertex = 0;
4638
4639 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4640 {
4642 next_unused_cell = triangulation.begin_raw(level + 1);
4643
4644 for (const auto &cell :
4645 triangulation.active_cell_iterators_on_level(level))
4646 if (cell->refine_flag_set())
4647 {
4648 // clear refinement flag
4649 cell->clear_refine_flag();
4650
4651 // search for next unused
4652 // vertex
4653 while (triangulation.vertices_used[next_unused_vertex] ==
4654 true)
4655 ++next_unused_vertex;
4656 Assert(
4657 next_unused_vertex < triangulation.vertices.size(),
4658 ExcMessage(
4659 "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."));
4660
4661 // Now we always ask the cell itself where to put
4662 // the new point. The cell in turn will query the
4663 // manifold object internally.
4664 triangulation.vertices[next_unused_vertex] =
4665 cell->center(true);
4666
4667 triangulation.vertices_used[next_unused_vertex] = true;
4668
4669 // search for next two unused cell (++ takes care of
4670 // the end of the vector)
4672 first_child,
4673 second_child;
4674 while (next_unused_cell->used() == true)
4675 ++next_unused_cell;
4676 first_child = next_unused_cell;
4677 first_child->set_used_flag();
4678 first_child->clear_user_data();
4679 ++next_unused_cell;
4680 Assert(
4681 next_unused_cell->used() == false,
4682 ExcMessage(
4683 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
4684 second_child = next_unused_cell;
4685 second_child->set_used_flag();
4686 second_child->clear_user_data();
4687
4688 types::subdomain_id subdomainid = cell->subdomain_id();
4689
4690 // insert first child
4691 cell->set_children(0, first_child->index());
4692 first_child->clear_children();
4693 first_child->set_bounding_object_indices(
4694 {cell->vertex_index(0), next_unused_vertex});
4695 first_child->set_material_id(cell->material_id());
4696 first_child->set_manifold_id(cell->manifold_id());
4697 first_child->set_subdomain_id(subdomainid);
4698 first_child->set_direction_flag(cell->direction_flag());
4699
4700 first_child->set_parent(cell->index());
4701
4702 // Set manifold id of the right face. Only do this
4703 // on the first child.
4704 first_child->face(1)->set_manifold_id(cell->manifold_id());
4705
4706 // reset neighborship info (refer to
4707 // internal::TriangulationImplementation::TriaLevel<0> for
4708 // details)
4709 first_child->set_neighbor(1, second_child);
4710 if (cell->neighbor(0).state() != IteratorState::valid)
4711 first_child->set_neighbor(0, cell->neighbor(0));
4712 else if (cell->neighbor(0)->is_active())
4713 {
4714 // since the neighbors level is always <=level,
4715 // if the cell is active, then there are no
4716 // cells to the left which may want to know
4717 // about this new child cell.
4718 Assert(cell->neighbor(0)->level() <= cell->level(),
4720 first_child->set_neighbor(0, cell->neighbor(0));
4721 }
4722 else
4723 // left neighbor is refined
4724 {
4725 // set neighbor to cell on same level
4726 const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4727 first_child->set_neighbor(0,
4728 cell->neighbor(0)->child(nbnb));
4729
4730 // reset neighbor info of all right descendant
4731 // of the left neighbor of cell
4733 left_neighbor = cell->neighbor(0);
4734 while (left_neighbor->has_children())
4735 {
4736 left_neighbor = left_neighbor->child(nbnb);
4737 left_neighbor->set_neighbor(nbnb, first_child);
4738 }
4739 }
4740
4741 // insert second child
4742 second_child->clear_children();
4743 second_child->set_bounding_object_indices(
4744 {next_unused_vertex, cell->vertex_index(1)});
4745 second_child->set_neighbor(0, first_child);
4746 second_child->set_material_id(cell->material_id());
4747 second_child->set_manifold_id(cell->manifold_id());
4748 second_child->set_subdomain_id(subdomainid);
4749 second_child->set_direction_flag(cell->direction_flag());
4750
4751 if (cell->neighbor(1).state() != IteratorState::valid)
4752 second_child->set_neighbor(1, cell->neighbor(1));
4753 else if (cell->neighbor(1)->is_active())
4754 {
4755 Assert(cell->neighbor(1)->level() <= cell->level(),
4757 second_child->set_neighbor(1, cell->neighbor(1));
4758 }
4759 else
4760 // right neighbor is refined same as above
4761 {
4762 const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4763 second_child->set_neighbor(
4764 1, cell->neighbor(1)->child(nbnb));
4765
4767 right_neighbor = cell->neighbor(1);
4768 while (right_neighbor->has_children())
4769 {
4770 right_neighbor = right_neighbor->child(nbnb);
4771 right_neighbor->set_neighbor(nbnb, second_child);
4772 }
4773 }
4774 // inform all listeners that cell refinement is done
4775 triangulation.signals.post_refinement_on_cell(cell);
4776 }
4777 }
4778
4779 // in 1d, we can not have distorted children unless the parent
4780 // was already distorted (that is because we don't use
4781 // boundary information for 1d triangulations). so return an
4782 // empty list
4784 }
4785
4786
4791 template <int spacedim>
4794 const bool check_for_distorted_cells)
4795 {
4796 const unsigned int dim = 2;
4797
4798
4799 // First check whether we can get away with isotropic refinement, or
4800 // whether we need to run through the full anisotropic algorithm
4801 {
4802 bool do_isotropic_refinement = true;
4803 for (const auto &cell : triangulation.active_cell_iterators())
4804 if (cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
4805 cell->refine_flag_set() == RefinementCase<dim>::cut_y)
4806 {
4807 do_isotropic_refinement = false;
4808 break;
4809 }
4810
4811 if (do_isotropic_refinement)
4813 check_for_distorted_cells);
4814 }
4815
4816 // Check whether a new level is needed. We have to check for
4817 // this on the highest level only
4818 for (const auto &cell : triangulation.active_cell_iterators_on_level(
4819 triangulation.levels.size() - 1))
4820 if (cell->refine_flag_set())
4821 {
4822 triangulation.levels.push_back(
4823 std::make_unique<
4825 break;
4826 }
4827
4828 // TODO[WB]: we clear user flags and pointers of lines; we're going
4829 // to use them to flag which lines need refinement
4831 triangulation.begin_line();
4832 line != triangulation.end_line();
4833 ++line)
4834 {
4835 line->clear_user_flag();
4836 line->clear_user_data();
4837 }
4838 // running over all cells and lines count the number
4839 // n_single_lines of lines which can be stored as single
4840 // lines, e.g. inner lines
4841 unsigned int n_single_lines = 0;
4842
4843 // New lines to be created: number lines which are stored in
4844 // pairs (the children of lines must be stored in pairs)
4845 unsigned int n_lines_in_pairs = 0;
4846
4847 // check how much space is needed on every level we need not
4848 // check the highest level since either - on the highest level
4849 // no cells are flagged for refinement - there are, but
4850 // prepare_refinement added another empty level
4851 unsigned int needed_vertices = 0;
4852 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4853 {
4854 // count number of flagged cells on this level and compute
4855 // how many new vertices and new lines will be needed
4856 unsigned int needed_cells = 0;
4857
4858 for (const auto &cell :
4859 triangulation.active_cell_iterators_on_level(level))
4860 if (cell->refine_flag_set())
4861 {
4862 if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
4863 {
4864 needed_cells += 4;
4865
4866 // new vertex at center of cell is needed in any
4867 // case
4868 ++needed_vertices;
4869
4870 // the four inner lines can be stored as singles
4871 n_single_lines += 4;
4872 }
4873 else // cut_x || cut_y
4874 {
4875 // set the flag showing that anisotropic
4876 // refinement is used for at least one cell
4877 triangulation.anisotropic_refinement = true;
4878
4879 needed_cells += 2;
4880 // no vertex at center
4881
4882 // the inner line can be stored as single
4883 n_single_lines += 1;
4884 }
4885
4886 // mark all faces (lines) for refinement; checking
4887 // locally whether the neighbor would also like to
4888 // refine them is rather difficult for lines so we
4889 // only flag them and after visiting all cells, we
4890 // decide which lines need refinement;
4891 for (const unsigned int line_no :
4893 {
4895 cell->refine_flag_set(), line_no) ==
4897 {
4899 line = cell->line(line_no);
4900 if (line->has_children() == false)
4901 line->set_user_flag();
4902 }
4903 }
4904 }
4905
4906
4907 // count number of used cells on the next higher level
4908 const unsigned int used_cells =
4909 std::count(triangulation.levels[level + 1]->cells.used.begin(),
4910 triangulation.levels[level + 1]->cells.used.end(),
4911 true);
4912
4913
4914 // reserve space for the used_cells cells already existing
4915 // on the next higher level as well as for the
4916 // needed_cells that will be created on that level
4917 reserve_space(*triangulation.levels[level + 1],
4918 used_cells + needed_cells,
4919 2,
4920 spacedim);
4921
4922 // reserve space for needed_cells new quads on the next
4923 // higher level
4924 reserve_space(triangulation.levels[level + 1]->cells,
4925 needed_cells,
4926 0);
4927 }
4928
4929 // now count the lines which were flagged for refinement
4931 triangulation.begin_line();
4932 line != triangulation.end_line();
4933 ++line)
4934 if (line->user_flag_set())
4935 {
4936 Assert(line->has_children() == false, ExcInternalError());
4937 n_lines_in_pairs += 2;
4938 needed_vertices += 1;
4939 }
4940 // reserve space for n_lines_in_pairs new lines. note, that
4941 // we can't reserve space for the single lines here as well,
4942 // as all the space reserved for lines in pairs would be
4943 // counted as unused and we would end up with too little space
4944 // to store all lines. memory reservation for n_single_lines
4945 // can only be done AFTER we refined the lines of the current
4946 // cells
4947 reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4948
4949 // add to needed vertices how many vertices are already in use
4950 needed_vertices += std::count(triangulation.vertices_used.begin(),
4951 triangulation.vertices_used.end(),
4952 true);
4953 // if we need more vertices: create them, if not: leave the
4954 // array as is, since shrinking is not really possible because
4955 // some of the vertices at the end may be in use
4956 if (needed_vertices > triangulation.vertices.size())
4957 {
4958 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4959 triangulation.vertices_used.resize(needed_vertices, false);
4960 }
4961
4962
4963 // Do REFINEMENT on every level; exclude highest level as
4964 // above
4965
4966 // index of next unused vertex
4967 unsigned int next_unused_vertex = 0;
4968
4969 // first the refinement of lines. children are stored
4970 // pairwise
4971 {
4972 // only active objects can be refined further
4974 line = triangulation.begin_active_line(),
4975 endl = triangulation.end_line();
4977 next_unused_line = triangulation.begin_raw_line();
4978
4979 for (; line != endl; ++line)
4980 if (line->user_flag_set())
4981 {
4982 // this line needs to be refined
4983
4984 // find the next unused vertex and set it
4985 // appropriately
4986 while (triangulation.vertices_used[next_unused_vertex] == true)
4987 ++next_unused_vertex;
4988 Assert(
4989 next_unused_vertex < triangulation.vertices.size(),
4990 ExcMessage(
4991 "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."));
4992 triangulation.vertices_used[next_unused_vertex] = true;
4993
4994 triangulation.vertices[next_unused_vertex] = line->center(true);
4995
4996 // now that we created the right point, make up the
4997 // two child lines. To this end, find a pair of
4998 // unused lines
4999 bool pair_found = false;
5000 (void)pair_found;
5001 for (; next_unused_line != endl; ++next_unused_line)
5002 if (!next_unused_line->used() &&
5003 !(++next_unused_line)->used())
5004 {
5005 // go back to the first of the two unused
5006 // lines
5007 --next_unused_line;
5008 pair_found = true;
5009 break;
5010 }
5011 Assert(pair_found, ExcInternalError());
5012
5013 // there are now two consecutive unused lines, such
5014 // that the children of a line will be consecutive.
5015 // then set the child pointer of the present line
5016 line->set_children(0, next_unused_line->index());
5017
5018 // set the two new lines
5020 children[2] = {next_unused_line, ++next_unused_line};
5021 // some tests; if any of the iterators should be
5022 // invalid, then already dereferencing will fail
5023 Assert(
5024 children[0]->used() == false,
5025 ExcMessage(
5026 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5027 Assert(
5028 children[1]->used() == false,
5029 ExcMessage(
5030 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5031
5032 children[0]->set_bounding_object_indices(
5033 {line->vertex_index(0), next_unused_vertex});
5034 children[1]->set_bounding_object_indices(
5035 {next_unused_vertex, line->vertex_index(1)});
5036
5037 children[0]->set_used_flag();
5038 children[1]->set_used_flag();
5039 children[0]->clear_children();
5040 children[1]->clear_children();
5041 children[0]->clear_user_data();
5042 children[1]->clear_user_data();
5043 children[0]->clear_user_flag();
5044 children[1]->clear_user_flag();
5045
5046
5047 children[0]->set_boundary_id_internal(line->boundary_id());
5048 children[1]->set_boundary_id_internal(line->boundary_id());
5049
5050 children[0]->set_manifold_id(line->manifold_id());
5051 children[1]->set_manifold_id(line->manifold_id());
5052
5053 // finally clear flag indicating the need for
5054 // refinement
5055 line->clear_user_flag();
5056 }
5057 }
5058
5059
5060 // Now set up the new cells
5061
5062 // reserve space for inner lines (can be stored as single
5063 // lines)
5064 reserve_space(triangulation.faces->lines, 0, n_single_lines);
5065
5067 cells_with_distorted_children;
5068
5069 // reset next_unused_line, as now also single empty places in
5070 // the vector can be used
5072 next_unused_line = triangulation.begin_raw_line();
5073
5074 for (int level = 0;
5075 level < static_cast<int>(triangulation.levels.size()) - 1;
5076 ++level)
5077 {
5079 next_unused_cell = triangulation.begin_raw(level + 1);
5080
5081 for (const auto &cell :
5082 triangulation.active_cell_iterators_on_level(level))
5083 if (cell->refine_flag_set())
5084 {
5085 // set the user flag to indicate, that at least one
5086 // line is at the boundary
5087
5088 // TODO[Tobias Leicht] find a better place to set
5089 // this flag, so that we do not need so much time to
5090 // check each cell here
5091 if (cell->at_boundary())
5092 cell->set_user_flag();
5093
5094 // actually set up the children and update neighbor
5095 // information
5097 next_unused_vertex,
5098 next_unused_line,
5099 next_unused_cell,
5100 cell);
5101
5102 if (check_for_distorted_cells &&
5103 has_distorted_children<dim, spacedim>(cell))
5104 cells_with_distorted_children.distorted_cells.push_back(
5105 cell);
5106 // inform all listeners that cell refinement is done
5107 triangulation.signals.post_refinement_on_cell(cell);
5108 }
5109 }
5110
5111 return cells_with_distorted_children;
5112 }
5113
5114
5119 template <int spacedim>
5122 const bool check_for_distorted_cells)
5123 {
5124 const unsigned int dim = 3;
5125
5126 // this function probably also works for spacedim>3 but it
5127 // isn't tested. it will probably be necessary to pull new
5128 // vertices onto the manifold just as we do for the other
5129 // functions above.
5130 Assert(spacedim == 3, ExcNotImplemented());
5131
5132 // Check whether a new level is needed. We have to check for
5133 // this on the highest level only
5134 for (const auto &cell : triangulation.active_cell_iterators_on_level(
5135 triangulation.levels.size() - 1))
5136 if (cell->refine_flag_set())
5137 {
5138 triangulation.levels.push_back(
5139 std::make_unique<
5141 break;
5142 }
5143
5144
5145 // first clear user flags for quads and lines; we're going to
5146 // use them to flag which lines and quads need refinement
5147 triangulation.faces->quads.clear_user_data();
5148
5150 triangulation.begin_line();
5151 line != triangulation.end_line();
5152 ++line)
5153 line->clear_user_flag();
5155 triangulation.begin_quad();
5156 quad != triangulation.end_quad();
5157 ++quad)
5158 quad->clear_user_flag();
5159
5160 // create an array of face refine cases. User indices of faces
5161 // will be set to values corresponding with indices in this
5162 // array.
5163 const RefinementCase<dim - 1> face_refinement_cases[4] = {
5164 RefinementCase<dim - 1>::no_refinement,
5165 RefinementCase<dim - 1>::cut_x,
5166 RefinementCase<dim - 1>::cut_y,
5167 RefinementCase<dim - 1>::cut_xy};
5168
5169 // check how much space is needed on every level we need not
5170 // check the highest level since either
5171 // - on the highest level no cells are flagged for refinement
5172 // - there are, but prepare_refinement added another empty
5173 // level which then is the highest level
5174
5175 // variables to hold the number of newly to be created
5176 // vertices, lines and quads. as these are stored globally,
5177 // declare them outside the loop over al levels. we need lines
5178 // and quads in pairs for refinement of old ones and lines and
5179 // quads, that can be stored as single ones, as they are newly
5180 // created in the inside of an existing cell
5181 unsigned int needed_vertices = 0;
5182 unsigned int needed_lines_single = 0;
5183 unsigned int needed_quads_single = 0;
5184 unsigned int needed_lines_pair = 0;
5185 unsigned int needed_quads_pair = 0;
5186 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5187 {
5188 // count number of flagged cells on this level and compute
5189 // how many new vertices and new lines will be needed
5190 unsigned int new_cells = 0;
5191
5192 for (const auto &acell :
5193 triangulation.active_cell_iterators_on_level(level))
5194 if (acell->refine_flag_set())
5195 {
5196 RefinementCase<dim> ref_case = acell->refine_flag_set();
5197
5198 // now for interior vertices, lines and quads, which
5199 // are needed in any case
5200 if (ref_case == RefinementCase<dim>::cut_x ||
5201 ref_case == RefinementCase<dim>::cut_y ||
5202 ref_case == RefinementCase<dim>::cut_z)
5203 {
5204 ++needed_quads_single;
5205 new_cells += 2;
5206 triangulation.anisotropic_refinement = true;
5207 }
5208 else if (ref_case == RefinementCase<dim>::cut_xy ||
5209 ref_case == RefinementCase<dim>::cut_xz ||
5210 ref_case == RefinementCase<dim>::cut_yz)
5211 {
5212 ++needed_lines_single;
5213 needed_quads_single += 4;
5214 new_cells += 4;
5215 triangulation.anisotropic_refinement = true;
5216 }
5217 else if (ref_case == RefinementCase<dim>::cut_xyz)
5218 {
5219 ++needed_vertices;
5220 needed_lines_single += 6;
5221 needed_quads_single += 12;
5222 new_cells += 8;
5223 }
5224 else
5225 {
5226 // we should never get here
5227 Assert(false, ExcInternalError());
5228 }
5229
5230 // mark all faces for refinement; checking locally
5231 // if and how the neighbor would like to refine
5232 // these is difficult so we only flag them and after
5233 // visiting all cells, we decide which faces need
5234 // which refinement;
5235 for (const unsigned int face :
5237 {
5239 aface = acell->face(face);
5240 // get the RefineCase this faces has for the
5241 // given RefineCase of the cell
5242 RefinementCase<dim - 1> face_ref_case =
5244 ref_case,
5245 face,
5246 acell->face_orientation(face),
5247 acell->face_flip(face),
5248 acell->face_rotation(face));
5249 // only do something, if this face has to be
5250 // refined
5251 if (face_ref_case)
5252 {
5253 if (face_ref_case ==
5255 {
5256 if (aface->n_active_descendants() < 4)
5257 // we use user_flags to denote needed
5258 // isotropic refinement
5259 aface->set_user_flag();
5260 }
5261 else if (aface->refinement_case() != face_ref_case)
5262 // we use user_indices to denote needed
5263 // anisotropic refinement. note, that we
5264 // can have at most one anisotropic
5265 // refinement case for this face, as
5266 // otherwise prepare_refinement() would
5267 // have changed one of the cells to yield
5268 // isotropic refinement at this
5269 // face. therefore we set the user_index
5270 // uniquely
5271 {
5272 Assert(aface->refinement_case() ==
5274 dim - 1>::isotropic_refinement ||
5275 aface->refinement_case() ==
5278 aface->set_user_index(face_ref_case);
5279 }
5280 }
5281 } // for all faces
5282
5283 // flag all lines, that have to be refined
5284 for (unsigned int line = 0;
5285 line < GeometryInfo<dim>::lines_per_cell;
5286 ++line)
5288 line) &&
5289 !acell->line(line)->has_children())
5290 acell->line(line)->set_user_flag();
5291
5292 } // if refine_flag set and for all cells on this level
5293
5294
5295 // count number of used cells on the next higher level
5296 const unsigned int used_cells =
5297 std::count(triangulation.levels[level + 1]->cells.used.begin(),
5298 triangulation.levels[level + 1]->cells.used.end(),
5299 true);
5300
5301
5302 // reserve space for the used_cells cells already existing
5303 // on the next higher level as well as for the
5304 // 8*flagged_cells that will be created on that level
5305 reserve_space(*triangulation.levels[level + 1],
5306 used_cells + new_cells,
5307 3,
5308 spacedim);
5309 // reserve space for 8*flagged_cells new hexes on the next
5310 // higher level
5311 reserve_space(triangulation.levels[level + 1]->cells, new_cells);
5312 } // for all levels
5313 // now count the quads and lines which were flagged for
5314 // refinement
5316 triangulation.begin_quad();
5317 quad != triangulation.end_quad();
5318 ++quad)
5319 {
5320 if (quad->user_flag_set())
5321 {
5322 // isotropic refinement: 1 interior vertex, 4 quads
5323 // and 4 interior lines. we store the interior lines
5324 // in pairs in case the face is already or will be
5325 // refined anisotropically
5326 needed_quads_pair += 4;
5327 needed_lines_pair += 4;
5328 needed_vertices += 1;
5329 }
5330 if (quad->user_index())
5331 {
5332 // anisotropic refinement: 1 interior
5333 // line and two quads
5334 needed_quads_pair += 2;
5335 needed_lines_single += 1;
5336 // there is a kind of complicated situation here which
5337 // requires our attention. if the quad is refined
5338 // isotropcally, two of the interior lines will get a
5339 // new mother line - the interior line of our
5340 // anisotropically refined quad. if those two lines
5341 // are not consecutive, we cannot do so and have to
5342 // replace them by two lines that are consecutive. we
5343 // try to avoid that situation, but it may happen
5344 // nevertheless through repeated refinement and
5345 // coarsening. thus we have to check here, as we will
5346 // need some additional space to store those new lines
5347 // in case we need them...
5348 if (quad->has_children())
5349 {
5350 Assert(quad->refinement_case() ==
5353 if ((face_refinement_cases[quad->user_index()] ==
5355 (quad->child(0)->line_index(1) + 1 !=
5356 quad->child(2)->line_index(1))) ||
5357 (face_refinement_cases[quad->user_index()] ==
5359 (quad->child(0)->line_index(3) + 1 !=
5360 quad->child(1)->line_index(3))))
5361 needed_lines_pair += 2;
5362 }
5363 }
5364 }
5365
5367 triangulation.begin_line();
5368 line != triangulation.end_line();
5369 ++line)
5370 if (line->user_flag_set())
5371 {
5372 needed_lines_pair += 2;
5373 needed_vertices += 1;
5374 }
5375
5376 // reserve space for needed_lines new lines stored in pairs
5377 reserve_space(triangulation.faces->lines,
5378 needed_lines_pair,
5379 needed_lines_single);
5380 // reserve space for needed_quads new quads stored in pairs
5382 needed_quads_pair,
5383 needed_quads_single);
5384 reserve_space(triangulation.faces->quads,
5385 needed_quads_pair,
5386 needed_quads_single);
5387
5388
5389 // add to needed vertices how many vertices are already in use
5390 needed_vertices += std::count(triangulation.vertices_used.begin(),
5391 triangulation.vertices_used.end(),
5392 true);
5393 // if we need more vertices: create them, if not: leave the
5394 // array as is, since shrinking is not really possible because
5395 // some of the vertices at the end may be in use
5396 if (needed_vertices > triangulation.vertices.size())
5397 {
5398 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5399 triangulation.vertices_used.resize(needed_vertices, false);
5400 }
5401
5402
5404 // Before we start with the actual refinement, we do some
5405 // sanity checks if in debug mode. especially, we try to catch
5406 // the notorious problem with lines being twice refined,
5407 // i.e. there are cells adjacent at one line ("around the
5408 // edge", but not at a face), with two cells differing by more
5409 // than one refinement level
5410 //
5411 // this check is very simple to implement here, since we have
5412 // all lines flagged if they shall be refined
5413#ifdef DEBUG
5414 for (const auto &cell : triangulation.active_cell_iterators())
5415 if (!cell->refine_flag_set())
5416 for (unsigned int line = 0;
5417 line < GeometryInfo<dim>::lines_per_cell;
5418 ++line)
5419 if (cell->line(line)->has_children())
5420 for (unsigned int c = 0; c < 2; ++c)
5421 Assert(cell->line(line)->child(c)->user_flag_set() == false,
5423#endif
5424
5426 // Do refinement on every level
5427 //
5428 // To make life a bit easier, we first refine those lines and
5429 // quads that were flagged for refinement and then compose the
5430 // newly to be created cells.
5431 //
5432 // index of next unused vertex
5433 unsigned int next_unused_vertex = 0;
5434
5435 // first for lines
5436 {
5437 // only active objects can be refined further
5439 line = triangulation.begin_active_line(),
5440 endl = triangulation.end_line();
5442 next_unused_line = triangulation.begin_raw_line();
5443
5444 for (; line != endl; ++line)
5445 if (line->user_flag_set())
5446 {
5447 // this line needs to be refined
5448
5449 // find the next unused vertex and set it
5450 // appropriately
5451 while (triangulation.vertices_used[next_unused_vertex] == true)
5452 ++next_unused_vertex;
5453 Assert(
5454 next_unused_vertex < triangulation.vertices.size(),
5455 ExcMessage(
5456 "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."));
5457 triangulation.vertices_used[next_unused_vertex] = true;
5458
5459 triangulation.vertices[next_unused_vertex] = line->center(true);
5460
5461 // now that we created the right point, make up the
5462 // two child lines (++ takes care of the end of the
5463 // vector)
5464 next_unused_line =
5465 triangulation.faces->lines.template next_free_pair_object<1>(
5467 Assert(next_unused_line.state() == IteratorState::valid,
5469
5470 // now we found two consecutive unused lines, such
5471 // that the children of a line will be consecutive.
5472 // then set the child pointer of the present line
5473 line->set_children(0, next_unused_line->index());
5474
5475 // set the two new lines
5477 children[2] = {next_unused_line, ++next_unused_line};
5478
5479 // some tests; if any of the iterators should be
5480 // invalid, then already dereferencing will fail
5481 Assert(
5482 children[0]->used() == false,
5483 ExcMessage(
5484 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5485 Assert(
5486 children[1]->used() == false,
5487 ExcMessage(
5488 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5489
5490 children[0]->set_bounding_object_indices(
5491 {line->vertex_index(0), next_unused_vertex});
5492 children[1]->set_bounding_object_indices(
5493 {next_unused_vertex, line->vertex_index(1)});
5494
5495 children[0]->set_used_flag();
5496 children[1]->set_used_flag();
5497 children[0]->clear_children();
5498 children[1]->clear_children();
5499 children[0]->clear_user_data();
5500 children[1]->clear_user_data();
5501 children[0]->clear_user_flag();
5502 children[1]->clear_user_flag();
5503
5504 children[0]->set_boundary_id_internal(line->boundary_id());
5505 children[1]->set_boundary_id_internal(line->boundary_id());
5506
5507 children[0]->set_manifold_id(line->manifold_id());
5508 children[1]->set_manifold_id(line->manifold_id());
5509
5510 // finally clear flag
5511 // indicating the need
5512 // for refinement
5513 line->clear_user_flag();
5514 }
5515 }
5516
5517
5519 // now refine marked quads
5521
5522 // here we encounter several cases:
5523
5524 // a) the quad is unrefined and shall be refined isotropically
5525
5526 // b) the quad is unrefined and shall be refined
5527 // anisotropically
5528
5529 // c) the quad is unrefined and shall be refined both
5530 // anisotropically and isotropically (this is reduced to case
5531 // b) and then case b) for the children again)
5532
5533 // d) the quad is refined anisotropically and shall be refined
5534 // isotropically (this is reduced to case b) for the
5535 // anisotropic children)
5536
5537 // e) the quad is refined isotropically and shall be refined
5538 // anisotropically (this is transformed to case c), however we
5539 // might have to renumber/rename children...)
5540
5541 // we need a loop in cases c) and d), as the anisotropic
5542 // children migt have a lower index than the mother quad
5543 for (unsigned int loop = 0; loop < 2; ++loop)
5544 {
5545 // usually, only active objects can be refined
5546 // further. however, in cases d) and e) that is not true,
5547 // so we have to use 'normal' iterators here
5549 quad = triangulation.begin_quad(),
5550 endq = triangulation.end_quad();
5552 next_unused_line = triangulation.begin_raw_line();
5554 next_unused_quad = triangulation.begin_raw_quad();
5555
5556 for (; quad != endq; ++quad)
5557 {
5558 if (quad->user_index())
5559 {
5560 RefinementCase<dim - 1> aniso_quad_ref_case =
5561 face_refinement_cases[quad->user_index()];
5562 // there is one unlikely event here, where we
5563 // already have refind the face: if the face was
5564 // refined anisotropically and we want to refine
5565 // it isotropically, both children are flagged for
5566 // anisotropic refinement. however, if those
5567 // children were already flagged for anisotropic
5568 // refinement, they might already be processed and
5569 // refined.
5570 if (aniso_quad_ref_case == quad->refinement_case())
5571 continue;
5572
5573 Assert(quad->refinement_case() ==
5575 quad->refinement_case() ==
5578
5579 // this quad needs to be refined anisotropically
5580 Assert(quad->user_index() ==
5582 quad->user_index() ==
5585
5586 // make the new line interior to the quad
5588 new_line;
5589
5590 new_line =
5591 triangulation.faces->lines
5592 .template next_free_single_object<1>(triangulation);
5593 Assert(
5594 new_line->used() == false,
5595 ExcMessage(
5596 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5597
5598 // first collect the
5599 // indices of the vertices:
5600 // *--1--*
5601 // | | |
5602 // | | | cut_x
5603 // | | |
5604 // *--0--*
5605 //
5606 // *-----*
5607 // | |
5608 // 0-----1 cut_y
5609 // | |
5610 // *-----*
5611 unsigned int vertex_indices[2];
5612 if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5613 {
5614 vertex_indices[0] =
5615 quad->line(2)->child(0)->vertex_index(1);
5616 vertex_indices[1] =
5617 quad->line(3)->child(0)->vertex_index(1);
5618 }
5619 else
5620 {
5621 vertex_indices[0] =
5622 quad->line(0)->child(0)->vertex_index(1);
5623 vertex_indices[1] =
5624 quad->line(1)->child(0)->vertex_index(1);
5625 }
5626
5627 new_line->set_bounding_object_indices(
5629 new_line->set_used_flag();
5630 new_line->clear_user_flag();
5631 new_line->clear_user_data();
5632 new_line->clear_children();
5633 new_line->set_boundary_id_internal(quad->boundary_id());
5634 new_line->set_manifold_id(quad->manifold_id());
5635
5636 // child 0 and 1 of a line are switched if the
5637 // line orientation is false. set up a miniature
5638 // table, indicating which child to take for line
5639 // orientations false and true. first index: child
5640 // index in standard orientation, second index:
5641 // line orientation
5642 const unsigned int index[2][2] = {
5643 {1, 0}, // child 0, line_orientation=false and true
5644 {0, 1}}; // child 1, line_orientation=false and true
5645
5646 // find some space (consecutive) for the two newly
5647 // to be created quads.
5649 new_quads[2];
5650
5651 next_unused_quad =
5652 triangulation.faces->quads
5653 .template next_free_pair_object<2>(triangulation);
5654 new_quads[0] = next_unused_quad;
5655 Assert(
5656 new_quads[0]->used() == false,
5657 ExcMessage(
5658 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5659
5660 ++next_unused_quad;
5661 new_quads[1] = next_unused_quad;
5662 Assert(
5663 new_quads[1]->used() == false,
5664 ExcMessage(
5665 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
5666
5667
5668 if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5669 {
5670 new_quads[0]->set_bounding_object_indices(
5671 {static_cast<int>(quad->line_index(0)),
5672 new_line->index(),
5673 quad->line(2)
5674 ->child(index[0][quad->line_orientation(2)])
5675 ->index(),
5676 quad->line(3)
5677 ->child(index[0][quad->line_orientation(3)])
5678 ->index()});
5679 new_quads[1]->set_bounding_object_indices(
5680 {new_line->index(),
5681 static_cast<int>(quad->line_index(1)),
5682 quad->line(2)
5683 ->child(index[1][quad->line_orientation(2)])
5684 ->index(),
5685 quad->line(3)
5686 ->child(index[1][quad->line_orientation(3)])
5687 ->index()});
5688 }
5689 else
5690 {
5691 new_quads[0]->set_bounding_object_indices(
5692 {quad->line(0)
5693 ->child(index[0][quad->line_orientation(0)])
5694 ->index(),
5695 quad->line(1)
5696 ->child(index[0][quad->line_orientation(1)])
5697 ->index(),
5698 static_cast<int>(quad->line_index(2)),
5699 new_line->index()});
5700 new_quads[1]->set_bounding_object_indices(
5701 {quad->line(0)
5702 ->child(index[1][quad->line_orientation(0)])
5703 ->index(),
5704 quad->line(1)
5705 ->child(index[1][quad->line_orientation(1)])
5706 ->index(),
5707 new_line->index(),
5708 static_cast<int>(quad->line_index(3))});
5709 }
5710
5711 for (const auto &new_quad : new_quads)
5712 {
5713 new_quad->set_used_flag();
5714 new_quad->clear_user_flag();
5715 new_quad->clear_user_data();
5716 new_quad->clear_children();
5717 new_quad->set_boundary_id_internal(quad->boundary_id());
5718 new_quad->set_manifold_id(quad->manifold_id());
5719 // set all line orientations to true, change
5720 // this after the loop, as we have to consider
5721 // different lines for each child
5722 for (unsigned int j = 0;
5723 j < GeometryInfo<dim>::lines_per_face;
5724 ++j)
5725 new_quad->set_line_orientation(j, true);
5726 }
5727 // now set the line orientation of children of
5728 // outer lines correctly, the lines in the
5729 // interior of the refined quad are automatically
5730 // oriented conforming to the standard
5731 new_quads[0]->set_line_orientation(
5732 0, quad->line_orientation(0));
5733 new_quads[0]->set_line_orientation(
5734 2, quad->line_orientation(2));
5735 new_quads[1]->set_line_orientation(
5736 1, quad->line_orientation(1));
5737 new_quads[1]->set_line_orientation(
5738 3, quad->line_orientation(3));
5739 if (aniso_quad_ref_case == RefinementCase<dim - 1>::cut_x)
5740 {
5741 new_quads[0]->set_line_orientation(
5742 3, quad->line_orientation(3));
5743 new_quads[1]->set_line_orientation(
5744 2, quad->line_orientation(2));
5745 }
5746 else
5747 {
5748 new_quads[0]->set_line_orientation(
5749 1, quad->line_orientation(1));
5750 new_quads[1]->set_line_orientation(
5751 0, quad->line_orientation(0));
5752 }
5753
5754 // test, whether this face is refined
5755 // isotropically already. if so, set the correct
5756 // children pointers.
5757 if (quad->refinement_case() ==
5758 RefinementCase<dim - 1>::cut_xy)
5759 {
5760 // we will put a new refinemnt level of
5761 // anisotropic refinement between the
5762 // unrefined and isotropically refined quad
5763 // ending up with the same fine quads but
5764 // introducing anisotropically refined ones as
5765 // children of the unrefined quad and mother
5766 // cells of the original fine ones.
5767
5768 // this process includes the creation of a new
5769 // middle line which we will assign as the
5770 // mother line of two of the existing inner
5771 // lines. If those inner lines are not
5772 // consecutive in memory, we won't find them
5773 // later on, so we have to create new ones
5774 // instead and replace all occurrences of the
5775 // old ones with those new ones. As this is
5776 // kind of ugly, we hope we don't have to do
5777 // it often...
5779 old_child[2];
5780 if (aniso_quad_ref_case ==
5782 {
5783 old_child[0] = quad->child(0)->line(1);
5784 old_child[1] = quad->child(2)->line(1);
5785 }
5786 else
5787 {
5788 Assert(aniso_quad_ref_case ==
5791
5792 old_child[0] = quad->child(0)->line(3);
5793 old_child[1] = quad->child(1)->line(3);
5794 }
5795
5796 if (old_child[0]->index() + 1 != old_child[1]->index())
5797 {
5798 // this is exactly the ugly case we taked
5799 // about. so, no coimplaining, lets get
5800 // two new lines and copy all info
5801 typename Triangulation<dim,
5802 spacedim>::raw_line_iterator
5803 new_child[2];
5804
5805 new_child[0] = new_child[1] =
5806 triangulation.faces->lines
5807 .template next_free_pair_object<1>(
5809 ++new_child[1];
5810
5811 new_child[0]->set_used_flag();
5812 new_child[1]->set_used_flag();
5813
5814 const int old_index_0 = old_child[0]->index(),
5815 old_index_1 = old_child[1]->index(),
5816 new_index_0 = new_child[0]->index(),
5817 new_index_1 = new_child[1]->index();
5818
5819 // loop over all quads and replace the old
5820 // lines
5821 for (unsigned int q = 0;
5822 q < triangulation.faces->quads.n_objects();
5823 ++q)
5824 for (unsigned int l = 0;
5825 l < GeometryInfo<dim>::lines_per_face;
5826 ++l)
5827 {
5828 const int this_index =
5829 triangulation.faces->quads
5830 .get_bounding_object_indices(q)[l];
5831 if (this_index == old_index_0)
5832 triangulation.faces->quads
5833 .get_bounding_object_indices(q)[l] =
5834 new_index_0;
5835 else if (this_index == old_index_1)
5836 triangulation.faces->quads
5837 .get_bounding_object_indices(q)[l] =
5838 new_index_1;
5839 }
5840 // now we have to copy all information of
5841 // the two lines
5842 for (unsigned int i = 0; i < 2; ++i)
5843 {
5844 Assert(!old_child[i]->has_children(),
5846
5847 new_child[i]->set_bounding_object_indices(
5848 {old_child[i]->vertex_index(0),
5849 old_child[i]->vertex_index(1)});
5850 new_child[i]->set_boundary_id_internal(
5851 old_child[i]->boundary_id());
5852 new_child[i]->set_manifold_id(
5853 old_child[i]->manifold_id());
5854 new_child[i]->set_user_index(
5855 old_child[i]->user_index());
5856 if (old_child[i]->user_flag_set())
5857 new_child[i]->set_user_flag();
5858 else
5859 new_child[i]->clear_user_flag();
5860
5861 new_child[i]->clear_children();
5862
5863 old_child[i]->clear_user_flag();
5864 old_child[i]->clear_user_index();
5865 old_child[i]->clear_used_flag();
5866 }
5867 }
5868 // now that we cared about the lines, go on
5869 // with the quads themselves, where we might
5870 // encounter similar situations...
5871 if (aniso_quad_ref_case ==
5873 {
5874 new_line->set_children(
5875 0, quad->child(0)->line_index(1));
5876 Assert(new_line->child(1) ==
5877 quad->child(2)->line(1),
5879 // now evereything is quite
5880 // complicated. we have the children
5881 // numbered according to
5882 //
5883 // *---*---*
5884 // |n+2|n+3|
5885 // *---*---*
5886 // | n |n+1|
5887 // *---*---*
5888 //
5889 // from the original isotropic
5890 // refinement. we have to reorder them as
5891 //
5892 // *---*---*
5893 // |n+1|n+3|
5894 // *---*---*
5895 // | n |n+2|
5896 // *---*---*
5897 //
5898 // such that n and n+1 are consecutive
5899 // children of m and n+2 and n+3 are
5900 // consecutive children of m+1, where m
5901 // and m+1 are given as in
5902 //
5903 // *---*---*
5904 // | | |
5905 // | m |m+1|
5906 // | | |
5907 // *---*---*
5908 //
5909 // this is a bit ugly, of course: loop
5910 // over all cells on all levels and look
5911 // for faces n+1 (switch_1) and n+2
5912 // (switch_2).
5913 const typename Triangulation<dim, spacedim>::
5914 quad_iterator switch_1 = quad->child(1),
5915 switch_2 = quad->child(2);
5916 const int switch_1_index = switch_1->index();
5917 const int switch_2_index = switch_2->index();
5918 for (unsigned int l = 0;
5919 l < triangulation.levels.size();
5920 ++l)
5921 for (unsigned int h = 0;
5922 h <
5923 triangulation.levels[l]->cells.n_objects();
5924 ++h)
5925 for (const unsigned int q :
5927 {
5928 const int face_index =
5929 triangulation.levels[l]
5930 ->cells.get_bounding_object_indices(
5931 h)[q];
5932 if (face_index == switch_1_index)
5933 triangulation.levels[l]
5934 ->cells.get_bounding_object_indices(
5935 h)[q] = switch_2_index;
5936 else if (face_index == switch_2_index)
5937 triangulation.levels[l]
5938 ->cells.get_bounding_object_indices(
5939 h)[q] = switch_1_index;
5940 }
5941 // now we have to copy all information of
5942 // the two quads
5943 const unsigned int switch_1_lines[4] = {
5944 switch_1->line_index(0),
5945 switch_1->line_index(1),
5946 switch_1->line_index(2),
5947 switch_1->line_index(3)};
5948 const bool switch_1_line_orientations[4] = {
5949 switch_1->line_orientation(0),
5950 switch_1->line_orientation(1),
5951 switch_1->line_orientation(2),
5952 switch_1->line_orientation(3)};
5953 const types::boundary_id switch_1_boundary_id =
5954 switch_1->boundary_id();
5955 const unsigned int switch_1_user_index =
5956 switch_1->user_index();
5957 const bool switch_1_user_flag =
5958 switch_1->user_flag_set();
5959 const RefinementCase<dim - 1>
5960 switch_1_refinement_case =
5961 switch_1->refinement_case();
5962 const int switch_1_first_child_pair =
5963 (switch_1_refinement_case ?
5964 switch_1->child_index(0) :
5965 -1);
5966 const int switch_1_second_child_pair =
5967 (switch_1_refinement_case ==
5968 RefinementCase<dim - 1>::cut_xy ?
5969 switch_1->child_index(2) :
5970 -1);
5971
5972 switch_1->set_bounding_object_indices(
5973 {switch_2->line_index(0),
5974 switch_2->line_index(1),
5975 switch_2->line_index(2),
5976 switch_2->line_index(3)});
5977 switch_1->set_line_orientation(
5978 0, switch_2->line_orientation(0));
5979 switch_1->set_line_orientation(
5980 1, switch_2->line_orientation(1));
5981 switch_1->set_line_orientation(
5982 2, switch_2->line_orientation(2));
5983 switch_1->set_line_orientation(
5984 3, switch_2->line_orientation(3));
5985 switch_1->set_boundary_id_internal(
5986 switch_2->boundary_id());
5987 switch_1->set_manifold_id(switch_2->manifold_id());
5988 switch_1->set_user_index(switch_2->user_index());
5989 if (switch_2->user_flag_set())
5990 switch_1->set_user_flag();
5991 else
5992 switch_1->clear_user_flag();
5993 switch_1->clear_refinement_case();
5994 switch_1->set_refinement_case(
5995 switch_2->refinement_case());
5996 switch_1->clear_children();
5997 if (switch_2->refinement_case())
5998 switch_1->set_children(0,
5999 switch_2->child_index(0));
6000 if (switch_2->refinement_case() ==
6001 RefinementCase<dim - 1>::cut_xy)
6002 switch_1->set_children(2,
6003 switch_2->child_index(2));
6004
6005 switch_2->set_bounding_object_indices(
6006 {switch_1_lines[0],
6007 switch_1_lines[1],
6008 switch_1_lines[2],
6009 switch_1_lines[3]});
6010 switch_2->set_line_orientation(
6011 0, switch_1_line_orientations[0]);
6012 switch_2->set_line_orientation(
6013 1, switch_1_line_orientations[1]);
6014 switch_2->set_line_orientation(
6015 2, switch_1_line_orientations[2]);
6016 switch_2->set_line_orientation(
6017 3, switch_1_line_orientations[3]);
6018 switch_2->set_boundary_id_internal(
6019 switch_1_boundary_id);
6020 switch_2->set_manifold_id(switch_1->manifold_id());
6021 switch_2->set_user_index(switch_1_user_index);
6022 if (switch_1_user_flag)
6023 switch_2->set_user_flag();
6024 else
6025 switch_2->clear_user_flag();
6026 switch_2->clear_refinement_case();
6027 switch_2->set_refinement_case(
6028 switch_1_refinement_case);
6029 switch_2->clear_children();
6030 switch_2->set_children(0,
6031 switch_1_first_child_pair);
6032 switch_2->set_children(2,
6033 switch_1_second_child_pair);
6034
6035 new_quads[0]->set_refinement_case(
6037 new_quads[0]->set_children(0, quad->child_index(0));
6038 new_quads[1]->set_refinement_case(
6040 new_quads[1]->set_children(0, quad->child_index(2));
6041 }
6042 else
6043 {
6044 new_quads[0]->set_refinement_case(
6046 new_quads[0]->set_children(0, quad->child_index(0));
6047 new_quads[1]->set_refinement_case(
6049 new_quads[1]->set_children(0, quad->child_index(2));
6050 new_line->set_children(
6051 0, quad->child(0)->line_index(3));
6052 Assert(new_line->child(1) ==
6053 quad->child(1)->line(3),
6055 }
6056 quad->clear_children();
6057 }
6058
6059 // note these quads as children to the present one
6060 quad->set_children(0, new_quads[0]->index());
6061
6062 quad->set_refinement_case(aniso_quad_ref_case);
6063
6064 // finally clear flag indicating the need for
6065 // refinement
6066 quad->clear_user_data();
6067 } // if (anisotropic refinement)
6068
6069 if (quad->user_flag_set())
6070 {
6071 // this quad needs to be refined isotropically
6072
6073 // first of all: we only get here in the first run
6074 // of the loop
6075 Assert(loop == 0, ExcInternalError());
6076
6077 // find the next unused vertex. we'll need this in
6078 // any case
6079 while (triangulation.vertices_used[next_unused_vertex] ==
6080 true)
6081 ++next_unused_vertex;
6082 Assert(
6083 next_unused_vertex < triangulation.vertices.size(),
6084 ExcMessage(
6085 "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."));
6086
6087 // now: if the quad is refined anisotropically
6088 // already, set the anisotropic refinement flag
6089 // for both children. Additionally, we have to
6090 // refine the inner line, as it is an outer line
6091 // of the two (anisotropic) children
6092 const RefinementCase<dim - 1> quad_ref_case =
6093 quad->refinement_case();
6094
6095 if (quad_ref_case == RefinementCase<dim - 1>::cut_x ||
6096 quad_ref_case == RefinementCase<dim - 1>::cut_y)
6097 {
6098 // set the 'opposite' refine case for children
6099 quad->child(0)->set_user_index(
6100 RefinementCase<dim - 1>::cut_xy - quad_ref_case);
6101 quad->child(1)->set_user_index(
6102 RefinementCase<dim - 1>::cut_xy - quad_ref_case);
6103 // refine the inner line
6105 middle_line;
6106 if (quad_ref_case == RefinementCase<dim - 1>::cut_x)
6107 middle_line = quad->child(0)->line(1);
6108 else
6109 middle_line = quad->child(0)->line(3);
6110
6111 // if the face has been refined
6112 // anisotropically in the last refinement step
6113 // it might be, that it is flagged already and
6114 // that the middle line is thus refined
6115 // already. if not create children.
6116 if (!middle_line->has_children())
6117 {
6118 // set the middle vertex
6119 // appropriately. double refinement of
6120 // quads can only happen in the interior
6121 // of the domain, so we need not care
6122 // about boundary quads here
6123 triangulation.vertices[next_unused_vertex] =
6124 middle_line->center(true);
6125 triangulation.vertices_used[next_unused_vertex] =
6126 true;
6127
6128 // now search a slot for the two
6129 // child lines
6130 next_unused_line =
6131 triangulation.faces->lines
6132 .template next_free_pair_object<1>(
6134
6135 // set the child pointer of the present
6136 // line
6137 middle_line->set_children(
6138 0, next_unused_line->index());
6139
6140 // set the two new lines
6141 const typename Triangulation<dim, spacedim>::
6142 raw_line_iterator children[2] = {
6143 next_unused_line, ++next_unused_line};
6144
6145 // some tests; if any of the iterators
6146 // should be invalid, then already
6147 // dereferencing will fail
6148 Assert(
6149 children[0]->used() == false,
6150 ExcMessage(
6151 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6152 Assert(
6153 children[1]->used() == false,
6154 ExcMessage(
6155 "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
6156
6157 children[0]->set_bounding_object_indices(
6158 {middle_line->vertex_index(0),
6159 next_unused_vertex});
6160 children[1]->set_bounding_object_indices(
6161 {next_unused_vertex,
6162 middle_line->vertex_index(1)});
6163
6164 children[0]->set_used_flag();
6165 children[1]->set_used_flag();
6166 children[0]->clear_children();
6167 children[1]->clear_children();
6168 children[0]->clear_user_data();
6169 children[1]->clear_user_data();
6170 children[0]->clear_user_flag();
6171 children[1]->clear_user_flag();
6172
6173 children[0]->set_boundary_id_internal(
6174 middle_line->boundary_id());
6175 children[1]->set_boundary_id_internal(
6176 middle_line->boundary_id());
6177
6178 children[0]->set_manifold_id(
6179 middle_line->manifold_id());
6180 children[1]->set_manifold_id(
6181 middle_line->manifold_id());
6182 }
6183 // now remove the flag from the quad and go to
6184 // the next quad, the actual refinement of the
6185 // quad takes place later on in this pass of
6186 // the loop or in the next one
6187 quad->clear_user_flag();
6188 continue;
6189 } // if (several refinement cases)
6190
6191 // if we got here, we have an unrefined quad and
6192 // have to do the usual work like in an purely
6193 // isotropic refinement
6194 Assert(quad_ref_case ==
6197
6198 // set the middle vertex appropriately: it might be that
6199 // the quad itself is not at the boundary, but that one of
6200 // its lines actually is. in this case, the newly created
6201 // vertices at the centers of the lines are not