deal.II version GIT relicensing-2594-g94a7da5c31 2025-02-10 10:20:00+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
tria_accessor.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tria_accessor_h
16#define dealii_tria_accessor_h
17
18
19#include <deal.II/base/config.h>
20
24#include <deal.II/base/point.h>
26
34
35#include <boost/container/small_vector.hpp>
36
37#include <cmath>
38#include <limits>
39#include <utility>
40
41
43
44// Forward declarations
45#ifndef DOXYGEN
46template <int dim, int spacedim>
48class Triangulation;
49template <typename Accessor>
50class TriaRawIterator;
51template <typename Accessor>
52class TriaIterator;
53template <typename Accessor>
55
56namespace parallel
57{
58 template <int dim, int spacedim>
60 class TriangulationBase;
61}
62
63template <int dim, int spacedim>
65class DoFHandler;
66template <int dim, int spacedim, bool lda>
67class DoFCellAccessor;
68
69
70template <int dim, int spacedim>
71class Manifold;
72
73template <int dim, int spacedim>
74class Mapping;
75#endif
76
77namespace internal
78{
79 namespace TriangulationImplementation
80 {
81 class TriaObjects;
82 struct Implementation;
83 struct ImplementationMixedMesh;
84 } // namespace TriangulationImplementation
85
86 namespace TriaAccessorImplementation
87 {
88 struct Implementation;
89
95 template <int structdim, int dim>
97 {
98 struct type
99 {
103 type() = default;
104
108 type(const int level)
109 {
111 (void)level; // removes -Wunused-parameter warning in optimized mode
112 }
113
117 operator int() const
118 {
119 return 0;
120 }
121
122 void
124 {
126 }
127
128 void
130 {
132 }
133 };
134 };
135
136
142 template <int dim>
143 struct PresentLevelType<dim, dim>
144 {
145 using type = int;
146 };
147 } // namespace TriaAccessorImplementation
148} // namespace internal
149template <int structdim, int dim, int spacedim>
150class TriaAccessor;
151template <int dim, int spacedim>
152class TriaAccessor<0, dim, spacedim>;
153template <int spacedim>
154class TriaAccessor<0, 1, spacedim>;
155
160{
165 "The operation you are attempting can only be performed for "
166 "(cell, face, or edge) iterators that point to valid "
167 "objects. These objects need not necessarily be active, "
168 "i.e., have no children, but they need to be part of a "
169 "triangulation. (The objects pointed to by an iterator "
170 "may -- after coarsening -- also be objects that used "
171 "to be part of a triangulation, but are now no longer "
172 "used. Their memory location may have been retained "
173 "for re-use upon the next mesh refinement, but is "
174 "currently unused.)");
185 "The operation you are attempting can only be performed for "
186 "(cell, face, or edge) iterators that point to 'active' "
187 "objects. 'Active' objects are those that do not have "
188 "children (in the case of cells), or that are part of "
189 "an active cell (in the case of faces or edges). However, "
190 "the object on which you are trying the current "
191 "operation is not 'active' in this sense.");
198 "The operation you are attempting can only be performed for "
199 "(cell, face, or edge) iterators that have children, "
200 "but the object on which you are trying the current "
201 "operation does not have any.");
209 "The operation you are attempting can only be performed for "
210 "(cell, face, or edge) iterators that have a parent object, "
211 "but the object on which you are trying the current "
212 "operation does not have one -- i.e., it is on the "
213 "coarsest level of the triangulation.");
218 int,
219 << "You can only set the child index if the cell does not "
220 << "currently have children registered; or you can clear it. "
221 << "The given index was " << arg1
222 << " (-1 means: clear children).");
226 template <typename AccessorType>
228 AccessorType,
229 << "You tried to dereference an iterator for which this "
230 << "is not possible. More information on this iterator: "
231 << "index=" << arg1.index() << ", state="
232 << (arg1.state() == IteratorState::valid ?
233 "valid" :
234 (arg1.state() == IteratorState::past_the_end ?
235 "past_the_end" :
236 "invalid")));
241 "Iterators can only be compared if they point to the same "
242 "triangulation, or if neither of them are associated "
243 "with a triangulation.");
244 // TODO: Write documentation!
249 // TODO: Write documentation!
269 // TODO: Write documentation!
275 int,
276 << "You can only set the child index of an even numbered child."
277 << "The number of the child given was " << arg1 << '.');
278} // namespace TriaAccessorExceptions
279
280
306template <int structdim, int dim, int spacedim = dim>
308{
309public:
316 static constexpr unsigned int space_dimension = spacedim;
317
323 static constexpr unsigned int dimension = dim;
324
330 static const unsigned int structure_dimension = structdim;
331
341 void
342 operator=(const TriaAccessorBase *) = delete;
343
344protected:
350 using AccessorData = void;
351
356 const int level = -1,
357 const int index = -1,
358 const AccessorData * = nullptr);
359
364
372 void
374
380
391 bool
392 operator<(const TriaAccessorBase &other) const;
393
394protected:
398 bool
400
404 bool
406
420 void
422
430 void
440 objects() const;
441
442public:
448 using LocalData = void *;
449
473 int
474 level() const;
475
502 int
503 index() const;
504
510 state() const;
511
518
522protected:
527 typename ::internal::TriaAccessorImplementation::
528 PresentLevelType<structdim, dim>::type present_level;
529
535
540
541private:
542 template <typename Accessor>
543 friend class TriaRawIterator;
544 template <typename Accessor>
545 friend class TriaIterator;
546 template <typename Accessor>
547 friend class TriaActiveIterator;
548};
549
550
551
572template <int structdim, int dim, int spacedim = dim>
574{
575public:
582 static constexpr unsigned int space_dimension = spacedim;
583
589 static constexpr unsigned int dimension = dim;
590
596 static const unsigned int structure_dimension = structdim;
597
603 using AccessorData = void;
604
612 InvalidAccessor(const void *parent = nullptr,
613 const int level = -1,
614 const int index = -1,
615 const AccessorData *local_data = nullptr);
616
625
630 template <typename OtherAccessor>
631 InvalidAccessor(const OtherAccessor &);
632
636 void
638
642 bool
644 bool
646
650 void
651 operator++() const;
652 void
653 operator--() const;
654
660 state();
661
662
667 static int
668 level();
669
674 static int
675 index();
676
681 bool
682 used() const;
683
688 bool
690
695 manifold_id() const;
696
700 unsigned int
701 user_index() const;
702
706 void
707 set_user_index(const unsigned int p) const;
708
712 void
714
719 vertex(const unsigned int i) const;
720
725 void *
726 line(const unsigned int i) const;
727
732 void *
733 quad(const unsigned int i) const;
734};
735
736
737
755template <int structdim, int dim, int spacedim>
756class TriaAccessor : public TriaAccessorBase<structdim, dim, spacedim>
757{
758public:
764
769 const int level = -1,
770 const int index = -1,
771 const AccessorData *local_data = nullptr);
772
777 TriaAccessor(const TriaAccessor &) = default;
778
782 TriaAccessor(TriaAccessor &&) = default; // NOLINT
783
796 template <int structdim2, int dim2, int spacedim2>
798
803 template <int structdim2, int dim2, int spacedim2>
805
816 operator=(const TriaAccessor &) = delete;
817
822 operator=(TriaAccessor &&) = default; // NOLINT
823
827 ~TriaAccessor() = default;
828
835 bool
836 used() const;
837
850 vertex_iterator(const unsigned int i) const;
851
867 unsigned int
868 vertex_index(const unsigned int i) const;
869
908 vertex(const unsigned int i) const;
909
913 typename ::internal::TriangulationImplementation::
914 Iterators<dim, spacedim>::line_iterator
915 line(const unsigned int i) const;
916
923 unsigned int
924 line_index(const unsigned int i) const;
925
929 typename ::internal::TriangulationImplementation::
930 Iterators<dim, spacedim>::quad_iterator
931 quad(const unsigned int i) const;
932
939 unsigned int
940 quad_index(const unsigned int i) const;
959 combined_face_orientation(const unsigned int face) const;
960
972 bool
973 face_orientation(const unsigned int face) const;
974
984 bool
985 face_flip(const unsigned int face) const;
986
996 bool
997 face_rotation(const unsigned int face) const;
998
1009 line_orientation(const unsigned int line) const;
1024 bool
1026
1031 unsigned int
1032 n_children() const;
1033
1047 unsigned int
1049
1063 unsigned int
1065
1070 child(const unsigned int i) const;
1071
1076 unsigned int
1079
1089 isotropic_child(const unsigned int i) const;
1090
1096
1102 int
1103 child_index(const unsigned int i) const;
1104
1110 int
1111 isotropic_child_index(const unsigned int i) const;
1135
1163 void
1165
1194 void
1196
1204 bool
1206
1218
1241
1259 void
1261
1275 void
1277
1294 bool
1296
1302 void
1304
1310 void
1312
1318 void
1320
1326 void
1328
1334 void
1336
1348 void
1349 set_user_pointer(void *p) const;
1350
1356 void
1358
1376 void *
1378
1400 void
1402
1409 void
1411
1421 void
1422 set_user_index(const unsigned int p) const;
1423
1429 void
1431
1443 unsigned int
1444 user_index() const;
1445
1463 void
1464 recursively_set_user_index(const unsigned int p) const;
1465
1474 void
1510 double
1511 diameter() const;
1512
1539 std::pair<Point<spacedim>, double>
1541
1552
1562 double
1563 extent_in_direction(const unsigned int axis) const;
1564
1568 double
1570
1585 intermediate_point(const Point<structdim> &coordinates) const;
1586
1611
1647 center(const bool respect_manifold = false,
1648 const bool interpolate_from_surrounding = false) const;
1649
1668 barycenter() const;
1669
1700 double
1701 measure() const;
1702
1717 bool
1720
1726
1730 unsigned int
1731 n_vertices() const;
1732
1736 unsigned int
1737 n_lines() const;
1738
1748 unsigned int
1749 n_faces() const;
1750
1757
1764
1773
1778private:
1783 void
1785
1793 void
1795 const std::initializer_list<int> &new_indices) const;
1796
1800 void
1802 const std::initializer_list<unsigned int> &new_indices) const;
1803
1811 void
1812 set_line_orientation(const unsigned int line,
1813 const types::geometric_orientation orientation) const;
1814
1822 void
1824 const unsigned int face,
1825 const types::geometric_orientation combined_orientation) const;
1826
1830 void
1832
1836 void
1838
1847 void
1849
1857 void
1859
1866 void
1867 set_children(const unsigned int i, const int index) const;
1868
1873 void
1875
1876private:
1877 friend class Triangulation<dim, spacedim>;
1878
1879 friend struct ::internal::TriangulationImplementation::Implementation;
1880 friend struct ::internal::TriangulationImplementation::
1881 ImplementationMixedMesh;
1882 friend struct ::internal::TriaAccessorImplementation::Implementation;
1883};
1884
1885
1886
1905template <int dim, int spacedim>
1906class TriaAccessor<0, dim, spacedim>
1907{
1908public:
1914 static constexpr unsigned int space_dimension = spacedim;
1915
1921 static constexpr unsigned int dimension = dim;
1922
1928 static const unsigned int structure_dimension = 0;
1929
1933 using AccessorData = void;
1934
1940 const unsigned int vertex_index);
1941
1948 const int level = 0,
1949 const int index = 0,
1950 const AccessorData * = nullptr);
1951
1955 template <int structdim2, int dim2, int spacedim2>
1957
1961 template <int structdim2, int dim2, int spacedim2>
1963
1968 state() const;
1969
1974 static int
1976
1981 int
1982 index() const;
1983
1990
2000 void
2002
2006 void
2011 bool
2012 operator==(const TriaAccessor &) const;
2013
2017 bool
2018 operator!=(const TriaAccessor &) const;
2019
2047 unsigned int
2048 vertex_index(const unsigned int i = 0) const;
2049
2056 vertex(const unsigned int i = 0) const;
2057
2062 typename ::internal::TriangulationImplementation::
2063 Iterators<dim, spacedim>::line_iterator static line(const unsigned int);
2064
2068 static unsigned int
2069 line_index(const unsigned int i);
2070
2074 static typename ::internal::TriangulationImplementation::
2075 Iterators<dim, spacedim>::quad_iterator
2076 quad(const unsigned int i);
2077
2081 static unsigned int
2082 quad_index(const unsigned int i);
2083
2099 double
2100 diameter() const;
2101
2109 double
2110 extent_in_direction(const unsigned int axis) const;
2111
2120 center(const bool respect_manifold = false,
2121 const bool interpolate_from_surrounding = false) const;
2122
2131 double
2132 measure() const;
2150 combined_face_orientation(const unsigned int face);
2151
2155 static bool
2156 face_orientation(const unsigned int face);
2157
2161 static bool
2162 face_flip(const unsigned int face);
2163
2167 static bool
2168 face_rotation(const unsigned int face);
2169
2174 line_orientation(const unsigned int line);
2175
2190 static bool
2192
2197 static unsigned int
2199
2204 static unsigned int
2206
2210 static unsigned int
2212
2216 static unsigned int
2218
2223 child(const unsigned int);
2224
2229 isotropic_child(const unsigned int);
2230
2234 static RefinementCase<0>
2236
2240 static int
2241 child_index(const unsigned int i);
2242
2246 static int
2247 isotropic_child_index(const unsigned int i);
2255 bool
2256 used() const;
2257
2258protected:
2266 void
2268
2277 bool
2278 operator<(const TriaAccessor &other) const;
2279
2284
2289
2290private:
2291 template <typename Accessor>
2292 friend class TriaRawIterator;
2293 template <typename Accessor>
2294 friend class TriaIterator;
2295 template <typename Accessor>
2297};
2298
2299
2300
2317template <int spacedim>
2318class TriaAccessor<0, 1, spacedim>
2319{
2320public:
2326 static constexpr unsigned int space_dimension = spacedim;
2327
2333 static constexpr unsigned int dimension = 1;
2334
2340 static const unsigned int structure_dimension = 0;
2341
2345 using AccessorData = void;
2346
2352 {
2364 right_vertex
2366
2379 const VertexKind vertex_kind,
2380 const unsigned int vertex_index);
2381
2388 const int = 0,
2389 const int = 0,
2390 const AccessorData * = nullptr);
2391
2395 template <int structdim2, int dim2, int spacedim2>
2397
2401 template <int structdim2, int dim2, int spacedim2>
2403
2408 void
2410
2416 void
2418
2426
2431 static int
2433
2438 int
2439 index() const;
2440
2447
2458 void
2459 operator++() const;
2460
2465 void
2466 operator--() const;
2470 bool
2471 operator==(const TriaAccessor &) const;
2472
2476 bool
2477 operator!=(const TriaAccessor &) const;
2478
2487 bool
2488 operator<(const TriaAccessor &other) const;
2489
2516 unsigned int
2517 vertex_index(const unsigned int i = 0) const;
2518
2525 vertex(const unsigned int i = 0) const;
2526
2532 center() const;
2533
2538 typename ::internal::TriangulationImplementation::
2539 Iterators<1, spacedim>::line_iterator static line(const unsigned int);
2540
2547 static unsigned int
2548 line_index(const unsigned int i);
2549
2553 static typename ::internal::TriangulationImplementation::
2554 Iterators<1, spacedim>::quad_iterator
2555 quad(const unsigned int i);
2556
2563 static unsigned int
2564 quad_index(const unsigned int i);
2565
2575 bool
2577
2594
2598 const Manifold<1, spacedim> &
2600
2609
2610
2622 bool
2623 user_flag_set() const;
2624
2630 void
2631 set_user_flag() const;
2632
2638 void
2639 clear_user_flag() const;
2640
2646 void
2648
2654 void
2656
2662 void
2663 clear_user_data() const;
2664
2676 void
2677 set_user_pointer(void *p) const;
2678
2684 void
2685 clear_user_pointer() const;
2686
2702 void *
2703 user_pointer() const;
2704
2726 void
2727 recursively_set_user_pointer(void *p) const;
2728
2735 void
2737
2747 void
2748 set_user_index(const unsigned int p) const;
2749
2755 void
2756 clear_user_index() const;
2757
2769 unsigned int
2770 user_index() const;
2771
2789 void
2790 recursively_set_user_index(const unsigned int p) const;
2791
2800 void
2819 combined_face_orientation(const unsigned int face);
2820
2824 static bool
2825 face_orientation(const unsigned int face);
2826
2830 static bool
2831 face_flip(const unsigned int face);
2832
2836 static bool
2837 face_rotation(const unsigned int face);
2838
2843 line_orientation(const unsigned int line);
2844
2859 static bool
2861
2866 static unsigned int
2868
2873 static unsigned int
2875
2879 static unsigned int
2881
2885 static unsigned int
2887
2892 child(const unsigned int);
2893
2898 isotropic_child(const unsigned int);
2899
2903 static RefinementCase<0>
2905
2909 static int
2910 child_index(const unsigned int i);
2911
2915 static int
2916 isotropic_child_index(const unsigned int i);
2947 void
2949
2956 void
2958
2968 void
2970
2982 void
2991 bool
2992 used() const;
2993
2999
3003 unsigned int
3004 n_vertices() const;
3005
3009 unsigned int
3010 n_lines() const;
3011
3018
3025
3026protected:
3031
3037
3042};
3043
3044
3045
3061template <int dim, int spacedim = dim>
3062class CellAccessor : public TriaAccessor<dim, dim, spacedim>
3063{
3064public:
3069
3074
3086 const int level = -1,
3087 const int index = -1,
3088 const AccessorData *local_data = nullptr);
3089
3094
3107 template <int structdim2, int dim2, int spacedim2>
3109
3114 template <int structdim2, int dim2, int spacedim2>
3116
3121
3125 // NOLINTNEXTLINE OSX does not compile with noexcept
3127
3131 ~CellAccessor() = default;
3132
3144
3148 // NOLINTNEXTLINE OSX does not compile with noexcept
3151
3175 as_dof_handler_iterator(const DoFHandler<dim, spacedim> &dof_handler) const;
3176
3187 const DoFHandler<dim, spacedim> &dof_handler) const;
3188
3189
3206 child(const unsigned int i) const;
3207
3211 boost::container::small_vector<TriaIterator<CellAccessor<dim, spacedim>>,
3214
3218 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
3219 face(const unsigned int i) const;
3220
3225 unsigned int
3228
3232 boost::container::small_vector<
3233 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
3234#ifndef _MSC_VER // MSVC prior to 2022 cannot use a constexpr function this way
3235 ReferenceCells::max_n_faces<dim>()
3236#else
3238#endif
3239 >
3241
3251 unsigned int
3252 face_index(const unsigned int i) const;
3253
3302 neighbor_child_on_subface(const unsigned int face_no,
3303 const unsigned int subface_no) const;
3304
3353 neighbor(const unsigned int face_no) const;
3354
3362 int
3363 neighbor_index(const unsigned int face_no) const;
3364
3372 int
3373 neighbor_level(const unsigned int face_no) const;
3374
3386 unsigned int
3387 neighbor_of_neighbor(const unsigned int face_no) const;
3388
3399 bool
3400 neighbor_is_coarser(const unsigned int face_no) const;
3401
3416 std::pair<unsigned int, unsigned int>
3417 neighbor_of_coarser_neighbor(const unsigned int neighbor) const;
3418
3425 unsigned int
3426 neighbor_face_no(const unsigned int neighbor) const;
3427
3431 static bool
3433
3447 bool
3448 has_periodic_neighbor(const unsigned int i) const;
3449
3467 periodic_neighbor(const unsigned int i) const;
3468
3477 neighbor_or_periodic_neighbor(const unsigned int i) const;
3478
3494 periodic_neighbor_child_on_subface(const unsigned int face_no,
3495 const unsigned int subface_no) const;
3496
3507 std::pair<unsigned int, unsigned int>
3508 periodic_neighbor_of_coarser_periodic_neighbor(const unsigned face_no) const;
3509
3515 int
3516 periodic_neighbor_index(const unsigned int i) const;
3517
3523 int
3524 periodic_neighbor_level(const unsigned int i) const;
3525
3540 unsigned int
3541 periodic_neighbor_of_periodic_neighbor(const unsigned int i) const;
3542
3548 unsigned int
3549 periodic_neighbor_face_no(const unsigned int i) const;
3550
3557 bool
3558 periodic_neighbor_is_coarser(const unsigned int i) const;
3559
3576 bool
3577 at_boundary(const unsigned int i) const;
3578
3587 bool
3588 at_boundary() const;
3589
3597 bool
3598 has_boundary_lines() const;
3626
3644 void
3647
3651 void
3653
3658 std::uint8_t
3660
3665 void
3666 set_refine_choice(const std::uint8_t refinement_choice = static_cast<char>(
3668
3672 void
3674
3682 bool
3684 const unsigned int face_no,
3685 const RefinementCase<dim - 1> &face_refinement_case =
3687
3693 bool
3694 flag_for_line_refinement(const unsigned int line_no) const;
3695
3705 subface_case(const unsigned int face_no) const;
3706
3710 bool
3712
3717 void
3719
3723 void
3748 material_id() const;
3749
3761 void
3762 set_material_id(const types::material_id new_material_id) const;
3763
3772 void
3773 recursively_set_material_id(const types::material_id new_material_id) const;
3801
3817 void
3818 set_subdomain_id(const types::subdomain_id new_subdomain_id) const;
3819
3830
3835 void
3837 const types::subdomain_id new_level_subdomain_id) const;
3838
3839
3855 void
3857 const types::subdomain_id new_subdomain_id) const;
3881
3891
3908 bool
3909 direction_flag() const;
3910
3936 unsigned int
3938
3946 int
3947 parent_index() const;
3948
3955 parent() const;
3956
3976 bool
3977 is_active() const;
3978
3998 bool
4000
4005 bool
4007
4041 bool
4042 is_ghost() const;
4043
4049 bool
4051
4078 bool
4080
4087 bool
4089
4103 bool
4105
4114 void
4115 set_neighbor(const unsigned int i,
4116 const TriaIterator<CellAccessor<dim, spacedim>> &pointer) const;
4117
4131 CellId
4132 id() const;
4133
4134 using TriaAccessor<dim, dim, spacedim>::diameter;
4135
4139 double
4140 diameter(const Mapping<dim, spacedim> &mapping) const;
4141
4159
4160protected:
4176 unsigned int
4177 neighbor_of_neighbor_internal(const unsigned int neighbor) const;
4178
4184 template <int dim_, int spacedim_>
4185 bool
4186 point_inside_codim(const Point<spacedim_> &p) const;
4187
4188
4189
4190private:
4195 void
4196 set_active_cell_index(const unsigned int active_cell_index) const;
4197
4201 void
4203
4207 void
4209
4213 void
4214 set_parent(const unsigned int parent_index);
4215
4226 void
4227 set_direction_flag(const bool new_direction_flag) const;
4228
4229 friend class Triangulation<dim, spacedim>;
4230
4231 friend class parallel::TriangulationBase<dim, spacedim>;
4232
4233 friend struct ::internal::TriangulationImplementation::Implementation;
4234 friend struct ::internal::TriangulationImplementation::
4235 ImplementationMixedMesh;
4236};
4237
4238
4239
4240/* ----- declaration of explicit specializations and general templates ----- */
4241
4242
4243template <int structdim, int dim, int spacedim>
4244template <typename OtherAccessor>
4246 const OtherAccessor &)
4247{
4248 Assert(false,
4249 ExcMessage("You are attempting an illegal conversion between "
4250 "iterator/accessor types. The constructor you call "
4251 "only exists to make certain template constructs "
4252 "easier to write as dimension independent code but "
4253 "the conversion is not valid in the current context."));
4254}
4255
4256
4257
4258template <int structdim, int dim, int spacedim>
4259template <int structdim2, int dim2, int spacedim2>
4262{
4263 Assert(false,
4264 ExcMessage("You are attempting an illegal conversion between "
4265 "iterator/accessor types. The constructor you call "
4266 "only exists to make certain template constructs "
4267 "easier to write as dimension independent code but "
4268 "the conversion is not valid in the current context."));
4269}
4270
4271
4272
4273template <int dim, int spacedim>
4274template <int structdim2, int dim2, int spacedim2>
4277{
4278 Assert(false,
4279 ExcMessage("You are attempting an illegal conversion between "
4280 "iterator/accessor types. The constructor you call "
4281 "only exists to make certain template constructs "
4282 "easier to write as dimension independent code but "
4283 "the conversion is not valid in the current context."));
4284}
4285
4286
4287
4288template <int structdim, int dim, int spacedim>
4289template <int structdim2, int dim2, int spacedim2>
4292{
4293 Assert(false,
4294 ExcMessage("You are attempting an illegal conversion between "
4295 "iterator/accessor types. The constructor you call "
4296 "only exists to make certain template constructs "
4297 "easier to write as dimension independent code but "
4298 "the conversion is not valid in the current context."));
4299}
4300
4301
4302
4303template <int dim, int spacedim>
4304template <int structdim2, int dim2, int spacedim2>
4307{
4308 Assert(false,
4309 ExcMessage("You are attempting an illegal conversion between "
4310 "iterator/accessor types. The constructor you call "
4311 "only exists to make certain template constructs "
4312 "easier to write as dimension independent code but "
4313 "the conversion is not valid in the current context."));
4314}
4315
4316
4317#ifndef DOXYGEN
4318
4319template <>
4320bool
4322template <>
4323bool
4325template <>
4326bool
4328template <>
4329bool
4331template <>
4332bool
4334template <>
4335bool
4337// -------------------------------------------------------------------
4338
4339template <>
4340void
4342
4343
4344
4345namespace internal
4346{
4347 namespace TriaAccessorImplementation
4348 {
4354 template <int dim, int spacedim>
4355 inline double
4356 diameter(
4357 const boost::container::small_vector<Point<spacedim>,
4359 vertices)
4360 {
4361 const ReferenceCell reference_cell =
4362 ReferenceCell::n_vertices_to_type(dim, vertices.size());
4363
4364 if (reference_cell == ReferenceCells::Line)
4365 // Return the distance between the two vertices
4366 return (vertices[1] - vertices[0]).norm();
4367 else if (reference_cell == ReferenceCells::Triangle)
4368 // Return the longest of the three edges
4369 return std::max({(vertices[1] - vertices[0]).norm(),
4370 (vertices[2] - vertices[1]).norm(),
4371 (vertices[2] - vertices[0]).norm()});
4372 else if (reference_cell == ReferenceCells::Quadrilateral)
4373 // Return the longer one of the two diagonals of the quadrilateral
4374 return std::max({(vertices[3] - vertices[0]).norm(),
4375 (vertices[2] - vertices[1]).norm()});
4376 else if (reference_cell == ReferenceCells::Tetrahedron)
4377 // Return the longest of the six edges of the tetrahedron
4378 return std::max({(vertices[1] - vertices[0]).norm(),
4379 (vertices[2] - vertices[0]).norm(),
4380 (vertices[2] - vertices[1]).norm(),
4381 (vertices[3] - vertices[0]).norm(),
4382 (vertices[3] - vertices[1]).norm(),
4383 (vertices[3] - vertices[2]).norm()});
4384 else if (reference_cell == ReferenceCells::Pyramid)
4385 // Return ...
4386 return std::max({// the longest diagonal of the quadrilateral base
4387 // of the pyramid or ...
4388 (vertices[3] - vertices[0]).norm(),
4389 (vertices[2] - vertices[1]).norm(),
4390 // the longest edge connected with the apex of the
4391 // pyramid
4392 (vertices[4] - vertices[0]).norm(),
4393 (vertices[4] - vertices[1]).norm(),
4394 (vertices[4] - vertices[2]).norm(),
4395 (vertices[4] - vertices[3]).norm()});
4396 else if (reference_cell == ReferenceCells::Wedge)
4397 // Return ...
4398 return std::max({// the longest of the 2*3=6 diagonals of the three
4399 // quadrilateral sides of the wedge or ...
4400 (vertices[4] - vertices[0]).norm(),
4401 (vertices[3] - vertices[1]).norm(),
4402 (vertices[5] - vertices[1]).norm(),
4403 (vertices[4] - vertices[2]).norm(),
4404 (vertices[5] - vertices[0]).norm(),
4405 (vertices[3] - vertices[2]).norm(),
4406 // the longest of the 3*2=6 edges of the two
4407 // triangular faces of the wedge
4408 (vertices[1] - vertices[0]).norm(),
4409 (vertices[2] - vertices[1]).norm(),
4410 (vertices[2] - vertices[0]).norm(),
4411 (vertices[4] - vertices[3]).norm(),
4412 (vertices[5] - vertices[4]).norm(),
4413 (vertices[5] - vertices[3]).norm()});
4414 else if (reference_cell == ReferenceCells::Hexahedron)
4415 // Return the longest of the four diagonals of the hexahedron
4416 return std::max({(vertices[7] - vertices[0]).norm(),
4417 (vertices[6] - vertices[1]).norm(),
4418 (vertices[2] - vertices[5]).norm(),
4419 (vertices[3] - vertices[4]).norm()});
4420
4422 return -1e10;
4423 }
4424 } // namespace TriaAccessorImplementation
4425} // namespace internal
4426
4427
4428/*--------------------- Functions: TriaAccessorBase -------------------------*/
4429
4430template <int structdim, int dim, int spacedim>
4432 const Triangulation<dim, spacedim> *tria,
4433 const int level,
4434 const int index,
4435 const AccessorData *)
4436 : present_level((structdim == dim) ? level : 0)
4437 , present_index(index)
4438 , tria(tria)
4439{
4440 // non-cells have no level, so a 0
4441 // should have been passed, or a -1
4442 // for an end-iterator, or -2 for
4443 // an invalid (default constructed)
4444 // iterator
4445 if (structdim != dim)
4446 {
4447 Assert((level == 0) || (level == -1) || (level == -2),
4449 }
4450}
4451
4452
4453template <int structdim, int dim, int spacedim>
4456 : present_level(a.present_level)
4457 , present_index(a.present_index)
4458 , tria(a.tria)
4459{}
4460
4461
4462template <int structdim, int dim, int spacedim>
4463inline void
4466{
4467 present_level = a.present_level;
4468 present_index = a.present_index;
4469 tria = a.tria;
4470
4471 if (structdim != dim)
4472 {
4473 Assert((present_level == 0) || (present_level == -1) ||
4474 (present_level == -2),
4476 }
4477}
4478
4479
4480
4481template <int structdim, int dim, int spacedim>
4485{
4486 present_level = a.present_level;
4487 present_index = a.present_index;
4488 tria = a.tria;
4489
4490 if (structdim != dim)
4491 {
4492 Assert((present_level == 0) || (present_level == -1) ||
4493 (present_level == -2),
4495 }
4496 return *this;
4497}
4498
4499
4500
4501template <int structdim, int dim, int spacedim>
4502inline bool
4505{
4506 Assert(tria == a.tria || tria == nullptr || a.tria == nullptr,
4508 return ((tria == a.tria) && (present_level == a.present_level) &&
4509 (present_index == a.present_index));
4510}
4511
4512
4513
4514template <int structdim, int dim, int spacedim>
4515inline bool
4518{
4519 Assert(tria == a.tria || tria == nullptr || a.tria == nullptr,
4521 return ((tria != a.tria) || (present_level != a.present_level) ||
4522 (present_index != a.present_index));
4523}
4524
4525
4526
4527template <int structdim, int dim, int spacedim>
4528inline bool
4531{
4533
4534 if (present_level != other.present_level)
4535 return (present_level < other.present_level);
4536
4537 return (present_index < other.present_index);
4538}
4539
4540
4541
4542template <int structdim, int dim, int spacedim>
4543inline int
4545{
4546 // This is always zero or invalid
4547 // if the object is not a cell
4548 return present_level;
4549}
4550
4551
4552
4553template <int structdim, int dim, int spacedim>
4554inline int
4556{
4557 return present_index;
4558}
4559
4560
4561
4562template <int structdim, int dim, int spacedim>
4565{
4566 if ((present_level >= 0) && (present_index >= 0))
4567 return IteratorState::valid;
4568 else if (present_index == -1)
4570 else
4572}
4573
4574
4575
4576template <int structdim, int dim, int spacedim>
4577inline const Triangulation<dim, spacedim> &
4579{
4580 return *tria;
4581}
4582
4583
4584
4585template <int structdim, int dim, int spacedim>
4586inline void
4588{
4589 // this iterator is used for
4590 // objects without level
4591 ++this->present_index;
4592
4593 if (structdim != dim)
4594 {
4595 // is index still in the range of
4596 // the vector? (note that we don't
4597 // have to set the level, since
4598 // dim!=1 and the object therefore
4599 // has no level)
4600 if (this->present_index >= static_cast<int>(objects().n_objects()))
4601 this->present_index = -1;
4602 }
4603 else
4604 {
4605 while (this->present_index >=
4606 static_cast<int>(
4607 this->tria->levels[this->present_level]->cells.n_objects()))
4608 {
4609 // no -> go one level up until we find
4610 // one with more than zero cells
4611 ++this->present_level;
4612 this->present_index = 0;
4613 // highest level reached?
4614 if (this->present_level >=
4615 static_cast<int>(this->tria->levels.size()))
4616 {
4617 // return with past the end pointer
4618 this->present_level = this->present_index = -1;
4619 return;
4620 }
4621 }
4622 }
4623}
4624
4625
4626template <int structdim, int dim, int spacedim>
4627inline void
4629{
4630 // same as operator++
4631 --this->present_index;
4632
4633 if (structdim != dim)
4634 {
4635 if (this->present_index < 0)
4636 this->present_index = -1;
4637 }
4638 else
4639 {
4640 while (this->present_index < 0)
4641 {
4642 // no -> go one level down
4643 --this->present_level;
4644 // lowest level reached?
4645 if (this->present_level == -1)
4646 {
4647 // return with past the end pointer
4648 this->present_level = this->present_index = -1;
4649 return;
4650 }
4651 // else
4652 this->present_index =
4653 this->tria->levels[this->present_level]->cells.n_objects() - 1;
4654 }
4655 }
4656}
4657
4658
4659
4660template <int structdim, int dim, int spacedim>
4661inline ::internal::TriangulationImplementation::TriaObjects &
4663{
4664 if (structdim == dim)
4665 return this->tria->levels[this->present_level]->cells;
4666
4667 if (structdim == 1 && dim > 1)
4668 return this->tria->faces->lines;
4669
4670 if (structdim == 2 && dim > 2)
4671 return this->tria->faces->quads;
4672
4674
4675 return this->tria->levels[this->present_level]->cells;
4676}
4677
4678
4679
4680/*---------------------- Functions: InvalidAccessor -------------------------*/
4681
4682template <int structdim, int dim, int spacedim>
4684 const int,
4685 const int,
4686 const AccessorData *)
4687{
4688 Assert(false,
4689 ExcMessage("You are attempting an invalid conversion between "
4690 "iterator/accessor types. The constructor you call "
4691 "only exists to make certain template constructs "
4692 "easier to write as dimension independent code but "
4693 "the conversion is not valid in the current context."));
4694}
4695
4696
4697
4698template <int structdim, int dim, int spacedim>
4700 const InvalidAccessor &)
4701{
4702 Assert(false,
4703 ExcMessage("You are attempting an invalid conversion between "
4704 "iterator/accessor types. The constructor you call "
4705 "only exists to make certain template constructs "
4706 "easier to write as dimension independent code but "
4707 "the conversion is not valid in the current context."));
4708}
4709
4710
4711
4712template <int structdim, int dim, int spacedim>
4713void
4715{
4716 // nothing to do here. we could
4717 // throw an exception but we can't
4718 // get here without first creating
4719 // an object which would have
4720 // already thrown
4721}
4722
4723
4724
4725template <int structdim, int dim, int spacedim>
4726bool
4728 const InvalidAccessor &) const
4729{
4730 // nothing to do here. we could
4731 // throw an exception but we can't
4732 // get here without first creating
4733 // an object which would have
4734 // already thrown
4735 return false;
4736}
4737
4738
4739
4740template <int structdim, int dim, int spacedim>
4741bool
4743 const InvalidAccessor &) const
4744{
4745 // nothing to do here. we could
4746 // throw an exception but we can't
4747 // get here without first creating
4748 // an object which would have
4749 // already thrown
4750 return true;
4751}
4752
4753
4754
4755template <int structdim, int dim, int spacedim>
4756bool
4758{
4759 // nothing to do here. we could
4760 // throw an exception but we can't
4761 // get here without first creating
4762 // an object which would have
4763 // already thrown
4764 return false;
4765}
4766
4767
4768
4769template <int structdim, int dim, int spacedim>
4770bool
4772{
4773 // nothing to do here. we could
4774 // throw an exception but we can't
4775 // get here without first creating
4776 // an object which would have
4777 // already thrown
4778 return false;
4779}
4780
4781
4782
4783template <int structdim, int dim, int spacedim>
4784void
4786{}
4787
4788
4789
4790template <int structdim, int dim, int spacedim>
4791void
4793{}
4794
4795
4796
4797template <int structdim, int dim, int spacedim>
4800{
4802}
4803
4804
4805
4806template <int structdim, int dim, int spacedim>
4807unsigned int
4809{
4811}
4812
4813
4814
4815template <int structdim, int dim, int spacedim>
4816void
4818 const unsigned int) const
4819{
4820 Assert(false,
4821 ExcMessage("You are trying to set the user index of an "
4822 "invalid object."));
4823}
4824
4825
4826
4827template <int structdim, int dim, int spacedim>
4828void
4830 const types::manifold_id) const
4831{
4832 Assert(false,
4833 ExcMessage("You are trying to set the manifold id of an "
4834 "invalid object."));
4835}
4836
4837
4838
4839template <int structdim, int dim, int spacedim>
4840inline Point<spacedim> &
4842{
4843 // nothing to do here. we could throw an exception but we can't get here
4844 // without first creating an object which would have already thrown
4845 static Point<spacedim> invalid_vertex;
4846 return invalid_vertex;
4847}
4848
4849
4850template <int structdim, int dim, int spacedim>
4851inline void *
4853{
4854 // nothing to do here. we could throw an exception but we can't get here
4855 // without first creating an object which would have already thrown
4856 return nullptr;
4857}
4858
4859
4860
4861template <int structdim, int dim, int spacedim>
4862inline void *
4864{
4865 // nothing to do here. we could throw an exception but we can't get here
4866 // without first creating an object which would have already thrown
4867 return nullptr;
4868}
4869
4870
4871/*------------------------ Functions: TriaAccessor ---------------------------*/
4872
4873
4874namespace internal
4875{
4876 namespace TriaAccessorImplementation
4877 {
4878 // make sure that if in the following we
4879 // write TriaAccessor
4880 // we mean the *class*
4881 // ::TriaAccessor, not the
4882 // enclosing namespace
4883 // ::internal::TriaAccessor
4884 using ::TriaAccessor;
4885
4890 struct Implementation
4891 {
4895 template <int structdim, int dim, int spacedim>
4896 inline static void
4897 set_combined_face_orientation(
4899 const unsigned int face,
4900 const types::geometric_orientation combined_orientation)
4901 {
4902 Assert(structdim == dim,
4903 ExcMessage("This function can only be used on objects that are "
4904 "cells and not on objects which bound cells."));
4905 AssertIndexRange(face, accessor.n_faces());
4906
4907 if (dim == 1)
4908 Assert(combined_orientation == numbers::default_geometric_orientation,
4909 ExcMessage("In 1d, faces do not have an orientation, so the "
4910 "only valid value is the default."));
4911 else if (dim == 2)
4912 Assert(combined_orientation ==
4914 combined_orientation == numbers::reverse_line_orientation,
4915 ExcMessage(
4916 "In 2d, the only valid values of the combined orientation "
4917 "are the standard orientation or the reversed line "
4918 "orientation."));
4919
4920 // face_orientations is not set up in 1d
4921 if (dim != 1)
4922 accessor.tria->levels[accessor.present_level]
4923 ->face_orientations.set_combined_orientation(
4924 accessor.present_index * ReferenceCells::max_n_faces<dim>() +
4925 face,
4926 combined_orientation);
4927 }
4928
4929
4930
4931 template <int dim, int spacedim>
4932 static std::array<unsigned int, 1>
4933 get_line_indices_of_cell(const TriaAccessor<1, dim, spacedim> &)
4934 {
4936 return {};
4937 }
4938
4939
4940
4941 template <int structdim, int dim, int spacedim>
4942 static std::array<unsigned int, 4>
4943 get_line_indices_of_cell(const TriaAccessor<2, dim, spacedim> &cell)
4944 {
4945 // For 2d cells the access cell->line_orientation() is already
4946 // efficient
4947 std::array<unsigned int, 4> line_indices = {};
4948 for (const unsigned int line : cell.line_indices())
4949 line_indices[line] = cell.line_index(line);
4950 return line_indices;
4951 }
4952
4957 template <int structdim, int dim, int spacedim>
4958 static std::array<unsigned int, 12>
4959 get_line_indices_of_cell(
4961 {
4962 std::array<unsigned int, 12> line_indices = {};
4963
4964 // For hexahedra, the classical access via quads -> lines is too
4965 // inefficient. Unroll this code here to allow the compiler to inline
4966 // the necessary functions.
4967 const auto ref_cell = cell.reference_cell();
4968 if (ref_cell == ReferenceCells::Hexahedron)
4969 {
4970 for (unsigned int f = 4; f < 6; ++f)
4971 {
4972 const auto orientation =
4973 cell.get_triangulation()
4974 .levels[cell.level()]
4975 ->face_orientations.get_combined_orientation(
4976 cell.index() * ReferenceCells::max_n_faces<dim>() + f);
4977
4978 // It might seem superfluous to spell out the four indices
4979 // that get later consumed by a for loop over these four
4980 // elements; however, for the compiler it is easier to inline
4981 // the statement of standard_to_real_face_line() when next to
4982 // each other, as opposed to be interleaved with a
4983 // line_index() call.
4984 const std::array<unsigned int, 4> my_indices{
4985 {ref_cell.standard_to_real_face_line(0, f, orientation),
4986 ref_cell.standard_to_real_face_line(1, f, orientation),
4987 ref_cell.standard_to_real_face_line(2, f, orientation),
4988 ref_cell.standard_to_real_face_line(3, f, orientation)}};
4989 const auto quad = cell.quad(f);
4990 for (unsigned int l = 0; l < 4; ++l)
4991 line_indices[4 * (f - 4) + l] =
4992 quad->line_index(my_indices[l]);
4993 }
4994 for (unsigned int f = 0; f < 2; ++f)
4995 {
4996 const auto orientation =
4997 cell.get_triangulation()
4998 .levels[cell.level()]
4999 ->face_orientations.get_combined_orientation(
5000 cell.index() * ReferenceCells::max_n_faces<dim>() + f);
5001 const std::array<unsigned int, 2> my_indices{
5002 {ref_cell.standard_to_real_face_line(0, f, orientation),
5003 ref_cell.standard_to_real_face_line(1, f, orientation)}};
5004 const auto quad = cell.quad(f);
5005 line_indices[8 + f] = quad->line_index(my_indices[0]);
5006 line_indices[10 + f] = quad->line_index(my_indices[1]);
5007 }
5008 }
5009 else if (ref_cell == ReferenceCells::Tetrahedron)
5010 {
5011 std::array<unsigned int, 3> orientations{
5014 cell.combined_face_orientation(2)}};
5015 const std::array<unsigned int, 6> my_indices{
5016 {ref_cell.standard_to_real_face_line(0, 0, orientations[0]),
5017 ref_cell.standard_to_real_face_line(1, 0, orientations[0]),
5018 ref_cell.standard_to_real_face_line(2, 0, orientations[0]),
5019 ref_cell.standard_to_real_face_line(1, 1, orientations[1]),
5020 ref_cell.standard_to_real_face_line(2, 1, orientations[1]),
5021 ref_cell.standard_to_real_face_line(1, 2, orientations[2])}};
5022 line_indices[0] = cell.quad(0)->line_index(my_indices[0]);
5023 line_indices[1] = cell.quad(0)->line_index(my_indices[1]);
5024 line_indices[2] = cell.quad(0)->line_index(my_indices[2]);
5025 line_indices[3] = cell.quad(1)->line_index(my_indices[3]);
5026 line_indices[4] = cell.quad(1)->line_index(my_indices[4]);
5027 line_indices[5] = cell.quad(2)->line_index(my_indices[5]);
5028 }
5029 else
5030 // For other shapes (wedges, pyramids), we do not currently
5031 // implement an optimized function.
5032 for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l)
5033 line_indices[l] = cell.line_index(l);
5034
5035 return line_indices;
5036 }
5037
5038
5039
5044 template <int dim, int spacedim>
5045 static std::array<types::geometric_orientation, 1>
5046 get_line_orientations_of_cell(const TriaAccessor<1, dim, spacedim> &)
5047 {
5049 return {};
5050 }
5051
5052
5053
5058 template <int dim, int spacedim>
5059 static std::array<types::geometric_orientation, 4>
5060 get_line_orientations_of_cell(const TriaAccessor<2, dim, spacedim> &cell)
5061 {
5062 // For 2d cells the access cell->line_orientation() is already
5063 // efficient
5064 std::array<types::geometric_orientation, 4> line_orientations = {};
5065 for (const unsigned int line : cell.line_indices())
5066 line_orientations[line] = cell.line_orientation(line);
5067 return line_orientations;
5068 }
5069
5070
5071
5076 template <int dim, int spacedim>
5077 static std::array<types::geometric_orientation, 12>
5078 get_line_orientations_of_cell(const TriaAccessor<3, dim, spacedim> &cell)
5079 {
5080 std::array<types::geometric_orientation, 12> line_orientations = {};
5081
5082 // For hexahedra, the classical access via quads -> lines is too
5083 // inefficient. Unroll this code here to allow the compiler to inline
5084 // the necessary functions.
5085 const auto ref_cell = cell.reference_cell();
5086 if (ref_cell == ReferenceCells::Hexahedron)
5087 {
5088 for (unsigned int f = 4; f < 6; ++f)
5089 {
5090 const auto orientation =
5091 cell.get_triangulation()
5092 .levels[cell.level()]
5093 ->face_orientations.get_combined_orientation(
5094 cell.index() * ReferenceCells::max_n_faces<dim>() + f);
5095
5096 // It might seem superfluous to spell out the four indices and
5097 // orientations that get later consumed by a for loop over
5098 // these four elements; however, for the compiler it is easier
5099 // to inline the statement of standard_to_real_face_line()
5100 // when next to each other, as opposed to be interleaved with
5101 // a line_index() call.
5102 const std::array<unsigned int, 4> my_indices{
5103 {ref_cell.standard_to_real_face_line(0, f, orientation),
5104 ref_cell.standard_to_real_face_line(1, f, orientation),
5105 ref_cell.standard_to_real_face_line(2, f, orientation),
5106 ref_cell.standard_to_real_face_line(3, f, orientation)}};
5107 const auto quad = cell.quad(f);
5108 const std::array<types::geometric_orientation, 4>
5109 my_orientations{{ref_cell.face_to_cell_line_orientation(
5110 0,
5111 f,
5112 orientation,
5113 quad->line_orientation(my_indices[0])),
5114 ref_cell.face_to_cell_line_orientation(
5115 1,
5116 f,
5117 orientation,
5118 quad->line_orientation(my_indices[1])),
5119 ref_cell.face_to_cell_line_orientation(
5120 2,
5121 f,
5122 orientation,
5123 quad->line_orientation(my_indices[2])),
5124 ref_cell.face_to_cell_line_orientation(
5125 3,
5126 f,
5127 orientation,
5128 quad->line_orientation(my_indices[3]))}};
5129 for (unsigned int l = 0; l < 4; ++l)
5130 line_orientations[4 * (f - 4) + l] = my_orientations[l];
5131 }
5132 for (unsigned int f = 0; f < 2; ++f)
5133 {
5134 const auto orientation =
5135 cell.get_triangulation()
5136 .levels[cell.level()]
5137 ->face_orientations.get_combined_orientation(
5138 cell.index() * ReferenceCells::max_n_faces<3>() + f);
5139 const std::array<unsigned int, 2> my_indices{
5140 {ref_cell.standard_to_real_face_line(0, f, orientation),
5141 ref_cell.standard_to_real_face_line(1, f, orientation)}};
5142 const auto quad = cell.quad(f);
5143 const std::array<types::geometric_orientation, 2>
5144 my_orientations{{ref_cell.face_to_cell_line_orientation(
5145 0,
5146 f,
5147 orientation,
5148 quad->line_orientation(my_indices[0])),
5149 ref_cell.face_to_cell_line_orientation(
5150 1,
5151 f,
5152 orientation,
5153 quad->line_orientation(my_indices[1]))}};
5154 line_orientations[8 + f] = my_orientations[0];
5155 line_orientations[10 + f] = my_orientations[1];
5156 }
5157 }
5158 else if (ref_cell == ReferenceCells::Tetrahedron)
5159 {
5160 std::array<unsigned int, 3> orientations{
5163 cell.combined_face_orientation(2)}};
5164 const std::array<unsigned int, 6> my_indices{
5165 {ref_cell.standard_to_real_face_line(0, 0, orientations[0]),
5166 ref_cell.standard_to_real_face_line(1, 0, orientations[0]),
5167 ref_cell.standard_to_real_face_line(2, 0, orientations[0]),
5168 ref_cell.standard_to_real_face_line(1, 1, orientations[1]),
5169 ref_cell.standard_to_real_face_line(2, 1, orientations[1]),
5170 ref_cell.standard_to_real_face_line(1, 2, orientations[2])}};
5171 line_orientations[0] = ref_cell.face_to_cell_line_orientation(
5172 0,
5173 0,
5174 orientations[0],
5175 cell.quad(0)->line_orientation(my_indices[0]));
5176 line_orientations[1] = ref_cell.face_to_cell_line_orientation(
5177 1,
5178 0,
5179 orientations[0],
5180 cell.quad(0)->line_orientation(my_indices[1]));
5181 line_orientations[2] = ref_cell.face_to_cell_line_orientation(
5182 2,
5183 0,
5184 orientations[0],
5185 cell.quad(0)->line_orientation(my_indices[2]));
5186 line_orientations[3] = ref_cell.face_to_cell_line_orientation(
5187 1,
5188 1,
5189 orientations[1],
5190 cell.quad(1)->line_orientation(my_indices[3]));
5191 line_orientations[4] = ref_cell.face_to_cell_line_orientation(
5192 2,
5193 1,
5194 orientations[1],
5195 cell.quad(1)->line_orientation(my_indices[4]));
5196 line_orientations[5] = ref_cell.face_to_cell_line_orientation(
5197 1,
5198 2,
5199 orientations[2],
5200 cell.quad(2)->line_orientation(my_indices[5]));
5201 }
5202 else
5203 // For other shapes (wedges, pyramids), we do not currently
5204 // implement an optimized function
5205 for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l)
5206 line_orientations[l] = cell.line_orientation(l);
5207
5208 return line_orientations;
5209 }
5210 };
5211 } // namespace TriaAccessorImplementation
5212} // namespace internal
5213
5214
5215
5216template <int structdim, int dim, int spacedim>
5218 const Triangulation<dim, spacedim> *parent,
5219 const int level,
5220 const int index,
5221 const AccessorData *local_data)
5222 : TriaAccessorBase<structdim, dim, spacedim>(parent, level, index, local_data)
5223{}
5224
5225
5226
5227template <int structdim, int dim, int spacedim>
5228inline bool
5230{
5231 Assert(this->state() == IteratorState::valid,
5232 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
5233 *this));
5234 return this->objects().used[this->present_index];
5235}
5236
5237
5238
5239template <int structdim, int dim, int spacedim>
5242 const unsigned int i) const
5243{
5245 0,
5246 vertex_index(i));
5247}
5248
5249
5250
5251template <int structdim, int dim, int spacedim>
5252inline ReferenceCell
5254{
5255 if (structdim == 0)
5257 else if (structdim == 1)
5258 return ReferenceCells::Line;
5259 else if (structdim == dim)
5260 return this->tria->levels[this->present_level]
5261 ->reference_cell[this->present_index];
5262 else
5263 return this->tria->faces->get_quad_type(this->present_index);
5264}
5265
5266
5267
5268template <int structdim, int dim, int spacedim>
5269inline unsigned int
5271 const unsigned int corner) const
5272{
5273 AssertIndexRange(corner, this->n_vertices());
5274
5275 if constexpr (structdim == 1)
5276 {
5277 // This branch needs to be first (and not combined with the structdim ==
5278 // dim branch) so that we can get line vertex indices when setting up the
5279 // cell vertex index cache
5280 return this->objects()
5281 .cells[this->present_index * ReferenceCells::max_n_faces<1>() + corner];
5282 }
5283 else if constexpr (structdim == dim)
5284 {
5285 // This branch should only be used after the cell vertex index cache is
5286 // set up
5287 const auto my_index = static_cast<std::size_t>(this->present_index) *
5288 ReferenceCells::max_n_vertices<dim>();
5289 AssertIndexRange(my_index + corner,
5290 this->tria->levels[this->present_level]
5291 ->cell_vertex_indices_cache.size());
5292 const unsigned int vertex_index =
5293 this->tria->levels[this->present_level]
5294 ->cell_vertex_indices_cache[my_index + corner];
5296 return vertex_index;
5297 }
5298 else if constexpr (structdim == 2)
5299 {
5300 const auto [line_index, vertex_index] =
5301 this->reference_cell().standard_vertex_to_face_and_vertex_index(corner);
5302 const auto vertex_within_line_index =
5303 this->reference_cell().standard_to_real_face_vertex(
5304 vertex_index, line_index, this->line_orientation(line_index));
5305
5306 return this->line(line_index)->vertex_index(vertex_within_line_index);
5307 }
5308 else
5309 {
5312 }
5313}
5314
5315
5316
5317template <int structdim, int dim, int spacedim>
5318inline Point<spacedim> &
5319TriaAccessor<structdim, dim, spacedim>::vertex(const unsigned int i) const
5320{
5321 return const_cast<Point<spacedim> &>(this->tria->vertices[vertex_index(i)]);
5322}
5323
5324
5325
5326template <int structdim, int dim, int spacedim>
5327inline typename ::internal::TriangulationImplementation::
5328 Iterators<dim, spacedim>::line_iterator
5329 TriaAccessor<structdim, dim, spacedim>::line(const unsigned int i) const
5330{
5331 // checks happen in line_index
5332 return typename ::internal::TriangulationImplementation::
5333 Iterators<dim, spacedim>::line_iterator(this->tria, 0, line_index(i));
5334}
5335
5336
5337
5338template <int structdim, int dim, int spacedim>
5339inline unsigned int
5341{
5342 (void)i;
5343 AssertIndexRange(i, this->n_lines());
5344 Assert(structdim != 1,
5345 ExcMessage("You can't ask for the index of a line bounding a "
5346 "one-dimensional cell because it is not bounded by "
5347 "lines."));
5348
5349 if constexpr (structdim == 2)
5350 {
5351 return this->objects()
5352 .cells[this->present_index * ReferenceCells::max_n_faces<2>() + i];
5353 }
5354 else if constexpr (structdim == 3)
5355 {
5356 const auto [face_index, line_index] =
5357 this->reference_cell().standard_line_to_face_and_line_index(i);
5358 const auto line_within_face_index =
5359 this->reference_cell().standard_to_real_face_line(
5360 line_index, face_index, this->combined_face_orientation(face_index));
5361
5362 return this->quad(face_index)->line_index(line_within_face_index);
5363 }
5364
5367}
5368
5369
5370
5371template <int structdim, int dim, int spacedim>
5372inline typename ::internal::TriangulationImplementation::
5373 Iterators<dim, spacedim>::quad_iterator
5374 TriaAccessor<structdim, dim, spacedim>::quad(const unsigned int i) const
5375{
5376 // checks happen in quad_index
5377 return typename ::internal::TriangulationImplementation::
5378 Iterators<dim, spacedim>::quad_iterator(this->tria, 0, quad_index(i));
5379}
5380
5381
5382
5383template <int structdim, int dim, int spacedim>
5384inline unsigned int
5386{
5387 Assert(structdim == 3,
5388 ExcMessage("You can't ask for the index of a quad bounding "
5389 "a one- or two-dimensional cell because it is not "
5390 "bounded by quads."));
5391 // work around a bogus GCC-9 warning which considers i unused except in 3d
5392 (void)i;
5393 if constexpr (structdim == 3)
5394 {
5395 AssertIndexRange(i, n_faces());
5396 return this->tria->levels[this->present_level]
5397 ->cells
5398 .cells[this->present_index * ReferenceCells::max_n_faces<3>() + i];
5399 }
5400 else
5402}
5403
5404
5405
5406template <int structdim, int dim, int spacedim>
5409 const unsigned int face) const
5410{
5412 AssertIndexRange(face, n_faces());
5413 Assert(structdim == dim,
5414 ExcMessage("This function can only be used on objects "
5415 "that are cells, but not on faces or edges "
5416 "that bound cells."));
5417 // work around a bogus GCC-9 warning which considers face unused except in 3d
5418 (void)face;
5419
5420 if constexpr (structdim == 1)
5422 else if constexpr (structdim == 2)
5423 {
5424 // if all elements are quads (or if we have a very special consistently
5425 // oriented triangular mesh) then we do not store this array
5426 if (this->tria->levels[this->present_level]
5427 ->face_orientations.n_objects() == 0)
5429 else
5430 return this->tria->levels[this->present_level]
5431 ->face_orientations.get_orientation(
5432 this->present_index * ReferenceCells::max_n_faces<dim>() + face);
5433 }
5434 else
5435 return this->tria->levels[this->present_level]
5436 ->face_orientations.get_combined_orientation(
5437 this->present_index * ReferenceCells::max_n_faces<dim>() + face);
5438}
5439
5440
5441
5442template <int structdim, int dim, int spacedim>
5443inline bool
5445 const unsigned int face) const
5446{
5448 AssertIndexRange(face, n_faces());
5449 Assert(structdim == dim,
5450 ExcMessage("This function can only be used on objects "
5451 "that are cells, but not on faces or edges "
5452 "that bound cells."));
5453 // work around a bogus GCC-9 warning which considers face unused in 1d
5454 (void)face;
5455
5456 if constexpr (structdim == 1)
5457 // in 1d 'faces' are vertices and those are always consistently oriented
5458 return true;
5459 else if constexpr (structdim == 2)
5460 return this->line_orientation(face) ==
5462 else
5463 return this->tria->levels[this->present_level]
5464 ->face_orientations.get_orientation(
5465 this->present_index * ReferenceCells::max_n_faces<structdim>() + face);
5466}
5467
5468
5469
5470template <int structdim, int dim, int spacedim>
5471inline bool
5472TriaAccessor<structdim, dim, spacedim>::face_flip(const unsigned int face) const
5473{
5475 Assert(structdim == dim,
5476 ExcMessage("This function can only be used on objects "
5477 "that are cells, but not on faces or edges "
5478 "that bound cells."));
5479 AssertIndexRange(face, n_faces());
5480 // work around a bogus GCC-9 warning which considers face unused except in 3d
5481 (void)face;
5482
5483 if constexpr (structdim == 3)
5484 return this->tria->levels[this->present_level]->face_orientations.get_flip(
5485 this->present_index * ReferenceCells::max_n_faces<structdim>() + face);
5486 else
5487 // In 1d and 2d, face_flip is always false as faces can only be
5488 // 'flipped' in 3d.
5489 return false;
5490}
5491
5492
5493template <int structdim, int dim, int spacedim>
5494inline bool
5496 const unsigned int face) const
5497{
5499 Assert(structdim == dim,
5500 ExcMessage("This function can only be used on objects "
5501 "that are cells, but not on faces or edges "
5502 "that bound cells."));
5503 AssertIndexRange(face, n_faces());
5504 // work around a bogus GCC-9 warning which considers face unused except in 3d
5505 (void)face;
5506
5507 if constexpr (structdim == 3)
5508 return this->tria->levels[this->present_level]
5509 ->face_orientations.get_rotation(
5510 this->present_index * ReferenceCells::max_n_faces<structdim>() + face);
5511 else
5512 // In 1d and 2d, face_rotation is always false as faces can only be
5513 // 'rotated' in 3d.
5514 return false;
5515}
5516
5517
5518
5519template <int structdim, int dim, int spacedim>
5522 const unsigned int line) const
5523{
5525 AssertIndexRange(line, this->n_lines());
5526 // work around a bogus GCC-9 warning which considers line unused in 1d
5527 (void)line;
5528
5529 if constexpr (structdim == 1)
5531 else if constexpr (structdim == 2 && dim == 2)
5532 // lines in 2d are faces
5533 {
5534 const auto combined_orientation = combined_face_orientation(line);
5535 Assert(combined_orientation == numbers::default_geometric_orientation ||
5536 combined_orientation == numbers::reverse_line_orientation,
5538 return combined_orientation;
5539 }
5540 else if constexpr (structdim == 2 && dim == 3)
5541 {
5542 // line orientations in 3d are stored in their own array as bools: here
5543 // 'true' is the default orientation and 'false' is the reversed one
5544 // (which matches set_line_orientation())
5545 const auto index =
5546 this->present_index * ReferenceCells::max_n_lines<2>() + line;
5547 Assert(index < this->tria->faces->quads_line_orientations.size(),
5549 return this->tria->faces->quads_line_orientations[index] ?
5552 }
5553 else if constexpr (structdim == 3 && dim == 3)
5554 {
5555 const auto reference_cell = this->reference_cell();
5556 // First pick a face on which this line is a part of, and the
5557 // index of the line within.
5558 const auto [face_index, line_index] =
5559 reference_cell.standard_line_to_face_and_line_index(line);
5560 const auto line_within_face_index =
5561 reference_cell.standard_to_real_face_line(
5562 line_index, face_index, this->combined_face_orientation(face_index));
5563
5564 // Then query how that line is oriented within that face:
5565 return reference_cell.face_to_cell_line_orientation(
5566 line_index,
5567 face_index,
5568 this->combined_face_orientation(face_index),
5569 this->quad(face_index)->line_orientation(line_within_face_index));
5570 }
5571 else
5572 {
5574 return false;
5575 }
5576}
5577
5578
5579
5580template <int structdim, int dim, int spacedim>
5581inline void
5583 const unsigned int line,
5584 const types::geometric_orientation value) const
5585{
5587 AssertIndexRange(line, this->n_lines());
5588 Assert(dim != 1,
5589 ExcMessage("In 1d lines are cells and thus do not need to have their "
5590 "orientations set."));
5591 Assert(dim != 2,
5592 ExcMessage("In 2d lines are faces, and, for compatibility with other "
5593 "dimensions, their orientations should be set via "
5594 "set_combined_orientation()."));
5595 // work around a bogus GCC-9 warning which considers line and value unused
5596 // except in 3d
5597 (void)line;
5598 (void)value;
5599
5600 if constexpr (dim == 3)
5601 {
5602 // We set line orientations per face, not per cell, so this only works for
5603 // faces in 3d.
5604 Assert(structdim == 2, ExcNotImplemented());
5605 const auto index =
5606 this->present_index * ReferenceCells::max_n_lines<2>() + line;
5607 Assert(index < this->tria->faces->quads_line_orientations.size(),
5609 this->tria->faces->quads_line_orientations[index] =
5611 }
5612}
5613
5614
5615
5616template <int structdim, int dim, int spacedim>
5617inline void
5619 const unsigned int face,
5620 const types::geometric_orientation combined_orientation) const
5621{
5623 AssertIndexRange(face, this->n_faces());
5624
5625 ::internal::TriaAccessorImplementation::Implementation::
5626 set_combined_face_orientation(*this, face, combined_orientation);
5627}
5628
5629
5630
5631template <int structdim, int dim, int spacedim>
5632void
5634{
5635 Assert(this->state() == IteratorState::valid,
5636 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
5637 *this));
5638 this->objects().used[this->present_index] = true;
5639}
5640
5641
5642
5643template <int structdim, int dim, int spacedim>
5644void
5646{
5647 Assert(this->state() == IteratorState::valid,
5648 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
5649 *this));
5650 this->objects().used[this->present_index] = false;
5651}
5652
5653
5654template <int structdim, int dim, int spacedim>
5655int
5657{
5659 AssertIndexRange(i, n_children());
5660
5661 // each set of two children are stored
5662 // consecutively, so we only have to find
5663 // the location of the set of children
5664 const unsigned int n_sets_of_two =
5666 return this->objects().children[n_sets_of_two * this->present_index + i / 2] +
5667 i % 2;
5668}
5669
5670
5671
5672template <int structdim, int dim, int spacedim>
5673int
5675 const unsigned int i) const
5676{
5678
5679 switch (structdim)
5680 {
5681 case 1:
5682 return child_index(i);
5683 case 2:
5684 {
5685 const RefinementCase<2> this_refinement_case(
5686 static_cast<std::uint8_t>(refinement_case()));
5687
5688 Assert(this_refinement_case != RefinementCase<2>::no_refinement,
5690
5691 if (this_refinement_case == RefinementCase<2>::cut_xy)
5692 return child_index(i);
5693 else if ((this_refinement_case == RefinementCase<2>::cut_x) &&
5694 (child(i % 2)->refinement_case() ==
5696 return child(i % 2)->child_index(i / 2);
5697 else if ((this_refinement_case == RefinementCase<2>::cut_y) &&
5698 (child(i / 2)->refinement_case() ==
5700 return child(i / 2)->child_index(i % 2);
5701 else
5702 Assert(
5703 false,
5704 ExcMessage(
5705 "This cell has no grandchildren equivalent to isotropic refinement"));
5706 break;
5707 }
5708
5709 case 3:
5711 }
5712 return -1;
5713}
5714
5715
5716
5717template <int structdim, int dim, int spacedim>
5720{
5721 Assert(this->state() == IteratorState::valid,
5722 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
5723 *this));
5724
5725 switch (structdim)
5726 {
5727 case 1:
5729 this->objects().children[this->present_index] != -1 ?
5730 // cast the branches
5731 // here first to uchar
5732 // and then (above) to
5733 // RefinementCase<structdim>
5734 // so that the
5735 // conversion is valid
5736 // even for the case
5737 // structdim>1 (for
5738 // which this part of
5739 // the code is dead
5740 // anyway)
5741 static_cast<std::uint8_t>(RefinementCase<1>::cut_x) :
5742 static_cast<std::uint8_t>(RefinementCase<1>::no_refinement)));
5743
5744 default:
5745 Assert(static_cast<unsigned int>(this->present_index) <
5746 this->objects().refinement_cases.size(),
5747 ExcIndexRange(this->present_index,
5748 0,
5749 this->objects().refinement_cases.size()));
5750
5751 return (static_cast<RefinementCase<structdim>>(
5752 this->objects().refinement_cases[this->present_index]));
5753 }
5754}
5755
5756
5757
5758template <int structdim, int dim, int spacedim>
5760TriaAccessor<structdim, dim, spacedim>::child(const unsigned int i) const
5761
5762{
5763 // checking of 'i' happens in child_index
5765 this->tria, (dim == structdim ? this->level() + 1 : 0), child_index(i));
5766
5767 Assert((q.state() == IteratorState::past_the_end) || q->used(),
5769
5770 return q;
5771}
5772
5773
5774
5775template <int structdim, int dim, int spacedim>
5776inline unsigned int
5779{
5780 const auto n_children = this->n_children();
5781 for (unsigned int child_n = 0; child_n < n_children; ++child_n)
5782 if (this->child(child_n) == child)
5783 return child_n;
5784
5785 Assert(false,
5786 ExcMessage("The given child is not a child of the current object."));
5788}
5789
5790
5791
5792template <int structdim, int dim, int spacedim>
5795 const unsigned int i) const
5796{
5797 // checking of 'i' happens in child() or
5798 // child_index() called below
5799 switch (structdim)
5800 {
5801 case 1:
5802 // no anisotropic refinement in 1d
5803 return child(i);
5804
5805 case 2:
5806 {
5807 const RefinementCase<2> this_refinement_case(
5808 static_cast<std::uint8_t>(refinement_case()));
5809
5810 Assert(this_refinement_case != RefinementCase<2>::no_refinement,
5812
5813 if (this_refinement_case == RefinementCase<2>::cut_xy)
5814 return child(i);
5815 else if ((this_refinement_case == RefinementCase<2>::cut_x) &&
5816 (child(i % 2)->refinement_case() ==
5818 return child(i % 2)->child(i / 2);
5819 else if ((this_refinement_case == RefinementCase<2>::cut_y) &&
5820 (child(i / 2)->refinement_case() ==
5822 return child(i / 2)->child(i % 2);
5823 else
5824 Assert(
5825 false,
5826 ExcMessage(
5827 "This cell has no grandchildren equivalent to isotropic refinement"));
5828 break;
5829 }
5830
5831 default:
5833 }
5834 // we don't get here but have to return
5835 // something...
5836 return child(0);
5837}
5838
5839
5840
5841template <int structdim, int dim, int spacedim>
5842inline bool
5844{
5845 Assert(this->state() == IteratorState::valid,
5846 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
5847 *this));
5848
5849 // each set of two children are stored
5850 // consecutively, so we only have to find
5851 // the location of the set of children
5852 const unsigned int n_sets_of_two =
5854 return (this->objects().children[n_sets_of_two * this->present_index] != -1);
5855}
5856
5857
5858
5859template <int structdim, int dim, int spacedim>
5860inline unsigned int
5862{
5863 Assert(this->state() == IteratorState::valid,
5864 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
5865 *this));
5869 else
5870 return GeometryInfo<structdim>::n_children(refinement_case());
5871}
5872
5873
5874
5875template <int structdim, int dim, int spacedim>
5876inline void
5878 const RefinementCase<structdim> &refinement_case) const
5879{
5880 Assert(this->state() == IteratorState::valid,
5881 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
5882 *this));
5883 Assert(static_cast<unsigned int>(this->present_index) <
5884 this->objects().refinement_cases.size(),
5885 ExcIndexRange(this->present_index,
5886 0,
5887 this->objects().refinement_cases.size()));
5888
5889 this->objects().refinement_cases[this->present_index] = refinement_case;
5890}
5891
5892
5893template <int structdim, int dim, int spacedim>
5894inline void
5896{
5897 Assert(this->state() == IteratorState::valid,
5898 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
5899 *this));
5900 Assert(static_cast<unsigned int>(this->present_index) <
5901 this->objects().refinement_cases.size(),
5902 ExcIndexRange(this->present_index,
5903 0,
5904 this->objects().refinement_cases.size()));
5905
5906 this->objects().refinement_cases[this->present_index] =
5908}
5909
5910
5911
5912template <int structdim, int dim, int spacedim>
5913void
5915 const int index) const
5916{
5919
5920 // each set of two children are stored
5921 // consecutively, so we only have to find
5922 // the location of the set of children
5923 const unsigned int n_sets_of_two =
5925
5926 Assert(
5927 // clearing the child index for a cell
5928 (index == -1) ||
5929 // if setting the child index for the i'th child (with i==0),
5930 // then the index must be a non-negative number
5931 (i == 0 && !this->has_children() && (index >= 0)) ||
5932 // if setting the child index for the i'th child (with i>0),
5933 // then the previously stored index must be the invalid
5934 // index
5935 (i > 0 && this->has_children() && (index >= 0) &&
5936 this->objects().children[n_sets_of_two * this->present_index + i / 2] ==
5937 -1),
5939
5940 this->objects().children[n_sets_of_two * this->present_index + i / 2] = index;
5941}
5942
5943
5944
5945template <int structdim, int dim, int spacedim>
5946void
5948{
5949 // each set of two children are stored
5950 // consecutively, so we only have to find
5951 // the location of the set of children
5952 const unsigned int n_sets_of_two =
5954
5955 for (unsigned int i = 0; i < n_sets_of_two; ++i)
5956 set_children(2 * i, -1);
5957}
5958
5959
5960
5961template <int structdim, int dim, int spacedim>
5962inline bool
5964{
5966 return this->objects().user_flags[this->present_index];
5967}
5968
5969
5970
5971template <int structdim, int dim, int spacedim>
5972inline void
5974{
5976 this->objects().user_flags[this->present_index] = true;
5977}
5978
5979
5980
5981template <int structdim, int dim, int spacedim>
5982inline void
5984{
5986 this->objects().user_flags[this->present_index] = false;
5987}
5988
5989
5990
5991template <int structdim, int dim, int spacedim>
5992void
5994{
5995 set_user_flag();
5996
5997 if (this->has_children())
5998 for (unsigned int c = 0; c < this->n_children(); ++c)
5999 this->child(c)->recursively_set_user_flag();
6000}
6001
6002
6003
6004template <int structdim, int dim, int spacedim>
6005void
6007{
6008 clear_user_flag();
6009
6010 if (this->has_children())
6011 for (unsigned int c = 0; c < this->n_children(); ++c)
6012 this->child(c)->recursively_clear_user_flag();
6013}
6014
6015
6016
6017template <int structdim, int dim, int spacedim>
6018void
6020{
6022 this->objects().clear_user_data(this->present_index);
6023}
6024
6025
6026
6027template <int structdim, int dim, int spacedim>
6028void
6030{
6032 this->objects().user_pointer(this->present_index) = p;
6033}
6034
6035
6036
6037template <int structdim, int dim, int spacedim>
6038void
6040{
6042 this->objects().user_pointer(this->present_index) = nullptr;
6043}
6044
6045
6046
6047template <int structdim, int dim, int spacedim>
6048void *
6050{
6052 return this->objects().user_pointer(this->present_index);
6053}
6054
6055
6056
6057template <int structdim, int dim, int spacedim>
6058void
6060 void *p) const
6061{
6062 set_user_pointer(p);
6063
6064 if (this->has_children())
6065 for (unsigned int c = 0; c < this->n_children(); ++c)
6066 this->child(c)->recursively_set_user_pointer(p);
6067}
6068
6069
6070
6071template <int structdim, int dim, int spacedim>
6072void
6074{
6075 clear_user_pointer();
6076
6077 if (this->has_children())
6078 for (unsigned int c = 0; c < this->n_children(); ++c)
6079 this->child(c)->recursively_clear_user_pointer();
6080}
6081
6082
6083
6084template <int structdim, int dim, int spacedim>
6085void
6087 const unsigned int p) const
6088{
6090 this->objects().user_index(this->present_index) = p;
6091}
6092
6093
6094
6095template <int structdim, int dim, int spacedim>
6096void
6098{
6100 this->objects().user_index(this->present_index) = 0;
6101}
6102
6103
6104
6105template <int structdim, int dim, int spacedim>
6106unsigned int
6108{
6110 return this->objects().user_index(this->present_index);
6111}
6112
6113
6114
6115template <int structdim, int dim, int spacedim>
6116void
6118 const unsigned int p) const
6119{
6120 set_user_index(p);
6121
6122 if (this->has_children())
6123 for (unsigned int c = 0; c < this->n_children(); ++c)
6124 this->child(c)->recursively_set_user_index(p);
6125}
6126
6127
6128
6129template <int structdim, int dim, int spacedim>
6130void
6132{
6133 clear_user_index();
6134
6135 if (this->has_children())
6136 for (unsigned int c = 0; c < this->n_children(); ++c)
6137 this->child(c)->recursively_clear_user_index();
6138}
6139
6140
6141
6142template <int structdim, int dim, int spacedim>
6143inline unsigned int
6145{
6146 if (!this->has_children())
6147 return 0;
6148
6149 unsigned int max_depth = 1;
6150 for (unsigned int c = 0; c < n_children(); ++c)
6151 max_depth = std::max(max_depth, child(c)->max_refinement_depth() + 1);
6152 return max_depth;
6153}
6154
6155
6156
6157template <int structdim, int dim, int spacedim>
6158unsigned int
6160{
6161 if (!this->has_children())
6162 return 1;
6163 else
6164 {
6165 unsigned int sum = 0;
6166 for (unsigned int c = 0; c < n_children(); ++c)
6167 sum += this->child(c)->n_active_descendants();
6168 return sum;
6169 }
6170}
6171
6172
6173
6174template <int structdim, int dim, int spacedim>
6177{
6178 Assert(structdim < dim, ExcImpossibleInDim(dim));
6180
6181 return this->objects()
6182 .boundary_or_material_id[this->present_index]
6183 .boundary_id;
6184}
6185
6186
6187
6188template <int structdim, int dim, int spacedim>
6189void
6191 const types::boundary_id boundary_ind) const
6192{
6193 Assert(structdim < dim, ExcImpossibleInDim(dim));
6196 ExcMessage("You are trying to set the boundary_id to an invalid "
6197 "value (numbers::internal_face_boundary_id is reserved)."));
6198 Assert(this->at_boundary(),
6199 ExcMessage("You are trying to set the boundary_id of an "
6200 "internal object, which is not allowed!"));
6201
6202 this->objects().boundary_or_material_id[this->present_index].boundary_id =
6203 boundary_ind;
6204}
6205
6206
6207
6208template <int structdim, int dim, int spacedim>
6209void
6211 const types::boundary_id boundary_ind) const
6212{
6213 Assert(structdim < dim, ExcImpossibleInDim(dim));
6215
6216 this->objects().boundary_or_material_id[this->present_index].boundary_id =
6217 boundary_ind;
6218}
6219
6220
6221
6222template <int structdim, int dim, int spacedim>
6223void
6225 const types::boundary_id boundary_ind) const
6226{
6227 set_boundary_id(boundary_ind);
6228
6229 switch (structdim)
6230 {
6231 case 1:
6232 // 1d objects have no sub-objects
6233 // where we have to do anything
6234 break;
6235
6236 case 2:
6237 // for boundary quads also set
6238 // boundary_id of bounding lines
6239 for (unsigned int i = 0; i < this->n_lines(); ++i)
6240 this->line(i)->set_boundary_id(boundary_ind);
6241 break;
6242
6243 default:
6245 }
6246}
6247
6248
6249
6250template <int structdim, int dim, int spacedim>
6251bool
6253{
6254 // error checking is done
6255 // in boundary_id()
6257}
6258
6259
6260
6261template <int structdim, int dim, int spacedim>
6264{
6266 return this->tria->get_manifold(this->manifold_id());
6267}
6268
6269
6270template <int structdim, int dim, int spacedim>
6273{
6275
6276 return this->objects().manifold_id[this->present_index];
6277}
6278
6279
6280
6281template <int structdim, int dim, int spacedim>
6282void
6284 const types::manifold_id manifold_ind) const
6285{
6287
6288 this->objects().manifold_id[this->present_index] = manifold_ind;
6289}
6290
6291
6292template <int structdim, int dim, int spacedim>
6293void
6295 const types::manifold_id manifold_ind) const
6296{
6297 set_manifold_id(manifold_ind);
6298
6299 if (this->has_children())
6300 for (unsigned int c = 0; c < this->n_children(); ++c)
6301 this->child(c)->set_all_manifold_ids(manifold_ind);
6302
6303 switch (structdim)
6304 {
6305 case 1:
6306 if (dim == 1)
6307 {
6308 (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(0)] =
6309 manifold_ind;
6310 (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(1)] =
6311 manifold_ind;
6312 }
6313 break;
6314
6315 case 2:
6316 // for quads/simplices also set manifold_id of bounding lines
6317 for (unsigned int i = 0; i < this->n_lines(); ++i)
6318 this->line(i)->set_manifold_id(manifold_ind);
6319 break;
6320 default:
6322 }
6323}
6324
6325
6326
6327template <int structdim, int dim, int spacedim>
6328double
6330{
6331 boost::container::small_vector<Point<spacedim>,
6332# ifndef _MSC_VER
6333 ReferenceCells::max_n_vertices<structdim>()
6334# else
6336# endif
6337 >
6338 vertices(this->n_vertices());
6339
6340 for (unsigned int v = 0; v < vertices.size(); ++v)
6341 vertices[v] = this->vertex(v);
6342
6343 return internal::TriaAccessorImplementation::diameter<structdim, spacedim>(
6344 vertices);
6345}
6346
6347
6348
6349template <int dim, int spacedim>
6350double
6352 const Mapping<dim, spacedim> &mapping) const
6353{
6354 return internal::TriaAccessorImplementation::diameter<dim, spacedim>(
6356 this->tria, this->level(), this->index())));
6357}
6358
6359
6360
6361template <int structdim, int dim, int spacedim>
6362std::pair<Point<spacedim>, double>
6364{
6365 // If the object is one dimensional,
6366 // the enclosing ball is the initial iterate
6367 // i.e., the ball's center and diameter are
6368 // the center and the diameter of the object.
6369 if (structdim == 1)
6370 return std::make_pair((this->vertex(1) + this->vertex(0)) * 0.5,
6371 (this->vertex(1) - this->vertex(0)).norm() * 0.5);
6372
6373 // The list is_initial_guess_vertex contains bool values and has
6374 // the same size as the number of vertices per object.
6375 // The entries of is_initial_guess_vertex are set true only for those
6376 // two vertices corresponding to the largest diagonal which is being used
6377 // to construct the initial ball.
6378 // We employ this mask to skip these two vertices while enlarging the ball.
6379 std::vector<bool> is_initial_guess_vertex(this->n_vertices());
6380
6381 // First let all the vertices be outside
6382 std::fill(is_initial_guess_vertex.begin(),
6383 is_initial_guess_vertex.end(),
6384 false);
6385
6386 // Get an initial guess by looking at the largest diagonal
6387 Point<spacedim> center;
6388 double radius = 0;
6389
6390 switch (structdim)
6391 {
6392 case 2:
6393 {
6394 const Point<spacedim> p30(this->vertex(3) - this->vertex(0));
6395 const Point<spacedim> p21(this->vertex(2) - this->vertex(1));
6396 if (p30.norm() > p21.norm())
6397 {
6398 center = this->vertex(0) + 0.5 * p30;
6399 radius = p30.norm() / 2.;
6400 is_initial_guess_vertex[3] = true;
6401 is_initial_guess_vertex[0] = true;
6402 }
6403 else
6404 {
6405 center = this->vertex(1) + 0.5 * p21;
6406 radius = p21.norm() / 2.;
6407 is_initial_guess_vertex[2] = true;
6408 is_initial_guess_vertex[1] = true;
6409 }
6410 break;
6411 }
6412 case 3:
6413 {
6414 const Point<spacedim> p70(this->vertex(7) - this->vertex(0));
6415 const Point<spacedim> p61(this->vertex(6) - this->vertex(1));
6416 const Point<spacedim> p25(this->vertex(2) - this->vertex(5));
6417 const Point<spacedim> p34(this->vertex(3) - this->vertex(4));
6418 const std::vector<double> diagonals = {p70.norm(),
6419 p61.norm(),
6420 p25.norm(),
6421 p34.norm()};
6422 const std::vector<double>::const_iterator it =
6423 std::max_element(diagonals.begin(), diagonals.end());
6424 if (it == diagonals.begin())
6425 {
6426 center = this->vertex(0) + 0.5 * p70;
6427 is_initial_guess_vertex[7] = true;
6428 is_initial_guess_vertex[0] = true;
6429 }
6430 else if (it == diagonals.begin() + 1)
6431 {
6432 center = this->vertex(1) + 0.5 * p61;
6433 is_initial_guess_vertex[6] = true;
6434 is_initial_guess_vertex[1] = true;
6435 }
6436 else if (it == diagonals.begin() + 2)
6437 {
6438 center = this->vertex(5) + 0.5 * p25;
6439 is_initial_guess_vertex[2] = true;
6440 is_initial_guess_vertex[5] = true;
6441 }
6442 else
6443 {
6444 center = this->vertex(4) + 0.5 * p34;
6445 is_initial_guess_vertex[3] = true;
6446 is_initial_guess_vertex[4] = true;
6447 }
6448 radius = *it * 0.5;
6449 break;
6450 }
6451 default:
6453 return std::pair<Point<spacedim>, double>();
6454 }
6455
6456 // For each vertex that is found to be geometrically outside the ball
6457 // enlarge the ball so that the new ball contains both the previous ball
6458 // and the given vertex.
6459 for (const unsigned int v : this->vertex_indices())
6460 if (!is_initial_guess_vertex[v])
6461 {
6462 const double distance = center.distance(this->vertex(v));
6463 if (distance > radius)
6464 {
6465 // we found a vertex which is outside of the ball
6466 // extend it (move center and change radius)
6467 const Point<spacedim> pCV(center - this->vertex(v));
6468 radius = (distance + radius) * 0.5;
6469 center = this->vertex(v) + pCV * (radius / distance);
6470
6471 // Now the new ball constructed in this block
6472 // encloses the vertex (v) that was found to be geometrically
6473 // outside the old ball.
6474 }
6475 }
6476# ifdef DEBUG
6477 bool all_vertices_within_ball = true;
6478
6479 // Set all_vertices_within_ball false if any of the vertices of the object
6480 // are geometrically outside the ball
6481 for (const unsigned int v : this->vertex_indices())
6482 if (center.distance(this->vertex(v)) >
6483 radius + 100. * std::numeric_limits<double>::epsilon())
6484 {
6485 all_vertices_within_ball = false;
6486 break;
6487 }
6488 // If all the vertices are not within the ball throw error
6489 Assert(all_vertices_within_ball, ExcInternalError());
6490# endif
6491 return std::make_pair(center, radius);
6492}
6493
6494
6495template <int structdim, int dim, int spacedim>
6496double
6498{
6499 switch (structdim)
6500 {
6501 case 1:
6502 return (this->vertex(1) - this->vertex(0)).norm();
6503 case 2:
6504 case 3:
6505 {
6506 double min = std::numeric_limits<double>::max();
6507 for (const unsigned int i : this->vertex_indices())
6508 for (unsigned int j = i + 1; j < this->n_vertices(); ++j)
6509 min = std::min(min,
6510 (this->vertex(i) - this->vertex(j)) *
6511 (this->vertex(i) - this->vertex(j)));
6512 return std::sqrt(min);
6513 }
6514 default:
6516 return -1e10;
6517 }
6518}
6519
6520
6521template <int structdim, int dim, int spacedim>
6522bool
6525{
6526 // go through the vertices and check... The
6527 // cell is a translation of the previous
6528 // one in case the distance between the
6529 // individual vertices in the two cell is
6530 // the same for all the vertices. So do the
6531 // check by first getting the distance on
6532 // the first vertex, and then checking
6533 // whether all others have the same down to
6534 // rounding errors (we have to be careful
6535 // here because the calculation of the
6536 // displacement between one cell and the
6537 // next can already result in the loss of
6538 // one or two digits), so we choose 1e-12
6539 // times the distance between the zeroth
6540 // vertices here.
6541 bool is_translation = true;
6542 const Tensor<1, spacedim> dist = o->vertex(0) - this->vertex(0);
6543 const double tol_square = 1e-24 * dist.norm_square();
6544 for (unsigned int i = 1; i < this->n_vertices(); ++i)
6545 {
6546 const Tensor<1, spacedim> dist_new =
6547 (o->vertex(i) - this->vertex(i)) - dist;
6548 if (dist_new.norm_square() > tol_square)
6549 {
6550 is_translation = false;
6551 break;
6552 }
6553 }
6554 return is_translation;
6555}
6556
6557
6558
6559template <int structdim, int dim, int spacedim>
6560unsigned int
6562{
6563 return this->reference_cell().n_vertices();
6564}
6565
6566
6567
6568template <int structdim, int dim, int spacedim>
6569unsigned int
6571{
6572 return this->reference_cell().n_lines();
6573}
6574
6575
6576
6577template <int structdim, int dim, int spacedim>
6578unsigned int
6580{
6581 Assert(structdim == dim,
6582 ExcMessage("This function can only be used on objects "
6583 "that are cells, but not on faces or edges "
6584 "that bound cells."));
6585
6586 return this->reference_cell().n_faces();
6587}
6588
6589
6590
6591template <int structdim, int dim, int spacedim>
6594{
6595 return {0U, n_vertices()};
6596}
6597
6598
6599
6600template <int structdim, int dim, int spacedim>
6603{
6604 return {0U, n_lines()};
6605}
6606
6607
6608
6609template <int structdim, int dim, int spacedim>
6612{
6613 return {0U, n_faces()};
6614}
6615
6616
6617
6618/*----------------- Functions: TriaAccessor<0,dim,spacedim> -----------------*/
6619
6620template <int dim, int spacedim>
6622 const Triangulation<dim, spacedim> *tria,
6623 const unsigned int vertex_index)
6624 : tria(tria)
6625 , global_vertex_index(vertex_index)
6626{}
6627
6628
6629
6630template <int dim, int spacedim>
6632 const Triangulation<dim, spacedim> *tria,
6633 const int /*level*/,
6634 const int index,
6635 const AccessorData *)
6636 : tria(tria)
6638{}
6639
6640
6641
6642template <int dim, int spacedim>
6643template <int structdim2, int dim2, int spacedim2>
6646 : tria(nullptr)
6648{
6649 Assert(false, ExcImpossibleInDim(0));
6650}
6651
6652
6653
6654template <int dim, int spacedim>
6655template <int structdim2, int dim2, int spacedim2>
6658 : tria(nullptr)
6660{
6661 Assert(false, ExcImpossibleInDim(0));
6662}
6663
6664
6665
6666template <int dim, int spacedim>
6667inline void
6669{
6670 tria = t.tria;
6671 global_vertex_index = t.global_vertex_index;
6672}
6673
6674
6675
6676template <int dim, int spacedim>
6677inline bool
6679 const TriaAccessor<0, dim, spacedim> &other) const
6680{
6682
6683 return (global_vertex_index < other.global_vertex_index);
6684}
6685
6686
6687
6688template <int dim, int spacedim>
6691{
6692 if (global_vertex_index != numbers::invalid_unsigned_int)
6693 return IteratorState::valid;
6694 else
6696}
6697
6698
6699
6700template <int dim, int spacedim>
6701inline int
6703{
6704 return 0;
6705}
6706
6707
6708
6709template <int dim, int spacedim>
6710inline int
6712{
6713 return global_vertex_index;
6714}
6715
6716
6717
6718template <int dim, int spacedim>
6719inline const Triangulation<dim, spacedim> &
6721{
6722 return *tria;
6723}
6724
6725
6726
6727template <int dim, int spacedim>
6728inline void
6730{
6732 if (global_vertex_index >= tria->n_vertices())
6734}
6735
6736
6737
6738template <int dim, int spacedim>
6739inline void
6741{
6742 if (global_vertex_index != numbers::invalid_unsigned_int)
6743 {
6744 if (global_vertex_index != 0)
6746 else
6748 }
6749}
6750
6751
6752
6753template <int dim, int spacedim>
6754inline bool
6756{
6757 const bool result =
6758 ((tria == t.tria) && (global_vertex_index == t.global_vertex_index));
6759
6760 return result;
6761}
6762
6763
6764
6765template <int dim, int spacedim>
6766inline bool
6768{
6769 return !(*this == t);
6770}
6771
6772
6773
6774template <int dim, int spacedim>
6775inline unsigned int
6776TriaAccessor<0, dim, spacedim>::vertex_index(const unsigned int) const
6777{
6778 return global_vertex_index;
6779}
6780
6781
6782
6783template <int dim, int spacedim>
6784inline Point<spacedim> &
6785TriaAccessor<0, dim, spacedim>::vertex(const unsigned int) const
6786{
6787 return const_cast<Point<spacedim> &>(
6788 this->tria->vertices[global_vertex_index]);
6789}
6790
6791
6792
6793template <int dim, int spacedim>
6794inline typename ::internal::TriangulationImplementation::
6795 Iterators<dim, spacedim>::line_iterator
6796 TriaAccessor<0, dim, spacedim>::line(const unsigned int)
6797{
6798 return typename ::internal::TriangulationImplementation::
6799 Iterators<dim, spacedim>::line_iterator();
6800}
6801
6802
6803
6804template <int dim, int spacedim>
6805inline unsigned int
6807{
6808 Assert(false, ExcImpossibleInDim(0));
6810}
6811
6812
6813
6814template <int dim, int spacedim>
6815inline typename ::internal::TriangulationImplementation::
6816 Iterators<dim, spacedim>::quad_iterator
6817 TriaAccessor<0, dim, spacedim>::quad(const unsigned int)
6818{
6819 return typename ::internal::TriangulationImplementation::
6820 Iterators<dim, spacedim>::quad_iterator();
6821}
6822
6823
6824
6825template <int dim, int spacedim>
6826inline unsigned int
6828{
6829 Assert(false, ExcImpossibleInDim(0));
6831}
6832
6833
6834
6835template <int dim, int spacedim>
6836inline double
6838{
6839 return 0.;
6840}
6841
6842
6843
6844template <int dim, int spacedim>
6845inline double
6847{
6848 return 0.;
6849}
6850
6851
6852
6853template <int dim, int spacedim>
6854inline Point<spacedim>
6855TriaAccessor<0, dim, spacedim>::center(const bool, const bool) const
6856{
6857 return this->tria->vertices[global_vertex_index];
6858}
6859
6860
6861
6862template <int dim, int spacedim>
6863inline double
6865{
6866 return 0.;
6867}
6868
6869
6870
6871template <int dim, int spacedim>
6874 const unsigned int /*face*/)
6875{
6877}
6878
6879
6880
6881template <int dim, int spacedim>
6882inline bool
6883TriaAccessor<0, dim, spacedim>::face_orientation(const unsigned int /*face*/)
6884{
6885 return false;
6886}
6887
6888
6889
6890template <int dim, int spacedim>
6891inline bool
6892TriaAccessor<0, dim, spacedim>::face_flip(const unsigned int /*face*/)
6893{
6894 return false;
6895}
6896
6897
6898
6899template <int dim, int spacedim>
6900inline bool
6901TriaAccessor<0, dim, spacedim>::face_rotation(const unsigned int /*face*/)
6902{
6903 return false;
6904}
6905
6906
6907
6908template <int dim, int spacedim>
6910TriaAccessor<0, dim, spacedim>::line_orientation(const unsigned int /*line*/)
6911{
6913}
6914
6915
6916
6917template <int dim, int spacedim>
6918inline bool
6920{
6921 return false;
6922}
6923
6924
6925
6926template <int dim, int spacedim>
6927inline unsigned int
6929{
6930 return 0;
6931}
6932
6933
6934
6935template <int dim, int spacedim>
6936inline unsigned int
6938{
6939 return 0;
6940}
6941
6942
6943
6944template <int dim, int spacedim>
6945inline unsigned int
6947{
6948 return 0;
6949}
6950
6951
6952
6953template <int dim, int spacedim>
6954inline unsigned int
6957{
6959}
6960
6961
6962
6963template <int dim, int spacedim>
6965TriaAccessor<0, dim, spacedim>::child(const unsigned int)
6966{
6968}
6969
6970
6971
6972template <int dim, int spacedim>
6975{
6977}
6978
6979
6980
6981template <int dim, int spacedim>
6982inline RefinementCase<0>
6984{
6986}
6987
6988
6989
6990template <int dim, int spacedim>
6991inline int
6993{
6994 return -1;
6995}
6996
6997
6998
6999template <int dim, int spacedim>
7000inline int
7002{
7003 return -1;
7004}
7005
7006
7007
7008template <int dim, int spacedim>
7009inline bool
7011{
7012 return tria->vertex_used(global_vertex_index);
7013}
7014
7015
7016
7017/*------------------- Functions: TriaAccessor<0,1,spacedim> -----------------*/
7018
7019template <int spacedim>
7021 const Triangulation<1, spacedim> *tria,
7022 const VertexKind vertex_kind,
7023 const unsigned int vertex_index)
7024 : tria(tria)
7025 , vertex_kind(vertex_kind)
7026 , global_vertex_index(vertex_index)
7027{}
7028
7029
7030
7031template <int spacedim>
7033 const Triangulation<1, spacedim> *tria,
7034 const int level,
7035 const int index,
7036 const AccessorData *)
7037 : tria(tria)
7038 , vertex_kind(interior_vertex)
7040{
7041 // in general, calling this constructor should yield an error -- users should
7042 // instead call the one immediately above. however, if you create something
7043 // like Triangulation<1>::face_iterator() then this calls the default
7044 // constructor of the iterator which calls the accessor with argument list
7045 // (0,-2,-2,0), so in this particular case accept this call and create an
7046 // object that corresponds to the default constructed (invalid) vertex
7047 // accessor
7048 (void)level;
7049 (void)index;
7050 Assert((level == -2) && (index == -2),
7051 ExcMessage(
7052 "This constructor can not be called for face iterators in 1d, "
7053 "except to default-construct iterator objects."));
7054}
7055
7056
7057
7058template <int spacedim>
7059template <int structdim2, int dim2, int spacedim2>
7062 : tria(nullptr)
7063 , vertex_kind(interior_vertex)
7065{
7066 Assert(false, ExcImpossibleInDim(0));
7067}
7068
7069
7070
7071template <int spacedim>
7072template <int structdim2, int dim2, int spacedim2>
7075 : tria(nullptr)
7076 , vertex_kind(interior_vertex)
7078{
7079 Assert(false, ExcImpossibleInDim(0));
7080}
7081
7082
7083
7084template <int spacedim>
7085inline void
7087{
7088 tria = t.tria;
7089 vertex_kind = t.vertex_kind;
7090 global_vertex_index = t.global_vertex_index;
7091}
7092
7093
7094
7095template <int spacedim>
7096inline void
7099{
7100 // We cannot convert from TriaAccessorBase to
7101 // TriaAccessor<0,1,spacedim> because the latter is not derived from
7102 // the former. We should never get here.
7104}
7105
7106
7107
7108template <int spacedim>
7109inline bool
7111 const TriaAccessor<0, 1, spacedim> &other) const
7112{
7114
7115 return (global_vertex_index < other.global_vertex_index);
7116}
7117
7118
7119
7120template <int spacedim>
7123{
7124 return IteratorState::valid;
7125}
7126
7127
7128template <int spacedim>
7129inline int
7131{
7132 return 0;
7133}
7134
7135
7136
7137template <int spacedim>
7138inline int
7140{
7141 return global_vertex_index;
7142}
7143
7144
7145
7146template <int spacedim>
7147inline const Triangulation<1, spacedim> &
7149{
7150 return *tria;
7151}
7152
7153
7154
7155template <int spacedim>
7156inline void
7158{
7160}
7161
7162
7163template <int spacedim>
7164inline void
7166{
7168}
7169
7170
7171
7172template <int spacedim>
7173inline bool
7175{
7176 const bool result =
7177 ((tria == t.tria) && (global_vertex_index == t.global_vertex_index));
7178 // if we point to the same vertex,
7179 // make sure we know the same about
7180 // it
7181 if (result == true)
7182 Assert(vertex_kind == t.vertex_kind, ExcInternalError());
7183
7184 return result;
7185}
7186
7187
7188
7189template <int spacedim>
7190inline bool
7192{
7193 return !(*this == t);
7194}
7195
7196
7197
7198template <int spacedim>
7199inline unsigned int
7200TriaAccessor<0, 1, spacedim>::vertex_index(const unsigned int i) const
7201{
7202 AssertIndexRange(i, 1);
7203 (void)i;
7204 return global_vertex_index;
7205}
7206
7207
7208
7209template <int spacedim>
7210inline Point<spacedim> &
7211TriaAccessor<0, 1, spacedim>::vertex(const unsigned int i) const
7212{
7213 AssertIndexRange(i, 1);
7214 (void)i;
7215 return const_cast<Point<spacedim> &>(
7216 this->tria->vertices[global_vertex_index]);
7217}
7218
7219
7220
7221template <int spacedim>
7222inline Point<spacedim>
7224{
7225 return this->tria->vertices[global_vertex_index];
7226}
7227
7228
7229
7230template <int spacedim>
7231inline typename ::internal::TriangulationImplementation::
7232 Iterators<1, spacedim>::line_iterator
7233 TriaAccessor<0, 1, spacedim>::line(const unsigned int)
7234{
7235 return {};
7236}
7237
7238
7239template <int spacedim>
7240inline unsigned int
7242{
7243 Assert(false, ExcImpossibleInDim(0));
7245}
7246
7247
7248template <int spacedim>
7249inline typename ::internal::TriangulationImplementation::
7250 Iterators<1, spacedim>::quad_iterator
7251 TriaAccessor<0, 1, spacedim>::quad(const unsigned int)
7252{
7253 return {};
7254}
7255
7256
7257
7258template <int spacedim>
7259inline unsigned int
7261{
7262 Assert(false, ExcImpossibleInDim(0));
7264}
7265
7266
7267template <int spacedim>
7268inline bool
7270{
7271 return vertex_kind != interior_vertex;
7272}
7273
7274
7275template <int spacedim>
7276inline types::boundary_id
7278{
7279 switch (vertex_kind)
7280 {
7281 case left_vertex:
7282 case right_vertex:
7283 {
7285 this->vertex_index()) !=
7286 tria->vertex_to_boundary_id_map_1d->end(),
7288
7289 return (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()];
7290 }
7291
7292 default:
7294 }
7295}
7296
7297
7298
7299template <int spacedim>
7300inline const Manifold<1, spacedim> &
7302{
7303 return this->tria->get_manifold(this->manifold_id());
7304}
7305
7306
7307
7308template <int spacedim>
7309inline types::manifold_id
7311{
7312 if (tria->vertex_to_manifold_id_map_1d->find(this->vertex_index()) !=
7313 tria->vertex_to_manifold_id_map_1d->end())
7314 return (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()];
7315 else
7317}
7318
7319
7320template <int spacedim>
7323 const unsigned int /*face*/)
7324{
7325 return 0;
7326}
7327
7328
7329template <int spacedim>
7330inline bool
7331TriaAccessor<0, 1, spacedim>::face_orientation(const unsigned int /*face*/)
7332{
7333 return false;
7334}
7335
7336
7337
7338template <int spacedim>
7339inline bool
7340TriaAccessor<0, 1, spacedim>::face_flip(const unsigned int /*face*/)
7341{
7342 return false;
7343}
7344
7345
7346
7347template <int spacedim>
7348inline bool
7349TriaAccessor<0, 1, spacedim>::face_rotation(const unsigned int /*face*/)
7350{
7351 return false;
7352}
7353
7354
7355
7356template <int spacedim>
7358TriaAccessor<0, 1, spacedim>::line_orientation(const unsigned int /*line*/)
7359{
7361}
7362
7363
7364
7365template <int spacedim>
7366inline bool
7368{
7369 return false;
7370}
7371
7372
7373
7374template <int spacedim>
7375inline unsigned int
7377{
7378 return 0;
7379}
7380
7381
7382
7383template <int spacedim>
7384inline unsigned int
7386{
7387 return 0;
7388}
7389
7390
7391
7392template <int spacedim>
7393inline unsigned int
7395{
7396 return 0;
7397}
7398
7399
7400
7401template <int spacedim>
7402inline unsigned int
7405{
7407}
7408
7409
7410
7411template <int spacedim>
7413TriaAccessor<0, 1, spacedim>::child(const unsigned int)
7414{
7416}
7417
7418
7419template <int spacedim>
7422{
7424}
7425
7426
7427template <int spacedim>
7428inline RefinementCase<0>
7430{
7432}
7433
7434template <int spacedim>
7435inline int
7437{
7438 return -1;
7439}
7440
7441
7442template <int spacedim>
7443inline int
7445{
7446 return -1;
7447}
7448
7449
7450
7451template <int spacedim>
7452inline void
7454{
7455 Assert(tria->vertex_to_boundary_id_map_1d->find(this->vertex_index()) !=
7456 tria->vertex_to_boundary_id_map_1d->end(),
7457 ExcMessage("You can't set the boundary_id of a face of a cell that is "
7458 "not actually at the boundary."));
7459
7460 (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()] = b;
7461}
7462
7463
7464
7465template <int spacedim>
7466inline void
7468{
7469 (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()] = b;
7470}
7471
7472
7473
7474template <int spacedim>
7475inline void
7477 const types::boundary_id b) const
7478{
7479 set_boundary_id(b);
7480}
7481
7482
7483
7484template <int spacedim>
7485inline void
7487{
7488 set_manifold_id(b);
7489}
7490
7491
7492
7493template <int spacedim>
7494inline bool
7496{
7497 return tria->vertex_used(global_vertex_index);
7498}
7499
7500
7501
7502template <int spacedim>
7503inline ReferenceCell
7505{
7507}
7508
7509
7510
7511template <int spacedim>
7512unsigned int
7514{
7515 return 1;
7516}
7517
7518
7519
7520template <int spacedim>
7521unsigned int
7523{
7524 return 0;
7525}
7526
7527
7528
7529template <int spacedim>
7532{
7533 return {0U, n_vertices()};
7534}
7535
7536
7537
7538template <int spacedim>
7541{
7542 return {0U, n_lines()};
7543}
7544
7545/*------------------ Functions: CellAccessor<dim,spacedim> ------------------*/
7546
7547
7548template <int dim, int spacedim>
7550 const Triangulation<dim, spacedim> *parent,
7551 const int level,
7552 const int index,
7553 const AccessorData *local_data)
7554 : TriaAccessor<dim, dim, spacedim>(parent, level, index, local_data)
7555{}
7556
7557
7558
7559template <int dim, int spacedim>
7561 const TriaAccessor<dim, dim, spacedim> &cell_accessor)
7562 : TriaAccessor<dim, dim, spacedim>(
7563 static_cast<const TriaAccessor<dim, dim, spacedim> &>(cell_accessor))
7564{}
7565
7566
7567
7568template <int dim, int spacedim>
7570CellAccessor<dim, spacedim>::child(const unsigned int i) const
7571{
7573 this->present_level + 1,
7574 this->child_index(i));
7575
7576 Assert((q.state() == IteratorState::past_the_end) || q->used(),
7578
7579 return q;
7580}
7581
7582
7583
7584template <int dim, int spacedim>
7585inline boost::container::small_vector<TriaIterator<CellAccessor<dim, spacedim>>,
7588{
7589 boost::container::small_vector<TriaIterator<CellAccessor<dim, spacedim>>,
7591 child_iterators(this->n_children());
7592
7593 for (unsigned int i = 0; i < this->n_children(); ++i)
7594 child_iterators[i] = this->child(i);
7595
7596 return child_iterators;
7597}
7598
7599
7600
7601template <int dim, int spacedim>
7602inline TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
7603CellAccessor<dim, spacedim>::face(const unsigned int i) const
7604{
7605 AssertIndexRange(i, this->n_faces());
7606 if constexpr (dim == 1)
7607 {
7608 using VertexKind = typename TriaAccessor<0, 1, spacedim>::VertexKind;
7609 VertexKind vertex_kind = VertexKind::interior_vertex;
7610 if (i == 0 && at_boundary(0))
7611 vertex_kind = VertexKind::left_vertex;
7612 if (i == 1 && at_boundary(1))
7613 vertex_kind = VertexKind::right_vertex;
7614 TriaAccessor<0, 1, spacedim> a(&this->get_triangulation(),
7615 vertex_kind,
7616 this->vertex_index(i));
7618 }
7619 else if constexpr (dim == 2)
7620 return this->line(i);
7621 else if constexpr (dim == 3)
7622 return this->quad(i);
7623 else
7624 {
7625 Assert(false, ExcNotImplemented());
7626 return {};
7627 }
7628}
7629
7630
7631
7632template <int dim, int spacedim>
7633inline unsigned int
7636{
7637 for (const unsigned int face_n : this->face_indices())
7638 if (this->face(face_n) == face)
7639 return face_n;
7640
7641 Assert(false,
7642 ExcMessage("The given face is not a face of the current cell."));
7644}
7645
7646
7647
7648template <int dim, int spacedim>
7649inline boost::container::small_vector<
7650 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
7651# ifndef _MSC_VER
7652 ReferenceCells::max_n_faces<dim>()
7653# else
7655# endif
7656 >
7658{
7659 boost::container::small_vector<
7660 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
7661# ifndef _MSC_VER
7662 ReferenceCells::max_n_faces<dim>()
7663# else
7665# endif
7666 >
7667 face_iterators(this->n_faces());
7668
7669 for (const unsigned int i : this->face_indices())
7670 face_iterators[i] = this->face(i);
7671
7672 return face_iterators;
7673}
7674
7675
7676
7677template <int dim, int spacedim>
7678inline unsigned int
7679CellAccessor<dim, spacedim>::face_index(const unsigned int i) const
7680{
7681 switch (dim)
7682 {
7683 case 1:
7684 return this->vertex_index(i);
7685
7686 case 2:
7687 return this->line_index(i);
7688
7689 case 3:
7690 return this->quad_index(i);
7691
7692 default:
7694 }
7695}
7696
7697
7698
7699template <int dim, int spacedim>
7700inline int
7701CellAccessor<dim, spacedim>::neighbor_index(const unsigned int face_no) const
7702{
7703 AssertIndexRange(face_no, this->n_faces());
7704 return this->tria->levels[this->present_level]
7705 ->neighbors[this->present_index * ReferenceCells::max_n_faces<dim>() +
7706 face_no]
7707 .second;
7708}
7709
7710
7711
7712template <int dim, int spacedim>
7713inline int
7714CellAccessor<dim, spacedim>::neighbor_level(const unsigned int face_no) const
7715{
7716 AssertIndexRange(face_no, this->n_faces());
7717 return this->tria->levels[this->present_level]
7718 ->neighbors[this->present_index * ReferenceCells::max_n_faces<dim>() +
7719 face_no]
7720 .first;
7721}
7722
7723
7724
7725template <int dim, int spacedim>
7728{
7730 // cells flagged for refinement must be active
7731 // (the @p set_refine_flag function checks this,
7732 // but activity may change when refinement is
7733 // executed and for some reason the refine
7734 // flag is not cleared).
7735 Assert(this->is_active() || !this->tria->levels[this->present_level]
7736 ->refine_flags[this->present_index],
7737 ExcRefineCellNotActive());
7738 return RefinementCase<dim>(
7739 this->tria->levels[this->present_level]->refine_flags[this->present_index]);
7740}
7741
7742
7743
7744template <int dim, int spacedim>
7745inline void
7747 const RefinementCase<dim> refinement_case) const
7748{
7749 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
7750 Assert(!coarsen_flag_set(), ExcCellFlaggedForCoarsening());
7751
7752 this->tria->levels[this->present_level]->refine_flags[this->present_index] =
7753 refinement_case;
7754}
7755
7756
7757
7758template <int dim, int spacedim>
7759inline void
7761{
7762 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
7763 this->tria->levels[this->present_level]->refine_flags[this->present_index] =
7765}
7766
7767
7768template <int dim, int spacedim>
7769inline std::uint8_t
7771{
7773 if (this->tria->levels[this->present_level]->refine_choice.size() == 0)
7774 return 0U;
7775 return this->tria->levels[this->present_level]
7776 ->refine_choice[this->present_index];
7777}
7778
7779
7780template <int dim, int spacedim>
7781inline void
7783 const std::uint8_t refinement_choice) const
7784{
7785 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
7786 if (this->tria->levels[this->present_level]->refine_choice.size() != 0)
7787 this->tria->levels[this->present_level]
7788 ->refine_choice[this->present_index] = refinement_choice;
7789}
7790
7791
7792template <int dim, int spacedim>
7793inline void
7795{
7796 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
7797 if (this->tria->levels[this->present_level]->refine_choice.size() != 0)
7798 this->tria->levels[this->present_level]
7799 ->refine_choice[this->present_index] =
7801}
7802
7803
7804template <int dim, int spacedim>
7805inline bool
7807 const unsigned int face_no,
7808 const RefinementCase<dim - 1> &face_refinement_case) const
7809{
7810 Assert(dim > 1, ExcImpossibleInDim(dim));
7811 AssertIndexRange(face_no, this->n_faces());
7812 AssertIndexRange(face_refinement_case,
7814
7815 // the new refinement case is a combination
7816 // of the minimum required one for the given
7817 // face refinement and the already existing
7818 // flagged refinement case
7819 RefinementCase<dim> old_ref_case = refine_flag_set();
7820 RefinementCase<dim> new_ref_case =
7821 (old_ref_case |
7823 face_refinement_case,
7824 face_no,
7825 this->face_orientation(face_no),
7826 this->face_flip(face_no),
7827 this->face_rotation(face_no)));
7828 set_refine_flag(new_ref_case);
7829 // return, whether we had to change the
7830 // refinement flag
7831 return new_ref_case != old_ref_case;
7832}
7833
7834
7835
7836template <int dim, int spacedim>
7837inline bool
7839 const unsigned int line_no) const
7840{
7841 Assert(dim > 1, ExcImpossibleInDim(dim));
7842 AssertIndexRange(line_no, this->n_lines());
7843
7844 // the new refinement case is a combination
7845 // of the minimum required one for the given
7846 // line refinement and the already existing
7847 // flagged refinement case
7849 old_ref_case = refine_flag_set(),
7850 new_ref_case =
7851 old_ref_case |
7853 set_refine_flag(new_ref_case);
7854 // return, whether we had to change the
7855 // refinement flag
7856 return new_ref_case != old_ref_case;
7857}
7858
7859
7860
7861template <>
7862inline ::internal::SubfaceCase<1>
7863CellAccessor<1>::subface_case(const unsigned int) const
7864{
7865 return ::internal::SubfaceCase<1>::case_none;
7866}
7867
7868template <>
7869inline ::internal::SubfaceCase<1>
7870CellAccessor<1, 2>::subface_case(const unsigned int) const
7871{
7872 return ::internal::SubfaceCase<1>::case_none;
7873}
7874
7875
7876template <>
7877inline ::internal::SubfaceCase<1>
7878CellAccessor<1, 3>::subface_case(const unsigned int) const
7879{
7880 return ::internal::SubfaceCase<1>::case_none;
7881}
7882
7883
7884template <>
7885inline ::internal::SubfaceCase<2>
7886CellAccessor<2>::subface_case(const unsigned int face_no) const
7887{
7889 AssertIndexRange(face_no, this->n_faces());
7890 return ((face(face_no)->has_children()) ?
7892 ::internal::SubfaceCase<2>::case_none);
7893}
7894
7895template <>
7896inline ::internal::SubfaceCase<2>
7897CellAccessor<2, 3>::subface_case(const unsigned int face_no) const
7898{
7900 AssertIndexRange(face_no, this->n_faces());
7901 return ((face(face_no)->has_children()) ?
7903 ::internal::SubfaceCase<2>::case_none);
7904}
7905
7906
7907template <>
7908inline ::internal::SubfaceCase<3>
7909CellAccessor<3>::subface_case(const unsigned int face_no) const
7910{
7912 AssertIndexRange(face_no, this->n_faces());
7913 switch (static_cast<std::uint8_t>(face(face_no)->refinement_case()))
7914 {
7916 return ::internal::SubfaceCase<3>::case_none;
7918 if (face(face_no)->child(0)->has_children())
7919 {
7920 Assert(face(face_no)->child(0)->refinement_case() ==
7923 if (face(face_no)->child(1)->has_children())
7924 {
7925 Assert(face(face_no)->child(1)->refinement_case() ==
7928 return ::internal::SubfaceCase<3>::case_x1y2y;
7929 }
7930 else
7931 return ::internal::SubfaceCase<3>::case_x1y;
7932 }
7933 else
7934 {
7935 if (face(face_no)->child(1)->has_children())
7936 {
7937 Assert(face(face_no)->child(1)->refinement_case() ==
7940 return ::internal::SubfaceCase<3>::case_x2y;
7941 }
7942 else
7943 return ::internal::SubfaceCase<3>::case_x;
7944 }
7946 if (face(face_no)->child(0)->has_children())
7947 {
7948 Assert(face(face_no)->child(0)->refinement_case() ==
7951 if (face(face_no)->child(1)->has_children())
7952 {
7953 Assert(face(face_no)->child(1)->refinement_case() ==
7956 return ::internal::SubfaceCase<3>::case_y1x2x;
7957 }
7958 else
7959 return ::internal::SubfaceCase<3>::case_y1x;
7960 }
7961 else
7962 {
7963 if (face(face_no)->child(1)->has_children())
7964 {
7965 Assert(face(face_no)->child(1)->refinement_case() ==
7968 return ::internal::SubfaceCase<3>::case_y2x;
7969 }
7970 else
7971 return ::internal::SubfaceCase<3>::case_y;
7972 }
7974 return ::internal::SubfaceCase<3>::case_xy;
7975 default:
7977 }
7978 // we should never get here
7979 return ::internal::SubfaceCase<3>::case_none;
7980}
7981
7982
7983
7984template <int dim, int spacedim>
7985inline bool
7987{
7989 // cells flagged for coarsening must be active
7990 // (the @p set_refine_flag function checks this,
7991 // but activity may change when refinement is
7992 // executed and for some reason the refine
7993 // flag is not cleared).
7994 Assert(this->is_active() || !this->tria->levels[this->present_level]
7995 ->coarsen_flags[this->present_index],
7996 ExcRefineCellNotActive());
7997 return this->tria->levels[this->present_level]
7998 ->coarsen_flags[this->present_index];
7999}
8000
8001
8002
8003template <int dim, int spacedim>
8004inline void
8006{
8007 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
8008 Assert(!refine_flag_set(), ExcCellFlaggedForRefinement());
8009
8010 this->tria->levels[this->present_level]->coarsen_flags[this->present_index] =
8011 true;
8012}
8013
8014
8015
8016template <int dim, int spacedim>
8017inline void
8019{
8020 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
8021 this->tria->levels[this->present_level]->coarsen_flags[this->present_index] =
8022 false;
8023}
8024
8025
8026
8027template <int dim, int spacedim>
8029CellAccessor<dim, spacedim>::neighbor(const unsigned int face_no) const
8030{
8032 neighbor_level(face_no),
8033 neighbor_index(face_no));
8034
8035 Assert((q.state() == IteratorState::past_the_end) || q->used(),
8037
8038 return q;
8039}
8040
8041
8042
8043template <int dim, int spacedim>
8044inline bool
8046{
8047 return !this->has_children();
8048}
8049
8050
8051
8052template <int dim, int spacedim>
8053inline bool
8055{
8056 Assert(this->is_active(),
8057 ExcMessage("is_locally_owned() can only be called on active cells!"));
8058# ifndef DEAL_II_WITH_MPI
8059 return true;
8060# else
8061
8062 // Serial triangulations report invalid_subdomain_id as their locally owned
8063 // subdomain, so the first condition checks whether we have a serial
8064 // triangulation, in which case all cells are locally owned. The second
8065 // condition compares the subdomain id in the parallel case.
8066 const types::subdomain_id locally_owned_subdomain =
8067 this->tria->locally_owned_subdomain();
8068 return (locally_owned_subdomain == numbers::invalid_subdomain_id ||
8069 this->subdomain_id() == locally_owned_subdomain);
8070
8071# endif
8072}
8073
8074
8075template <int dim, int spacedim>
8076inline bool
8078{
8079# ifndef DEAL_II_WITH_MPI
8080 return true;
8081# else
8082
8083 // Serial triangulations report invalid_subdomain_id as their locally owned
8084 // subdomain, so the first condition checks whether we have a serial
8085 // triangulation, in which case all cells are locally owned. The second
8086 // condition compares the subdomain id in the parallel case.
8087 const types::subdomain_id locally_owned_subdomain =
8088 this->tria->locally_owned_subdomain();
8089 return (locally_owned_subdomain == numbers::invalid_subdomain_id ||
8090 this->level_subdomain_id() == locally_owned_subdomain);
8091
8092# endif
8093}
8094
8095
8096template <int dim, int spacedim>
8097inline bool
8099{
8100 Assert(this->is_active(),
8101 ExcMessage("is_ghost() can only be called on active cells!"));
8102 if (this->has_children())
8103 return false;
8104
8105# ifndef DEAL_II_WITH_MPI
8106 return false;
8107# else
8108
8109 // Serial triangulations report invalid_subdomain_id as their locally owned
8110 // subdomain, so the first condition rules out that case as all cells to a
8111 // serial triangulation are locally owned and none is ghosted. The second
8112 // and third conditions check whether the cell's subdomain is not the
8113 // locally owned one and not artificial.
8114 const types::subdomain_id locally_owned_subdomain =
8115 this->tria->locally_owned_subdomain();
8116 const types::subdomain_id subdomain_id = this->subdomain_id();
8117 return (locally_owned_subdomain != numbers::invalid_subdomain_id &&
8118 subdomain_id != locally_owned_subdomain &&
8119 subdomain_id != numbers::artificial_subdomain_id);
8120
8121# endif
8122}
8123
8124
8125template <int dim, int spacedim>
8126inline bool
8128{
8129# ifndef DEAL_II_WITH_MPI
8130 return false;
8131# else
8132
8133 // Serial triangulations report invalid_subdomain_id as their locally owned
8134 // subdomain, so the first condition checks whether we have a serial
8135 // triangulation, in which case all cells are locally owned. The second
8136 // condition compares the subdomain id in the parallel case.
8137 const types::subdomain_id locally_owned_subdomain =
8138 this->tria->locally_owned_subdomain();
8139 const types::subdomain_id subdomain_id = this->level_subdomain_id();
8140 return (locally_owned_subdomain != numbers::invalid_subdomain_id &&
8141 subdomain_id != locally_owned_subdomain &&
8142 subdomain_id != numbers::artificial_subdomain_id);
8143
8144# endif
8145}
8146
8147
8148
8149template <int dim, int spacedim>
8150inline bool
8152{
8153 Assert(this->is_active(),
8154 ExcMessage("is_artificial() can only be called on active cells!"));
8155# ifndef DEAL_II_WITH_MPI
8156 return false;
8157# else
8158
8159 // Serial triangulations report invalid_subdomain_id as their locally owned
8160 // subdomain, so the first condition rules out that case as all cells to a
8161 // serial triangulation are locally owned and none is artificial.
8162 return (this->tria->locally_owned_subdomain() !=
8164 this->subdomain_id() == numbers::artificial_subdomain_id);
8165
8166# endif
8167}
8168
8169
8170
8171template <int dim, int spacedim>
8172inline bool
8174{
8175# ifndef DEAL_II_WITH_MPI
8176 return false;
8177# else
8178 return (this->tria->locally_owned_subdomain() !=
8180 this->level_subdomain_id() == numbers::artificial_subdomain_id);
8181# endif
8182}
8183
8184
8185
8186template <int dim, int spacedim>
8189{
8191 Assert(this->is_active(),
8192 ExcMessage("subdomain_id() can only be called on active cells!"));
8193 return this->tria->levels[this->present_level]
8194 ->subdomain_ids[this->present_index];
8195}
8196
8197
8198
8199template <int dim, int spacedim>
8202{
8204 return this->tria->levels[this->present_level]
8205 ->level_subdomain_ids[this->present_index];
8206}
8207
8208
8209
8210template <int dim, int spacedim>
8211inline unsigned int
8212CellAccessor<dim, spacedim>::neighbor_face_no(const unsigned int neighbor) const
8213{
8214 const unsigned int n2 = neighbor_of_neighbor_internal(neighbor);
8216 // return this value as the
8217 // neighbor is not coarser
8218 return n2;
8219 else
8220 // the neighbor is coarser
8221 return neighbor_of_coarser_neighbor(neighbor).first;
8222}
8223
8224
8225
8226template <int dim, int spacedim>
8227inline bool
8229{
8230 return false;
8231}
8232
8233
8234
8235template <int dim, int spacedim>
8236inline unsigned int
8238{
8240 return this->tria->levels[this->present_level]
8241 ->active_cell_indices[this->present_index];
8242}
8243
8244
8245
8246template <int dim, int spacedim>
8249{
8251 Assert(this->is_active(),
8252 ExcMessage(
8253 "global_active_cell_index() can only be called on active cells!"));
8254
8255 return this->tria->levels[this->present_level]
8256 ->global_active_cell_indices[this->present_index];
8257}
8258
8259
8260
8261template <int dim, int spacedim>
8264{
8265 return this->tria->levels[this->present_level]
8266 ->global_level_cell_indices[this->present_index];
8267}
8268
8269#endif // DOXYGEN
8270
8272
8273
8274#endif
void recursively_set_subdomain_id(const types::subdomain_id new_subdomain_id) const
CellAccessor(const Triangulation< dim, spacedim > *parent=nullptr, const int level=-1, const int index=-1, const AccessorData *local_data=nullptr)
TriaIterator< CellAccessor< dim, spacedim > > parent() const
unsigned int neighbor_face_no(const unsigned int neighbor) const
void set_active_cell_index(const unsigned int active_cell_index) const
unsigned int neighbor_of_neighbor_internal(const unsigned int neighbor) const
TriaIterator< CellAccessor< dim, spacedim > > periodic_neighbor(const unsigned int i) const
types::global_cell_index global_active_cell_index() const
void set_neighbor(const unsigned int i, const TriaIterator< CellAccessor< dim, spacedim > > &pointer) const
TriaIterator< CellAccessor< dim, spacedim > > neighbor(const unsigned int face_no) const
boost::container::small_vector< TriaIterator< TriaAccessor< dim - 1, dim, spacedim > >, ReferenceCells::max_n_faces< dim >() > face_iterators() const
bool is_artificial_on_level() const
void set_direction_flag(const bool new_direction_flag) const
void recursively_set_material_id(const types::material_id new_material_id) const
void set_level_subdomain_id(const types::subdomain_id new_level_subdomain_id) const
types::subdomain_id level_subdomain_id() const
std::uint8_t refine_choice() const
RefinementCase< dim > refine_flag_set() const
bool flag_for_face_refinement(const unsigned int face_no, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement) const
unsigned int face_index(const unsigned int i) const
double diameter(const Mapping< dim, spacedim > &mapping) const
TriaActiveIterator< DoFCellAccessor< dim, spacedim, false > > as_dof_handler_iterator(const DoFHandler< dim, spacedim > &dof_handler) const
CellAccessor< dim, spacedim > & operator=(CellAccessor< dim, spacedim > &&)=default
bool is_active() const
bool is_ghost() const
TriaIterator< CellAccessor< dim, spacedim > > neighbor_child_on_subface(const unsigned int face_no, const unsigned int subface_no) const
void set_subdomain_id(const types::subdomain_id new_subdomain_id) const
bool neighbor_is_coarser(const unsigned int face_no) const
TriaIterator< CellAccessor< dim, spacedim > > child(const unsigned int i) const
void set_global_level_cell_index(const types::global_cell_index index) const
bool has_periodic_neighbor(const unsigned int i) const
void clear_refine_choice() const
int periodic_neighbor_level(const unsigned int i) const
std::pair< unsigned int, unsigned int > neighbor_of_coarser_neighbor(const unsigned int neighbor) const
~CellAccessor()=default
CellAccessor(const CellAccessor< dim, spacedim > &)=default
void set_coarsen_flag() const
TriaIterator< TriaAccessor< dim - 1, dim, spacedim > > face(const unsigned int i) const
CellAccessor< dim, spacedim > & operator=(const CellAccessor< dim, spacedim > &)=delete
unsigned int neighbor_of_neighbor(const unsigned int face_no) const
unsigned int face_iterator_to_index(const TriaIterator< TriaAccessor< dim - 1, dim, spacedim > > &face) const
void set_material_id(const types::material_id new_material_id) const
bool is_locally_owned() const
void set_refine_flag(const RefinementCase< dim > ref_case=RefinementCase< dim >::isotropic_refinement) const
CellAccessor(CellAccessor< dim, spacedim > &&)=default
bool point_inside_codim(const Point< spacedim_ > &p) const
typename TriaAccessor< dim, dim, spacedim >::AccessorData AccessorData
bool is_locally_owned_on_level() const
bool has_boundary_lines() const
TriaIterator< CellAccessor< dim, spacedim > > periodic_neighbor_child_on_subface(const unsigned int face_no, const unsigned int subface_no) const
int neighbor_level(const unsigned int face_no) const
int periodic_neighbor_index(const unsigned int i) const
bool periodic_neighbor_is_coarser(const unsigned int i) const
void set_global_active_cell_index(const types::global_cell_index index) const
void clear_coarsen_flag() const
void set_refine_choice(const std::uint8_t refinement_choice=static_cast< char >(IsotropicRefinementChoice::isotropic_refinement)) const
boost::container::small_vector< TriaIterator< CellAccessor< dim, spacedim > >, GeometryInfo< dim >::max_children_per_cell > child_iterators() const
void set_parent(const unsigned int parent_index)
std::pair< unsigned int, unsigned int > periodic_neighbor_of_coarser_periodic_neighbor(const unsigned face_no) const
CellAccessor(const TriaAccessor< dim, dim, spacedim > &cell_accessor)
bool at_boundary() const
unsigned int active_cell_index() const
static bool is_level_cell()
void clear_refine_flag() const
int neighbor_index(const unsigned int face_no) const
bool point_inside(const Point< spacedim > &p) const
types::subdomain_id subdomain_id() const
bool direction_flag() const
types::material_id material_id() const
bool coarsen_flag_set() const
types::global_cell_index global_level_cell_index() const
bool flag_for_line_refinement(const unsigned int line_no) const
CellId id() const
::internal::SubfaceCase< dim > subface_case(const unsigned int face_no) const
bool is_artificial() const
TriaIterator< DoFCellAccessor< dim, spacedim, true > > as_dof_handler_level_iterator(const DoFHandler< dim, spacedim > &dof_handler) const
bool is_ghost_on_level() const
TriaIterator< CellAccessor< dim, spacedim > > neighbor_or_periodic_neighbor(const unsigned int i) const
int parent_index() const
unsigned int periodic_neighbor_of_periodic_neighbor(const unsigned int i) const
unsigned int periodic_neighbor_face_no(const unsigned int i) const
void set_manifold_id(const types::manifold_id) const
InvalidAccessor(const InvalidAccessor &)
static constexpr unsigned int space_dimension
static int level()
bool operator==(const InvalidAccessor &) const
Point< spacedim > & vertex(const unsigned int i) const
void * quad(const unsigned int i) const
void operator++() const
void * line(const unsigned int i) const
static IteratorState::IteratorStates state()
bool operator!=(const InvalidAccessor &) const
bool has_children() const
unsigned int user_index() const
static const unsigned int structure_dimension
void copy_from(const InvalidAccessor &)
void set_user_index(const unsigned int p) const
static int index()
types::manifold_id manifold_id() const
InvalidAccessor(const void *parent=nullptr, const int level=-1, const int index=-1, const AccessorData *local_data=nullptr)
static constexpr unsigned int dimension
void operator--() const
bool used() const
InvalidAccessor(const OtherAccessor &)
Abstract base class for mapping classes.
Definition mapping.h:318
virtual boost::container::small_vector< Point< spacedim >, ReferenceCells::max_n_vertices< dim >() > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition point.h:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
unsigned int n_vertices() const
unsigned int standard_to_real_face_vertex(const unsigned int vertex, const unsigned int face, const types::geometric_orientation face_orientation) const
std::array< unsigned int, 2 > standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const
std::array< unsigned int, 2 > standard_line_to_face_and_line_index(const unsigned int line) const
unsigned int n_faces() const
unsigned int n_lines() const
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
unsigned int standard_to_real_face_line(const unsigned int line, const unsigned int face, const types::geometric_orientation face_orientation) const
numbers::NumberTraits< Number >::real_type norm() const
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
bool operator!=(const TriaAccessorBase &) const
static constexpr unsigned int space_dimension
static constexpr unsigned int dimension
TriaAccessorBase(const TriaAccessorBase &)
void operator=(const TriaAccessorBase *)=delete
static const unsigned int structure_dimension
TriaAccessorBase(const Triangulation< dim, spacedim > *parent=nullptr, const int level=-1, const int index=-1, const AccessorData *=nullptr)
void copy_from(const TriaAccessorBase &)
const Triangulation< dim, spacedim > & get_triangulation() const
IteratorState::IteratorStates state() const
int index() const
bool operator<(const TriaAccessorBase &other) const
typename::internal::TriaAccessorImplementation::PresentLevelType< structdim, dim >::type present_level
::internal::TriangulationImplementation::TriaObjects & objects() const
TriaAccessorBase & operator=(const TriaAccessorBase &)
int level() const
const Triangulation< dim, spacedim > * tria
bool operator==(const TriaAccessorBase &) const
static TriaIterator< TriaAccessor< 0, 1, spacedim > > isotropic_child(const unsigned int)
Return an invalid object.
static unsigned int line_index(const unsigned int i)
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
static unsigned int child_iterator_to_index(const TriaIterator< TriaAccessor< 0, 1, spacedim > > &)
Return an invalid unsigned integer.
static unsigned int max_refinement_depth()
static unsigned int quad_index(const unsigned int i)
static typename::internal::TriangulationImplementation::Iterators< 1, spacedim >::quad_iterator quad(const unsigned int i)
const Manifold< 1, spacedim > & get_manifold() const
TriaAccessor(const InvalidAccessor< structdim2, dim2, spacedim2 > &)
TriaAccessor(const Triangulation< 1, spacedim > *tria, const VertexKind vertex_kind, const unsigned int vertex_index)
static TriaIterator< TriaAccessor< 0, 1, spacedim > > child(const unsigned int)
Return an invalid object.
unsigned int n_lines() const
const Triangulation< 1, spacedim > & get_triangulation() const
bool operator!=(const TriaAccessor &) const
TriaAccessor(const Triangulation< 1, spacedim > *tria=nullptr, const int=0, const int=0, const AccessorData *=nullptr)
static unsigned int n_active_descendants()
static int isotropic_child_index(const unsigned int i)
Returns -1.
types::boundary_id boundary_id() const
static typename::internal::TriangulationImplementation::Iterators< 1, spacedim >::line_iterator line(const unsigned int)
static RefinementCase< 0 > refinement_case()
ReferenceCell reference_cell() const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > line_indices() const
unsigned int vertex_index(const unsigned int i=0) const
void copy_from(const TriaAccessor &)
void set_all_boundary_ids(const types::boundary_id) const
static IteratorState::IteratorStates state()
Point< spacedim > center() const
static int child_index(const unsigned int i)
Returns -1.
types::manifold_id manifold_id() const
static unsigned int n_children()
TriaAccessor(const TriaAccessor< structdim2, dim2, spacedim2 > &)
void set_boundary_id(const types::boundary_id) const
<