Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2022 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
481 reorder_compatibility(std::vector<CellData<2>> &cells, 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
490 reorder_compatibility(std::vector<CellData<3>> &cells,
491 SubCellData & subcelldata)
492 {
493 unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
494 for (auto &cell : cells)
495 if (cell.vertices.size() == GeometryInfo<3>::vertices_per_cell)
496 {
497 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
498 tmp[i] = cell.vertices[i];
499 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
500 cell.vertices[GeometryInfo<3>::ucd_to_deal[i]] = tmp[i];
501 }
502
503 // now points in boundary quads
504 for (auto &boundary_quad : subcelldata.boundary_quads)
505 if (boundary_quad.vertices.size() == GeometryInfo<2>::vertices_per_cell)
506 std::swap(boundary_quad.vertices[2], boundary_quad.vertices[3]);
507 }
508
509
510
528 template <int dim, int spacedim>
529 unsigned int
530 middle_vertex_index(
532 {
533 if (line->has_children())
534 return line->child(0)->vertex_index(1);
536 }
537
538
539 template <int dim, int spacedim>
540 unsigned int
541 middle_vertex_index(
543 {
544 switch (static_cast<unsigned char>(quad->refinement_case()))
545 {
547 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
548 break;
550 return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
551 break;
553 return quad->child(0)->vertex_index(3);
554 break;
555 default:
556 break;
557 }
559 }
560
561
562 template <int dim, int spacedim>
563 unsigned int
564 middle_vertex_index(
566 {
567 switch (static_cast<unsigned char>(hex->refinement_case()))
568 {
570 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
571 break;
573 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
574 break;
576 return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
577 break;
579 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
580 break;
582 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
583 break;
585 return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
586 break;
588 return hex->child(0)->vertex_index(7);
589 break;
590 default:
591 break;
592 }
594 }
595
596
609 template <class TRIANGULATION>
610 inline typename TRIANGULATION::DistortedCellList
611 collect_distorted_coarse_cells(const TRIANGULATION &)
612 {
613 return typename TRIANGULATION::DistortedCellList();
614 }
615
616
617
626 template <int dim>
628 collect_distorted_coarse_cells(const Triangulation<dim, dim> &triangulation)
629 {
630 typename Triangulation<dim, dim>::DistortedCellList distorted_cells;
631 for (const auto &cell : triangulation.cell_iterators_on_level(0))
632 {
634 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
635 vertices[i] = cell->vertex(i);
636
639
640 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
641 if (determinants[i] <= 1e-9 * std::pow(cell->diameter(), 1. * dim))
642 {
643 distorted_cells.distorted_cells.push_back(cell);
644 break;
645 }
646 }
647
648 return distorted_cells;
649 }
650
651
658 template <int dim>
659 bool
660 has_distorted_children(
661 const typename Triangulation<dim, dim>::cell_iterator &cell)
662 {
663 Assert(cell->has_children(), ExcInternalError());
664
665 for (unsigned int c = 0; c < cell->n_children(); ++c)
666 {
668 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
669 vertices[i] = cell->child(c)->vertex(i);
670
673
674 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
675 if (determinants[i] <=
676 1e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
677 return true;
678 }
679
680 return false;
681 }
682
683
691 template <int dim, int spacedim>
692 bool
693 has_distorted_children(
695 {
696 return false;
697 }
698
699
700 template <int dim, int spacedim>
701 void
702 update_periodic_face_map_recursively(
703 const typename Triangulation<dim, spacedim>::cell_iterator &cell_1,
704 const typename Triangulation<dim, spacedim>::cell_iterator &cell_2,
705 unsigned int n_face_1,
706 unsigned int n_face_2,
707 const std::bitset<3> & orientation,
708 typename std::map<
710 unsigned int>,
711 std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
712 unsigned int>,
713 std::bitset<3>>> &periodic_face_map)
714 {
715 using FaceIterator = typename Triangulation<dim, spacedim>::face_iterator;
716 const FaceIterator face_1 = cell_1->face(n_face_1);
717 const FaceIterator face_2 = cell_2->face(n_face_2);
718
719 const bool face_orientation = orientation[0];
720 const bool face_flip = orientation[1];
721 const bool face_rotation = orientation[2];
722
723 Assert((dim != 1) || (face_orientation == true && face_flip == false &&
724 face_rotation == false),
725 ExcMessage("The supplied orientation "
726 "(face_orientation, face_flip, face_rotation) "
727 "is invalid for 1D"));
728
729 Assert((dim != 2) || (face_orientation == true && face_rotation == false),
730 ExcMessage("The supplied orientation "
731 "(face_orientation, face_flip, face_rotation) "
732 "is invalid for 2D"));
733
734 Assert(face_1 != face_2, ExcMessage("face_1 and face_2 are equal!"));
735
736 Assert(face_1->at_boundary() && face_2->at_boundary(),
737 ExcMessage("Periodic faces must be on the boundary"));
738
739 // Check if the requirement that each edge can only have at most one hanging
740 // node, and as a consequence neighboring cells can differ by at most
741 // one refinement level is enforced. In 1d, there are no hanging nodes and
742 // so neighboring cells can differ by more than one refinement level.
743 Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
745
746 // insert periodic face pair for both cells
747 using CellFace =
748 std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
749 unsigned int>;
750 const CellFace cell_face_1(cell_1, n_face_1);
751 const CellFace cell_face_2(cell_2, n_face_2);
752 const std::pair<CellFace, std::bitset<3>> cell_face_orientation_2(
753 cell_face_2, orientation);
754
755 const std::pair<CellFace, std::pair<CellFace, std::bitset<3>>>
756 periodic_faces(cell_face_1, cell_face_orientation_2);
757
758 // Only one periodic neighbor is allowed
759 Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
760 periodic_face_map.insert(periodic_faces);
761
762 if (dim == 1)
763 {
764 if (cell_1->has_children())
765 {
766 if (cell_2->has_children())
767 {
768 update_periodic_face_map_recursively<dim, spacedim>(
769 cell_1->child(n_face_1),
770 cell_2->child(n_face_2),
771 n_face_1,
772 n_face_2,
773 orientation,
774 periodic_face_map);
775 }
776 else // only face_1 has children
777 {
778 update_periodic_face_map_recursively<dim, spacedim>(
779 cell_1->child(n_face_1),
780 cell_2,
781 n_face_1,
782 n_face_2,
783 orientation,
784 periodic_face_map);
785 }
786 }
787 }
788 else // dim == 2 || dim == 3
789 {
790 // A lookup table on how to go through the child cells depending on the
791 // orientation:
792 // see Documentation of GeometryInfo for details
793
794 static const int lookup_table_2d[2][2] =
795 // flip:
796 {
797 {0, 1}, // false
798 {1, 0} // true
799 };
800
801 static const int lookup_table_3d[2][2][2][4] =
802 // orientation flip rotation
803 {{{
804 {0, 2, 1, 3}, // false false false
805 {2, 3, 0, 1} // false false true
806 },
807 {
808 {3, 1, 2, 0}, // false true false
809 {1, 0, 3, 2} // false true true
810 }},
811 {{
812 {0, 1, 2, 3}, // true false false
813 {1, 3, 0, 2} // true false true
814 },
815 {
816 {3, 2, 1, 0}, // true true false
817 {2, 0, 3, 1} // true true true
818 }}};
819
820 if (cell_1->has_children())
821 {
822 if (cell_2->has_children())
823 {
824 // In the case that both faces have children, we loop over all
825 // children and apply update_periodic_face_map_recursively
826 // recursively:
827
828 Assert(face_1->n_children() ==
830 face_2->n_children() ==
833
834 for (unsigned int i = 0;
835 i < GeometryInfo<dim>::max_children_per_face;
836 ++i)
837 {
838 // Lookup the index for the second face
839 unsigned int j = 0;
840 switch (dim)
841 {
842 case 2:
843 j = lookup_table_2d[face_flip][i];
844 break;
845 case 3:
846 j = lookup_table_3d[face_orientation][face_flip]
847 [face_rotation][i];
848 break;
849 default:
851 }
852
853 // find subcell ids that belong to the subface indices
854 unsigned int child_cell_1 =
856 cell_1->refinement_case(),
857 n_face_1,
858 i,
859 cell_1->face_orientation(n_face_1),
860 cell_1->face_flip(n_face_1),
861 cell_1->face_rotation(n_face_1),
862 face_1->refinement_case());
863 unsigned int child_cell_2 =
865 cell_2->refinement_case(),
866 n_face_2,
867 j,
868 cell_2->face_orientation(n_face_2),
869 cell_2->face_flip(n_face_2),
870 cell_2->face_rotation(n_face_2),
871 face_2->refinement_case());
872
873 Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
874 face_1->child(i),
876 Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
877 face_2->child(j),
879
880 // precondition: subcell has the same orientation as cell
881 // (so that the face numbers coincide) recursive call
882 update_periodic_face_map_recursively<dim, spacedim>(
883 cell_1->child(child_cell_1),
884 cell_2->child(child_cell_2),
885 n_face_1,
886 n_face_2,
887 orientation,
888 periodic_face_map);
889 }
890 }
891 else // only face_1 has children
892 {
893 for (unsigned int i = 0;
894 i < GeometryInfo<dim>::max_children_per_face;
895 ++i)
896 {
897 // find subcell ids that belong to the subface indices
898 unsigned int child_cell_1 =
900 cell_1->refinement_case(),
901 n_face_1,
902 i,
903 cell_1->face_orientation(n_face_1),
904 cell_1->face_flip(n_face_1),
905 cell_1->face_rotation(n_face_1),
906 face_1->refinement_case());
907
908 // recursive call
909 update_periodic_face_map_recursively<dim, spacedim>(
910 cell_1->child(child_cell_1),
911 cell_2,
912 n_face_1,
913 n_face_2,
914 orientation,
915 periodic_face_map);
916 }
917 }
918 }
919 }
920 }
921
922
923} // end of anonymous namespace
924
925
926namespace internal
927{
928 namespace TriangulationImplementation
929 {
930 // make sure that if in the following we
931 // write Triangulation<dim,spacedim>
932 // we mean the *class*
933 // ::Triangulation, not the
934 // enclosing namespace
935 // internal::TriangulationImplementation
936 using ::Triangulation;
937
943 int,
944 << "Something went wrong when making cell " << arg1
945 << ". Read the docs and the source code "
946 << "for more information.");
952 int,
953 << "Something went wrong upon construction of cell "
954 << arg1);
965 int,
966 << "Cell " << arg1
967 << " has negative measure. This typically "
968 << "indicates some distortion in the cell, or a mistakenly "
969 << "swapped pair of vertices in the input to "
970 << "Triangulation::create_triangulation().");
979 int,
980 int,
981 int,
982 << "Error while creating cell " << arg1
983 << ": the vertex index " << arg2 << " must be between 0 and "
984 << arg3 << '.');
990 int,
991 int,
992 << "While trying to assign a boundary indicator to a line: "
993 << "the line with end vertices " << arg1 << " and " << arg2
994 << " does not exist.");
1000 int,
1001 int,
1002 int,
1003 int,
1004 << "While trying to assign a boundary indicator to a quad: "
1005 << "the quad with bounding lines " << arg1 << ", " << arg2
1006 << ", " << arg3 << ", " << arg4 << " does not exist.");
1013 int,
1014 int,
1016 << "The input data for creating a triangulation contained "
1017 << "information about a line with indices " << arg1 << " and " << arg2
1018 << " that is described to have boundary indicator "
1019 << static_cast<int>(arg3)
1020 << ". However, this is an internal line not located on the "
1021 << "boundary. You cannot assign a boundary indicator to it." << std::endl
1022 << std::endl
1023 << "If this happened at a place where you call "
1024 << "Triangulation::create_triangulation() yourself, you need "
1025 << "to check the SubCellData object you pass to this function."
1026 << std::endl
1027 << std::endl
1028 << "If this happened in a place where you are reading a mesh "
1029 << "from a file, then you need to investigate why such a line "
1030 << "ended up in the input file. A typical case is a geometry "
1031 << "that consisted of multiple parts and for which the mesh "
1032 << "generator program assumes that the interface between "
1033 << "two parts is a boundary when that isn't supposed to be "
1034 << "the case, or where the mesh generator simply assigns "
1035 << "'geometry indicators' to lines at the perimeter of "
1036 << "a part that are not supposed to be interpreted as "
1037 << "'boundary indicators'.");
1044 int,
1045 int,
1046 int,
1047 int,
1049 << "The input data for creating a triangulation contained "
1050 << "information about a quad with indices " << arg1 << ", " << arg2
1051 << ", " << arg3 << ", and " << arg4
1052 << " that is described to have boundary indicator "
1053 << static_cast<int>(arg5)
1054 << ". However, this is an internal quad not located on the "
1055 << "boundary. You cannot assign a boundary indicator to it." << std::endl
1056 << std::endl
1057 << "If this happened at a place where you call "
1058 << "Triangulation::create_triangulation() yourself, you need "
1059 << "to check the SubCellData object you pass to this function."
1060 << std::endl
1061 << std::endl
1062 << "If this happened in a place where you are reading a mesh "
1063 << "from a file, then you need to investigate why such a quad "
1064 << "ended up in the input file. A typical case is a geometry "
1065 << "that consisted of multiple parts and for which the mesh "
1066 << "generator program assumes that the interface between "
1067 << "two parts is a boundary when that isn't supposed to be "
1068 << "the case, or where the mesh generator simply assigns "
1069 << "'geometry indicators' to quads at the surface of "
1070 << "a part that are not supposed to be interpreted as "
1071 << "'boundary indicators'.");
1078 int,
1079 int,
1080 << "In SubCellData the line info of the line with vertex indices " << arg1
1081 << " and " << arg2 << " appears more than once. "
1082 << "This is not allowed.");
1089 int,
1090 int,
1091 std::string,
1092 << "In SubCellData the line info of the line with vertex indices " << arg1
1093 << " and " << arg2 << " appears multiple times with different (valid) "
1094 << arg3 << ". This is not allowed.");
1101 int,
1102 int,
1103 int,
1104 int,
1105 std::string,
1106 << "In SubCellData the quad info of the quad with line indices " << arg1
1107 << ", " << arg2 << ", " << arg3 << " and " << arg4
1108 << " appears multiple times with different (valid) " << arg5
1109 << ". This is not allowed.");
1110
1111 /*
1112 * Reserve space for TriaFaces. Details:
1113 *
1114 * Reserve space for line_orientations.
1115 *
1116 * @note Used only for dim=3.
1117 */
1118 void
1120 const unsigned int new_quads_in_pairs,
1121 const unsigned int new_quads_single)
1122 {
1123 AssertDimension(tria_faces.dim, 3);
1124
1125 Assert(new_quads_in_pairs % 2 == 0, ExcInternalError());
1126
1127 unsigned int next_free_single = 0;
1128 unsigned int next_free_pair = 0;
1129
1130 // count the number of objects, of unused single objects and of
1131 // unused pairs of objects
1132 unsigned int n_quads = 0;
1133 unsigned int n_unused_pairs = 0;
1134 unsigned int n_unused_singles = 0;
1135 for (unsigned int i = 0; i < tria_faces.quads.used.size(); ++i)
1136 {
1137 if (tria_faces.quads.used[i])
1138 ++n_quads;
1139 else if (i + 1 < tria_faces.quads.used.size())
1140 {
1141 if (tria_faces.quads.used[i + 1])
1142 {
1143 ++n_unused_singles;
1144 if (next_free_single == 0)
1145 next_free_single = i;
1146 }
1147 else
1148 {
1149 ++n_unused_pairs;
1150 if (next_free_pair == 0)
1151 next_free_pair = i;
1152 ++i;
1153 }
1154 }
1155 else
1156 ++n_unused_singles;
1157 }
1158 Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
1159 tria_faces.quads.used.size(),
1161
1162 // how many single quads are needed in addition to n_unused_quads?
1163 const int additional_single_quads = new_quads_single - n_unused_singles;
1164
1165 unsigned int new_size =
1166 tria_faces.quads.used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
1167 if (additional_single_quads > 0)
1168 new_size += additional_single_quads;
1169
1170 // see above...
1171 if (new_size > tria_faces.quads.n_objects())
1172 {
1173 // reserve the field of the derived class
1174 tria_faces.quads_line_orientations.reserve(
1176 tria_faces.quads_line_orientations.insert(
1177 tria_faces.quads_line_orientations.end(),
1179 tria_faces.quads_line_orientations.size(),
1180 1u);
1181
1182 tria_faces.quad_reference_cell.reserve(new_size);
1183 tria_faces.quad_reference_cell.insert(
1184 tria_faces.quad_reference_cell.end(),
1185 new_size - tria_faces.quad_reference_cell.size(),
1187 }
1188 }
1189
1190
1191
1205 void
1207 const unsigned int total_cells,
1208 const unsigned int dimension,
1209 const unsigned int space_dimension)
1210 {
1211 // we need space for total_cells cells. Maybe we have more already
1212 // with those cells which are unused, so only allocate new space if
1213 // needed.
1214 //
1215 // note that all arrays should have equal sizes (checked by
1216 // @p{monitor_memory}
1217 if (total_cells > tria_level.refine_flags.size())
1218 {
1219 tria_level.refine_flags.reserve(total_cells);
1220 tria_level.refine_flags.insert(tria_level.refine_flags.end(),
1221 total_cells -
1222 tria_level.refine_flags.size(),
1223 /*RefinementCase::no_refinement=*/0);
1224
1225 tria_level.coarsen_flags.reserve(total_cells);
1226 tria_level.coarsen_flags.insert(tria_level.coarsen_flags.end(),
1227 total_cells -
1228 tria_level.coarsen_flags.size(),
1229 false);
1230
1231 tria_level.active_cell_indices.reserve(total_cells);
1232 tria_level.active_cell_indices.insert(
1233 tria_level.active_cell_indices.end(),
1234 total_cells - tria_level.active_cell_indices.size(),
1236
1237 tria_level.subdomain_ids.reserve(total_cells);
1238 tria_level.subdomain_ids.insert(tria_level.subdomain_ids.end(),
1239 total_cells -
1240 tria_level.subdomain_ids.size(),
1241 0);
1242
1243 tria_level.level_subdomain_ids.reserve(total_cells);
1244 tria_level.level_subdomain_ids.insert(
1245 tria_level.level_subdomain_ids.end(),
1246 total_cells - tria_level.level_subdomain_ids.size(),
1247 0);
1248
1249 tria_level.global_active_cell_indices.reserve(total_cells);
1250 tria_level.global_active_cell_indices.insert(
1251 tria_level.global_active_cell_indices.end(),
1252 total_cells - tria_level.global_active_cell_indices.size(),
1254
1255 tria_level.global_level_cell_indices.reserve(total_cells);
1256 tria_level.global_level_cell_indices.insert(
1257 tria_level.global_level_cell_indices.end(),
1258 total_cells - tria_level.global_level_cell_indices.size(),
1260
1261 if (dimension < space_dimension)
1262 {
1263 tria_level.direction_flags.reserve(total_cells);
1264 tria_level.direction_flags.insert(
1265 tria_level.direction_flags.end(),
1266 total_cells - tria_level.direction_flags.size(),
1267 true);
1268 }
1269 else
1270 tria_level.direction_flags.clear();
1271
1272 tria_level.parents.reserve((total_cells + 1) / 2);
1273 tria_level.parents.insert(tria_level.parents.end(),
1274 (total_cells + 1) / 2 -
1275 tria_level.parents.size(),
1276 -1);
1277
1278 tria_level.neighbors.reserve(total_cells * (2 * dimension));
1279 tria_level.neighbors.insert(tria_level.neighbors.end(),
1280 total_cells * (2 * dimension) -
1281 tria_level.neighbors.size(),
1282 std::make_pair(-1, -1));
1283
1284 if (tria_level.dim == 2 || tria_level.dim == 3)
1285 {
1286 const unsigned int max_faces_per_cell = 2 * dimension;
1287 tria_level.face_orientations.reserve(total_cells *
1288 max_faces_per_cell);
1289 tria_level.face_orientations.insert(
1290 tria_level.face_orientations.end(),
1291 total_cells * max_faces_per_cell -
1292 tria_level.face_orientations.size(),
1293 1u);
1294
1295 tria_level.reference_cell.reserve(total_cells);
1296 tria_level.reference_cell.insert(
1297 tria_level.reference_cell.end(),
1298 total_cells - tria_level.reference_cell.size(),
1299 tria_level.dim == 2 ? ::ReferenceCells::Quadrilateral :
1301 }
1302 }
1303 }
1304
1305
1306
1311 int,
1312 int,
1313 << "The containers have sizes " << arg1 << " and " << arg2
1314 << ", which is not as expected.");
1315
1321 void
1322 monitor_memory(const TriaLevel & tria_level,
1323 const unsigned int true_dimension)
1324 {
1325 (void)tria_level;
1326 (void)true_dimension;
1327 Assert(2 * true_dimension * tria_level.refine_flags.size() ==
1328 tria_level.neighbors.size(),
1329 ExcMemoryInexact(tria_level.refine_flags.size(),
1330 tria_level.neighbors.size()));
1331 Assert(2 * true_dimension * tria_level.coarsen_flags.size() ==
1332 tria_level.neighbors.size(),
1333 ExcMemoryInexact(tria_level.coarsen_flags.size(),
1334 tria_level.neighbors.size()));
1335 }
1336
1337
1338
1351 void
1353 const unsigned int new_objects_in_pairs,
1354 const unsigned int new_objects_single = 0)
1355 {
1356 if (tria_objects.structdim <= 2)
1357 {
1358 Assert(new_objects_in_pairs % 2 == 0, ExcInternalError());
1359
1360 tria_objects.next_free_single = 0;
1361 tria_objects.next_free_pair = 0;
1362 tria_objects.reverse_order_next_free_single = false;
1363
1364 // count the number of objects, of unused single objects and of
1365 // unused pairs of objects
1366 unsigned int n_objects = 0;
1367 unsigned int n_unused_pairs = 0;
1368 unsigned int n_unused_singles = 0;
1369 for (unsigned int i = 0; i < tria_objects.used.size(); ++i)
1370 {
1371 if (tria_objects.used[i])
1372 ++n_objects;
1373 else if (i + 1 < tria_objects.used.size())
1374 {
1375 if (tria_objects.used[i + 1])
1376 {
1377 ++n_unused_singles;
1378 if (tria_objects.next_free_single == 0)
1379 tria_objects.next_free_single = i;
1380 }
1381 else
1382 {
1383 ++n_unused_pairs;
1384 if (tria_objects.next_free_pair == 0)
1385 tria_objects.next_free_pair = i;
1386 ++i;
1387 }
1388 }
1389 else
1390 ++n_unused_singles;
1391 }
1392 Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
1393 tria_objects.used.size(),
1395
1396 // how many single objects are needed in addition to
1397 // n_unused_objects?
1398 const int additional_single_objects =
1399 new_objects_single - n_unused_singles;
1400
1401 unsigned int new_size = tria_objects.used.size() +
1402 new_objects_in_pairs - 2 * n_unused_pairs;
1403 if (additional_single_objects > 0)
1404 new_size += additional_single_objects;
1405
1406 // only allocate space if necessary
1407 if (new_size > tria_objects.n_objects())
1408 {
1409 const unsigned int max_faces_per_cell =
1410 2 * tria_objects.structdim;
1411 const unsigned int max_children_per_cell =
1412 1 << tria_objects.structdim;
1413
1414 tria_objects.cells.reserve(new_size * max_faces_per_cell);
1415 tria_objects.cells.insert(tria_objects.cells.end(),
1416 (new_size - tria_objects.n_objects()) *
1417 max_faces_per_cell,
1418 -1);
1419
1420 tria_objects.used.reserve(new_size);
1421 tria_objects.used.insert(tria_objects.used.end(),
1422 new_size - tria_objects.used.size(),
1423 false);
1424
1425 tria_objects.user_flags.reserve(new_size);
1426 tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1427 new_size -
1428 tria_objects.user_flags.size(),
1429 false);
1430
1431 const unsigned int factor = max_children_per_cell / 2;
1432 tria_objects.children.reserve(factor * new_size);
1433 tria_objects.children.insert(tria_objects.children.end(),
1434 factor * new_size -
1435 tria_objects.children.size(),
1436 -1);
1437
1438 if (tria_objects.structdim > 1)
1439 {
1440 tria_objects.refinement_cases.reserve(new_size);
1441 tria_objects.refinement_cases.insert(
1442 tria_objects.refinement_cases.end(),
1443 new_size - tria_objects.refinement_cases.size(),
1444 /*RefinementCase::no_refinement=*/0);
1445 }
1446
1447 // first reserve, then resize. Otherwise the std library can
1448 // decide to allocate more entries.
1449 tria_objects.boundary_or_material_id.reserve(new_size);
1450 tria_objects.boundary_or_material_id.resize(new_size);
1451
1452 tria_objects.user_data.reserve(new_size);
1453 tria_objects.user_data.resize(new_size);
1454
1455 tria_objects.manifold_id.reserve(new_size);
1456 tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1457 new_size -
1458 tria_objects.manifold_id.size(),
1460 }
1461
1462 if (n_unused_singles == 0)
1463 {
1464 tria_objects.next_free_single = new_size - 1;
1465 tria_objects.reverse_order_next_free_single = true;
1466 }
1467 }
1468 else
1469 {
1470 const unsigned int new_hexes = new_objects_in_pairs;
1471
1472 const unsigned int new_size =
1473 new_hexes + std::count(tria_objects.used.begin(),
1474 tria_objects.used.end(),
1475 true);
1476
1477 // see above...
1478 if (new_size > tria_objects.n_objects())
1479 {
1480 const unsigned int max_faces_per_cell =
1481 2 * tria_objects.structdim;
1482
1483 tria_objects.cells.reserve(new_size * max_faces_per_cell);
1484 tria_objects.cells.insert(tria_objects.cells.end(),
1485 (new_size - tria_objects.n_objects()) *
1486 max_faces_per_cell,
1487 -1);
1488
1489 tria_objects.used.reserve(new_size);
1490 tria_objects.used.insert(tria_objects.used.end(),
1491 new_size - tria_objects.used.size(),
1492 false);
1493
1494 tria_objects.user_flags.reserve(new_size);
1495 tria_objects.user_flags.insert(tria_objects.user_flags.end(),
1496 new_size -
1497 tria_objects.user_flags.size(),
1498 false);
1499
1500 tria_objects.children.reserve(4 * new_size);
1501 tria_objects.children.insert(tria_objects.children.end(),
1502 4 * new_size -
1503 tria_objects.children.size(),
1504 -1);
1505
1506 // for the following fields, we know exactly how many elements
1507 // we need, so first reserve then resize (resize itself, at least
1508 // with some compiler libraries, appears to round up the size it
1509 // actually reserves)
1510 tria_objects.boundary_or_material_id.reserve(new_size);
1511 tria_objects.boundary_or_material_id.resize(new_size);
1512
1513 tria_objects.manifold_id.reserve(new_size);
1514 tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
1515 new_size -
1516 tria_objects.manifold_id.size(),
1518
1519 tria_objects.user_data.reserve(new_size);
1520 tria_objects.user_data.resize(new_size);
1521
1522 tria_objects.refinement_cases.reserve(new_size);
1523 tria_objects.refinement_cases.insert(
1524 tria_objects.refinement_cases.end(),
1525 new_size - tria_objects.refinement_cases.size(),
1526 /*RefinementCase::no_refinement=*/0);
1527 }
1528 tria_objects.next_free_single = tria_objects.next_free_pair = 0;
1529 }
1530 }
1531
1532
1533
1539 void
1540 monitor_memory(const TriaObjects &tria_object, const unsigned int)
1541 {
1542 Assert(tria_object.n_objects() == tria_object.used.size(),
1543 ExcMemoryInexact(tria_object.n_objects(),
1544 tria_object.used.size()));
1545 Assert(tria_object.n_objects() == tria_object.user_flags.size(),
1546 ExcMemoryInexact(tria_object.n_objects(),
1547 tria_object.user_flags.size()));
1548 Assert(tria_object.n_objects() ==
1549 tria_object.boundary_or_material_id.size(),
1550 ExcMemoryInexact(tria_object.n_objects(),
1551 tria_object.boundary_or_material_id.size()));
1552 Assert(tria_object.n_objects() == tria_object.manifold_id.size(),
1553 ExcMemoryInexact(tria_object.n_objects(),
1554 tria_object.manifold_id.size()));
1555 Assert(tria_object.n_objects() == tria_object.user_data.size(),
1556 ExcMemoryInexact(tria_object.n_objects(),
1557 tria_object.user_data.size()));
1558
1559 if (tria_object.structdim == 1)
1560 {
1561 Assert(1 * tria_object.n_objects() == tria_object.children.size(),
1562 ExcMemoryInexact(tria_object.n_objects(),
1563 tria_object.children.size()));
1564 }
1565 else if (tria_object.structdim == 2)
1566 {
1567 Assert(2 * tria_object.n_objects() == tria_object.children.size(),
1568 ExcMemoryInexact(tria_object.n_objects(),
1569 tria_object.children.size()));
1570 }
1571 else if (tria_object.structdim == 3)
1572 {
1573 Assert(4 * tria_object.n_objects() == tria_object.children.size(),
1574 ExcMemoryInexact(tria_object.n_objects(),
1575 tria_object.children.size()));
1576 }
1577 }
1578
1579
1580
1585 template <int dim, int spacedim>
1587 {
1588 public:
1592 virtual ~Policy() = default;
1593
1597 virtual void
1599
1603 virtual void
1607 std::vector<unsigned int> & line_cell_count,
1608 std::vector<unsigned int> &quad_cell_count) = 0;
1609
1615 const bool check_for_distorted_cells) = 0;
1616
1620 virtual void
1623
1627 virtual void
1630
1634 virtual bool
1636 const typename Triangulation<dim, spacedim>::cell_iterator &cell) = 0;
1637
1644 virtual std::unique_ptr<Policy<dim, spacedim>>
1645 clone() = 0;
1646 };
1647
1648
1649
1655 template <int dim, int spacedim, typename T>
1656 class PolicyWrapper : public Policy<dim, spacedim>
1657 {
1658 public:
1659 void
1661 {
1662 T::update_neighbors(tria);
1663 }
1664
1665 void
1669 std::vector<unsigned int> & line_cell_count,
1670 std::vector<unsigned int> &quad_cell_count) override
1671 {
1672 T::delete_children(tria, cell, line_cell_count, quad_cell_count);
1673 }
1674
1677 const bool check_for_distorted_cells) override
1678 {
1679 return T::execute_refinement(triangulation, check_for_distorted_cells);
1680 }
1681
1682 void
1685 {
1686 T::prevent_distorted_boundary_cells(triangulation);
1687 }
1688
1689 void
1692 {
1693 T::prepare_refinement_dim_dependent(triangulation);
1694 }
1695
1696 bool
1699 override
1700 {
1701 return T::template coarsening_allowed<dim, spacedim>(cell);
1702 }
1703
1704 std::unique_ptr<Policy<dim, spacedim>>
1705 clone() override
1706 {
1707 return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
1708 }
1709 };
1710
1711
1712
1809 {
1821 template <int dim, int spacedim>
1822 static void
1825 const unsigned int level_objects,
1827 {
1828 using line_iterator =
1830
1831 number_cache.n_levels = 0;
1832 if (level_objects > 0)
1833 // find the last level on which there are used cells
1834 for (unsigned int level = 0; level < level_objects; ++level)
1835 if (triangulation.begin(level) != triangulation.end(level))
1836 number_cache.n_levels = level + 1;
1837
1838 // no cells at all?
1839 Assert(number_cache.n_levels > 0, ExcInternalError());
1840
1841 //---------------------------------
1842 // update the number of lines on the different levels in the
1843 // cache
1844 number_cache.n_lines = 0;
1845 number_cache.n_active_lines = 0;
1846
1847 // for 1d, lines have levels so take count the objects per
1848 // level and globally
1849 if (dim == 1)
1850 {
1851 number_cache.n_lines_level.resize(number_cache.n_levels);
1852 number_cache.n_active_lines_level.resize(number_cache.n_levels);
1853
1854 for (unsigned int level = 0; level < number_cache.n_levels; ++level)
1855 {
1856 // count lines on this level
1857 number_cache.n_lines_level[level] = 0;
1858 number_cache.n_active_lines_level[level] = 0;
1859
1860 line_iterator line = triangulation.begin_line(level),
1861 endc =
1862 (level == number_cache.n_levels - 1 ?
1863 line_iterator(triangulation.end_line()) :
1864 triangulation.begin_line(level + 1));
1865 for (; line != endc; ++line)
1866 {
1867 ++number_cache.n_lines_level[level];
1868 if (line->has_children() == false)
1869 ++number_cache.n_active_lines_level[level];
1870 }
1871
1872 // update total number of lines
1873 number_cache.n_lines += number_cache.n_lines_level[level];
1874 number_cache.n_active_lines +=
1875 number_cache.n_active_lines_level[level];
1876 }
1877 }
1878 else
1879 {
1880 // for dim>1, there are no levels for lines
1881 number_cache.n_lines_level.clear();
1882 number_cache.n_active_lines_level.clear();
1883
1884 line_iterator line = triangulation.begin_line(),
1885 endc = triangulation.end_line();
1886 for (; line != endc; ++line)
1887 {
1888 ++number_cache.n_lines;
1889 if (line->has_children() == false)
1890 ++number_cache.n_active_lines;
1891 }
1892 }
1893 }
1894
1909 template <int dim, int spacedim>
1910 static void
1913 const unsigned int level_objects,
1915 {
1916 // update lines and n_levels in number_cache. since we don't
1917 // access any of these numbers, we can do this in the
1918 // background
1920 static_cast<
1921 void (*)(const Triangulation<dim, spacedim> &,
1922 const unsigned int,
1924 &compute_number_cache<dim, spacedim>),
1926 level_objects,
1928 number_cache));
1929
1930 using quad_iterator =
1932
1933 //---------------------------------
1934 // update the number of quads on the different levels in the
1935 // cache
1936 number_cache.n_quads = 0;
1937 number_cache.n_active_quads = 0;
1938
1939 // for 2d, quads have levels so take count the objects per
1940 // level and globally
1941 if (dim == 2)
1942 {
1943 // count the number of levels; the function we called above
1944 // on a separate Task for lines also does this and puts it into
1945 // number_cache.n_levels, but this datum may not yet be
1946 // available as we call the function on a separate task
1947 unsigned int n_levels = 0;
1948 if (level_objects > 0)
1949 // find the last level on which there are used cells
1950 for (unsigned int level = 0; level < level_objects; ++level)
1951 if (triangulation.begin(level) != triangulation.end(level))
1952 n_levels = level + 1;
1953
1954 number_cache.n_quads_level.resize(n_levels);
1955 number_cache.n_active_quads_level.resize(n_levels);
1956
1957 for (unsigned int level = 0; level < n_levels; ++level)
1958 {
1959 // count quads on this level
1960 number_cache.n_quads_level[level] = 0;
1961 number_cache.n_active_quads_level[level] = 0;
1962
1963 quad_iterator quad = triangulation.begin_quad(level),
1964 endc =
1965 (level == n_levels - 1 ?
1966 quad_iterator(triangulation.end_quad()) :
1967 triangulation.begin_quad(level + 1));
1968 for (; quad != endc; ++quad)
1969 {
1970 ++number_cache.n_quads_level[level];
1971 if (quad->has_children() == false)
1972 ++number_cache.n_active_quads_level[level];
1973 }
1974
1975 // update total number of quads
1976 number_cache.n_quads += number_cache.n_quads_level[level];
1977 number_cache.n_active_quads +=
1978 number_cache.n_active_quads_level[level];
1979 }
1980 }
1981 else
1982 {
1983 // for dim>2, there are no levels for quads
1984 number_cache.n_quads_level.clear();
1985 number_cache.n_active_quads_level.clear();
1986
1987 quad_iterator quad = triangulation.begin_quad(),
1988 endc = triangulation.end_quad();
1989 for (; quad != endc; ++quad)
1990 {
1991 ++number_cache.n_quads;
1992 if (quad->has_children() == false)
1993 ++number_cache.n_active_quads;
1994 }
1995 }
1996
1997 // wait for the background computation for lines
1998 update_lines.join();
1999 }
2000
2016 template <int dim, int spacedim>
2017 static void
2020 const unsigned int level_objects,
2022 {
2023 // update quads, lines and n_levels in number_cache. since we
2024 // don't access any of these numbers, we can do this in the
2025 // background
2026 Threads::Task<void> update_quads_and_lines = Threads::new_task(
2027 static_cast<
2028 void (*)(const Triangulation<dim, spacedim> &,
2029 const unsigned int,
2031 &compute_number_cache<dim, spacedim>),
2033 level_objects,
2035 number_cache));
2036
2037 using hex_iterator =
2039
2040 //---------------------------------
2041 // update the number of hexes on the different levels in the
2042 // cache
2043 number_cache.n_hexes = 0;
2044 number_cache.n_active_hexes = 0;
2045
2046 // for 3d, hexes have levels so take count the objects per
2047 // level and globally
2048 if (dim == 3)
2049 {
2050 // count the number of levels; the function we called
2051 // above on a separate Task for quads (recursively, via
2052 // the lines function) also does this and puts it into
2053 // number_cache.n_levels, but this datum may not yet be
2054 // available as we call the function on a separate task
2055 unsigned int n_levels = 0;
2056 if (level_objects > 0)
2057 // find the last level on which there are used cells
2058 for (unsigned int level = 0; level < level_objects; ++level)
2059 if (triangulation.begin(level) != triangulation.end(level))
2060 n_levels = level + 1;
2061
2062 number_cache.n_hexes_level.resize(n_levels);
2063 number_cache.n_active_hexes_level.resize(n_levels);
2064
2065 for (unsigned int level = 0; level < n_levels; ++level)
2066 {
2067 // count hexes on this level
2068 number_cache.n_hexes_level[level] = 0;
2069 number_cache.n_active_hexes_level[level] = 0;
2070
2071 hex_iterator hex = triangulation.begin_hex(level),
2072 endc = (level == n_levels - 1 ?
2073 hex_iterator(triangulation.end_hex()) :
2074 triangulation.begin_hex(level + 1));
2075 for (; hex != endc; ++hex)
2076 {
2077 ++number_cache.n_hexes_level[level];
2078 if (hex->has_children() == false)
2079 ++number_cache.n_active_hexes_level[level];
2080 }
2081
2082 // update total number of hexes
2083 number_cache.n_hexes += number_cache.n_hexes_level[level];
2084 number_cache.n_active_hexes +=
2085 number_cache.n_active_hexes_level[level];
2086 }
2087 }
2088 else
2089 {
2090 // for dim>3, there are no levels for hexes
2091 number_cache.n_hexes_level.clear();
2092 number_cache.n_active_hexes_level.clear();
2093
2094 hex_iterator hex = triangulation.begin_hex(),
2095 endc = triangulation.end_hex();
2096 for (; hex != endc; ++hex)
2097 {
2098 ++number_cache.n_hexes;
2099 if (hex->has_children() == false)
2100 ++number_cache.n_active_hexes;
2101 }
2102 }
2103
2104 // wait for the background computation for quads
2105 update_quads_and_lines.join();
2106 }
2107
2108
2109
2110 template <int spacedim>
2111 static void
2113 {}
2114
2115
2116 template <int dim, int spacedim>
2117 static void
2119 {
2120 // each face can be neighbored on two sides
2121 // by cells. according to the face's
2122 // intrinsic normal we define the left
2123 // neighbor as the one for which the face
2124 // normal points outward, and store that
2125 // one first; the second one is then
2126 // the right neighbor for which the
2127 // face normal points inward. This
2128 // information depends on the type of cell
2129 // and local number of face for the
2130 // 'standard ordering and orientation' of
2131 // faces and then on the face_orientation
2132 // information for the real mesh. Set up a
2133 // table to have fast access to those
2134 // offsets (0 for left and 1 for
2135 // right). Some of the values are invalid
2136 // as they reference too large face
2137 // numbers, but we just leave them at a
2138 // zero value.
2139 //
2140 // Note, that in 2d for lines as faces the
2141 // normal direction given in the
2142 // GeometryInfo class is not consistent. We
2143 // thus define here that the normal for a
2144 // line points to the right if the line
2145 // points upwards.
2146 //
2147 // There is one more point to
2148 // consider, however: if we have
2149 // dim<spacedim, then we may have
2150 // cases where cells are
2151 // inverted. In effect, both
2152 // cells think they are the left
2153 // neighbor of an edge, for
2154 // example, which leads us to
2155 // forget neighborship
2156 // information (a case that shows
2157 // this is
2158 // codim_one/hanging_nodes_02). We
2159 // store whether a cell is
2160 // inverted using the
2161 // direction_flag, so if a cell
2162 // has a false direction_flag,
2163 // then we need to invert our
2164 // selection whether we are a
2165 // left or right neighbor in all
2166 // following computations.
2167 //
2168 // first index: dimension (minus 2)
2169 // second index: local face index
2170 // third index: face_orientation (false and true)
2171 static const unsigned int left_right_offset[2][6][2] = {
2172 // quadrilateral
2173 {{0, 1}, // face 0, face_orientation = false and true
2174 {1, 0}, // face 1, face_orientation = false and true
2175 {1, 0}, // face 2, face_orientation = false and true
2176 {0, 1}, // face 3, face_orientation = false and true
2177 {0, 0}, // face 4, invalid face
2178 {0, 0}}, // face 5, invalid face
2179 // hexahedron
2180 {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
2181
2182 // now create a vector of the two active
2183 // neighbors (left and right) for each face
2184 // and fill it by looping over all cells. For
2185 // cases with anisotropic refinement and more
2186 // then one cell neighboring at a given side
2187 // of the face we will automatically get the
2188 // active one on the highest level as we loop
2189 // over cells from lower levels first.
2191 std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
2192 adjacent_cells(2 * triangulation.n_raw_faces(), dummy);
2193
2194 for (const auto &cell : triangulation.cell_iterators())
2195 for (auto f : cell->face_indices())
2196 {
2198 cell->face(f);
2199
2200 const unsigned int offset =
2201 (cell->direction_flag() ?
2202 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2203 1 -
2204 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2205
2206 adjacent_cells[2 * face->index() + offset] = cell;
2207
2208 // if this cell is not refined, but the
2209 // face is, then we'll have to set our
2210 // cell as neighbor for the child faces
2211 // as well. Fortunately the normal
2212 // orientation of children will be just
2213 // the same.
2214 if (dim == 2)
2215 {
2216 if (cell->is_active() && face->has_children())
2217 {
2218 adjacent_cells[2 * face->child(0)->index() + offset] =
2219 cell;
2220 adjacent_cells[2 * face->child(1)->index() + offset] =
2221 cell;
2222 }
2223 }
2224 else // -> dim == 3
2225 {
2226 // We need the same as in 2d
2227 // here. Furthermore, if the face is
2228 // refined with cut_x or cut_y then
2229 // those children again in the other
2230 // direction, and if this cell is
2231 // refined isotropically (along the
2232 // face) then the neighbor will
2233 // (probably) be refined as cut_x or
2234 // cut_y along the face. For those
2235 // neighboring children cells, their
2236 // neighbor will be the current,
2237 // inactive cell, as our children are
2238 // too fine to be neighbors. Catch that
2239 // case by also acting on inactive
2240 // cells with isotropic refinement
2241 // along the face. If the situation
2242 // described is not present, the data
2243 // will be overwritten later on when we
2244 // visit cells on finer levels, so no
2245 // harm will be done.
2246 if (face->has_children() &&
2247 (cell->is_active() ||
2249 cell->refinement_case(), f) ==
2250 RefinementCase<dim - 1>::isotropic_refinement))
2251 {
2252 for (unsigned int c = 0; c < face->n_children(); ++c)
2253 adjacent_cells[2 * face->child(c)->index() + offset] =
2254 cell;
2255 if (face->child(0)->has_children())
2256 {
2257 adjacent_cells[2 * face->child(0)->child(0)->index() +
2258 offset] = cell;
2259 adjacent_cells[2 * face->child(0)->child(1)->index() +
2260 offset] = cell;
2261 }
2262 if (face->child(1)->has_children())
2263 {
2264 adjacent_cells[2 * face->child(1)->child(0)->index() +
2265 offset] = cell;
2266 adjacent_cells[2 * face->child(1)->child(1)->index() +
2267 offset] = cell;
2268 }
2269 } // if cell active and face refined
2270 } // else -> dim==3
2271 } // for all faces of all cells
2272
2273 // now loop again over all cells and set the
2274 // corresponding neighbor cell. Note, that we
2275 // have to use the opposite of the
2276 // left_right_offset in this case as we want
2277 // the offset of the neighbor, not our own.
2278 for (const auto &cell : triangulation.cell_iterators())
2279 for (auto f : cell->face_indices())
2280 {
2281 const unsigned int offset =
2282 (cell->direction_flag() ?
2283 left_right_offset[dim - 2][f][cell->face_orientation(f)] :
2284 1 -
2285 left_right_offset[dim - 2][f][cell->face_orientation(f)]);
2286 cell->set_neighbor(
2287 f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
2288 }
2289 }
2290
2291
2295 template <int dim, int spacedim>
2296 static void
2298 const std::vector<CellData<dim>> & cells,
2299 const SubCellData & subcelldata,
2301 {
2302 AssertThrow(vertices.size() > 0, ExcMessage("No vertices given"));
2303 AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
2304
2305 // Check that all cells have positive volume.
2306#ifndef _MSC_VER
2307 // TODO: The following code does not compile with MSVC. Find a way
2308 // around it
2309 if (dim == spacedim)
2310 for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
2311 {
2312 // If we should check for distorted cells, then we permit them
2313 // to exist. If a cell has negative measure, then it must be
2314 // distorted (the converse is not necessarily true); hence
2315 // throw an exception if no such cells should exist.
2317 {
2318 const double cell_measure = GridTools::cell_measure<spacedim>(
2319 vertices,
2320 ArrayView<const unsigned int>(cells[cell_no].vertices));
2322 }
2323 }
2324#endif
2325
2326 // clear old content
2327 tria.levels.clear();
2328 tria.levels.push_back(
2329 std::make_unique<
2331
2332 if (dim > 1)
2333 tria.faces = std::make_unique<
2335
2336 // copy vertices
2338 tria.vertices_used.assign(vertices.size(), true);
2339
2340 // compute connectivity
2341 const auto connectivity = build_connectivity<unsigned int>(cells);
2342 const unsigned int n_cell = cells.size();
2343
2344 // TriaObjects: lines
2345 if (dim >= 2)
2346 {
2347 auto &lines_0 = tria.faces->lines; // data structure to be filled
2348
2349 // get connectivity between quads and lines
2350 const auto & crs = connectivity.entity_to_entities(1, 0);
2351 const unsigned int n_lines = crs.ptr.size() - 1;
2352
2353 // allocate memory
2354 reserve_space_(lines_0, n_lines);
2355
2356 // loop over lines
2357 for (unsigned int line = 0; line < n_lines; ++line)
2358 for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
2359 ++i, ++j)
2360 lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
2361 crs.col[i]; // set vertex indices
2362 }
2363
2364 // TriaObjects: quads
2365 if (dim == 3)
2366 {
2367 auto &quads_0 = tria.faces->quads; // data structures to be filled
2368 auto &faces = *tria.faces;
2369
2370 // get connectivity between quads and lines
2371 const auto & crs = connectivity.entity_to_entities(2, 1);
2372 const unsigned int n_quads = crs.ptr.size() - 1;
2373
2374 // allocate memory
2375 reserve_space_(quads_0, n_quads);
2376 reserve_space_(faces, 2 /*structdim*/, n_quads);
2377
2378 // loop over all quads -> entity type, line indices/orientations
2379 for (unsigned int q = 0, k = 0; q < n_quads; ++q)
2380 {
2381 // set entity type of quads
2382 faces.quad_reference_cell[q] = connectivity.entity_types(2)[q];
2383
2384 // loop over all its lines
2385 for (unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
2386 ++i, ++j, ++k)
2387 {
2388 // set line index
2389 quads_0.cells[q * GeometryInfo<2>::faces_per_cell + j] =
2390 crs.col[i];
2391
2392 // set line orientations
2393 faces.quads_line_orientations
2395 connectivity.entity_orientations(1)[k];
2396 }
2397 }
2398 }
2399
2400 // TriaObjects/TriaLevel: cell
2401 {
2402 auto &cells_0 = tria.levels[0]->cells; // data structure to be filled
2403 auto &level = *tria.levels[0];
2404
2405 // get connectivity between cells/faces and cells/cells
2406 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2407 const auto &nei = connectivity.entity_to_entities(dim, dim);
2408
2409 // in 2D optional: since in in pure QUAD meshes same line
2410 // orientations can be guaranteed
2411 const bool orientation_needed =
2412 dim == 3 ||
2413 (dim == 2 &&
2414 std::any_of(connectivity.entity_orientations(1).begin(),
2415 connectivity.entity_orientations(1).end(),
2416 [](const auto &i) { return i == 0; }));
2417
2418 // allocate memory
2419 reserve_space_(cells_0, n_cell);
2420 reserve_space_(level, spacedim, n_cell, orientation_needed);
2421
2422 // loop over all cells
2423 for (unsigned int cell = 0; cell < n_cell; ++cell)
2424 {
2425 // set material ids
2426 cells_0.boundary_or_material_id[cell].material_id =
2427 cells[cell].material_id;
2428
2429 // set manifold ids
2430 cells_0.manifold_id[cell] = cells[cell].manifold_id;
2431
2432 // set entity types
2433 level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
2434
2435 // loop over faces
2436 for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2437 ++i, ++j)
2438 {
2439 // set neighbor if not at boundary
2440 if (nei.col[i] != static_cast<unsigned int>(-1))
2441 level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
2442 j] = {0, nei.col[i]};
2443
2444 // set face indices
2445 cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
2446 crs.col[i];
2447
2448 // set face orientation if needed
2449 if (orientation_needed)
2450 {
2451 level.face_orientations
2453 connectivity.entity_orientations(dim - 1)[i];
2454 }
2455 }
2456 }
2457 }
2458
2459 // TriaFaces: boundary id of boundary faces
2460 if (dim > 1)
2461 {
2462 auto &bids_face = dim == 3 ?
2463 tria.faces->quads.boundary_or_material_id :
2464 tria.faces->lines.boundary_or_material_id;
2465
2466 // count number of cells a face is belonging to
2467 std::vector<unsigned int> count(bids_face.size(), 0);
2468
2469 // get connectivity between cells/faces
2470 const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
2471
2472 // count how many cells are adjacent to the same face
2473 for (unsigned int cell = 0; cell < cells.size(); ++cell)
2474 for (unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
2475 count[crs.col[i]]++;
2476
2477 // loop over all faces
2478 for (unsigned int face = 0; face < count.size(); ++face)
2479 {
2480 if (count[face] != 1) // inner face
2481 continue;
2482
2483 // boundary faces ...
2484 bids_face[face].boundary_id = 0;
2485
2486 if (dim != 3)
2487 continue;
2488
2489 // ... and the lines of quads in 3D
2490 const auto &crs = connectivity.entity_to_entities(2, 1);
2491 for (unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
2492 tria.faces->lines.boundary_or_material_id[crs.col[i]]
2493 .boundary_id = 0;
2494 }
2495 }
2496 else // 1D
2497 {
2498 static const unsigned int t_tba = static_cast<unsigned int>(-1);
2499 static const unsigned int t_inner = static_cast<unsigned int>(-2);
2500
2501 std::vector<unsigned int> type(vertices.size(), t_tba);
2502
2503 const auto &crs = connectivity.entity_to_entities(1, 0);
2504
2505 for (unsigned int cell = 0; cell < cells.size(); ++cell)
2506 for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
2507 ++i, ++j)
2508 if (type[crs.col[i]] != t_inner)
2509 type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
2510
2511 for (unsigned int face = 0; face < type.size(); ++face)
2512 {
2513 // note: we also treat manifolds here!?
2516 if (type[face] != t_inner && type[face] != t_tba)
2517 (*tria.vertex_to_boundary_id_map_1d)[face] = type[face];
2518 }
2519 }
2520
2521 // SubCellData: line
2522 if (dim >= 2)
2523 process_subcelldata(connectivity.entity_to_entities(1, 0),
2524 tria.faces->lines,
2525 subcelldata.boundary_lines,
2526 vertices);
2527
2528 // SubCellData: quad
2529 if (dim == 3)
2530 process_subcelldata(connectivity.entity_to_entities(2, 0),
2531 tria.faces->quads,
2532 subcelldata.boundary_quads,
2533 vertices);
2534 }
2535
2536
2537 template <int structdim, int spacedim, typename T>
2538 static void
2540 const CRS<T> & crs,
2541 TriaObjects & obj,
2542 const std::vector<CellData<structdim>> &boundary_objects_in,
2543 const std::vector<Point<spacedim>> & vertex_locations)
2544 {
2545 AssertDimension(obj.structdim, structdim);
2546
2547 if (boundary_objects_in.size() == 0)
2548 return; // empty subcelldata -> nothing to do
2549
2550 // pre-sort subcelldata
2551 auto boundary_objects = boundary_objects_in;
2552
2553 // ... sort vertices
2554 for (auto &boundary_object : boundary_objects)
2555 std::sort(boundary_object.vertices.begin(),
2556 boundary_object.vertices.end());
2557
2558 // ... sort cells
2559 std::sort(boundary_objects.begin(),
2560 boundary_objects.end(),
2561 [](const auto &a, const auto &b) {
2562 return a.vertices < b.vertices;
2563 });
2564
2565 unsigned int counter = 0;
2566
2567 std::vector<unsigned int> key;
2569
2570 for (unsigned int o = 0; o < obj.n_objects(); ++o)
2571 {
2572 auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
2573 auto &manifold_id = obj.manifold_id[o];
2574
2575 // assert that object has not been visited yet and its value
2576 // has not been modified yet
2577 AssertThrow(boundary_id == 0 ||
2582
2583 // create key
2584 key.assign(crs.col.data() + crs.ptr[o],
2585 crs.col.data() + crs.ptr[o + 1]);
2586 std::sort(key.begin(), key.end());
2587
2588 // is subcelldata provided? -> binary search
2589 const auto subcell_object =
2590 std::lower_bound(boundary_objects.begin(),
2591 boundary_objects.end(),
2592 key,
2593 [&](const auto &cell, const auto &key) {
2594 return cell.vertices < key;
2595 });
2596
2597 // no subcelldata provided for this object
2598 if (subcell_object == boundary_objects.end() ||
2599 subcell_object->vertices != key)
2600 continue;
2601
2602 counter++;
2603
2604 // set manifold id
2605 manifold_id = subcell_object->manifold_id;
2606
2607 // set boundary id
2608 if (subcell_object->boundary_id !=
2610 {
2611 (void)vertex_locations;
2614 ExcMessage(
2615 "The input arguments for creating a triangulation "
2616 "specified a boundary id for an internal face. This "
2617 "is not allowed."
2618 "\n\n"
2619 "The object in question has vertex indices " +
2620 [subcell_object]() {
2621 std::string s;
2622 for (const auto v : subcell_object->vertices)
2623 s += std::to_string(v) + ',';
2624 return s;
2625 }() +
2626 " which are located at positions " +
2627 [vertex_locations, subcell_object]() {
2628 std::ostringstream s;
2629 for (const auto v : subcell_object->vertices)
2630 s << '(' << vertex_locations[v] << ')';
2631 return s.str();
2632 }() +
2633 "."));
2634 boundary_id = subcell_object->boundary_id;
2635 }
2636 }
2637
2638 // make sure that all subcelldata entries have been processed
2639 // TODO: this is not guaranteed, why?
2640 // AssertDimension(counter, boundary_objects_in.size());
2641 }
2642
2643
2644
2645 static void
2647 const unsigned structdim,
2648 const unsigned int size)
2649 {
2650 const unsigned int dim = faces.dim;
2651
2652 const unsigned int max_faces_per_cell = 2 * structdim;
2653
2654 if (dim == 3 && structdim == 2)
2655 {
2656 // quad entity types
2657 faces.quad_reference_cell.assign(size,
2659
2660 // quad line orientations
2661 faces.quads_line_orientations.assign(size * max_faces_per_cell, -1);
2662 }
2663 }
2664
2665
2666
2667 static void
2669 const unsigned int spacedim,
2670 const unsigned int size,
2671 const bool orientation_needed)
2672 {
2673 const unsigned int dim = level.dim;
2674
2675 const unsigned int max_faces_per_cell = 2 * dim;
2676
2677 level.active_cell_indices.assign(size, -1);
2678 level.subdomain_ids.assign(size, 0);
2679 level.level_subdomain_ids.assign(size, 0);
2680
2681 level.refine_flags.assign(size, 0u);
2682 level.coarsen_flags.assign(size, false);
2683
2684 level.parents.assign((size + 1) / 2, -1);
2685
2686 if (dim < spacedim)
2687 level.direction_flags.assign(size, true);
2688
2689 level.neighbors.assign(size * max_faces_per_cell, {-1, -1});
2690
2691 level.reference_cell.assign(size, ::ReferenceCells::Invalid);
2692
2693 if (orientation_needed)
2694 level.face_orientations.assign(size * max_faces_per_cell, -1);
2695
2696 level.global_active_cell_indices.assign(size,
2698 level.global_level_cell_indices.assign(size,
2700 }
2701
2702
2703
2704 static void
2705 reserve_space_(TriaObjects &obj, const unsigned int size)
2706 {
2707 const unsigned int structdim = obj.structdim;
2708
2709 const unsigned int max_children_per_cell = 1 << structdim;
2710 const unsigned int max_faces_per_cell = 2 * structdim;
2711
2712 obj.used.assign(size, true);
2713 obj.boundary_or_material_id.assign(
2714 size,
2716 BoundaryOrMaterialId());
2717 obj.manifold_id.assign(size, -1);
2718 obj.user_flags.assign(size, false);
2719 obj.user_data.resize(size);
2720
2721 if (structdim > 1) // TODO: why?
2722 obj.refinement_cases.assign(size, 0);
2723
2724 obj.children.assign(max_children_per_cell / 2 * size, -1);
2725
2726 obj.cells.assign(max_faces_per_cell * size, -1);
2727
2728 if (structdim <= 2)
2729 {
2730 obj.next_free_single = size - 1;
2731 obj.next_free_pair = 0;
2733 }
2734 else
2735 {
2736 obj.next_free_single = obj.next_free_pair = 0;
2737 }
2738 }
2739
2740
2756 template <int spacedim>
2757 static void
2760 std::vector<unsigned int> &,
2761 std::vector<unsigned int> &)
2762 {
2763 const unsigned int dim = 1;
2764
2765 // first we need to reset the
2766 // neighbor pointers of the
2767 // neighbors of this cell's
2768 // children to this cell. This is
2769 // different for one dimension,
2770 // since there neighbors can have a
2771 // refinement level differing from
2772 // that of this cell's children by
2773 // more than one level.
2774
2775 Assert(!cell->child(0)->has_children() &&
2776 !cell->child(1)->has_children(),
2778
2779 // first do it for the cells to the
2780 // left
2781 if (cell->neighbor(0).state() == IteratorState::valid)
2782 if (cell->neighbor(0)->has_children())
2783 {
2785 cell->neighbor(0);
2786 Assert(neighbor->level() == cell->level(), ExcInternalError());
2787
2788 // right child
2789 neighbor = neighbor->child(1);
2790 while (true)
2791 {
2792 Assert(neighbor->neighbor(1) == cell->child(0),
2794 neighbor->set_neighbor(1, 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(1);
2802 else
2803 break;
2804 }
2805 }
2806
2807 // now do it for the cells to the
2808 // left
2809 if (cell->neighbor(1).state() == IteratorState::valid)
2810 if (cell->neighbor(1)->has_children())
2811 {
2813 cell->neighbor(1);
2814 Assert(neighbor->level() == cell->level(), ExcInternalError());
2815
2816 // left child
2817 neighbor = neighbor->child(0);
2818 while (true)
2819 {
2820 Assert(neighbor->neighbor(0) == cell->child(1),
2822 neighbor->set_neighbor(0, cell);
2823
2824 // move on to further
2825 // children on the
2826 // boundary between this
2827 // cell and its neighbor
2828 if (neighbor->has_children())
2829 neighbor = neighbor->child(0);
2830 else
2831 break;
2832 }
2833 }
2834
2835
2836 // delete the vertex which will not
2837 // be needed anymore. This vertex
2838 // is the second of the first child
2839 triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
2840
2841 // invalidate children. clear user
2842 // pointers, to avoid that they may
2843 // appear at unwanted places later
2844 // on...
2845 for (unsigned int child = 0; child < cell->n_children(); ++child)
2846 {
2847 cell->child(child)->clear_user_data();
2848 cell->child(child)->clear_user_flag();
2849 cell->child(child)->clear_used_flag();
2850 }
2851
2852
2853 // delete pointer to children
2854 cell->clear_children();
2855 cell->clear_user_flag();
2856 }
2857
2858
2859
2860 template <int spacedim>
2861 static void
2864 std::vector<unsigned int> &line_cell_count,
2865 std::vector<unsigned int> &)
2866 {
2867 const unsigned int dim = 2;
2868 const RefinementCase<dim> ref_case = cell->refinement_case();
2869
2870 Assert(line_cell_count.size() == triangulation.n_raw_lines(),
2872
2873 // vectors to hold all lines which
2874 // may be deleted
2875 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
2876 lines_to_delete(0);
2877
2878 lines_to_delete.reserve(4 * 2 + 4);
2879
2880 // now we decrease the counters for
2881 // lines contained in the child
2882 // cells
2883 for (unsigned int c = 0; c < cell->n_children(); ++c)
2884 {
2886 cell->child(c);
2887 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2888 --line_cell_count[child->line_index(l)];
2889 }
2890
2891
2892 // delete the vertex which will not
2893 // be needed anymore. This vertex
2894 // is the second of the second line
2895 // of the first child, if the cell
2896 // is refined with cut_xy, else there
2897 // is no inner vertex.
2898 // additionally delete unneeded inner
2899 // lines
2900 if (ref_case == RefinementCase<dim>::cut_xy)
2901 {
2903 .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
2904
2905 lines_to_delete.push_back(cell->child(0)->line(1));
2906 lines_to_delete.push_back(cell->child(0)->line(3));
2907 lines_to_delete.push_back(cell->child(3)->line(0));
2908 lines_to_delete.push_back(cell->child(3)->line(2));
2909 }
2910 else
2911 {
2912 unsigned int inner_face_no =
2913 ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
2914
2915 // the inner line will not be
2916 // used any more
2917 lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
2918 }
2919
2920 // invalidate children
2921 for (unsigned int child = 0; child < cell->n_children(); ++child)
2922 {
2923 cell->child(child)->clear_user_data();
2924 cell->child(child)->clear_user_flag();
2925 cell->child(child)->clear_used_flag();
2926 }
2927
2928
2929 // delete pointer to children
2930 cell->clear_children();
2931 cell->clear_refinement_case();
2932 cell->clear_user_flag();
2933
2934 // look at the refinement of outer
2935 // lines. if nobody needs those
2936 // anymore we can add them to the
2937 // list of lines to be deleted.
2938 for (unsigned int line_no = 0;
2939 line_no < GeometryInfo<dim>::lines_per_cell;
2940 ++line_no)
2941 {
2943 cell->line(line_no);
2944
2945 if (line->has_children())
2946 {
2947 // if one of the cell counters is
2948 // zero, the other has to be as well
2949
2950 Assert((line_cell_count[line->child_index(0)] == 0 &&
2951 line_cell_count[line->child_index(1)] == 0) ||
2952 (line_cell_count[line->child_index(0)] > 0 &&
2953 line_cell_count[line->child_index(1)] > 0),
2955
2956 if (line_cell_count[line->child_index(0)] == 0)
2957 {
2958 for (unsigned int c = 0; c < 2; ++c)
2959 Assert(!line->child(c)->has_children(),
2961
2962 // we may delete the line's
2963 // children and the middle vertex
2964 // as no cell references them
2965 // anymore
2967 .vertices_used[line->child(0)->vertex_index(1)] = false;
2968
2969 lines_to_delete.push_back(line->child(0));
2970 lines_to_delete.push_back(line->child(1));
2971
2972 line->clear_children();
2973 }
2974 }
2975 }
2976
2977 // finally, delete unneeded lines
2978
2979 // clear user pointers, to avoid that
2980 // they may appear at unwanted places
2981 // later on...
2982 // same for user flags, then finally
2983 // delete the lines
2984 typename std::vector<
2986 line = lines_to_delete.begin(),
2987 endline = lines_to_delete.end();
2988 for (; line != endline; ++line)
2989 {
2990 (*line)->clear_user_data();
2991 (*line)->clear_user_flag();
2992 (*line)->clear_used_flag();
2993 }
2994 }
2995
2996
2997
2998 template <int spacedim>
2999 static void
3002 std::vector<unsigned int> &line_cell_count,
3003 std::vector<unsigned int> &quad_cell_count)
3004 {
3005 const unsigned int dim = 3;
3006
3007 Assert(line_cell_count.size() == triangulation.n_raw_lines(),
3009 Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
3011
3012 // first of all, we store the RefineCase of
3013 // this cell
3014 const RefinementCase<dim> ref_case = cell->refinement_case();
3015 // vectors to hold all lines and quads which
3016 // may be deleted
3017 std::vector<typename Triangulation<dim, spacedim>::line_iterator>
3018 lines_to_delete(0);
3019 std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
3020 quads_to_delete(0);
3021
3022 lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
3023 quads_to_delete.reserve(6 * 4 + 12);
3024
3025 // now we decrease the counters for lines and
3026 // quads contained in the child cells
3027 for (unsigned int c = 0; c < cell->n_children(); ++c)
3028 {
3030 cell->child(c);
3031 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3032 --line_cell_count[child->line_index(l)];
3033 for (auto f : GeometryInfo<dim>::face_indices())
3034 --quad_cell_count[child->quad_index(f)];
3035 }
3036
3037 //-------------------------------------
3038 // delete interior quads and lines and the
3039 // interior vertex, depending on the
3040 // refinement case of the cell
3041 //
3042 // for append quads and lines: only append
3043 // them to the list of objects to be deleted
3044
3045 switch (ref_case)
3046 {
3048 quads_to_delete.push_back(cell->child(0)->face(1));
3049 break;
3051 quads_to_delete.push_back(cell->child(0)->face(3));
3052 break;
3054 quads_to_delete.push_back(cell->child(0)->face(5));
3055 break;
3057 quads_to_delete.push_back(cell->child(0)->face(1));
3058 quads_to_delete.push_back(cell->child(0)->face(3));
3059 quads_to_delete.push_back(cell->child(3)->face(0));
3060 quads_to_delete.push_back(cell->child(3)->face(2));
3061
3062 lines_to_delete.push_back(cell->child(0)->line(11));
3063 break;
3065 quads_to_delete.push_back(cell->child(0)->face(1));
3066 quads_to_delete.push_back(cell->child(0)->face(5));
3067 quads_to_delete.push_back(cell->child(3)->face(0));
3068 quads_to_delete.push_back(cell->child(3)->face(4));
3069
3070 lines_to_delete.push_back(cell->child(0)->line(5));
3071 break;
3073 quads_to_delete.push_back(cell->child(0)->face(3));
3074 quads_to_delete.push_back(cell->child(0)->face(5));
3075 quads_to_delete.push_back(cell->child(3)->face(2));
3076 quads_to_delete.push_back(cell->child(3)->face(4));
3077
3078 lines_to_delete.push_back(cell->child(0)->line(7));
3079 break;
3081 quads_to_delete.push_back(cell->child(0)->face(1));
3082 quads_to_delete.push_back(cell->child(2)->face(1));
3083 quads_to_delete.push_back(cell->child(4)->face(1));
3084 quads_to_delete.push_back(cell->child(6)->face(1));
3085
3086 quads_to_delete.push_back(cell->child(0)->face(3));
3087 quads_to_delete.push_back(cell->child(1)->face(3));
3088 quads_to_delete.push_back(cell->child(4)->face(3));
3089 quads_to_delete.push_back(cell->child(5)->face(3));
3090
3091 quads_to_delete.push_back(cell->child(0)->face(5));
3092 quads_to_delete.push_back(cell->child(1)->face(5));
3093 quads_to_delete.push_back(cell->child(2)->face(5));
3094 quads_to_delete.push_back(cell->child(3)->face(5));
3095
3096 lines_to_delete.push_back(cell->child(0)->line(5));
3097 lines_to_delete.push_back(cell->child(0)->line(7));
3098 lines_to_delete.push_back(cell->child(0)->line(11));
3099 lines_to_delete.push_back(cell->child(7)->line(0));
3100 lines_to_delete.push_back(cell->child(7)->line(2));
3101 lines_to_delete.push_back(cell->child(7)->line(8));
3102 // delete the vertex which will not
3103 // be needed anymore. This vertex
3104 // is the vertex at the heart of
3105 // this cell, which is the sixth of
3106 // the first child
3107 triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
3108 false;
3109 break;
3110 default:
3111 // only remaining case is
3112 // no_refinement, thus an error
3113 Assert(false, ExcInternalError());
3114 break;
3115 }
3116
3117
3118 // invalidate children
3119 for (unsigned int child = 0; child < cell->n_children(); ++child)
3120 {
3121 cell->child(child)->clear_user_data();
3122 cell->child(child)->clear_user_flag();
3123
3124 for (auto f : GeometryInfo<dim>::face_indices())
3125 {
3126 // set flags denoting deviations from
3127 // standard orientation of faces back
3128 // to initialization values
3129 cell->child(child)->set_face_orientation(f, true);
3130 cell->child(child)->set_face_flip(f, false);
3131 cell->child(child)->set_face_rotation(f, false);
3132 }
3133
3134 cell->child(child)->clear_used_flag();
3135 }
3136
3137
3138 // delete pointer to children
3139 cell->clear_children();
3140 cell->clear_refinement_case();
3141 cell->clear_user_flag();
3142
3143 // so far we only looked at inner quads,
3144 // lines and vertices. Now we have to
3145 // consider outer ones as well. here, we have
3146 // to check, whether there are other cells
3147 // still needing these objects. otherwise we
3148 // can delete them. first for quads (and
3149 // their inner lines).
3150
3151 for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
3152 {
3154 cell->face(quad_no);
3155
3156 Assert(
3157 (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
3158 quad->has_children()) ||
3159 GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
3162
3163 switch (quad->refinement_case())
3164 {
3165 case RefinementCase<dim - 1>::no_refinement:
3166 // nothing to do as the quad
3167 // is not refined
3168 break;
3169 case RefinementCase<dim - 1>::cut_x:
3170 case RefinementCase<dim - 1>::cut_y:
3171 {
3172 // if one of the cell counters is
3173 // zero, the other has to be as
3174 // well
3175 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3176 quad_cell_count[quad->child_index(1)] == 0) ||
3177 (quad_cell_count[quad->child_index(0)] > 0 &&
3178 quad_cell_count[quad->child_index(1)] > 0),
3180 // it might be, that the quad is
3181 // refined twice anisotropically,
3182 // first check, whether we may
3183 // delete possible grand_children
3184 unsigned int deleted_grandchildren = 0;
3185 unsigned int number_of_child_refinements = 0;
3186
3187 for (unsigned int c = 0; c < 2; ++c)
3188 if (quad->child(c)->has_children())
3189 {
3190 ++number_of_child_refinements;
3191 // if one of the cell counters is
3192 // zero, the other has to be as
3193 // well
3194 Assert(
3195 (quad_cell_count[quad->child(c)->child_index(0)] ==
3196 0 &&
3197 quad_cell_count[quad->child(c)->child_index(1)] ==
3198 0) ||
3199 (quad_cell_count[quad->child(c)->child_index(0)] >
3200 0 &&
3201 quad_cell_count[quad->child(c)->child_index(1)] >
3202 0),
3204 if (quad_cell_count[quad->child(c)->child_index(0)] ==
3205 0)
3206 {
3207 // Assert, that the two
3208 // anisotropic
3209 // refinements add up to
3210 // isotropic refinement
3211 Assert(quad->refinement_case() +
3212 quad->child(c)->refinement_case() ==
3215 // we may delete the
3216 // quad's children and
3217 // the inner line as no
3218 // cell references them
3219 // anymore
3220 quads_to_delete.push_back(
3221 quad->child(c)->child(0));
3222 quads_to_delete.push_back(
3223 quad->child(c)->child(1));
3224 if (quad->child(c)->refinement_case() ==
3226 lines_to_delete.push_back(
3227 quad->child(c)->child(0)->line(1));
3228 else
3229 lines_to_delete.push_back(
3230 quad->child(c)->child(0)->line(3));
3231 quad->child(c)->clear_children();
3232 quad->child(c)->clear_refinement_case();
3233 ++deleted_grandchildren;
3234 }
3235 }
3236 // if no grandchildren are left, we
3237 // may as well delete the
3238 // refinement of the inner line
3239 // between our children and the
3240 // corresponding vertex
3241 if (number_of_child_refinements > 0 &&
3242 deleted_grandchildren == number_of_child_refinements)
3243 {
3245 middle_line;
3246 if (quad->refinement_case() == RefinementCase<2>::cut_x)
3247 middle_line = quad->child(0)->line(1);
3248 else
3249 middle_line = quad->child(0)->line(3);
3250
3251 lines_to_delete.push_back(middle_line->child(0));
3252 lines_to_delete.push_back(middle_line->child(1));
3254 .vertices_used[middle_vertex_index<dim, spacedim>(
3255 middle_line)] = false;
3256 middle_line->clear_children();
3257 }
3258
3259 // now consider the direct children
3260 // of the given quad
3261 if (quad_cell_count[quad->child_index(0)] == 0)
3262 {
3263 // we may delete the quad's
3264 // children and the inner line
3265 // as no cell references them
3266 // anymore
3267 quads_to_delete.push_back(quad->child(0));
3268 quads_to_delete.push_back(quad->child(1));
3269 if (quad->refinement_case() == RefinementCase<2>::cut_x)
3270 lines_to_delete.push_back(quad->child(0)->line(1));
3271 else
3272 lines_to_delete.push_back(quad->child(0)->line(3));
3273
3274 // if the counters just dropped
3275 // to zero, otherwise the
3276 // children would have been
3277 // deleted earlier, then this
3278 // cell's children must have
3279 // contained the anisotropic
3280 // quad children. thus, if
3281 // those have again anisotropic
3282 // children, which are in
3283 // effect isotropic children of
3284 // the original quad, those are
3285 // still needed by a
3286 // neighboring cell and we
3287 // cannot delete them. instead,
3288 // we have to reset this quad's
3289 // refine case to isotropic and
3290 // set the children
3291 // accordingly.
3292 if (quad->child(0)->has_children())
3293 if (quad->refinement_case() ==
3295 {
3296 // now evereything is
3297 // quite complicated. we
3298 // have the children
3299 // numbered according to
3300 //
3301 // *---*---*
3302 // |n+1|m+1|
3303 // *---*---*
3304 // | n | m |
3305 // *---*---*
3306 //
3307 // from the original
3308 // anisotropic
3309 // refinement. we have to
3310 // reorder them as
3311 //
3312 // *---*---*
3313 // | m |m+1|
3314 // *---*---*
3315 // | n |n+1|
3316 // *---*---*
3317 //
3318 // for isotropic refinement.
3319 //
3320 // this is a bit ugly, of
3321 // course: loop over all
3322 // cells on all levels
3323 // and look for faces n+1
3324 // (switch_1) and m
3325 // (switch_2).
3326 const typename Triangulation<dim, spacedim>::
3327 quad_iterator switch_1 =
3328 quad->child(0)->child(1),
3329 switch_2 =
3330 quad->child(1)->child(0);
3331
3332 Assert(!switch_1->has_children(),
3334 Assert(!switch_2->has_children(),
3336
3337 const int switch_1_index = switch_1->index();
3338 const int switch_2_index = switch_2->index();
3339 for (unsigned int l = 0;
3340 l < triangulation.levels.size();
3341 ++l)
3342 for (unsigned int h = 0;
3343 h <
3344 triangulation.levels[l]->cells.n_objects();
3345 ++h)
3346 for (const unsigned int q :
3348 {
3349 const int index =
3350 triangulation.levels[l]
3351 ->cells.get_bounding_object_indices(
3352 h)[q];
3353 if (index == switch_1_index)
3354 triangulation.levels[l]
3355 ->cells.get_bounding_object_indices(
3356 h)[q] = switch_2_index;
3357 else if (index == switch_2_index)
3358 triangulation.levels[l]
3359 ->cells.get_bounding_object_indices(
3360 h)[q] = switch_1_index;
3361 }
3362 // now we have to copy
3363 // all information of the
3364 // two quads
3365 const int switch_1_lines[4] = {
3366 static_cast<signed int>(
3367 switch_1->line_index(0)),
3368 static_cast<signed int>(
3369 switch_1->line_index(1)),
3370 static_cast<signed int>(
3371 switch_1->line_index(2)),
3372 static_cast<signed int>(
3373 switch_1->line_index(3))};
3374 const bool switch_1_line_orientations[4] = {
3375 switch_1->line_orientation(0),
3376 switch_1->line_orientation(1),
3377 switch_1->line_orientation(2),
3378 switch_1->line_orientation(3)};
3379 const types::boundary_id switch_1_boundary_id =
3380 switch_1->boundary_id();
3381 const unsigned int switch_1_user_index =
3382 switch_1->user_index();
3383 const bool switch_1_user_flag =
3384 switch_1->user_flag_set();
3385
3386 switch_1->set_bounding_object_indices(
3387 {switch_2->line_index(0),
3388 switch_2->line_index(1),
3389 switch_2->line_index(2),
3390 switch_2->line_index(3)});
3391 switch_1->set_line_orientation(
3392 0, switch_2->line_orientation(0));
3393 switch_1->set_line_orientation(
3394 1, switch_2->line_orientation(1));
3395 switch_1->set_line_orientation(
3396 2, switch_2->line_orientation(2));
3397 switch_1->set_line_orientation(
3398 3, switch_2->line_orientation(3));
3399 switch_1->set_boundary_id_internal(
3400 switch_2->boundary_id());
3401 switch_1->set_manifold_id(
3402 switch_2->manifold_id());
3403 switch_1->set_user_index(switch_2->user_index());
3404 if (switch_2->user_flag_set())
3405 switch_1->set_user_flag();
3406 else
3407 switch_1->clear_user_flag();
3408
3409 switch_2->set_bounding_object_indices(
3410 {switch_1_lines[0],
3411 switch_1_lines[1],
3412 switch_1_lines[2],
3413 switch_1_lines[3]});
3414 switch_2->set_line_orientation(
3415 0, switch_1_line_orientations[0]);
3416 switch_2->set_line_orientation(
3417 1, switch_1_line_orientations[1]);
3418 switch_2->set_line_orientation(
3419 2, switch_1_line_orientations[2]);
3420 switch_2->set_line_orientation(
3421 3, switch_1_line_orientations[3]);
3422 switch_2->set_boundary_id_internal(
3423 switch_1_boundary_id);
3424 switch_2->set_manifold_id(
3425 switch_1->manifold_id());
3426 switch_2->set_user_index(switch_1_user_index);
3427 if (switch_1_user_flag)
3428 switch_2->set_user_flag();
3429 else
3430 switch_2->clear_user_flag();
3431
3432 const unsigned int child_0 =
3433 quad->child(0)->child_index(0);
3434 const unsigned int child_2 =
3435 quad->child(1)->child_index(0);
3436 quad->clear_children();
3437 quad->clear_refinement_case();
3438 quad->set_refinement_case(
3440 quad->set_children(0, child_0);
3441 quad->set_children(2, child_2);
3442 std::swap(quad_cell_count[child_0 + 1],
3443 quad_cell_count[child_2]);
3444 }
3445 else
3446 {
3447 // the face was refined
3448 // with cut_y, thus the
3449 // children are already
3450 // in correct order. we
3451 // only have to set them
3452 // correctly, deleting
3453 // the indirection of two
3454 // anisotropic refinement
3455 // and going directly
3456 // from the quad to
3457 // isotropic children
3458 const unsigned int child_0 =
3459 quad->child(0)->child_index(0);
3460 const unsigned int child_2 =
3461 quad->child(1)->child_index(0);
3462 quad->clear_children();
3463 quad->clear_refinement_case();
3464 quad->set_refinement_case(
3466 quad->set_children(0, child_0);
3467 quad->set_children(2, child_2);
3468 }
3469 else
3470 {
3471 quad->clear_children();
3472 quad->clear_refinement_case();
3473 }
3474 }
3475 break;
3476 }
3477 case RefinementCase<dim - 1>::cut_xy:
3478 {
3479 // if one of the cell counters is
3480 // zero, the others have to be as
3481 // well
3482
3483 Assert((quad_cell_count[quad->child_index(0)] == 0 &&
3484 quad_cell_count[quad->child_index(1)] == 0 &&
3485 quad_cell_count[quad->child_index(2)] == 0 &&
3486 quad_cell_count[quad->child_index(3)] == 0) ||
3487 (quad_cell_count[quad->child_index(0)] > 0 &&
3488 quad_cell_count[quad->child_index(1)] > 0 &&
3489 quad_cell_count[quad->child_index(2)] > 0 &&
3490 quad_cell_count[quad->child_index(3)] > 0),
3492
3493 if (quad_cell_count[quad->child_index(0)] == 0)
3494 {
3495 // we may delete the quad's
3496 // children, the inner lines
3497 // and the middle vertex as no
3498 // cell references them anymore
3499 lines_to_delete.push_back(quad->child(0)->line(1));
3500 lines_to_delete.push_back(quad->child(3)->line(0));
3501 lines_to_delete.push_back(quad->child(0)->line(3));
3502 lines_to_delete.push_back(quad->child(3)->line(2));
3503
3504 for (unsigned int child = 0; child < quad->n_children();
3505 ++child)
3506 quads_to_delete.push_back(quad->child(child));
3507
3509 .vertices_used[quad->child(0)->vertex_index(3)] =
3510 false;
3511
3512 quad->clear_children();
3513 quad->clear_refinement_case();
3514 }
3515 }
3516 break;
3517
3518 default:
3519 Assert(false, ExcInternalError());
3520 break;
3521 }
3522 }
3523
3524 // now we repeat a similar procedure
3525 // for the outer lines of this cell.
3526
3527 // if in debug mode: check that each
3528 // of the lines for which we consider
3529 // deleting the children in fact has
3530 // children (the bits/coarsening_3d
3531 // test tripped over this initially)
3532 for (unsigned int line_no = 0;
3533 line_no < GeometryInfo<dim>::lines_per_cell;
3534 ++line_no)
3535 {
3537 cell->line(line_no);
3538
3539 Assert(
3540 (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
3541 line->has_children()) ||
3542 GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
3545
3546 if (line->has_children())
3547 {
3548 // if one of the cell counters is
3549 // zero, the other has to be as well
3550
3551 Assert((line_cell_count[line->child_index(0)] == 0 &&
3552 line_cell_count[line->child_index(1)] == 0) ||
3553 (line_cell_count[line->child_index(0)] > 0 &&
3554 line_cell_count[line->child_index(1)] > 0),
3556
3557 if (line_cell_count[line->child_index(0)] == 0)
3558 {
3559 for (unsigned int c = 0; c < 2; ++c)
3560 Assert(!line->child(c)->has_children(),
3562
3563 // we may delete the line's
3564 // children and the middle vertex
3565 // as no cell references them
3566 // anymore
3568 .vertices_used[line->child(0)->vertex_index(1)] = false;
3569
3570 lines_to_delete.push_back(line->child(0));
3571 lines_to_delete.push_back(line->child(1));
3572
3573 line->clear_children();
3574 }
3575 }
3576 }
3577
3578 // finally, delete unneeded quads and lines
3579
3580 // clear user pointers, to avoid that
3581 // they may appear at unwanted places
3582 // later on...
3583 // same for user flags, then finally
3584 // delete the quads and lines
3585 typename std::vector<
3587 line = lines_to_delete.begin(),
3588 endline = lines_to_delete.end();
3589 for (; line != endline; ++line)
3590 {
3591 (*line)->clear_user_data();
3592 (*line)->clear_user_flag();
3593 (*line)->clear_used_flag();
3594 }
3595
3596 typename std::vector<
3598 quad = quads_to_delete.begin(),
3599 endquad = quads_to_delete.end();
3600 for (; quad != endquad; ++quad)
3601 {
3602 (*quad)->clear_user_data();
3603 (*quad)->clear_children();
3604 (*quad)->clear_refinement_case();
3605 (*quad)->clear_user_flag();
3606 (*quad)->clear_used_flag();
3607 }
3608 }
3609
3610
3628 template <int spacedim>
3629 static void
3632 unsigned int & next_unused_vertex,
3634 &next_unused_line,
3636 &next_unused_cell,
3637 const typename Triangulation<2, spacedim>::cell_iterator &cell)
3638 {
3639 const unsigned int dim = 2;
3640 // clear refinement flag
3641 const RefinementCase<dim> ref_case = cell->refine_flag_set();
3642 cell->clear_refine_flag();
3643
3644 /* For the refinement process: since we go the levels up from the
3645 lowest, there are (unlike above) only two possibilities: a neighbor
3646 cell is on the same level or one level up (in both cases, it may or
3647 may not be refined later on, but we don't care here).
3648
3649 First:
3650 Set up an array of the 3x3 vertices, which are distributed on the
3651 cell (the array consists of indices into the @p{vertices} std::vector
3652
3653 2--7--3
3654 | | |
3655 4--8--5
3656 | | |
3657 0--6--1
3658
3659 note: in case of cut_x or cut_y not all these vertices are needed for
3660 the new cells
3661
3662 Second:
3663 Set up an array of the new lines (the array consists of iterator
3664 pointers into the lines arrays)
3665
3666 .-6-.-7-. The directions are: .->-.->-.
3667 1 9 3 ^ ^ ^
3668 .-10.11-. .->-.->-.
3669 0 8 2 ^ ^ ^
3670 .-4-.-5-. .->-.->-.
3671
3672 cut_x:
3673 .-4-.-5-.
3674 | | |
3675 0 6 1
3676 | | |
3677 .-2-.-3-.
3678
3679 cut_y:
3680 .---5---.
3681 1 3
3682 .---6---.
3683 0 2
3684 .---4---.
3685
3686
3687 Third:
3688 Set up an array of neighbors:
3689
3690 6 7
3691 .--.--.
3692 1| | |3
3693 .--.--.
3694 0| | |2
3695 .--.--.
3696 4 5
3697
3698 We need this array for two reasons: first to get the lines which will
3699 bound the four subcells (if the neighboring cell is refined, these
3700 lines already exist), and second to update neighborship information.
3701 Since if a neighbor is not refined, its neighborship record only
3702 points to the present, unrefined, cell rather than the children we
3703 are presently creating, we only need the neighborship information
3704 if the neighbor cells are refined. In all other cases, we store
3705 the unrefined neighbor address
3706
3707 We also need for every neighbor (if refined) which number among its
3708 neighbors the present (unrefined) cell has, since that number is to
3709 be replaced and because that also is the number of the subline which
3710 will be the interface between that neighbor and the to be created
3711 cell. We will store this number (between 0 and 3) in the field
3712 @p{neighbors_neighbor}.
3713
3714 It would be sufficient to use the children of the common line to the
3715 neighbor, if we only wanted to get the new sublines and the new
3716 vertex, but because we need to update the neighborship information of
3717 the two refined subcells of the neighbor, we need to search these
3718 anyway.
3719
3720 Convention:
3721 The created children are numbered like this:
3722
3723 .--.--.
3724 |2 . 3|
3725 .--.--.
3726 |0 | 1|
3727 .--.--.
3728 */
3729 // collect the indices of the eight surrounding vertices
3730 // 2--7--3
3731 // | | |
3732 // 4--8--5
3733 // | | |
3734 // 0--6--1
3735 int new_vertices[9];
3736 for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
3737 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
3738 for (unsigned int line_no = 0; line_no < 4; ++line_no)
3739 if (cell->line(line_no)->has_children())
3740 new_vertices[4 + line_no] =
3741 cell->line(line_no)->child(0)->vertex_index(1);
3742
3743 if (ref_case == RefinementCase<dim>::cut_xy)
3744 {
3745 // find the next
3746 // unused vertex and
3747 // allocate it for
3748 // the new vertex we
3749 // need here
3750 while (triangulation.vertices_used[next_unused_vertex] == true)
3751 ++next_unused_vertex;
3752 Assert(next_unused_vertex < triangulation.vertices.size(),
3753 ExcMessage(
3754 "Internal error: During refinement, the triangulation "
3755 "wants to access an element of the 'vertices' array "
3756 "but it turns out that the array is not large enough."));
3757 triangulation.vertices_used[next_unused_vertex] = true;
3758
3759 new_vertices[8] = next_unused_vertex;
3760
3761 // determine middle vertex by transfinite interpolation to be
3762 // consistent with what happens to quads in a
3763 // Triangulation<3,3> when they are refined
3764 triangulation.vertices[next_unused_vertex] =
3765 cell->center(true, true);
3766 }
3767
3768
3769 // Now the lines:
3771 unsigned int lmin = 8;
3772 unsigned int lmax = 12;
3773 if (ref_case != RefinementCase<dim>::cut_xy)
3774 {
3775 lmin = 6;
3776 lmax = 7;
3777 }
3778
3779 for (unsigned int l = lmin; l < lmax; ++l)
3780 {
3781 while (next_unused_line->used() == true)
3782 ++next_unused_line;
3783 new_lines[l] = next_unused_line;
3784 ++next_unused_line;
3785
3786 AssertIsNotUsed(new_lines[l]);
3787 }
3788
3789 if (ref_case == RefinementCase<dim>::cut_xy)
3790 {
3791 // .-6-.-7-.
3792 // 1 9 3
3793 // .-10.11-.
3794 // 0 8 2
3795 // .-4-.-5-.
3796
3797 // lines 0-7 already exist, create only the four interior
3798 // lines 8-11
3799 unsigned int l = 0;
3800 for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
3801 for (unsigned int c = 0; c < 2; ++c, ++l)
3802 new_lines[l] = cell->line(face_no)->child(c);
3803 Assert(l == 8, ExcInternalError());
3804
3805 new_lines[8]->set_bounding_object_indices(
3806 {new_vertices[6], new_vertices[8]});
3807 new_lines[9]->set_bounding_object_indices(
3808 {new_vertices[8], new_vertices[7]});
3809 new_lines[10]->set_bounding_object_indices(
3810 {new_vertices[4], new_vertices[8]});
3811 new_lines[11]->set_bounding_object_indices(
3812 {new_vertices[8], new_vertices[5]});
3813 }
3814 else if (ref_case == RefinementCase<dim>::cut_x)
3815 {
3816 // .-4-.-5-.
3817 // | | |
3818 // 0 6 1
3819 // | | |
3820 // .-2-.-3-.
3821 new_lines[0] = cell->line(0);
3822 new_lines[1] = cell->line(1);
3823 new_lines[2] = cell->line(2)->child(0);
3824 new_lines[3] = cell->line(2)->child(1);
3825 new_lines[4] = cell->line(3)->child(0);
3826 new_lines[5] = cell->line(3)->child(1);
3827 new_lines[6]->set_bounding_object_indices(
3828 {new_vertices[6], new_vertices[7]});
3829 }
3830 else
3831 {
3833 // .---5---.
3834 // 1 3
3835 // .---6---.
3836 // 0 2
3837 // .---4---.
3838 new_lines[0] = cell->line(0)->child(0);
3839 new_lines[1] = cell->line(0)->child(1);
3840 new_lines[2] = cell->line(1)->child(0);
3841 new_lines[3] = cell->line(1)->child(1);
3842 new_lines[4] = cell->line(2);
3843 new_lines[5] = cell->line(3);
3844 new_lines[6]->set_bounding_object_indices(
3845 {new_vertices[4], new_vertices[5]});
3846 }
3847
3848 for (unsigned int l = lmin; l < lmax; ++l)
3849 {
3850 new_lines[l]->set_used_flag();
3851 new_lines[l]->clear_user_flag();
3852 new_lines[l]->clear_user_data();
3853 new_lines[l]->clear_children();
3854 // interior line
3855 new_lines[l]->set_boundary_id_internal(
3857 new_lines[l]->set_manifold_id(cell->manifold_id());
3858 }
3859
3860 // Now add the four (two)
3861 // new cells!
3864 while (next_unused_cell->used() == true)
3865 ++next_unused_cell;
3866
3867 const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
3868 for (unsigned int i = 0; i < n_children; ++i)
3869 {
3870 AssertIsNotUsed(next_unused_cell);
3871 subcells[i] = next_unused_cell;
3872 ++next_unused_cell;
3873 if (i % 2 == 1 && i < n_children - 1)
3874 while (next_unused_cell->used() == true)
3875 ++next_unused_cell;
3876 }
3877
3878 if (ref_case == RefinementCase<dim>::cut_xy)
3879 {
3880 // children:
3881 // .--.--.
3882 // |2 . 3|
3883 // .--.--.
3884 // |0 | 1|
3885 // .--.--.
3886 // lines:
3887 // .-6-.-7-.
3888 // 1 9 3
3889 // .-10.11-.
3890 // 0 8 2
3891 // .-4-.-5-.
3892 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3893 new_lines[8]->index(),
3894 new_lines[4]->index(),
3895 new_lines[10]->index()});
3896 subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
3897 new_lines[2]->index(),
3898 new_lines[5]->index(),
3899 new_lines[11]->index()});
3900 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
3901 new_lines[9]->index(),
3902 new_lines[10]->index(),
3903 new_lines[6]->index()});
3904 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
3905 new_lines[3]->index(),
3906 new_lines[11]->index(),
3907 new_lines[7]->index()});
3908 }
3909 else if (ref_case == RefinementCase<dim>::cut_x)
3910 {
3911 // children:
3912 // .--.--.
3913 // | . |
3914 // .0 . 1.
3915 // | | |
3916 // .--.--.
3917 // lines:
3918 // .-4-.-5-.
3919 // | | |
3920 // 0 6 1
3921 // | | |
3922 // .-2-.-3-.
3923 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3924 new_lines[6]->index(),
3925 new_lines[2]->index(),
3926 new_lines[4]->index()});
3927 subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
3928 new_lines[1]->index(),
3929 new_lines[3]->index(),
3930 new_lines[5]->index()});
3931 }
3932 else
3933 {
3935 // children:
3936 // .-----.
3937 // | 1 |
3938 // .-----.
3939 // | 0 |
3940 // .-----.
3941 // lines:
3942 // .---5---.
3943 // 1 3
3944 // .---6---.
3945 // 0 2
3946 // .---4---.
3947 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
3948 new_lines[2]->index(),
3949 new_lines[4]->index(),
3950 new_lines[6]->index()});
3951 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
3952 new_lines[3]->index(),
3953 new_lines[6]->index(),
3954 new_lines[5]->index()});
3955 }
3956
3957 types::subdomain_id subdomainid = cell->subdomain_id();
3958
3959 for (unsigned int i = 0; i < n_children; ++i)
3960 {
3961 subcells[i]->set_used_flag();
3962 subcells[i]->clear_refine_flag();
3963 subcells[i]->clear_user_flag();
3964 subcells[i]->clear_user_data();
3965 subcells[i]->clear_children();
3966 // inherit material properties
3967 subcells[i]->set_material_id(cell->material_id());
3968 subcells[i]->set_manifold_id(cell->manifold_id());
3969 subcells[i]->set_subdomain_id(subdomainid);
3970
3971 if (i % 2 == 0)
3972 subcells[i]->set_parent(cell->index());
3973 }
3974
3975
3976
3977 // set child index for even children i=0,2 (0)
3978 for (unsigned int i = 0; i < n_children / 2; ++i)
3979 cell->set_children(2 * i, subcells[2 * i]->index());
3980 // set the refine case
3981 cell->set_refinement_case(ref_case);
3982
3983 // note that the
3984 // refinement flag was
3985 // already cleared at the
3986 // beginning of this function
3987
3988 if (dim < spacedim)
3989 for (unsigned int c = 0; c < n_children; ++c)
3990 cell->child(c)->set_direction_flag(cell->direction_flag());
3991 }
3992
3993
3994
3995 template <int dim, int spacedim>
3998 const bool check_for_distorted_cells)
3999 {
4000 AssertDimension(dim, 2);
4001
4002 // Check whether a new level is needed. We have to check for
4003 // this on the highest level only
4004 for (const auto &cell : triangulation.active_cell_iterators_on_level(
4005 triangulation.levels.size() - 1))
4006 if (cell->refine_flag_set())
4007 {
4008 triangulation.levels.push_back(
4009 std::make_unique<
4011 break;
4012 }
4013
4015 triangulation.begin_line();
4016 line != triangulation.end_line();
4017 ++line)
4018 {
4019 line->clear_user_flag();
4020 line->clear_user_data();
4021 }
4022
4023 unsigned int n_single_lines = 0;
4024 unsigned int n_lines_in_pairs = 0;
4025 unsigned int needed_vertices = 0;
4026
4027 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4028 {
4029 // count number of flagged cells on this level and compute
4030 // how many new vertices and new lines will be needed
4031 unsigned int needed_cells = 0;
4032
4033 for (const auto &cell :
4034 triangulation.active_cell_iterators_on_level(level))
4035 if (cell->refine_flag_set())
4036 {
4037 if (cell->reference_cell() ==
4039 {
4040 needed_cells += 4;
4041 needed_vertices += 0;
4042 n_single_lines += 3;
4043 }
4044 else if (cell->reference_cell() ==
4046 {
4047 needed_cells += 4;
4048 needed_vertices += 1;
4049 n_single_lines += 4;
4050 }
4051 else
4052 {
4054 }
4055
4056 for (const auto line_no : cell->face_indices())
4057 {
4058 auto line = cell->line(line_no);
4059 if (line->has_children() == false)
4060 line->set_user_flag();
4061 }
4062 }
4063
4064
4065 const unsigned int used_cells =
4066 std::count(triangulation.levels[level + 1]->cells.used.begin(),
4067 triangulation.levels[level + 1]->cells.used.end(),
4068 true);
4069
4070
4071 reserve_space(*triangulation.levels[level + 1],
4072 used_cells + needed_cells,
4073 2,
4074 spacedim);
4075
4076 reserve_space(triangulation.levels[level + 1]->cells,
4077 needed_cells,
4078 0);
4079 }
4080
4081 for (auto line = triangulation.begin_line();
4082 line != triangulation.end_line();
4083 ++line)
4084 if (line->user_flag_set())
4085 {
4086 Assert(line->has_children() == false, ExcInternalError());
4087 n_lines_in_pairs += 2;
4088 needed_vertices += 1;
4089 }
4090
4091 reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4092
4093 needed_vertices += std::count(triangulation.vertices_used.begin(),
4094 triangulation.vertices_used.end(),
4095 true);
4096
4097 if (needed_vertices > triangulation.vertices.size())
4098 {
4099 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4100 triangulation.vertices_used.resize(needed_vertices, false);
4101 }
4102
4103 unsigned int next_unused_vertex = 0;
4104
4105 {
4107 line = triangulation.begin_active_line(),
4108 endl = triangulation.end_line();
4110 next_unused_line = triangulation.begin_raw_line();
4111
4112 for (; line != endl; ++line)
4113 if (line->user_flag_set())
4114 {
4115 // this line needs to be refined
4116
4117 // find the next unused vertex and set it
4118 // appropriately
4119 while (triangulation.vertices_used[next_unused_vertex] == true)
4120 ++next_unused_vertex;
4121 Assert(
4122 next_unused_vertex < triangulation.vertices.size(),
4123 ExcMessage(
4124 "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."));
4125 triangulation.vertices_used[next_unused_vertex] = true;
4126
4127 triangulation.vertices[next_unused_vertex] = line->center(true);
4128
4129 bool pair_found = false;
4130 (void)pair_found;
4131 for (; next_unused_line != endl; ++next_unused_line)
4132 if (!next_unused_line->used() &&
4133 !(++next_unused_line)->used())
4134 {
4135 --next_unused_line;
4136 pair_found = true;
4137 break;
4138 }
4139 Assert(pair_found, ExcInternalError());
4140
4141 line->set_children(0, next_unused_line->index());
4142
4144 children[2] = {next_unused_line, ++next_unused_line};
4145
4146 AssertIsNotUsed(children[0]);
4147 AssertIsNotUsed(children[1]);
4148
4149 children[0]->set_bounding_object_indices(
4150 {line->vertex_index(0), next_unused_vertex});
4151 children[1]->set_bounding_object_indices(
4152 {next_unused_vertex, line->vertex_index(1)});
4153
4154 children[0]->set_used_flag();
4155 children[1]->set_used_flag();
4156 children[0]->clear_children();
4157 children[1]->clear_children();
4158 children[0]->clear_user_data();
4159 children[1]->clear_user_data();
4160 children[0]->clear_user_flag();
4161 children[1]->clear_user_flag();
4162
4163
4164 children[0]->set_boundary_id_internal(line->boundary_id());
4165 children[1]->set_boundary_id_internal(line->boundary_id());
4166
4167 children[0]->set_manifold_id(line->manifold_id());
4168 children[1]->set_manifold_id(line->manifold_id());
4169
4170 line->clear_user_flag();
4171 }
4172 }
4173
4174 reserve_space(triangulation.faces->lines, 0, n_single_lines);
4175
4177 cells_with_distorted_children;
4178
4180 next_unused_line = triangulation.begin_raw_line();
4181
4182 const auto create_children = [](auto & triangulation,
4183 unsigned int &next_unused_vertex,
4184 auto & next_unused_line,
4185 auto & next_unused_cell,
4186 const auto & cell) {
4187 const auto ref_case = cell->refine_flag_set();
4188 cell->clear_refine_flag();
4189
4190 unsigned int n_new_vertices = 0;
4191
4192 if (cell->reference_cell() == ::ReferenceCells::Triangle)
4193 n_new_vertices = 6;
4194 else if (cell->reference_cell() ==
4196 n_new_vertices = 9;
4197 else
4199
4200 std::vector<int> new_vertices(n_new_vertices);
4201 for (unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
4202 ++vertex_no)
4203 new_vertices[vertex_no] = cell->vertex_index(vertex_no);
4204 for (unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
4205 if (cell->line(line_no)->has_children())
4206 new_vertices[cell->n_vertices() + line_no] =
4207 cell->line(line_no)->child(0)->vertex_index(1);
4208
4209 if (cell->reference_cell() == ::ReferenceCells::Quadrilateral)
4210 {
4211 while (triangulation.vertices_used[next_unused_vertex] == true)
4212 ++next_unused_vertex;
4213 Assert(
4214 next_unused_vertex < triangulation.vertices.size(),
4215 ExcMessage(
4216 "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."));
4217 triangulation.vertices_used[next_unused_vertex] = true;
4218
4219 new_vertices[8] = next_unused_vertex;
4220
4221 triangulation.vertices[next_unused_vertex] =
4222 cell->center(true, true);
4223 }
4224
4225 std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
4226 12>
4227 new_lines;
4228 unsigned int lmin = 0;
4229 unsigned int lmax = 0;
4230
4231 if (cell->reference_cell() == ::ReferenceCells::Triangle)
4232 {
4233 lmin = 6;
4234 lmax = 9;
4235 }
4236 else if (cell->reference_cell() ==
4238 {
4239 lmin = 8;
4240 lmax = 12;
4241 }
4242 else
4243 {
4245 }
4246
4247 for (unsigned int l = lmin; l < lmax; ++l)
4248 {
4249 while (next_unused_line->used() == true)
4250 ++next_unused_line;
4251 new_lines[l] = next_unused_line;
4252 ++next_unused_line;
4253
4254 AssertIsNotUsed(new_lines[l]);
4255 }
4256
4257 if (true)
4258 {
4259 if (cell->reference_cell() == ::ReferenceCells::Triangle)
4260 {
4261 // add lines in the right order [TODO: clean up]
4262 const auto ref = [&](const unsigned int face_no,
4263 const unsigned int vertex_no) {
4264 if (cell->line(face_no)->child(0)->vertex_index(0) ==
4265 static_cast<unsigned int>(new_vertices[vertex_no]) ||
4266 cell->line(face_no)->child(0)->vertex_index(1) ==
4267 static_cast<unsigned int>(new_vertices[vertex_no]))
4268 {
4269 new_lines[2 * face_no + 0] =
4270 cell->line(face_no)->child(0);
4271 new_lines[2 * face_no + 1] =
4272 cell->line(face_no)->child(1);
4273 }
4274 else
4275 {
4276 new_lines[2 * face_no + 0] =
4277 cell->line(face_no)->child(1);
4278 new_lines[2 * face_no + 1] =
4279 cell->line(face_no)->child(0);
4280 }
4281 };
4282
4283 ref(0, 0);
4284 ref(1, 1);
4285 ref(2, 2);
4286
4287 new_lines[6]->set_bounding_object_indices(
4288 {new_vertices[3], new_vertices[4]});
4289 new_lines[7]->set_bounding_object_indices(
4290 {new_vertices[4], new_vertices[5]});
4291 new_lines[8]->set_bounding_object_indices(
4292 {new_vertices[5], new_vertices[3]});
4293 }
4294 else if (cell->reference_cell() ==
4296 {
4297 unsigned int l = 0;
4298 for (const unsigned int face_no : cell->face_indices())
4299 for (unsigned int c = 0; c < 2; ++c, ++l)
4300 new_lines[l] = cell->line(face_no)->child(c);
4301
4302 new_lines[8]->set_bounding_object_indices(
4303 {new_vertices[6], new_vertices[8]});
4304 new_lines[9]->set_bounding_object_indices(
4305 {new_vertices[8], new_vertices[7]});
4306 new_lines[10]->set_bounding_object_indices(
4307 {new_vertices[4], new_vertices[8]});
4308 new_lines[11]->set_bounding_object_indices(
4309 {new_vertices[8], new_vertices[5]});
4310 }
4311 else
4312 {
4314 }
4315 }
4316
4317
4318 for (unsigned int l = lmin; l < lmax; ++l)
4319 {
4320 new_lines[l]->set_used_flag();
4321 new_lines[l]->clear_user_flag();
4322 new_lines[l]->clear_user_data();
4323 new_lines[l]->clear_children();
4324 // interior line
4325 new_lines[l]->set_boundary_id_internal(
4327 new_lines[l]->set_manifold_id(cell->manifold_id());
4328 }
4329
4332 while (next_unused_cell->used() == true)
4333 ++next_unused_cell;
4334
4335 unsigned int n_children = 0;
4336
4337 if (cell->reference_cell() == ::ReferenceCells::Triangle)
4338 n_children = 4;
4339 else if (cell->reference_cell() ==
4341 n_children = 4;
4342 else
4344
4345 for (unsigned int i = 0; i < n_children; ++i)
4346 {
4347 AssertIsNotUsed(next_unused_cell);
4348 subcells[i] = next_unused_cell;
4349 ++next_unused_cell;
4350 if (i % 2 == 1 && i < n_children - 1)
4351 while (next_unused_cell->used() == true)
4352 ++next_unused_cell;
4353 }
4354
4355 if ((dim == 2) &&
4356 (cell->reference_cell() == ::ReferenceCells::Triangle))
4357 {
4358 subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
4359 new_lines[8]->index(),
4360 new_lines[5]->index()});
4361 subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
4362 new_lines[2]->index(),
4363 new_lines[6]->index()});
4364 subcells[2]->set_bounding_object_indices({new_lines[7]->index(),
4365 new_lines[3]->index(),
4366 new_lines[4]->index()});
4367 subcells[3]->set_bounding_object_indices({new_lines[6]->index(),
4368 new_lines[7]->index(),
4369 new_lines[8]->index()});
4370
4371 // subcell 0
4372
4373 const auto ref = [&](const unsigned int line_no,
4374 const unsigned int vertex_no,
4375 const unsigned int subcell_no,
4376 const unsigned int subcell_line_no) {
4377 if (new_lines[line_no]->vertex_index(1) !=
4378 static_cast<unsigned int>(new_vertices[vertex_no]))
4379 triangulation.levels[subcells[subcell_no]->level()]
4380 ->face_orientations[subcells[subcell_no]->index() *
4382 subcell_line_no] = 0;
4383 };
4384
4385 ref(0, 3, 0, 0);
4386 ref(8, 5, 0, 1);
4387 ref(5, 0, 0, 2);
4388
4389 ref(1, 1, 1, 0);
4390 ref(2, 4, 1, 1);
4391 ref(6, 3, 1, 2);
4392
4393 ref(7, 4, 2, 0);
4394 ref(3, 2, 2, 1);
4395 ref(4, 5, 2, 2);
4396
4397 ref(6, 4, 3, 0);
4398 ref(7, 5, 3, 1);
4399 ref(8, 3, 3, 2);
4400
4401 // triangulation.levels[subcells[1]->level()]->face_orientations[subcells[1]->index()
4402 // * GeometryInfo<2>::faces_per_cell + 2] = 0;
4403 // triangulation.levels[subcells[2]->level()]->face_orientations[subcells[2]->index()
4404 // * GeometryInfo<2>::faces_per_cell + 0] = 0;
4405 }
4406 else if ((dim == 2) && (cell->reference_cell() ==
4408 {
4409 subcells[0]->set_bounding_object_indices(
4410 {new_lines[0]->index(),
4411 new_lines[8]->index(),
4412 new_lines[4]->index(),
4413 new_lines[10]->index()});
4414 subcells[1]->set_bounding_object_indices(
4415 {new_lines[8]->index(),
4416 new_lines[2]->index(),
4417 new_lines[5]->index(),
4418 new_lines[11]->index()});
4419 subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
4420 new_lines[9]->index(),
4421 new_lines[10]->index(),
4422 new_lines[6]->index()});
4423 subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
4424 new_lines[3]->index(),
4425 new_lines[11]->index(),
4426 new_lines[7]->index()});
4427 }
4428 else
4429 {
4431 }
4432
4433 types::subdomain_id subdomainid = cell->subdomain_id();
4434
4435 for (unsigned int i = 0; i < n_children; ++i)
4436 {
4437 subcells[i]->set_used_flag();
4438 subcells[i]->clear_refine_flag();
4439 subcells[i]->clear_user_flag();
4440 subcells[i]->clear_user_data();
4441 subcells[i]->clear_children();
4442 // inherit material
4443 // properties
4444 subcells[i]->set_material_id(cell->material_id());
4445 subcells[i]->set_manifold_id(cell->manifold_id());
4446 subcells[i]->set_subdomain_id(subdomainid);
4447
4448 // TODO: here we assume that all children have the same reference
4449 // cell type as the parent! This is justified for 2D.
4450 triangulation.levels[subcells[i]->level()]
4451 ->reference_cell[subcells[i]->index()] = cell->reference_cell();
4452
4453 if (i % 2 == 0)
4454 subcells[i]->set_parent(cell->index());
4455 }
4456
4457 for (unsigned int i = 0; i < n_children / 2; ++i)
4458 cell->set_children(2 * i, subcells[2 * i]->index());
4459
4460 cell->set_refinement_case(ref_case);
4461
4462 if (dim < spacedim)
4463 for (unsigned int c = 0; c < n_children; ++c)
4464 cell->child(c)->set_direction_flag(cell->direction_flag());
4465 };
4466
4467 for (int level = 0;
4468 level < static_cast<int>(triangulation.levels.size()) - 1;
4469 ++level)
4470 {
4472 next_unused_cell = triangulation.begin_raw(level + 1);
4473
4474 for (const auto &cell :
4475 triangulation.active_cell_iterators_on_level(level))
4476 if (cell->refine_flag_set())
4477 {
4479 next_unused_vertex,
4480 next_unused_line,
4481 next_unused_cell,
4482 cell);
4483
4484 if (cell->reference_cell() ==
4486 check_for_distorted_cells &&
4487 has_distorted_children<dim, spacedim>(cell))
4488 cells_with_distorted_children.distorted_cells.push_back(
4489 cell);
4490
4491 triangulation.signals.post_refinement_on_cell(cell);
4492 }
4493 }
4494
4495 return cells_with_distorted_children;
4496 }
4497
4498
4499
4504 template <int spacedim>
4507 const bool /*check_for_distorted_cells*/)
4508 {
4509 const unsigned int dim = 1;
4510
4511 // Check whether a new level is needed. We have to check for
4512 // this on the highest level only
4513 for (const auto &cell : triangulation.active_cell_iterators_on_level(
4514 triangulation.levels.size() - 1))
4515 if (cell->refine_flag_set())
4516 {
4517 triangulation.levels.push_back(
4518 std::make_unique<
4520 break;
4521 }
4522
4523
4524 // check how much space is needed on every level. We need not
4525 // check the highest level since either - on the highest level
4526 // no cells are flagged for refinement - there are, but
4527 // prepare_refinement added another empty level
4528 unsigned int needed_vertices = 0;
4529 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4530 {
4531 // count number of flagged
4532 // cells on this level
4533 unsigned int flagged_cells = 0;
4534
4535 for (const auto &acell :
4536 triangulation.active_cell_iterators_on_level(level))
4537 if (acell->refine_flag_set())
4538 ++flagged_cells;
4539
4540 // count number of used cells
4541 // on the next higher level
4542 const unsigned int used_cells =
4543 std::count(triangulation.levels[level + 1]->cells.used.begin(),
4544 triangulation.levels[level + 1]->cells.used.end(),
4545 true);
4546
4547 // reserve space for the used_cells cells already existing
4548 // on the next higher level as well as for the
4549 // 2*flagged_cells that will be created on that level
4550 reserve_space(*triangulation.levels[level + 1],
4552 flagged_cells,
4553 1,
4554 spacedim);
4555 // reserve space for 2*flagged_cells new lines on the next
4556 // higher level
4557 reserve_space(triangulation.levels[level + 1]->cells,
4559 flagged_cells,
4560 0);
4561
4562 needed_vertices += flagged_cells;
4563 }
4564
4565 // add to needed vertices how many
4566 // vertices are already in use
4567 needed_vertices += std::count(triangulation.vertices_used.begin(),
4568 triangulation.vertices_used.end(),
4569 true);
4570 // if we need more vertices: create them, if not: leave the
4571 // array as is, since shrinking is not really possible because
4572 // some of the vertices at the end may be in use
4573 if (needed_vertices > triangulation.vertices.size())
4574 {
4575 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4576 triangulation.vertices_used.resize(needed_vertices, false);
4577 }
4578
4579
4580 // Do REFINEMENT on every level; exclude highest level as
4581 // above
4582
4583 // index of next unused vertex
4584 unsigned int next_unused_vertex = 0;
4585
4586 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4587 {
4589 next_unused_cell = triangulation.begin_raw(level + 1);
4590
4591 for (const auto &cell :
4592 triangulation.active_cell_iterators_on_level(level))
4593 if (cell->refine_flag_set())
4594 {
4595 // clear refinement flag
4596 cell->clear_refine_flag();
4597
4598 // search for next unused
4599 // vertex
4600 while (triangulation.vertices_used[next_unused_vertex] ==
4601 true)
4602 ++next_unused_vertex;
4603 Assert(
4604 next_unused_vertex < triangulation.vertices.size(),
4605 ExcMessage(
4606 "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."));
4607
4608 // Now we always ask the cell itself where to put
4609 // the new point. The cell in turn will query the
4610 // manifold object internally.
4611 triangulation.vertices[next_unused_vertex] =
4612 cell->center(true);
4613
4614 triangulation.vertices_used[next_unused_vertex] = true;
4615
4616 // search for next two unused cell (++ takes care of
4617 // the end of the vector)
4619 first_child,
4620 second_child;
4621 while (next_unused_cell->used() == true)
4622 ++next_unused_cell;
4623 first_child = next_unused_cell;
4624 first_child->set_used_flag();
4625 first_child->clear_user_data();
4626 ++next_unused_cell;
4627 AssertIsNotUsed(next_unused_cell);
4628 second_child = next_unused_cell;
4629 second_child->set_used_flag();
4630 second_child->clear_user_data();
4631
4632 types::subdomain_id subdomainid = cell->subdomain_id();
4633
4634 // insert first child
4635 cell->set_children(0, first_child->index());
4636 first_child->clear_children();
4637 first_child->set_bounding_object_indices(
4638 {cell->vertex_index(0), next_unused_vertex});
4639 first_child->set_material_id(cell->material_id());
4640 first_child->set_manifold_id(cell->manifold_id());
4641 first_child->set_subdomain_id(subdomainid);
4642 first_child->set_direction_flag(cell->direction_flag());
4643
4644 first_child->set_parent(cell->index());
4645
4646 // Set manifold id of the right face. Only do this
4647 // on the first child.
4648 first_child->face(1)->set_manifold_id(cell->manifold_id());
4649
4650 // reset neighborship info (refer to
4651 // internal::TriangulationImplementation::TriaLevel<0> for
4652 // details)
4653 first_child->set_neighbor(1, second_child);
4654 if (cell->neighbor(0).state() != IteratorState::valid)
4655 first_child->set_neighbor(0, cell->neighbor(0));
4656 else if (cell->neighbor(0)->is_active())
4657 {
4658 // since the neighbors level is always <=level,
4659 // if the cell is active, then there are no
4660 // cells to the left which may want to know
4661 // about this new child cell.
4662 Assert(cell->neighbor(0)->level() <= cell->level(),
4664 first_child->set_neighbor(0, cell->neighbor(0));
4665 }
4666 else
4667 // left neighbor is refined
4668 {
4669 // set neighbor to cell on same level
4670 const unsigned int nbnb = cell->neighbor_of_neighbor(0);
4671 first_child->set_neighbor(0,
4672 cell->neighbor(0)->child(nbnb));
4673
4674 // reset neighbor info of all right descendant
4675 // of the left neighbor of cell
4677 left_neighbor = cell->neighbor(0);
4678 while (left_neighbor->has_children())
4679 {
4680 left_neighbor = left_neighbor->child(nbnb);
4681 left_neighbor->set_neighbor(nbnb, first_child);
4682 }
4683 }
4684
4685 // insert second child
4686 second_child->clear_children();
4687 second_child->set_bounding_object_indices(
4688 {next_unused_vertex, cell->vertex_index(1)});
4689 second_child->set_neighbor(0, first_child);
4690 second_child->set_material_id(cell->material_id());
4691 second_child->set_manifold_id(cell->manifold_id());
4692 second_child->set_subdomain_id(subdomainid);
4693 second_child->set_direction_flag(cell->direction_flag());
4694
4695 if (cell->neighbor(1).state() != IteratorState::valid)
4696 second_child->set_neighbor(1, cell->neighbor(1));
4697 else if (cell->neighbor(1)->is_active())
4698 {
4699 Assert(cell->neighbor(1)->level() <= cell->level(),
4701 second_child->set_neighbor(1, cell->neighbor(1));
4702 }
4703 else
4704 // right neighbor is refined same as above
4705 {
4706 const unsigned int nbnb = cell->neighbor_of_neighbor(1);
4707 second_child->set_neighbor(
4708 1, cell->neighbor(1)->child(nbnb));
4709
4711 right_neighbor = cell->neighbor(1);
4712 while (right_neighbor->has_children())
4713 {
4714 right_neighbor = right_neighbor->child(nbnb);
4715 right_neighbor->set_neighbor(nbnb, second_child);
4716 }
4717 }
4718 // inform all listeners that cell refinement is done
4719 triangulation.signals.post_refinement_on_cell(cell);
4720 }
4721 }
4722
4723 // in 1d, we can not have distorted children unless the parent
4724 // was already distorted (that is because we don't use
4725 // boundary information for 1d triangulations). so return an
4726 // empty list
4728 }
4729
4730
4735 template <int spacedim>
4738 const bool check_for_distorted_cells)
4739 {
4740 const unsigned int dim = 2;
4741
4742
4743 // First check whether we can get away with isotropic refinement, or
4744 // whether we need to run through the full anisotropic algorithm
4745 {
4746 bool do_isotropic_refinement = true;
4747 for (const auto &cell : triangulation.active_cell_iterators())
4748 if (cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
4749 cell->refine_flag_set() == RefinementCase<dim>::cut_y)
4750 {
4751 do_isotropic_refinement = false;
4752 break;
4753 }
4754
4755 if (do_isotropic_refinement)
4757 check_for_distorted_cells);
4758 }
4759
4760 // Check whether a new level is needed. We have to check for
4761 // this on the highest level only
4762 for (const auto &cell : triangulation.active_cell_iterators_on_level(
4763 triangulation.levels.size() - 1))
4764 if (cell->refine_flag_set())
4765 {
4766 triangulation.levels.push_back(
4767 std::make_unique<
4769 break;
4770 }
4771
4772 // TODO[WB]: we clear user flags and pointers of lines; we're going
4773 // to use them to flag which lines need refinement
4775 triangulation.begin_line();
4776 line != triangulation.end_line();
4777 ++line)
4778 {
4779 line->clear_user_flag();
4780 line->clear_user_data();
4781 }
4782 // running over all cells and lines count the number
4783 // n_single_lines of lines which can be stored as single
4784 // lines, e.g. inner lines
4785 unsigned int n_single_lines = 0;
4786
4787 // New lines to be created: number lines which are stored in
4788 // pairs (the children of lines must be stored in pairs)
4789 unsigned int n_lines_in_pairs = 0;
4790
4791 // check how much space is needed on every level. We need not
4792 // check the highest level since either - on the highest level
4793 // no cells are flagged for refinement - there are, but
4794 // prepare_refinement added another empty level
4795 unsigned int needed_vertices = 0;
4796 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
4797 {
4798 // count number of flagged cells on this level and compute
4799 // how many new vertices and new lines will be needed
4800 unsigned int needed_cells = 0;
4801
4802 for (const auto &cell :
4803 triangulation.active_cell_iterators_on_level(level))
4804 if (cell->refine_flag_set())
4805 {
4806 if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
4807 {
4808 needed_cells += 4;
4809
4810 // new vertex at center of cell is needed in any
4811 // case
4812 ++needed_vertices;
4813
4814 // the four inner lines can be stored as singles
4815 n_single_lines += 4;
4816 }
4817 else // cut_x || cut_y
4818 {
4819 // set the flag showing that anisotropic
4820 // refinement is used for at least one cell
4821 triangulation.anisotropic_refinement = true;
4822
4823 needed_cells += 2;
4824 // no vertex at center
4825
4826 // the inner line can be stored as single
4827 n_single_lines += 1;
4828 }
4829
4830 // mark all faces (lines) for refinement; checking
4831 // locally whether the neighbor would also like to
4832 // refine them is rather difficult for lines so we
4833 // only flag them and after visiting all cells, we
4834 // decide which lines need refinement;
4835 for (const unsigned int line_no :
4837 {
4839 cell->refine_flag_set(), line_no) ==
4841 {
4843 line = cell->line(line_no);
4844 if (line->has_children() == false)
4845 line->set_user_flag();
4846 }
4847 }
4848 }
4849
4850
4851 // count number of used cells on the next higher level
4852 const unsigned int used_cells =
4853 std::count(triangulation.levels[level + 1]->cells.used.begin(),
4854 triangulation.levels[level + 1]->cells.used.end(),
4855 true);
4856
4857
4858 // reserve space for the used_cells cells already existing
4859 // on the next higher level as well as for the
4860 // needed_cells that will be created on that level
4861 reserve_space(*triangulation.levels[level + 1],
4862 used_cells + needed_cells,
4863 2,
4864 spacedim);
4865
4866 // reserve space for needed_cells new quads on the next
4867 // higher level
4868 reserve_space(triangulation.levels[level + 1]->cells,
4869 needed_cells,
4870 0);
4871 }
4872
4873 // now count the lines which were flagged for refinement
4875 triangulation.begin_line();
4876 line != triangulation.end_line();
4877 ++line)
4878 if (line->user_flag_set())
4879 {
4880 Assert(line->has_children() == false, ExcInternalError());
4881 n_lines_in_pairs += 2;
4882 needed_vertices += 1;
4883 }
4884 // reserve space for n_lines_in_pairs new lines. note, that
4885 // we can't reserve space for the single lines here as well,
4886 // as all the space reserved for lines in pairs would be
4887 // counted as unused and we would end up with too little space
4888 // to store all lines. memory reservation for n_single_lines
4889 // can only be done AFTER we refined the lines of the current
4890 // cells
4891 reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
4892
4893 // add to needed vertices how many vertices are already in use
4894 needed_vertices += std::count(triangulation.vertices_used.begin(),
4895 triangulation.vertices_used.end(),
4896 true);
4897 // if we need more vertices: create them, if not: leave the
4898 // array as is, since shrinking is not really possible because
4899 // some of the vertices at the end may be in use
4900 if (needed_vertices > triangulation.vertices.size())
4901 {
4902 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
4903 triangulation.vertices_used.resize(needed_vertices, false);
4904 }
4905
4906
4907 // Do REFINEMENT on every level; exclude highest level as
4908 // above
4909
4910 // index of next unused vertex
4911 unsigned int next_unused_vertex = 0;
4912
4913 // first the refinement of lines. children are stored
4914 // pairwise
4915 {
4916 // only active objects can be refined further
4918 line = triangulation.begin_active_line(),
4919 endl = triangulation.end_line();
4921 next_unused_line = triangulation.begin_raw_line();
4922
4923 for (; line != endl; ++line)
4924 if (line->user_flag_set())
4925 {
4926 // this line needs to be refined
4927
4928 // find the next unused vertex and set it
4929 // appropriately
4930 while (triangulation.vertices_used[next_unused_vertex] == true)
4931 ++next_unused_vertex;
4932 Assert(
4933 next_unused_vertex < triangulation.vertices.size(),
4934 ExcMessage(
4935 "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."));
4936 triangulation.vertices_used[next_unused_vertex] = true;
4937
4938 triangulation.vertices[next_unused_vertex] = line->center(true);
4939
4940 // now that we created the right point, make up the
4941 // two child lines. To this end, find a pair of
4942 // unused lines
4943 bool pair_found = false;
4944 (void)pair_found;
4945 for (; next_unused_line != endl; ++next_unused_line)
4946 if (!next_unused_line->used() &&
4947 !(++next_unused_line)->used())
4948 {
4949 // go back to the first of the two unused
4950 // lines
4951 --next_unused_line;
4952 pair_found = true;
4953 break;
4954 }
4955 Assert(pair_found, ExcInternalError());
4956
4957 // there are now two consecutive unused lines, such
4958 // that the children of a line will be consecutive.
4959 // then set the child pointer of the present line
4960 line->set_children(0, next_unused_line->index());
4961
4962 // set the two new lines
4964 children[2] = {next_unused_line, ++next_unused_line};
4965 // some tests; if any of the iterators should be
4966 // invalid, then already dereferencing will fail
4967 AssertIsNotUsed(children[0]);
4968 AssertIsNotUsed(children[1]);
4969
4970 children[0]->set_bounding_object_indices(
4971 {line->vertex_index(0), next_unused_vertex});
4972 children[1]->set_bounding_object_indices(
4973 {next_unused_vertex, line->vertex_index(1)});
4974
4975 children[0]->set_used_flag();
4976 children[1]->set_used_flag();
4977 children[0]->clear_children();
4978 children[1]->clear_children();
4979 children[0]->clear_user_data();
4980 children[1]->clear_user_data();
4981 children[0]->clear_user_flag();
4982 children[1]->clear_user_flag();
4983
4984
4985 children[0]->set_boundary_id_internal(line->boundary_id());
4986 children[1]->set_boundary_id_internal(line->boundary_id());
4987
4988 children[0]->set_manifold_id(line->manifold_id());
4989 children[1]->set_manifold_id(line->manifold_id());
4990
4991 // finally clear flag indicating the need for
4992 // refinement
4993 line->clear_user_flag();
4994 }
4995 }
4996
4997
4998 // Now set up the new cells
4999
5000 // reserve space for inner lines (can be stored as single
5001 // lines)
5002 reserve_space(triangulation.faces->lines, 0, n_single_lines);
5003
5005 cells_with_distorted_children;
5006
5007 // reset next_unused_line, as now also single empty places in
5008 // the vector can be used
5010 next_unused_line = triangulation.begin_raw_line();
5011
5012 for (int level = 0;
5013 level < static_cast<int>(triangulation.levels.size()) - 1;
5014 ++level)
5015 {
5017 next_unused_cell = triangulation.begin_raw(level + 1);
5018
5019 for (const auto &cell :
5020 triangulation.active_cell_iterators_on_level(level))
5021 if (cell->refine_flag_set())
5022 {
5023 // actually set up the children and update neighbor
5024 // information
5026 next_unused_vertex,
5027 next_unused_line,
5028 next_unused_cell,
5029 cell);
5030
5031 if (check_for_distorted_cells &&
5032 has_distorted_children<dim, spacedim>(cell))
5033 cells_with_distorted_children.distorted_cells.push_back(
5034 cell);
5035 // inform all listeners that cell refinement is done
5036 triangulation.signals.post_refinement_on_cell(cell);
5037 }
5038 }
5039
5040 return cells_with_distorted_children;
5041 }
5042
5043
5044 template <int spacedim>
5047 const bool check_for_distorted_cells)
5048 {
5049 static const int dim = 3;
5050 static const unsigned int X = numbers::invalid_unsigned_int;
5051
5052 Assert(spacedim == 3, ExcNotImplemented());
5053
5054 Assert(triangulation.vertices.size() ==
5055 triangulation.vertices_used.size(),
5057
5058 // Check whether a new level is needed. We have to check for
5059 // this on the highest level only
5060 for (const auto &cell : triangulation.active_cell_iterators_on_level(
5061 triangulation.levels.size() - 1))
5062 if (cell->refine_flag_set())
5063 {
5064 triangulation.levels.push_back(
5065 std::make_unique<
5067 break;
5068 }
5069
5070 // first clear user flags for quads and lines; we're going to
5071 // use them to flag which lines and quads need refinement
5072 triangulation.faces->quads.clear_user_data();
5073
5075 triangulation.begin_line();
5076 line != triangulation.end_line();
5077 ++line)
5078 line->clear_user_flag();
5079
5081 triangulation.begin_quad();
5082 quad != triangulation.end_quad();
5083 ++quad)
5084 quad->clear_user_flag();
5085
5086 // check how much space is needed on every level. We need not
5087 // check the highest level since either
5088 // - on the highest level no cells are flagged for refinement
5089 // - there are, but prepare_refinement added another empty
5090 // level which then is the highest level
5091
5092 // variables to hold the number of newly to be created
5093 // vertices, lines and quads. as these are stored globally,
5094 // declare them outside the loop over al levels. we need lines
5095 // and quads in pairs for refinement of old ones and lines and
5096 // quads, that can be stored as single ones, as they are newly
5097 // created in the inside of an existing cell
5098 unsigned int needed_vertices = 0;
5099 unsigned int needed_lines_single = 0;
5100 unsigned int needed_quads_single = 0;
5101 unsigned int needed_lines_pair = 0;
5102 unsigned int needed_quads_pair = 0;
5103 for (int level = triangulation.levels.size() - 2; level >= 0; --level)
5104 {
5105 unsigned int new_cells = 0;
5106
5107 for (const auto &cell :
5108 triangulation.active_cell_iterators_on_level(level))
5109 if (cell->refine_flag_set())
5110 {
5111 // Only support isotropic refinement
5112 Assert(cell->refine_flag_set() ==
5115
5116 // Now count up how many new cells, faces, edges, and vertices
5117 // we will need to allocate to do this refinement.
5118 new_cells += cell->reference_cell().n_isotropic_children();
5119
5120 if (cell->reference_cell() == ReferenceCells::Hexahedron)
5121 {
5122 ++needed_vertices;
5123 needed_lines_single += 6;
5124 needed_quads_single += 12;
5125 }
5126 else if (cell->reference_cell() ==
5128 {
5129 needed_lines_single += 1;
5130 needed_quads_single += 8;
5131 }
5132 else
5133 {
5134 Assert(false, ExcInternalError());
5135 }
5136
5137 // Also check whether we have to refine any of the faces and
5138 // edges that bound this cell. They may of course already be
5139 // refined, so we only *mark* them for refinement by setting
5140 // the user flags
5141 for (const auto face : cell->face_indices())
5142 if (cell->face(face)->n_children() == 0)
5143 cell->face(face)->set_user_flag();
5144 else
5145 Assert(cell->face(face)->n_children() ==
5146 cell->reference_cell()
5147 .face_reference_cell(face)
5148 .n_isotropic_children(),
5150
5151 for (const auto line : cell->line_indices())
5152 if (cell->line(line)->has_children() == false)
5153 cell->line(line)->set_user_flag();
5154 else
5155 Assert(cell->line(line)->n_children() == 2,
5157 }
5158
5159 const unsigned int used_cells =
5160 std::count(triangulation.levels[level + 1]->cells.used.begin(),
5161 triangulation.levels[level + 1]->cells.used.end(),
5162 true);
5163
5164 reserve_space(*triangulation.levels[level + 1],
5165 used_cells + new_cells,
5166 3,
5167 spacedim);
5168
5169 reserve_space(triangulation.levels[level + 1]->cells, new_cells);
5170 }
5171
5172 // now count the quads and lines which were flagged for
5173 // refinement
5175 triangulation.begin_quad();
5176 quad != triangulation.end_quad();
5177 ++quad)
5178 {
5179 if (quad->user_flag_set() == false)
5180 continue;
5181
5182 if (quad->reference_cell() == ReferenceCells::Quadrilateral)
5183 {
5184 needed_quads_pair += 4;
5185 needed_lines_pair += 4;
5186 needed_vertices += 1;
5187 }
5188 else if (quad->reference_cell() == ReferenceCells::Triangle)
5189 {
5190 needed_quads_pair += 4;
5191 needed_lines_single += 3;
5192 }
5193 else
5194 {
5195 Assert(false, ExcInternalError());
5196 }
5197 }
5198
5200 triangulation.begin_line();
5201 line != triangulation.end_line();
5202 ++line)
5203 {
5204 if (line->user_flag_set() == false)
5205 continue;
5206
5207 needed_lines_pair += 2;
5208 needed_vertices += 1;
5209 }
5210
5211 reserve_space(triangulation.faces->lines,
5212 needed_lines_pair,
5213 needed_lines_single);
5215 needed_quads_pair,
5216 needed_quads_single);
5217 reserve_space(triangulation.faces->quads,
5218 needed_quads_pair,
5219 needed_quads_single);
5220
5221
5222 // add to needed vertices how many vertices are already in use
5223 needed_vertices += std::count(triangulation.vertices_used.begin(),
5224 triangulation.vertices_used.end(),
5225 true);
5226
5227 if (needed_vertices > triangulation.vertices.size())
5228 {
5229 triangulation.vertices.resize(needed_vertices, Point<spacedim>());
5230 triangulation.vertices_used.resize(needed_vertices, false);
5231 }
5232
5233 //-----------------------------------------
5234 // Before we start with the actual refinement, we do some
5235 // sanity checks if in debug mode. especially, we try to catch
5236 // the notorious problem with lines being twice refined,
5237 // i.e. there are cells adjacent at one line ("around the
5238 // edge", but not at a face), with two cells differing by more
5239 // than one refinement level
5240 //
5241 // this check is very simple to implement here, since we have
5242 // all lines flagged if they shall be refined
5243#ifdef DEBUG
5244 for (const auto &cell : triangulation.active_cell_iterators())
5245 if (!cell->refine_flag_set())
5246 for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
5247 if (cell->line(line_n)->has_children())
5248 for (unsigned int c = 0; c < 2; ++c)
5249 Assert(cell->line(line_n)->child(c)->user_flag_set() == false,
5251#endif
5252
5253 unsigned int current_vertex = 0;
5254
5255 // helper function - find the next available vertex number and mark it
5256 // as used.
5257 auto get_next_unused_vertex = [](const unsigned int current_vertex,
5258 std::vector<bool> &vertices_used) {
5259 unsigned int next_vertex = current_vertex;
5260 while (next_vertex < vertices_used.size() &&
5261 vertices_used[next_vertex] == true)
5262 ++next_vertex;
5263 Assert(next_vertex < vertices_used.size(), ExcInternalError());
5264 vertices_used[next_vertex] = true;
5265
5266 return next_vertex;
5267 };
5268
5269 // LINES
5270 {
5272 line = triangulation.begin_active_line(),
5273 endl = triangulation.end_line();
5275 next_unused_line = triangulation.begin_raw_line();
5276
5277 for (; line != endl; ++line)
5278 {
5279 if (line->user_flag_set() == false)
5280 continue;
5281
5282 current_vertex =
5283 get_next_unused_vertex(current_vertex,
5284 triangulation.vertices_used);
5285 triangulation.vertices[current_vertex] = line->center(true);
5286
5287 next_unused_line =
5288 triangulation.faces->lines.template next_free_pair_object<1>(
5290 Assert(next_unused_line.state() == IteratorState::valid,
5292
5293 // now we found two consecutive unused lines, such
5294 // that the children of a line will be consecutive.
5295 // then set the child pointer of the present line
5296 line->set_children(0, next_unused_line->index());
5297
5299 children[2] = {next_unused_line, ++next_unused_line};
5300
5301 AssertIsNotUsed(children[0]);
5302 AssertIsNotUsed(children[1]);
5303
5304 children[0]->set_bounding_object_indices(
5305 {line->vertex_index(0), current_vertex});
5306 children[1]->set_bounding_object_indices(
5307 {current_vertex, line->vertex_index(1)});
5308
5309 children[0]->set_used_flag();
5310 children[1]->set_used_flag();
5311 children[0]->clear_children();
5312 children[1]->clear_children();
5313 children[0]->clear_user_data();
5314 children[1]->clear_user_data();
5315 children[0]->clear_user_flag();
5316 children[1]->clear_user_flag();
5317
5318 children[0]->set_boundary_id_internal(line->boundary_id());
5319 children[1]->set_boundary_id_internal(line->boundary_id());
5320
5321 children[0]->set_manifold_id(line->manifold_id());
5322 children[1]->set_manifold_id(line->manifold_id());
5323
5324 line->clear_user_flag();
5325 }
5326 }
5327
5328 // QUADS
5329 {
5331 quad = triangulation.begin_quad(),
5332 endq = triangulation.end_quad();
5334 next_unused_line = triangulation.begin_raw_line();
5336 next_unused_quad = triangulation.begin_raw_quad();
5337
5338 for (; quad != endq; ++quad)
5339 {
5340 if (quad->user_flag_set() == false)
5341 continue;
5342
5343 const auto reference_face_type = quad->reference_cell();
5344
5345 // 1) create new vertex (at the center of the face)
5346 if (reference_face_type == ReferenceCells::Quadrilateral)
5347 {
5348 current_vertex =
5349 get_next_unused_vertex(current_vertex,
5350 triangulation.vertices_used);
5351 triangulation.vertices[current_vertex] =
5352 quad->center(true, true);
5353 }
5354
5355 // 2) create new lines (property is set later)
5356 boost::container::small_vector<
5359 new_lines(quad->n_lines());
5360 {
5361 for (unsigned int i = 0; i < new_lines.size(); ++i)
5362 {
5363 if (reference_face_type == ReferenceCells::Quadrilateral)
5364 {
5365 if (i % 2 == 0)
5366 next_unused_line =
5367 triangulation.faces->lines
5368 .template next_free_pair_object<1>(triangulation);
5369 }
5370 else if (reference_face_type == ReferenceCells::Triangle)
5371 {
5372 next_unused_line =
5373 triangulation.faces->lines
5374 .template next_free_single_object<1>(triangulation);
5375 }
5376 else
5377 {
5378 Assert(false, ExcNotImplemented());
5379 }
5380
5381 new_lines[i] = next_unused_line;
5382 ++next_unused_line;
5383 AssertIsNotUsed(new_lines[i]);
5384 }
5385 }
5386
5387 // 3) create new quads (properties are set below). Both triangles
5388 // and quads are divided in four.
5389 std::array<
5391 4>
5392 new_quads;
5393 {
5394 next_unused_quad =
5395 triangulation.faces->quads.template next_free_pair_object<2>(
5397
5398 new_quads[0] = next_unused_quad;
5399 AssertIsNotUsed(new_quads[0]);
5400
5401 ++next_unused_quad;
5402 new_quads[1] = next_unused_quad;
5403 AssertIsNotUsed(new_quads[1]);
5404
5405 next_unused_quad =
5406 triangulation.faces->quads.template next_free_pair_object<2>(
5408 new_quads[2] = next_unused_quad;
5409 AssertIsNotUsed(new_quads[2]);
5410
5411 ++next_unused_quad;
5412 new_quads[3] = next_unused_quad;
5413 AssertIsNotUsed(new_quads[3]);
5414
5415 quad->set_children(0, new_quads[0]->index());
5416 quad->set_children(2, new_quads[2]->index());
5417 quad->set_refinement_case(RefinementCase<2>::cut_xy);
5418 }
5419
5420 // Maximum of 9 vertices per refined quad (9 for Quadrilateral, 6
5421 // for Triangle)
5422 std::array<unsigned int, 9> vertex_indices = {};
5423 {
5424 unsigned int k = 0;
5425 for (const auto i : quad->vertex_indices())
5426 vertex_indices[k++] = quad->vertex_index(i);
5427
5428 for (const auto i : quad->line_indices())
5429 vertex_indices[k++] =
5430 quad->line(i)->child(0)->vertex_index(1);
5431
5432 vertex_indices[k++] = current_vertex;
5433 }
5434
5435 boost::container::small_vector<
5437 12>
5438 lines(reference_face_type == ReferenceCells::Quadrilateral ?
5439 12 :
5440 9);
5441 {
5442 unsigned int k = 0;
5443
5444 for (unsigned int l = 0; l < quad->n_lines(); ++l)
5445 for (unsigned int c = 0; c < 2; ++c)
5446 {
5447 static constexpr std::array<std::array<unsigned int, 2>,
5448 2>
5449 index = {// child 0, line_orientation=false and true
5450 {{{1, 0}},
5451 // child 1, line_orientation=false and true
5452 {{0, 1}}}};
5453
5454 lines[k++] = quad->line(l)->child(
5455 index[c][quad->line_orientation(l)]);
5456 }
5457
5458 for (unsigned int l = 0; l < new_lines.size(); ++l)
5459 lines[k++] = new_lines[l];
5460 }
5461
5462 boost::container::small_vector<int, 12> line_indices(
5463 lines.size());
5464 for (unsigned int i = 0; i < line_indices.size(); ++i)
5465 line_indices[i] = lines[i]->index();
5466
5467 static constexpr std::array<std::array<unsigned int, 2>, 12>
5468 line_vertices_quad{{{{0, 4}},
5469 {{4, 2}},
5470 {{1, 5}},
5471 {{5, 3}},
5472 {{0, 6}},
5473 {{6, 1}},
5474 {{2, 7}},
5475 {{7, 3}},
5476 {{6, 8}},
5477 {{8, 7}},
5478 {{4, 8}},
5479 {{8, 5}}}};
5480
5481 static constexpr std::array<std::array<unsigned int, 4>, 4>
5482 quad_lines_quad{{{{0, 8, 4, 10}},
5483 {{8, 2, 5, 11}},
5484 {{1, 9, 10, 6}},
5485 {{9, 3, 11, 7}}}};
5486
5487 static constexpr std::
5488 array<std::array<std::array<unsigned int, 2>, 4>, 4>
5489 quad_line_vertices_quad{
5490 {{{{{0, 4}}, {{6, 8}}, {{0, 6}}, {{4, 8}}}},
5491 {{{{6, 8}}, {{1, 5}}, {{6, 1}}, {{8, 5}}}},
5492 {{{{4, 2}}, {{8, 7}}, {{4, 8}}, {{2, 7}}}},
5493 {{{{8, 7}}, {{5, 3}}, {{8, 5}}, {{7, 3}}}}}};
5494
5495 static constexpr std::array<std::array<unsigned int, 2>, 12>
5496 line_vertices_tri{{{{0, 3}},
5497 {{3, 1}},
5498 {{1, 4}},
5499 {{4, 2}},
5500 {{2, 5}},
5501 {{5, 0}},
5502 {{3, 4}},
5503 {{4, 5}},
5504 {{3, 5}},
5505 {{X, X}},
5506 {{X, X}},
5507 {{X, X}}}};
5508
5509 static constexpr std::array<std::array<unsigned int, 4>, 4>
5510 quad_lines_tri{{{{0, 8, 5, X}},
5511 {{1, 2, 6, X}},
5512 {{7, 3, 4, X}},
5513 {{6, 7, 8, X}}}};
5514
5515 static constexpr std::
5516 array<std::array<std::array<unsigned int, 2>, 4>, 4>
5517 quad_line_vertices_tri{
5518 {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}},
5519 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
5520 {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}},
5521 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}};
5522
5523 const auto &line_vertices =
5524 (reference_face_type == ReferenceCells::Quadrilateral) ?
5525 line_vertices_quad :
5526 line_vertices_tri;
5527 const auto &quad_lines =
5528 (reference_face_type == ReferenceCells::Quadrilateral) ?
5529 quad_lines_quad :
5530 quad_lines_tri;
5531 const auto &quad_line_vertices =
5532 (reference_face_type == ReferenceCells::Quadrilateral) ?
5533 quad_line_vertices_quad :
5534 quad_line_vertices_tri;
5535
5536 // 4) set properties of lines
5537 for (unsigned int i = 0, j = lines.size() - new_lines.size();
5538 i < new_lines.size();
5539 ++i, ++j)
5540 {
5541 auto &new_line = new_lines[i];
5542 new_line->set_bounding_object_indices(
5543 {vertex_indices[line_vertices[j][0]],
5544 vertex_indices[line_vertices[j][1]]});
5545 new_line->set_used_flag();
5546 new_line->clear_user_flag();
5547 new_line->clear_user_data();
5548 new_line->clear_children();
5549 new_line->set_boundary_id_internal(quad->boundary_id());
5550 new_line->set_manifold_id(quad->manifold_id());
5551 }
5552
5553 // 5) set properties of quads
5554 for (unsigned int i = 0; i < new_quads.size(); ++i)
5555 {
5556 auto &new_quad = new_quads[i];
5557
5558 // TODO: we assume here that all children have the same type
5559 // as the parent
5560 triangulation.faces->quad_reference_cell[new_quad->index()] =
5561 reference_face_type;
5562
5563 if (new_quad->n_lines() == 3)
5564 new_quad->set_bounding_object_indices(
5565 {line_indices[quad_lines[i][0]],
5566 line_indices[quad_lines[i][1]],
5567 line_indices[quad_lines[i][2]]});
5568 else if (new_quad->n_lines() == 4)
5569 new_quad->set_bounding_object_indices(
5570 {line_indices[quad_lines[i][0]],
5571 line_indices[quad_lines[i][1]],
5572 line_indices[quad_lines[i][2]],
5573 line_indices[quad_lines[i][3]]});
5574 else
5575 Assert(false, ExcNotImplemented());
5576
5577 new_quad->set_used_flag();
5578 new_quad->clear_user_flag();
5579 new_quad->clear_user_data();
5580 new_quad->clear_children();
5581 new_quad->set_boundary_id_internal(quad->boundary_id());
5582 new_quad->set_manifold_id(quad->manifold_id());
5583
5584#ifdef DEBUG
5585 std::set<unsigned int> s;
5586#endif
5587
5588 // ... and fix orientation of faces (lines) of quad
5589 for (const auto f : new_quad->line_indices())
5590 {
5591 std::array<unsigned int, 2> vertices_0, vertices_1;
5592
5593 for (unsigned int v = 0; v < 2; ++v)
5594 vertices_0[v] =
5595 lines[quad_lines[i][f]]->vertex_index(v);
5596
5597 for (unsigned int v = 0; v < 2; ++v)
5598 vertices_1[v] =
5599 vertex_indices[quad_line_vertices[i][f][v]];
5600
5601 const auto orientation =
5603 vertices_1);
5604
5605#ifdef DEBUG
5606 for (const auto i : vertices_0)
5607 s.insert(i);
5608 for (const auto i : vertices_1)
5609 s.insert(i);
5610#endif
5611
5612 new_quad->set_line_orientation(f, orientation);
5613 }
5614#ifdef DEBUG
5616 s.size(),
5617 (reference_face_type == ReferenceCells::Quadrilateral ? 4 :
5618 3));
5619#endif
5620 }
5621
5622 quad->clear_user_flag();
5623 }
5624 }
5625
5627 cells_with_distorted_children;
5628
5629 for (unsigned int level = 0; level != triangulation.levels.size() - 1;
5630 ++level)
5631 {
5633 hex = triangulation.begin_active_hex(level),
5634 endh = triangulation.begin_active_hex(level + 1);
5636 next_unused_hex = triangulation.begin_raw_hex(level + 1);
5637
5638 for (; hex != endh; ++hex)
5639 {
5640 if (hex->refine_flag_set() ==
5642 continue;
5643
5644 const auto &reference_cell_type = hex->reference_cell();
5645
5646 const RefinementCase<dim> ref_case = hex->refine_flag_set();
5647 hex->clear_refine_flag();
5648 hex->set_refinement_case(ref_case);
5649
5650 unsigned int n_new_lines = 0;
5651 unsigned int n_new_quads = 0;
5652 unsigned int n_new_hexes = 0;
5653
5654 if (reference_cell_type == ReferenceCells::Hexahedron)
5655 {
5656 n_new_lines = 6;
5657 n_new_quads = 12;
5658 n_new_hexes = 8;
5659 }
5660 else if (reference_cell_type == ReferenceCells::Tetrahedron)
5661 {
5662 n_new_lines = 1;
5663 n_new_quads = 8;
5664 n_new_hexes = 8;
5665 }
5666 else
5667 Assert(false, ExcNotImplemented());
5668
5669 // Hexes add a single new internal vertex
5670 if (reference_cell_type == ReferenceCells::Hexahedron)
5671 {
5672 current_vertex =
5673 get_next_unused_vertex(current_vertex,
5674 triangulation.vertices_used);
5675 triangulation.vertices[current_vertex] =
5676 hex->center(true, true);
5677 }
5678
5679 boost::container::small_vector<
5681 6>
5682 new_lines(n_new_lines);
5683 for (unsigned int i = 0; i < n_new_lines; ++i)
5684 {
5685 new_lines[i] =
5686 triangulation.faces->lines
5687 .template next_free_single_object<1>(triangulation);
5688
5689 AssertIsNotUsed(new_lines[i]);
5690 new_lines[i]->set_used_flag();
5691 new_lines[i]->clear_user_flag();
5692 new_lines[i]->clear_user_data();
5693 new_lines[i]->clear_children();
5694 new_lines[i]->set_boundary_id_internal(
5696 new_lines[i]->set_manifold_id(hex->manifold_id());
5697 }
5698
5699 boost::container::small_vector<
5701 12>
5702 new_quads(n_new_quads);
5703 for (unsigned int i = 0; i < n_new_quads; ++i)
5704 {
5705 new_quads[i] =
5706 triangulation.faces->quads
5707 .template next_free_single_object<2>(triangulation);
5708
5709 auto &new_quad = new_quads[i];
5710
5711 // TODO: faces of children have the same type as the faces
5712 // of the parent
5713 triangulation.faces
5714 ->quad_reference_cell[new_quad->index()] =
5715 (reference_cell_type == ReferenceCells::Hexahedron) ?
5718
5719 AssertIsNotUsed(new_quad);
5720 new_quad->set_used_flag();
5721 new_quad->clear_user_flag();
5722 new_quad->clear_user_data();
5723 new_quad->clear_children();
5724 new_quad->set_boundary_id_internal(
5726 new_quad->set_manifold_id(hex->manifold_id());
5727 for (const auto j : new_quads[i]->line_indices())
5728 new_quad->set_line_orientation(j, true);
5729 }
5730
5731 // we always get 8 children per refined cell
5732 std::array<
5734 8>
5735 new_hexes;
5736 {
5737 for (unsigned int i = 0; i < n_new_hexes; ++i)
5738 {
5739 if (i % 2 == 0)
5740 next_unused_hex =
5741 triangulation.levels[level + 1]->cells.next_free_hex(
5742 triangulation, level + 1);
5743 else
5744 ++next_unused_hex;
5745
5746 new_hexes[i] = next_unused_hex;
5747
5748 auto &new_hex = new_hexes[i];
5749
5750 // TODO: children have the same type as the parent
5751 triangulation.levels[new_hex->level()]
5752 ->reference_cell[new_hex->index()] =
5753 reference_cell_type;
5754
5755 AssertIsNotUsed(new_hex);
5756 new_hex->set_used_flag();
5757 new_hex->clear_user_flag();
5758 new_hex->clear_user_data();
5759 new_hex->clear_children();
5760 new_hex->set_material_id(hex->material_id());
5761 new_hex->set_manifold_id(hex->manifold_id());
5762 new_hex->set_subdomain_id(hex->subdomain_id());
5763
5764 if (i % 2)
5765 new_hex->set_parent(hex->index());
5766 // set the face_orientation flag to true for all
5767 // faces initially, as this is the default value
5768 // which is true for all faces interior to the
5769 // hex. later on go the other way round and
5770 // reset faces that are at the boundary of the
5771 // mother cube
5772 //
5773 // the same is true for the face_flip and
5774 // face_rotation flags. however, the latter two
5775 // are set to false by default as this is the
5776 // standard value
5777 for (const auto f : new_hex->face_indices())
5778 {
5779 new_hex->set_face_orientation(f, true);
5780 new_hex->set_face_flip(f, false);
5781 new_hex->set_face_rotation(f, false);
5782 }
5783 }
5784 for (unsigned int i = 0; i < n_new_hexes / 2; ++i)
5785 hex->set_children(2 * i, new_hexes[2 * i]->index());
5786 }
5787
5788 {
5789 // load vertex indices
5790 std::array<unsigned int, 27> vertex_indices = {};
5791
5792 {
5793 unsigned int k = 0;
5794
5795 for (const unsigned int i : hex->vertex_indices())
5796 vertex_indices[k++] = hex->vertex_index(i);
5797
5798 for (const unsigned int i : hex->line_indices())
5799 vertex_indices[k++] =
5800 hex->line(i)->child(0)->vertex_index(1);
5801
5802 if (reference_cell_type == ReferenceCells::Hexahedron)
5803 {
5804 for (const unsigned int i : hex->face_indices())
5805 vertex_indices[k++] =
5806 middle_vertex_index<dim, spacedim>(hex->face(i));
5807
5808 vertex_indices[k++] = current_vertex;
5809 }
5810 }
5811
5812 // set up new lines
5813 {
5814 static constexpr std::array<std::array<unsigned int, 2>, 6>
5815 new_line_vertices_hex = {{{{22, 26}},
5816 {{26, 23}},
5817 {{20, 26}},
5818 {{26, 21}},
5819 {{24, 26}},
5820 {{26, 25}}}};
5821
5822 static constexpr std::array<std::array<unsigned int, 2>, 6>
5823 new_line_vertices_tet = {{{{6, 8}},
5824 {{X, X}},
5825 {{X, X}},
5826 {{X, X}},
5827 {{X, X}},
5828 {{X, X}}}};
5829
5830 const auto &new_line_vertices =
5831 (reference_cell_type == ReferenceCells::Hexahedron) ?
5832 new_line_vertices_hex :
5833 new_line_vertices_tet;
5834
5835 for (unsigned int i = 0; i < new_lines.size(); ++i)
5836 new_lines[i]->set_bounding_object_indices(
5837 {vertex_indices[new_line_vertices[i][0]],
5838 vertex_indices[new_line_vertices[i][1]]});
5839 }
5840
5841 // set up new quads
5842 {
5843 boost::container::small_vector<
5845 30>
5846 relevant_lines(0);
5847
5848 if (reference_cell_type == ReferenceCells::Hexahedron)
5849 {
5850 relevant_lines.resize(30);
5851 for (unsigned int f = 0, k = 0; f < 6; ++f)
5852 for (unsigned int c = 0; c < 4; ++c, ++k)
5853 {
5854 static constexpr std::
5855 array<std::array<unsigned int, 2>, 4>
5856 temp = {
5857 {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}};
5858
5859 relevant_lines[k] =
5860 hex->face(f)
5861 ->isotropic_child(
5863 standard_to_real_face_vertex(
5864 temp[c][0],
5865 hex->face_orientation(f),
5866 hex->face_flip(f),
5867 hex->face_rotation(f)))
5868 ->line(GeometryInfo<dim>::
5869 standard_to_real_face_line(
5870 temp[c][1],
5871 hex->face_orientation(f),
5872 hex->face_flip(f),
5873 hex->face_rotation(f)));
5874 }
5875
5876 for (unsigned int i = 0, k = 24; i < 6; ++i, ++k)
5877 relevant_lines[k] = new_lines[i];
5878 }
5879 else if (reference_cell_type == ReferenceCells::Tetrahedron)
5880 {
5881 relevant_lines.resize(13);
5882
5883 unsigned int k = 0;
5884 for (unsigned int f = 0; f < 4; ++f)
5885 for (unsigned int l = 0; l < 3; ++l, ++k)
5886 {
5887 // TODO: add comment
5888 static const std::
5889 array<std::array<unsigned int, 3>, 6>
5890 table = {{{{1, 0, 2}}, // 0
5891 {{0, 1, 2}},
5892 {{0, 2, 1}}, // 2
5893 {{1, 2, 0}},
5894 {{2, 1, 0}}, // 4
5895 {{2, 0, 1}}}};
5896
5897 relevant_lines[k] =
5898 hex->face(f)
5899 ->child(3 /*center triangle*/)
5900 ->line(
5901 table[triangulation.levels[hex->level()]
5902 ->face_orientations
5903 [hex->index() *
5905 dim>::faces_per_cell +
5906 f]][l]);
5907 }
5908
5909 relevant_lines[k++] = new_lines[0];
5910
5911 AssertDimension(k, 13);
5912 }
5913 else
5914 Assert(false, ExcNotImplemented());
5915
5916 boost::container::small_vector<unsigned int, 30>
5917 relevant_line_indices(relevant_lines.size());
5918 for (unsigned int i = 0; i < relevant_line_indices.size();
5919 ++i)
5920 relevant_line_indices[i] = relevant_lines[i]->index();
5921
5922 static constexpr std::array<std::array<unsigned int, 4>, 12>
5923 new_quad_lines_hex = {{{{10, 28, 16, 24}},
5924 {{28, 14, 17, 25}},
5925 {{11, 29, 24, 20}},
5926 {{29, 15, 25, 21}},
5927 {{18, 26, 0, 28}},
5928 {{26, 22, 1, 29}},
5929 {{19, 27, 28, 4}},
5930 {{27, 23, 29, 5}},
5931 {{2, 24, 8, 26}},
5932 {{24, 6, 9, 27}},
5933 {{3, 25, 26, 12}},
5934 {{25, 7, 27, 13}}}};
5935
5936 static constexpr std::array<std::array<unsigned int, 4>, 12>
5937 new_quad_lines_tet = {{{{2, 3, 8, X}},
5938 {{0, 9, 5, X}},
5939 {{1, 6, 11, X}},
5940 {{4, 10, 7, X}},
5941 {{2, 12, 5, X}},
5942 {{1, 9, 12, X}},
5943 {{4, 8, 12, X}},
5944 {{6, 12, 10, X}},
5945 {{X, X, X, X}},
5946 {{X, X, X, X}},
5947 {{X, X, X, X}},
5948 {{X, X, X, X}}}};
5949
5950 static constexpr std::
5951 array<std::array<std::array<unsigned int, 2>, 4>, 12>
5952 table_hex = {
5953 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
5954 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
5955 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
5956 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
5957 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
5958 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
5959 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
5960 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
5961 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
5962 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
5963 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
5964 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
5965
5966 static constexpr std::
5967 array<std::array<std::array<unsigned int, 2>, 4>, 12>
5968 table_tet = {
5969 {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
5970 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
5971 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
5972 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
5973 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
5974 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
5975 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
5976 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
5977 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5978 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5979 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
5980 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
5981
5982 const auto &new_quad_lines =
5983 (reference_cell_type == ReferenceCells::Hexahedron) ?
5984 new_quad_lines_hex :
5985 new_quad_lines_tet;
5986
5987 const auto &table =
5988 (reference_cell_type == ReferenceCells::Hexahedron) ?
5989 table_hex :
5990 table_tet;
5991
5992 for (unsigned int q = 0; q < new_quads.size(); ++q)
5993 {
5994 for (unsigned int l = 0; l < 3; ++l)
5995 {
5996 std::array<unsigned int, 2> vertices_0, vertices_1;
5997
5998 for (unsigned int v = 0; v < 2; ++v)
5999 vertices_0[v] =
6000 relevant_lines[new_quad_lines[q][l]]
6001 ->vertex_index(v);
6002
6003 for (unsigned int v = 0; v < 2; ++v)
6004 vertices_1[v] = vertex_indices[table[q][l][v]];
6005 }
6006 }
6007
6008 for (unsigned int q = 0; q < new_quads.size(); ++q)
6009 {
6010 auto &new_quad = new_quads[q];
6011
6012 if (new_quad->n_lines() == 3)
6013 new_quad->set_bounding_object_indices(
6014 {relevant_line_indices[new_quad_lines[q][0]],
6015 relevant_line_indices[new_quad_lines[q][1]],
6016 relevant_line_indices[new_quad_lines[q][2]]});
6017 else if (new_quad->n_lines() == 4)
6018 new_quad->set_bounding_object_indices(
6019 {relevant_line_indices[new_quad_lines[q][0]],
6020 relevant_line_indices[new_quad_lines[q][1]],
6021 relevant_line_indices[new_quad_lines[q][2]],
6022 relevant_line_indices[new_quad_lines[q][3]]});
6023 else
6024 Assert(false, ExcNotImplemented());
6025
6026 for (const auto l : new_quad->line_indices())
6027 {
6028 std::array<unsigned int, 2> vertices_0, vertices_1;
6029
6030 for (unsigned int v = 0; v < 2; ++v)
6031 vertices_0[v] =
6032 relevant_lines[new_quad_lines[q][l]]
6033 ->vertex_index(v);
6034
6035 for (unsigned int v = 0; v < 2; ++v)
6036 vertices_1[v] = vertex_indices[table[q][l][v]];
6037
6038 const auto orientation =
6040 vertices_0, vertices_1);
6041
6042 new_quad->set_line_orientation(l, orientation);
6043 }
6044 }
6045 }
6046
6047 // set up new hex
6048 {
6049 std::array<int, 36> quad_indices;
6050
6051 if (reference_cell_type == ReferenceCells::Hexahedron)
6052 {
6053 for (unsigned int i = 0; i < new_quads.size(); ++i)
6054 quad_indices[i] = new_quads[i]->index();
6055
6056 for (unsigned int f = 0, k = new_quads.size(); f < 6;
6057 ++f)
6058 for (unsigned int c = 0; c < 4; ++c, ++k)
6059 quad_indices[k] =
6060 hex->face(f)->isotropic_child_index(
6062 c,
6063 hex->face_orientation(f),
6064 hex->face_flip(f),
6065 hex->face_rotation(f)));
6066 }
6067 else if (reference_cell_type == ReferenceCells::Tetrahedron)
6068 {
6069 for (unsigned int i = 0; i < new_quads.size(); ++i)
6070 quad_indices[i] = new_quads[i]->index();
6071
6072 for (unsigned int f = 0, k = new_quads.size(); f < 4;
6073 ++f)
6074 for (unsigned int c = 0; c < 4; ++c, ++k)
6075 {
6076 quad_indices[k] = hex->face(f)->child_index(
6077 (c == 3) ?
6078 3 :
6079 reference_cell_type
6080 .standard_to_real_face_vertex(
6081 c,
6082 f,
6083 triangulation.levels[hex->level()]
6084 ->face_orientations
6085 [hex->index() *
6087 f]));
6088 }
6089 }
6090 else
6091 {
6092 Assert(false, ExcNotImplemented());
6093 }
6094
6095 static constexpr std::array<std::array<unsigned int, 6>, 8>
6096 cell_quads_hex = {{
6097 {{12, 0, 20, 4, 28, 8}}, // bottom children
6098 {{0, 16, 22, 6, 29, 9}}, //
6099 {{13, 1, 4, 24, 30, 10}}, //
6100 {{1, 17, 6, 26, 31, 11}}, //
6101 {{14, 2, 21, 5, 8, 32}}, // top children
6102 {{2, 18, 23, 7, 9, 33}}, //
6103 {{15, 3, 5, 25, 10, 34}}, //
6104 {{3, 19, 7, 27, 11, 35}} //
6105 }};
6106
6107 static constexpr std::array<std::array<unsigned int, 6>, 8>
6108 cell_quads_tet{{{{8, 13, 16, 0, X, X}},
6109 {{9, 12, 1, 21, X, X}},
6110 {{10, 2, 17, 20, X, X}},
6111 {{3, 14, 18, 22, X, X}},
6112 {{11, 1, 4, 5, X, X}},
6113 {{15, 0, 4, 6, X, X}},
6114 {{19, 7, 6, 3, X, X}},
6115 {{23, 5, 2, 7, X, X}}}};
6116
6117 static constexpr std::
6118 array<std::array<std::array<unsigned int, 4>, 6>, 8>
6119 cell_face_vertices_hex{{{{{{0, 8, 16, 20}},
6120 {{10, 24, 22, 26}},
6121 {{0, 16, 10, 22}},
6122 {{8, 20, 24, 26}},
6123 {{0, 10, 8, 24}},
6124 {{16, 22, 20, 26}}}},
6125 {{{{10, 24, 22, 26}},
6126 {{1, 9, 17, 21}},
6127 {{10, 22, 1, 17}},
6128 {{24, 26, 9, 21}},
6129 {{10, 1, 24, 9}},
6130 {{22, 17, 26, 21}}}},
6131 {{{{8, 2, 20, 18}},
6132 {{24, 11, 26, 23}},
6133 {{8, 20, 24, 26}},
6134 {{2, 18, 11, 23}},
6135 {{8, 24, 2, 11}},
6136 {{20, 26, 18, 23}}}},
6137 {{{{24, 11, 26, 23}},
6138 {{9, 3, 21, 19}},
6139 {{24, 26, 9, 21}},
6140 {{11, 23, 3, 19}},
6141 {{24, 9, 11, 3}},
6142 {{26, 21, 23, 19}}}},
6143 {{{{16, 20, 4, 12}},
6144 {{22, 26, 14, 25}},
6145 {{16, 4, 22, 14}},
6146 {{20, 12, 26, 25}},
6147 {{16, 22, 20, 26}},
6148 {{4, 14, 12, 25}}}},
6149 {{{{22, 26, 14, 25}},
6150 {{17, 21, 5, 13}},
6151 {{22, 14, 17, 5}},
6152 {{26, 25, 21, 13}},
6153 {{22, 17, 26, 21}},
6154 {{14, 5, 25, 13}}}},
6155 {{{{20, 18, 12, 6}},
6156 {{26, 23, 25, 15}},