1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #ifndef dealii_tria_accessor_templates_h
17 #define dealii_tria_accessor_templates_h
18
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/base/geometry_info.h>
23 #include <deal.II/base/template_constraints.h>
24
25 #include <deal.II/grid/tria.h>
26 #include <deal.II/grid/tria_accessor.h>
27 #include <deal.II/grid/tria_faces.h>
28 #include <deal.II/grid/tria_iterator.h>
29 #include <deal.II/grid/tria_iterator.templates.h>
30 #include <deal.II/grid/tria_levels.h>
31
32 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
33 #include <boost/container/small_vector.hpp>
34 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
35
36 #include <cmath>
37
38 DEAL_II_NAMESPACE_OPEN
39
40
41 /*--------------------- Functions: TriaAccessorBase -------------------------*/
42
43 template <int structdim, int dim, int spacedim>
TriaAccessorBase(const Triangulation<dim,spacedim> * tria,const int level,const int index,const AccessorData *)44 inline TriaAccessorBase<structdim, dim, spacedim>::TriaAccessorBase(
45 const Triangulation<dim, spacedim> *tria,
46 const int level,
47 const int index,
48 const AccessorData *)
49 : present_level((structdim == dim) ? level : 0)
50 , present_index(index)
51 , tria(tria)
52 {
53 // non-cells have no level, so a 0
54 // should have been passed, or a -1
55 // for an end-iterator, or -2 for
56 // an invalid (default constructed)
57 // iterator
58 if (structdim != dim)
59 {
60 Assert((level == 0) || (level == -1) || (level == -2),
61 ExcInternalError());
62 }
63 }
64
65
66 template <int structdim, int dim, int spacedim>
TriaAccessorBase(const TriaAccessorBase<structdim,dim,spacedim> & a)67 inline TriaAccessorBase<structdim, dim, spacedim>::TriaAccessorBase(
68 const TriaAccessorBase<structdim, dim, spacedim> &a)
69 : present_level(a.present_level)
70 , present_index(a.present_index)
71 , tria(a.tria)
72 {}
73
74
75 template <int structdim, int dim, int spacedim>
76 inline void
copy_from(const TriaAccessorBase<structdim,dim,spacedim> & a)77 TriaAccessorBase<structdim, dim, spacedim>::copy_from(
78 const TriaAccessorBase<structdim, dim, spacedim> &a)
79 {
80 present_level = a.present_level;
81 present_index = a.present_index;
82 tria = a.tria;
83
84 if (structdim != dim)
85 {
86 Assert((present_level == 0) || (present_level == -1) ||
87 (present_level == -2),
88 ExcInternalError());
89 }
90 }
91
92
93
94 template <int structdim, int dim, int spacedim>
95 inline TriaAccessorBase<structdim, dim, spacedim> &
96 TriaAccessorBase<structdim, dim, spacedim>::
97 operator=(const TriaAccessorBase<structdim, dim, spacedim> &a)
98 {
99 present_level = a.present_level;
100 present_index = a.present_index;
101 tria = a.tria;
102
103 if (structdim != dim)
104 {
105 Assert((present_level == 0) || (present_level == -1) ||
106 (present_level == -2),
107 ExcInternalError());
108 }
109 return *this;
110 }
111
112
113
114 template <int structdim, int dim, int spacedim>
115 inline bool
116 TriaAccessorBase<structdim, dim, spacedim>::
117 operator==(const TriaAccessorBase<structdim, dim, spacedim> &a) const
118 {
119 Assert(tria == a.tria || tria == nullptr || a.tria == nullptr,
120 TriaAccessorExceptions::ExcCantCompareIterators());
121 return ((tria == a.tria) && (present_level == a.present_level) &&
122 (present_index == a.present_index));
123 }
124
125
126
127 template <int structdim, int dim, int spacedim>
128 inline bool
129 TriaAccessorBase<structdim, dim, spacedim>::
130 operator!=(const TriaAccessorBase<structdim, dim, spacedim> &a) const
131 {
132 Assert(tria == a.tria || tria == nullptr || a.tria == nullptr,
133 TriaAccessorExceptions::ExcCantCompareIterators());
134 return ((tria != a.tria) || (present_level != a.present_level) ||
135 (present_index != a.present_index));
136 }
137
138
139
140 template <int structdim, int dim, int spacedim>
141 inline bool
142 TriaAccessorBase<structdim, dim, spacedim>::
143 operator<(const TriaAccessorBase<structdim, dim, spacedim> &other) const
144 {
145 Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators());
146
147 if (present_level != other.present_level)
148 return (present_level < other.present_level);
149
150 return (present_index < other.present_index);
151 }
152
153
154
155 template <int structdim, int dim, int spacedim>
156 inline int
level()157 TriaAccessorBase<structdim, dim, spacedim>::level() const
158 {
159 // This is always zero or invalid
160 // if the object is not a cell
161 return present_level;
162 }
163
164
165
166 template <int structdim, int dim, int spacedim>
167 inline int
index()168 TriaAccessorBase<structdim, dim, spacedim>::index() const
169 {
170 return present_index;
171 }
172
173
174
175 template <int structdim, int dim, int spacedim>
176 inline IteratorState::IteratorStates
state()177 TriaAccessorBase<structdim, dim, spacedim>::state() const
178 {
179 if ((present_level >= 0) && (present_index >= 0))
180 return IteratorState::valid;
181 else if (present_index == -1)
182 return IteratorState::past_the_end;
183 else
184 return IteratorState::invalid;
185 }
186
187
188
189 template <int structdim, int dim, int spacedim>
190 inline const Triangulation<dim, spacedim> &
get_triangulation()191 TriaAccessorBase<structdim, dim, spacedim>::get_triangulation() const
192 {
193 return *tria;
194 }
195
196
197
198 template <int structdim, int dim, int spacedim>
199 inline void
200 TriaAccessorBase<structdim, dim, spacedim>::operator++()
201 {
202 // this iterator is used for
203 // objects without level
204 ++this->present_index;
205
206 if (structdim != dim)
207 {
208 // is index still in the range of
209 // the vector? (note that we don't
210 // have to set the level, since
211 // dim!=1 and the object therefore
212 // has no level)
213 if (this->present_index >= static_cast<int>(objects().n_objects()))
214 this->present_index = -1;
215 }
216 else
217 {
218 while (this->present_index >=
219 static_cast<int>(
220 this->tria->levels[this->present_level]->cells.n_objects()))
221 {
222 // no -> go one level up until we find
223 // one with more than zero cells
224 ++this->present_level;
225 this->present_index = 0;
226 // highest level reached?
227 if (this->present_level >=
228 static_cast<int>(this->tria->levels.size()))
229 {
230 // return with past the end pointer
231 this->present_level = this->present_index = -1;
232 return;
233 }
234 }
235 }
236 }
237
238
239 template <int structdim, int dim, int spacedim>
240 inline void
241 TriaAccessorBase<structdim, dim, spacedim>::operator--()
242 {
243 // same as operator++
244 --this->present_index;
245
246 if (structdim != dim)
247 {
248 if (this->present_index < 0)
249 this->present_index = -1;
250 }
251 else
252 {
253 while (this->present_index < 0)
254 {
255 // no -> go one level down
256 --this->present_level;
257 // lowest level reached?
258 if (this->present_level == -1)
259 {
260 // return with past the end pointer
261 this->present_level = this->present_index = -1;
262 return;
263 }
264 // else
265 this->present_index =
266 this->tria->levels[this->present_level]->cells.n_objects() - 1;
267 }
268 }
269 }
270
271
272
273 template <int structdim, int dim, int spacedim>
274 inline dealii::internal::TriangulationImplementation::TriaObjects &
objects()275 TriaAccessorBase<structdim, dim, spacedim>::objects() const
276 {
277 if (structdim == dim)
278 return this->tria->levels[this->present_level]->cells;
279
280 if (structdim == 1 && dim > 1)
281 return this->tria->faces->lines;
282
283 if (structdim == 2 && dim > 2)
284 return this->tria->faces->quads;
285
286 Assert(false, ExcInternalError());
287
288 return this->tria->levels[this->present_level]->cells;
289 }
290
291
292
293 /*---------------------- Functions: InvalidAccessor -------------------------*/
294
295 template <int structdim, int dim, int spacedim>
InvalidAccessor(const Triangulation<dim,spacedim> *,const int,const int,const AccessorData *)296 InvalidAccessor<structdim, dim, spacedim>::InvalidAccessor(
297 const Triangulation<dim, spacedim> *,
298 const int,
299 const int,
300 const AccessorData *)
301 {
302 Assert(false,
303 ExcMessage("You are attempting an illegal conversion between "
304 "iterator/accessor types. The constructor you call "
305 "only exists to make certain template constructs "
306 "easier to write as dimension independent code but "
307 "the conversion is not valid in the current context."));
308 }
309
310
311
312 template <int structdim, int dim, int spacedim>
InvalidAccessor(const InvalidAccessor & i)313 InvalidAccessor<structdim, dim, spacedim>::InvalidAccessor(
314 const InvalidAccessor &i)
315 : TriaAccessorBase<structdim, dim, spacedim>(
316 static_cast<const TriaAccessorBase<structdim, dim, spacedim> &>(i))
317 {
318 Assert(false,
319 ExcMessage("You are attempting an illegal conversion between "
320 "iterator/accessor types. The constructor you call "
321 "only exists to make certain template constructs "
322 "easier to write as dimension independent code but "
323 "the conversion is not valid in the current context."));
324 }
325
326
327
328 template <int structdim, int dim, int spacedim>
329 void
copy_from(const InvalidAccessor &)330 InvalidAccessor<structdim, dim, spacedim>::copy_from(const InvalidAccessor &)
331 {
332 // nothing to do here. we could
333 // throw an exception but we can't
334 // get here without first creating
335 // an object which would have
336 // already thrown
337 }
338
339
340
341 template <int structdim, int dim, int spacedim>
342 bool
343 InvalidAccessor<structdim, dim, spacedim>::
344 operator==(const InvalidAccessor &) const
345 {
346 // nothing to do here. we could
347 // throw an exception but we can't
348 // get here without first creating
349 // an object which would have
350 // already thrown
351 return false;
352 }
353
354
355
356 template <int structdim, int dim, int spacedim>
357 bool
358 InvalidAccessor<structdim, dim, spacedim>::
359 operator!=(const InvalidAccessor &) const
360 {
361 // nothing to do here. we could
362 // throw an exception but we can't
363 // get here without first creating
364 // an object which would have
365 // already thrown
366 return true;
367 }
368
369
370
371 template <int structdim, int dim, int spacedim>
372 bool
used()373 InvalidAccessor<structdim, dim, spacedim>::used() const
374 {
375 // nothing to do here. we could
376 // throw an exception but we can't
377 // get here without first creating
378 // an object which would have
379 // already thrown
380 return false;
381 }
382
383
384
385 template <int structdim, int dim, int spacedim>
386 bool
has_children()387 InvalidAccessor<structdim, dim, spacedim>::has_children() const
388 {
389 // nothing to do here. we could
390 // throw an exception but we can't
391 // get here without first creating
392 // an object which would have
393 // already thrown
394 return false;
395 }
396
397
398
399 template <int structdim, int dim, int spacedim>
400 void
401 InvalidAccessor<structdim, dim, spacedim>::operator++() const
402 {}
403
404
405
406 template <int structdim, int dim, int spacedim>
407 void
408 InvalidAccessor<structdim, dim, spacedim>::operator--() const
409 {}
410
411
412
413 template <int structdim, int dim, int spacedim>
414 types::manifold_id
manifold_id()415 InvalidAccessor<structdim, dim, spacedim>::manifold_id() const
416 {
417 return numbers::flat_manifold_id;
418 }
419
420
421
422 template <int structdim, int dim, int spacedim>
423 unsigned int
user_index()424 InvalidAccessor<structdim, dim, spacedim>::user_index() const
425 {
426 return numbers::invalid_unsigned_int;
427 }
428
429
430
431 template <int structdim, int dim, int spacedim>
432 void
set_user_index(const unsigned int)433 InvalidAccessor<structdim, dim, spacedim>::set_user_index(
434 const unsigned int) const
435 {
436 Assert(false,
437 ExcMessage("You are trying to set the user index of an "
438 "invalid object."));
439 }
440
441
442
443 template <int structdim, int dim, int spacedim>
444 void
set_manifold_id(const types::manifold_id)445 InvalidAccessor<structdim, dim, spacedim>::set_manifold_id(
446 const types::manifold_id) const
447 {
448 Assert(false,
449 ExcMessage("You are trying to set the manifold id of an "
450 "invalid object."));
451 }
452
453
454
455 template <int structdim, int dim, int spacedim>
456 inline Point<spacedim> &
vertex(const unsigned int)457 InvalidAccessor<structdim, dim, spacedim>::vertex(const unsigned int) const
458 {
459 // nothing to do here. we could throw an exception but we can't get here
460 // without first creating an object which would have already thrown
461 static Point<spacedim> invalid_vertex;
462 return invalid_vertex;
463 }
464
465
466 template <int structdim, int dim, int spacedim>
467 inline typename dealii::internal::TriangulationImplementation::
468 Iterators<dim, spacedim>::line_iterator
line(const unsigned int)469 InvalidAccessor<structdim, dim, spacedim>::line(const unsigned int) const
470 {
471 // nothing to do here. we could throw an exception but we can't get here
472 // without first creating an object which would have already thrown
473 return typename dealii::internal::TriangulationImplementation::
474 Iterators<dim, spacedim>::line_iterator();
475 }
476
477
478
479 template <int structdim, int dim, int spacedim>
480 inline typename dealii::internal::TriangulationImplementation::
481 Iterators<dim, spacedim>::quad_iterator
quad(const unsigned int)482 InvalidAccessor<structdim, dim, spacedim>::quad(const unsigned int) const
483 {
484 // nothing to do here. we could throw an exception but we can't get here
485 // without first creating an object which would have already thrown
486 return dealii::internal::TriangulationImplementation::
487 Iterators<dim, spacedim>::quad_iterator();
488 }
489
490
491 /*------------------------ Functions: TriaAccessor ---------------------------*/
492
493
494 namespace internal
495 {
496 namespace TriaAccessorImplementation
497 {
498 // make sure that if in the following we
499 // write TriaAccessor
500 // we mean the *class*
501 // dealii::TriaAccessor, not the
502 // enclosing namespace
503 // dealii::internal::TriaAccessor
504 using dealii::TriaAccessor;
505
506 /**
507 * A class with the same purpose as the similarly named class of the
508 * Triangulation class. See there for more information.
509 */
510 struct Implementation
511 {
512 /**
513 * Implementation of the function of some name in the mother class.
514 */
515 template <int dim, int spacedim>
516 inline static unsigned int
line_indexImplementation517 line_index(const TriaAccessor<1, dim, spacedim> &, const unsigned int)
518 {
519 Assert(false,
520 ExcMessage("You can't ask for the index of a line bounding "
521 "a one-dimensional cell because it is not "
522 "bounded by lines."));
523 return numbers::invalid_unsigned_int;
524 }
525
526
527 template <int dim, int spacedim>
528 inline static unsigned int
line_indexImplementation529 line_index(const TriaAccessor<2, dim, spacedim> &accessor,
530 const unsigned int i)
531 {
532 return accessor.objects().get_bounding_object_indices(
533 accessor.present_index)[i];
534 }
535
536
537 inline static unsigned int
line_indexImplementation538 line_index(const TriaAccessor<3, 3, 3> &accessor, const unsigned int i)
539 {
540 const auto pair =
541 accessor.reference_cell_info().standard_line_to_face_and_line_index(
542 i);
543 const auto quad_index = pair[0];
544 const auto line_index =
545 accessor.reference_cell_info().standard_to_real_face_line(
546 pair[1], pair[0], face_orientation_raw(accessor, quad_index));
547
548 return accessor.quad(quad_index)->line_index(line_index);
549 }
550
551
552
553 /**
554 * Implementation of the function of some name in the mother class.
555 */
556 template <int structdim, int dim, int spacedim>
557 inline static unsigned int
quad_indexImplementation558 quad_index(const TriaAccessor<structdim, dim, spacedim> &,
559 const unsigned int)
560 {
561 Assert(false,
562 ExcMessage("You can't ask for the index of a quad bounding "
563 "a one- or two-dimensional cell because it is not "
564 "bounded by quads."));
565 return numbers::invalid_unsigned_int;
566 }
567
568
569 inline static unsigned int
quad_indexImplementation570 quad_index(const TriaAccessor<3, 3, 3> &accessor, const unsigned int i)
571 {
572 return accessor.tria->levels[accessor.present_level]
573 ->cells.get_bounding_object_indices(accessor.present_index)[i];
574 }
575
576
577
578 /**
579 * Implementation of the function of some name in the mother class
580 */
581 template <int structdim, int dim, int spacedim>
582 inline static bool
face_orientationImplementation583 face_orientation(const TriaAccessor<structdim, dim, spacedim> &,
584 const unsigned int)
585 {
586 /*
587 * Default implementation used in 1d and 2d
588 *
589 * In 1d, face_orientation is always true
590 */
591
592 return true;
593 }
594
595
596 template <int spacedim>
597 inline static bool
face_orientationImplementation598 face_orientation(const TriaAccessor<2, 2, spacedim> &accessor,
599 const unsigned int face)
600 {
601 return line_orientation(accessor, face);
602 }
603
604
605 inline static bool
face_orientationImplementation606 face_orientation(const TriaAccessor<3, 3, 3> &accessor,
607 const unsigned int face)
608 {
609 return ReferenceCell::internal::get_bit(
610 accessor.tria->levels[accessor.present_level]->face_orientations
611 [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
612 0 /*=orientation_bit*/);
613 }
614
615
616 inline static unsigned int
face_orientation_rawImplementation617 face_orientation_raw(const TriaAccessor<3, 3, 3> &accessor,
618 const unsigned int face)
619 {
620 AssertIndexRange(face, accessor.n_faces());
621 return accessor.tria->levels[accessor.present_level]->face_orientations
622 [accessor.present_index * GeometryInfo<3>::faces_per_cell + face];
623 }
624
625
626 /**
627 * Implementation of the function of some name in the mother class.
628 */
629 template <int structdim, int dim, int spacedim>
630 inline static bool
face_flipImplementation631 face_flip(const TriaAccessor<structdim, dim, spacedim> &,
632 const unsigned int)
633 {
634 /*
635 * Default implementation used in 1d and 2d
636 *
637 * In 1d, face_flip is always false as there is no such concept as
638 * "flipped" faces in 1d.
639 *
640 * In 2d, we currently only support meshes where all faces are in
641 * standard orientation, so the result is also false. This also
642 * matches the fact that one can *always* orient faces in 2d in such a
643 * way that the don't need to be flipped
644 */
645 return false;
646 }
647
648
649
650 inline static bool
face_flipImplementation651 face_flip(const TriaAccessor<3, 3, 3> &accessor, const unsigned int face)
652 {
653 AssertIndexRange(face, accessor.n_faces());
654 Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face <
655 accessor.tria->levels[accessor.present_level]
656 ->face_orientations.size(),
657 ExcInternalError());
658
659 return ReferenceCell::internal::get_bit(
660 accessor.tria->levels[accessor.present_level]->face_orientations
661 [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
662 2 /*=flip_bit*/);
663 }
664
665
666
667 /**
668 * Implementation of the function of some name in the mother class.
669 */
670 template <int structdim, int dim, int spacedim>
671 inline static bool
face_rotationImplementation672 face_rotation(const TriaAccessor<structdim, dim, spacedim> &,
673 const unsigned int)
674 {
675 /*
676 * Default implementation used in 1d and 2d
677 *
678 * In 1d and 2d, face_rotation is always false as there is no such
679 * concept as "rotated" faces in 1d and 2d.
680 */
681 return false;
682 }
683
684
685 inline static bool
face_rotationImplementation686 face_rotation(const TriaAccessor<3, 3, 3> &accessor,
687 const unsigned int face)
688 {
689 AssertIndexRange(face, accessor.n_faces());
690 Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face <
691 accessor.tria->levels[accessor.present_level]
692 ->face_orientations.size(),
693 ExcInternalError());
694
695 return ReferenceCell::internal::get_bit(
696 accessor.tria->levels[accessor.present_level]->face_orientations
697 [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
698 1 /*=rotation_bit*/);
699 }
700
701 /**
702 * Implementation of the function of some name in the mother class.
703 */
704 template <int dim, int spacedim>
705 inline static bool
line_orientationImplementation706 line_orientation(const TriaAccessor<1, dim, spacedim> &,
707 const unsigned int)
708 {
709 return true;
710 }
711
712
713 template <int spacedim>
714 inline static bool
line_orientationImplementation715 line_orientation(const TriaAccessor<2, 2, spacedim> &accessor,
716 const unsigned int line)
717 {
718 if (accessor.tria->levels[accessor.present_level]
719 ->face_orientations.size() == 0)
720 return true; // quads in 2d have no non-standard orientation and
721 // the array TriaLevel::face_orientations is left empty
722 else
723 return accessor.tria->levels[accessor.present_level]
724 ->face_orientations[accessor.present_index *
725 GeometryInfo<2>::faces_per_cell +
726 line];
727 }
728
729
730 template <int spacedim>
731 inline static bool
line_orientationImplementation732 line_orientation(const TriaAccessor<2, 3, spacedim> &accessor,
733 const unsigned int line)
734 {
735 Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed());
736 Assert(accessor.present_index * GeometryInfo<3>::lines_per_face + line <
737 accessor.tria->faces->quads_line_orientations.size(),
738 ExcInternalError());
739
740 // quads as part of 3d hexes can have non-standard orientation
741 return accessor.tria->faces->quads_line_orientations
742 [accessor.present_index * GeometryInfo<3>::lines_per_face + line];
743 }
744
745
746 inline static bool
line_orientationImplementation747 line_orientation(const TriaAccessor<3, 3, 3> &accessor,
748 const unsigned int line)
749 {
750 Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed());
751 AssertIndexRange(line, accessor.n_lines());
752
753 const auto pair =
754 accessor.reference_cell_info().standard_line_to_face_and_line_index(
755 line);
756 const auto quad_index = pair[0];
757 const auto line_index =
758 accessor.reference_cell_info().standard_to_real_face_line(
759 pair[1], pair[0], face_orientation_raw(accessor, quad_index));
760
761 return accessor.reference_cell_info().combine_face_and_line_orientation(
762 pair[1],
763 face_orientation_raw(accessor, quad_index),
764 accessor.quad(quad_index)->line_orientation(line_index));
765 }
766
767
768
769 /**
770 * Implementation of the function of some name in the mother class.
771 */
772 template <int structdim, int dim, int spacedim>
773 inline static void
set_face_orientationImplementation774 set_face_orientation(const TriaAccessor<structdim, dim, spacedim> &,
775 const unsigned int,
776 const bool)
777 {
778 Assert(false, ExcInternalError());
779 }
780
781
782 inline static void
set_face_orientationImplementation783 set_face_orientation(const TriaAccessor<3, 3, 3> &accessor,
784 const unsigned int face,
785 const bool value)
786 {
787 Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed());
788 AssertIndexRange(face, accessor.n_faces());
789 Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face <
790 accessor.tria->levels[accessor.present_level]
791 ->face_orientations.size(),
792 ExcInternalError());
793 ReferenceCell::internal::set_bit(
794 accessor.tria->levels[accessor.present_level]->face_orientations
795 [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
796 0 /*=orientation_bit*/,
797 value);
798 }
799
800
801
802 /**
803 * Implementation of the function of some name in the mother class.
804 */
805 template <int structdim, int dim, int spacedim>
806 inline static void
set_face_flipImplementation807 set_face_flip(const TriaAccessor<structdim, dim, spacedim> &,
808 const unsigned int,
809 const bool)
810 {
811 Assert(false, ExcInternalError());
812 }
813
814
815 inline static void
set_face_flipImplementation816 set_face_flip(const TriaAccessor<3, 3, 3> &accessor,
817 const unsigned int face,
818 const bool value)
819 {
820 AssertIndexRange(face, accessor.n_faces());
821 Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face <
822 accessor.tria->levels[accessor.present_level]
823 ->face_orientations.size(),
824 ExcInternalError());
825
826 ReferenceCell::internal::set_bit(
827 accessor.tria->levels[accessor.present_level]->face_orientations
828 [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
829 2 /*=flip_bit*/,
830 value);
831 }
832
833
834
835 /**
836 * Implementation of the function of some name in the mother class.
837 */
838 template <int structdim, int dim, int spacedim>
839 inline static void
set_face_rotationImplementation840 set_face_rotation(const TriaAccessor<structdim, dim, spacedim> &,
841 const unsigned int,
842 const bool)
843 {
844 Assert(false, ExcInternalError());
845 }
846
847
848 inline static void
set_face_rotationImplementation849 set_face_rotation(const TriaAccessor<3, 3, 3> &accessor,
850 const unsigned int face,
851 const bool value)
852 {
853 AssertIndexRange(face, accessor.n_faces());
854 Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face <
855 accessor.tria->levels[accessor.present_level]
856 ->face_orientations.size(),
857 ExcInternalError());
858
859 ReferenceCell::internal::set_bit(
860 accessor.tria->levels[accessor.present_level]->face_orientations
861 [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
862 1 /*=rotation_bit*/,
863 value);
864 }
865
866 /**
867 * Implementation of the function of some name in the mother class.
868 */
869 template <int dim, int spacedim>
870 inline static void
set_line_orientationImplementation871 set_line_orientation(const TriaAccessor<1, dim, spacedim> &,
872 const unsigned int,
873 const bool)
874 {
875 Assert(false, ExcInternalError());
876 }
877
878
879 template <int spacedim>
880 inline static void
set_line_orientationImplementation881 set_line_orientation(const TriaAccessor<2, 2, spacedim> &,
882 const unsigned int,
883 const bool)
884 {
885 // quads in 2d have no
886 // non-standard orientation
887 Assert(false, ExcInternalError());
888 }
889
890
891 template <int spacedim>
892 inline static void
set_line_orientationImplementation893 set_line_orientation(const TriaAccessor<2, 3, spacedim> &accessor,
894 const unsigned int line,
895 const bool value)
896 {
897 Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed());
898 AssertIndexRange(line, accessor.n_lines());
899 Assert(accessor.present_index * GeometryInfo<2>::lines_per_cell + line <
900 accessor.tria->faces->quads_line_orientations.size(),
901 ExcInternalError());
902
903 // quads as part of 3d hexes can have non-standard orientation
904 accessor.tria->faces->quads_line_orientations
905 [accessor.present_index * GeometryInfo<2>::lines_per_cell + line] =
906 value;
907 }
908
909
910 inline static void
set_line_orientationImplementation911 set_line_orientation(const TriaAccessor<3, 3, 3> &,
912 const unsigned int,
913 const bool)
914 {
915 // it seems like we don't need this
916 // one
917 Assert(false, ExcNotImplemented());
918 }
919
920
921 /**
922 * Implementation of the function of same name in the enclosing class.
923 */
924 template <int dim, int spacedim>
925 inline static unsigned int
vertex_indexImplementation926 vertex_index(const TriaAccessor<1, dim, spacedim> &accessor,
927 const unsigned int corner)
928 {
929 return accessor.objects().get_bounding_object_indices(
930 accessor.present_index)[corner];
931 }
932
933
934 template <int dim, int spacedim>
935 inline static unsigned int
vertex_indexImplementation936 vertex_index(const TriaAccessor<2, dim, spacedim> &accessor,
937 const unsigned int corner)
938 {
939 const auto pair = accessor.reference_cell_info()
940 .standard_vertex_to_face_and_vertex_index(corner);
941 const auto line_index = pair[0];
942 const auto vertex_index =
943 accessor.reference_cell_info().standard_to_real_face_vertex(
944 pair[1], pair[0], accessor.line_orientation(line_index));
945
946 return accessor.line(line_index)->vertex_index(vertex_index);
947 }
948
949
950
951 inline static unsigned int
vertex_indexImplementation952 vertex_index(const TriaAccessor<3, 3, 3> &accessor,
953 const unsigned int corner)
954 {
955 const auto pair = accessor.reference_cell_info()
956 .standard_vertex_to_face_and_vertex_index(corner);
957 const auto face_index = pair[0];
958 const auto vertex_index =
959 accessor.reference_cell_info().standard_to_real_face_vertex(
960 pair[1], pair[0], face_orientation_raw(accessor, face_index));
961
962 return accessor.quad(face_index)->vertex_index(vertex_index);
963 }
964 };
965 } // namespace TriaAccessorImplementation
966 } // namespace internal
967
968
969
970 template <int structdim, int dim, int spacedim>
TriaAccessor(const Triangulation<dim,spacedim> * parent,const int level,const int index,const AccessorData * local_data)971 inline TriaAccessor<structdim, dim, spacedim>::TriaAccessor(
972 const Triangulation<dim, spacedim> *parent,
973 const int level,
974 const int index,
975 const AccessorData * local_data)
976 : TriaAccessorBase<structdim, dim, spacedim>(parent, level, index, local_data)
977 {}
978
979
980
981 template <int structdim, int dim, int spacedim>
982 inline bool
used()983 TriaAccessor<structdim, dim, spacedim>::used() const
984 {
985 Assert(this->state() == IteratorState::valid,
986 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
987 *this));
988 return this->objects().used[this->present_index];
989 }
990
991
992
993 template <int structdim, int dim, int spacedim>
994 inline TriaIterator<TriaAccessor<0, dim, spacedim>>
vertex_iterator(const unsigned int i)995 TriaAccessor<structdim, dim, spacedim>::vertex_iterator(
996 const unsigned int i) const
997 {
998 return TriaIterator<TriaAccessor<0, dim, spacedim>>(this->tria,
999 0,
1000 vertex_index(i));
1001 }
1002
1003
1004
1005 template <int structdim, int dim, int spacedim>
1006 inline ReferenceCell::Type
reference_cell_type()1007 TriaAccessor<structdim, dim, spacedim>::reference_cell_type() const
1008 {
1009 if (structdim == 0)
1010 return ReferenceCell::Type::Vertex;
1011 else if (structdim == 1)
1012 return ReferenceCell::Type::Line;
1013 else if (structdim == dim)
1014 return this->tria->levels[this->present_level]
1015 ->reference_cell_type[this->present_index];
1016 else
1017 return this->tria->faces->quad_reference_cell_type[this->present_index];
1018 }
1019
1020
1021
1022 template <int structdim, int dim, int spacedim>
1023 inline unsigned int
vertex_index(const unsigned int corner)1024 TriaAccessor<structdim, dim, spacedim>::vertex_index(
1025 const unsigned int corner) const
1026 {
1027 AssertIndexRange(corner, this->n_vertices());
1028
1029 return dealii::internal::TriaAccessorImplementation::Implementation::
1030 vertex_index(*this, corner);
1031 }
1032
1033
1034
1035 template <int structdim, int dim, int spacedim>
1036 inline Point<spacedim> &
vertex(const unsigned int i)1037 TriaAccessor<structdim, dim, spacedim>::vertex(const unsigned int i) const
1038 {
1039 return const_cast<Point<spacedim> &>(this->tria->vertices[vertex_index(i)]);
1040 }
1041
1042
1043
1044 template <int structdim, int dim, int spacedim>
1045 inline typename dealii::internal::TriangulationImplementation::
1046 Iterators<dim, spacedim>::line_iterator
line(const unsigned int i)1047 TriaAccessor<structdim, dim, spacedim>::line(const unsigned int i) const
1048 {
1049 // checks happen in line_index
1050 return typename dealii::internal::TriangulationImplementation::
1051 Iterators<dim, spacedim>::line_iterator(this->tria, 0, line_index(i));
1052 }
1053
1054
1055
1056 template <int structdim, int dim, int spacedim>
1057 inline unsigned int
line_index(const unsigned int i)1058 TriaAccessor<structdim, dim, spacedim>::line_index(const unsigned int i) const
1059 {
1060 AssertIndexRange(i, this->n_lines());
1061
1062 return dealii::internal::TriaAccessorImplementation::Implementation::
1063 line_index(*this, i);
1064 }
1065
1066
1067
1068 template <int structdim, int dim, int spacedim>
1069 inline typename dealii::internal::TriangulationImplementation::
1070 Iterators<dim, spacedim>::quad_iterator
quad(const unsigned int i)1071 TriaAccessor<structdim, dim, spacedim>::quad(const unsigned int i) const
1072 {
1073 // checks happen in quad_index
1074 return typename dealii::internal::TriangulationImplementation::
1075 Iterators<dim, spacedim>::quad_iterator(this->tria, 0, quad_index(i));
1076 }
1077
1078
1079
1080 template <int structdim, int dim, int spacedim>
1081 inline unsigned int
quad_index(const unsigned int i)1082 TriaAccessor<structdim, dim, spacedim>::quad_index(const unsigned int i) const
1083 {
1084 return dealii::internal::TriaAccessorImplementation::Implementation::
1085 quad_index(*this, i);
1086 }
1087
1088
1089
1090 template <int structdim, int dim, int spacedim>
1091 inline bool
face_orientation(const unsigned int face)1092 TriaAccessor<structdim, dim, spacedim>::face_orientation(
1093 const unsigned int face) const
1094 {
1095 Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
1096
1097 return dealii::internal::TriaAccessorImplementation::Implementation::
1098 face_orientation(*this, face);
1099 }
1100
1101
1102
1103 template <int structdim, int dim, int spacedim>
1104 inline bool
face_flip(const unsigned int face)1105 TriaAccessor<structdim, dim, spacedim>::face_flip(const unsigned int face) const
1106 {
1107 Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
1108
1109 return dealii::internal::TriaAccessorImplementation::Implementation::
1110 face_flip(*this, face);
1111 }
1112
1113
1114 template <int structdim, int dim, int spacedim>
1115 inline bool
face_rotation(const unsigned int face)1116 TriaAccessor<structdim, dim, spacedim>::face_rotation(
1117 const unsigned int face) const
1118 {
1119 Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
1120
1121 return dealii::internal::TriaAccessorImplementation::Implementation::
1122 face_rotation(*this, face);
1123 }
1124
1125
1126
1127 template <int structdim, int dim, int spacedim>
1128 inline bool
line_orientation(const unsigned int line)1129 TriaAccessor<structdim, dim, spacedim>::line_orientation(
1130 const unsigned int line) const
1131 {
1132 Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
1133 AssertIndexRange(line, this->n_lines());
1134
1135 return dealii::internal::TriaAccessorImplementation::Implementation::
1136 line_orientation(*this, line);
1137 }
1138
1139
1140
1141 template <int structdim, int dim, int spacedim>
1142 inline void
set_face_orientation(const unsigned int face,const bool value)1143 TriaAccessor<structdim, dim, spacedim>::set_face_orientation(
1144 const unsigned int face,
1145 const bool value) const
1146 {
1147 Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
1148
1149 dealii::internal::TriaAccessorImplementation::Implementation::
1150 set_face_orientation(*this, face, value);
1151 }
1152
1153
1154
1155 template <int structdim, int dim, int spacedim>
1156 inline void
set_face_flip(const unsigned int face,const bool value)1157 TriaAccessor<structdim, dim, spacedim>::set_face_flip(const unsigned int face,
1158 const bool value) const
1159 {
1160 Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
1161
1162 dealii::internal::TriaAccessorImplementation::Implementation::set_face_flip(
1163 *this, face, value);
1164 }
1165
1166
1167 template <int structdim, int dim, int spacedim>
1168 inline void
set_face_rotation(const unsigned int face,const bool value)1169 TriaAccessor<structdim, dim, spacedim>::set_face_rotation(
1170 const unsigned int face,
1171 const bool value) const
1172 {
1173 Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
1174
1175 dealii::internal::TriaAccessorImplementation::Implementation::
1176 set_face_rotation(*this, face, value);
1177 }
1178
1179
1180
1181 template <int structdim, int dim, int spacedim>
1182 inline void
set_line_orientation(const unsigned int line,const bool value)1183 TriaAccessor<structdim, dim, spacedim>::set_line_orientation(
1184 const unsigned int line,
1185 const bool value) const
1186 {
1187 Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
1188 AssertIndexRange(line, this->n_lines());
1189
1190 dealii::internal::TriaAccessorImplementation::Implementation::
1191 set_line_orientation(*this, line, value);
1192 }
1193
1194
1195
1196 template <int structdim, int dim, int spacedim>
1197 void
set_used_flag()1198 TriaAccessor<structdim, dim, spacedim>::set_used_flag() const
1199 {
1200 Assert(this->state() == IteratorState::valid,
1201 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
1202 *this));
1203 this->objects().used[this->present_index] = true;
1204 }
1205
1206
1207
1208 template <int structdim, int dim, int spacedim>
1209 void
clear_used_flag()1210 TriaAccessor<structdim, dim, spacedim>::clear_used_flag() const
1211 {
1212 Assert(this->state() == IteratorState::valid,
1213 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
1214 *this));
1215 this->objects().used[this->present_index] = false;
1216 }
1217
1218
1219 template <int structdim, int dim, int spacedim>
1220 int
child_index(const unsigned int i)1221 TriaAccessor<structdim, dim, spacedim>::child_index(const unsigned int i) const
1222 {
1223 Assert(has_children(), TriaAccessorExceptions::ExcCellHasNoChildren());
1224 AssertIndexRange(i, n_children());
1225
1226 // each set of two children are stored
1227 // consecutively, so we only have to find
1228 // the location of the set of children
1229 const unsigned int n_sets_of_two =
1230 GeometryInfo<structdim>::max_children_per_cell / 2;
1231 return this->objects().children[n_sets_of_two * this->present_index + i / 2] +
1232 i % 2;
1233 }
1234
1235
1236
1237 template <int structdim, int dim, int spacedim>
1238 int
isotropic_child_index(const unsigned int i)1239 TriaAccessor<structdim, dim, spacedim>::isotropic_child_index(
1240 const unsigned int i) const
1241 {
1242 AssertIndexRange(i, GeometryInfo<structdim>::max_children_per_cell);
1243
1244 switch (structdim)
1245 {
1246 case 1:
1247 return child_index(i);
1248 case 2:
1249 {
1250 const RefinementCase<2> this_refinement_case(
1251 static_cast<std::uint8_t>(refinement_case()));
1252
1253 Assert(this_refinement_case != RefinementCase<2>::no_refinement,
1254 TriaAccessorExceptions::ExcCellHasNoChildren());
1255
1256 if (this_refinement_case == RefinementCase<2>::cut_xy)
1257 return child_index(i);
1258 else if ((this_refinement_case == RefinementCase<2>::cut_x) &&
1259 (child(i % 2)->refinement_case() ==
1260 RefinementCase<2>::cut_y))
1261 return child(i % 2)->child_index(i / 2);
1262 else if ((this_refinement_case == RefinementCase<2>::cut_y) &&
1263 (child(i / 2)->refinement_case() ==
1264 RefinementCase<2>::cut_x))
1265 return child(i / 2)->child_index(i % 2);
1266 else
1267 Assert(
1268 false,
1269 ExcMessage(
1270 "This cell has no grandchildren equivalent to isotropic refinement"));
1271 break;
1272 }
1273
1274 case 3:
1275 Assert(false, ExcNotImplemented());
1276 }
1277 return -1;
1278 }
1279
1280
1281
1282 template <int structdim, int dim, int spacedim>
1283 RefinementCase<structdim>
refinement_case()1284 TriaAccessor<structdim, dim, spacedim>::refinement_case() const
1285 {
1286 Assert(this->state() == IteratorState::valid,
1287 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
1288 *this));
1289
1290 switch (structdim)
1291 {
1292 case 1:
1293 return (RefinementCase<structdim>(
1294 this->objects().children[this->present_index] != -1 ?
1295 // cast the branches
1296 // here first to uchar
1297 // and then (above) to
1298 // RefinementCase<structdim>
1299 // so that the
1300 // conversion is valid
1301 // even for the case
1302 // structdim>1 (for
1303 // which this part of
1304 // the code is dead
1305 // anyway)
1306 static_cast<std::uint8_t>(RefinementCase<1>::cut_x) :
1307 static_cast<std::uint8_t>(RefinementCase<1>::no_refinement)));
1308
1309 default:
1310 Assert(static_cast<unsigned int>(this->present_index) <
1311 this->objects().refinement_cases.size(),
1312 ExcIndexRange(this->present_index,
1313 0,
1314 this->objects().refinement_cases.size()));
1315
1316 return (static_cast<RefinementCase<structdim>>(
1317 this->objects().refinement_cases[this->present_index]));
1318 }
1319 }
1320
1321
1322
1323 template <int structdim, int dim, int spacedim>
1324 inline TriaIterator<TriaAccessor<structdim, dim, spacedim>>
child(const unsigned int i)1325 TriaAccessor<structdim, dim, spacedim>::child(const unsigned int i) const
1326
1327 {
1328 // checking of 'i' happens in child_index
1329 const TriaIterator<TriaAccessor<structdim, dim, spacedim>> q(
1330 this->tria, (dim == structdim ? this->level() + 1 : 0), child_index(i));
1331
1332 Assert((q.state() == IteratorState::past_the_end) || q->used(),
1333 ExcInternalError());
1334
1335 return q;
1336 }
1337
1338
1339
1340 template <int structdim, int dim, int spacedim>
1341 inline unsigned int
child_iterator_to_index(const TriaIterator<TriaAccessor<structdim,dim,spacedim>> & child)1342 TriaAccessor<structdim, dim, spacedim>::child_iterator_to_index(
1343 const TriaIterator<TriaAccessor<structdim, dim, spacedim>> &child) const
1344 {
1345 const auto n_children = this->n_children();
1346 for (unsigned int child_n = 0; child_n < n_children; ++child_n)
1347 if (this->child(child_n) == child)
1348 return child_n;
1349
1350 Assert(false,
1351 ExcMessage("The given child is not a child of the current object."));
1352 return numbers::invalid_unsigned_int;
1353 }
1354
1355
1356
1357 template <int structdim, int dim, int spacedim>
1358 inline TriaIterator<TriaAccessor<structdim, dim, spacedim>>
isotropic_child(const unsigned int i)1359 TriaAccessor<structdim, dim, spacedim>::isotropic_child(
1360 const unsigned int i) const
1361 {
1362 // checking of 'i' happens in child() or
1363 // child_index() called below
1364 switch (structdim)
1365 {
1366 case 1:
1367 // no anisotropic refinement in 1D
1368 return child(i);
1369
1370 case 2:
1371 {
1372 const RefinementCase<2> this_refinement_case(
1373 static_cast<std::uint8_t>(refinement_case()));
1374
1375 Assert(this_refinement_case != RefinementCase<2>::no_refinement,
1376 TriaAccessorExceptions::ExcCellHasNoChildren());
1377
1378 if (this_refinement_case == RefinementCase<2>::cut_xy)
1379 return child(i);
1380 else if ((this_refinement_case == RefinementCase<2>::cut_x) &&
1381 (child(i % 2)->refinement_case() ==
1382 RefinementCase<2>::cut_y))
1383 return child(i % 2)->child(i / 2);
1384 else if ((this_refinement_case == RefinementCase<2>::cut_y) &&
1385 (child(i / 2)->refinement_case() ==
1386 RefinementCase<2>::cut_x))
1387 return child(i / 2)->child(i % 2);
1388 else
1389 Assert(
1390 false,
1391 ExcMessage(
1392 "This cell has no grandchildren equivalent to isotropic refinement"));
1393 break;
1394 }
1395
1396 default:
1397 Assert(false, ExcNotImplemented());
1398 }
1399 // we don't get here but have to return
1400 // something...
1401 return child(0);
1402 }
1403
1404
1405
1406 template <int structdim, int dim, int spacedim>
1407 inline bool
has_children()1408 TriaAccessor<structdim, dim, spacedim>::has_children() const
1409 {
1410 Assert(this->state() == IteratorState::valid,
1411 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
1412 *this));
1413
1414 // each set of two children are stored
1415 // consecutively, so we only have to find
1416 // the location of the set of children
1417 const unsigned int n_sets_of_two =
1418 GeometryInfo<structdim>::max_children_per_cell / 2;
1419 return (this->objects().children[n_sets_of_two * this->present_index] != -1);
1420 }
1421
1422
1423
1424 template <int structdim, int dim, int spacedim>
1425 inline unsigned int
n_children()1426 TriaAccessor<structdim, dim, spacedim>::n_children() const
1427 {
1428 return GeometryInfo<structdim>::n_children(refinement_case());
1429 }
1430
1431
1432
1433 template <int structdim, int dim, int spacedim>
1434 inline void
set_refinement_case(const RefinementCase<structdim> & refinement_case)1435 TriaAccessor<structdim, dim, spacedim>::set_refinement_case(
1436 const RefinementCase<structdim> &refinement_case) const
1437 {
1438 Assert(this->state() == IteratorState::valid,
1439 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
1440 *this));
1441 Assert(static_cast<unsigned int>(this->present_index) <
1442 this->objects().refinement_cases.size(),
1443 ExcIndexRange(this->present_index,
1444 0,
1445 this->objects().refinement_cases.size()));
1446
1447 this->objects().refinement_cases[this->present_index] = refinement_case;
1448 }
1449
1450
1451 template <int structdim, int dim, int spacedim>
1452 inline void
clear_refinement_case()1453 TriaAccessor<structdim, dim, spacedim>::clear_refinement_case() const
1454 {
1455 Assert(this->state() == IteratorState::valid,
1456 TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
1457 *this));
1458 Assert(static_cast<unsigned int>(this->present_index) <
1459 this->objects().refinement_cases.size(),
1460 ExcIndexRange(this->present_index,
1461 0,
1462 this->objects().refinement_cases.size()));
1463
1464 this->objects().refinement_cases[this->present_index] =
1465 RefinementCase<structdim>::no_refinement;
1466 }
1467
1468
1469
1470 template <int structdim, int dim, int spacedim>
1471 void
set_children(const unsigned int i,const int index)1472 TriaAccessor<structdim, dim, spacedim>::set_children(const unsigned int i,
1473 const int index) const
1474 {
1475 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1476 Assert(i % 2 == 0, TriaAccessorExceptions::ExcSetOnlyEvenChildren(i));
1477
1478 // each set of two children are stored
1479 // consecutively, so we only have to find
1480 // the location of the set of children
1481 const unsigned int n_sets_of_two =
1482 GeometryInfo<structdim>::max_children_per_cell / 2;
1483
1484 Assert(
1485 // clearing the child index for a cell
1486 (index == -1) ||
1487 // if setting the child index for the i'th child (with i==0),
1488 // then the index must be a non-negative number
1489 (i == 0 && !this->has_children() && (index >= 0)) ||
1490 // if setting the child index for the i'th child (with i>0),
1491 // then the previously stored index must be the invalid
1492 // index
1493 (i > 0 && this->has_children() && (index >= 0) &&
1494 this->objects().children[n_sets_of_two * this->present_index + i / 2] ==
1495 -1),
1496 TriaAccessorExceptions::ExcCantSetChildren(index));
1497
1498 this->objects().children[n_sets_of_two * this->present_index + i / 2] = index;
1499 }
1500
1501
1502
1503 template <int structdim, int dim, int spacedim>
1504 void
clear_children()1505 TriaAccessor<structdim, dim, spacedim>::clear_children() const
1506 {
1507 // each set of two children are stored
1508 // consecutively, so we only have to find
1509 // the location of the set of children
1510 const unsigned int n_sets_of_two =
1511 GeometryInfo<structdim>::max_children_per_cell / 2;
1512
1513 for (unsigned int i = 0; i < n_sets_of_two; ++i)
1514 set_children(2 * i, -1);
1515 }
1516
1517
1518
1519 template <int structdim, int dim, int spacedim>
1520 inline bool
user_flag_set()1521 TriaAccessor<structdim, dim, spacedim>::user_flag_set() const
1522 {
1523 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1524 return this->objects().user_flags[this->present_index];
1525 }
1526
1527
1528
1529 template <int structdim, int dim, int spacedim>
1530 inline void
set_user_flag()1531 TriaAccessor<structdim, dim, spacedim>::set_user_flag() const
1532 {
1533 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1534 this->objects().user_flags[this->present_index] = true;
1535 }
1536
1537
1538
1539 template <int structdim, int dim, int spacedim>
1540 inline void
clear_user_flag()1541 TriaAccessor<structdim, dim, spacedim>::clear_user_flag() const
1542 {
1543 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1544 this->objects().user_flags[this->present_index] = false;
1545 }
1546
1547
1548
1549 template <int structdim, int dim, int spacedim>
1550 void
recursively_set_user_flag()1551 TriaAccessor<structdim, dim, spacedim>::recursively_set_user_flag() const
1552 {
1553 set_user_flag();
1554
1555 if (this->has_children())
1556 for (unsigned int c = 0; c < this->n_children(); ++c)
1557 this->child(c)->recursively_set_user_flag();
1558 }
1559
1560
1561
1562 template <int structdim, int dim, int spacedim>
1563 void
recursively_clear_user_flag()1564 TriaAccessor<structdim, dim, spacedim>::recursively_clear_user_flag() const
1565 {
1566 clear_user_flag();
1567
1568 if (this->has_children())
1569 for (unsigned int c = 0; c < this->n_children(); ++c)
1570 this->child(c)->recursively_clear_user_flag();
1571 }
1572
1573
1574
1575 template <int structdim, int dim, int spacedim>
1576 void
clear_user_data()1577 TriaAccessor<structdim, dim, spacedim>::clear_user_data() const
1578 {
1579 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1580 this->objects().clear_user_data(this->present_index);
1581 }
1582
1583
1584
1585 template <int structdim, int dim, int spacedim>
1586 void
set_user_pointer(void * p)1587 TriaAccessor<structdim, dim, spacedim>::set_user_pointer(void *p) const
1588 {
1589 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1590 this->objects().user_pointer(this->present_index) = p;
1591 }
1592
1593
1594
1595 template <int structdim, int dim, int spacedim>
1596 void
clear_user_pointer()1597 TriaAccessor<structdim, dim, spacedim>::clear_user_pointer() const
1598 {
1599 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1600 this->objects().user_pointer(this->present_index) = nullptr;
1601 }
1602
1603
1604
1605 template <int structdim, int dim, int spacedim>
1606 void *
user_pointer()1607 TriaAccessor<structdim, dim, spacedim>::user_pointer() const
1608 {
1609 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1610 return this->objects().user_pointer(this->present_index);
1611 }
1612
1613
1614
1615 template <int structdim, int dim, int spacedim>
1616 void
recursively_set_user_pointer(void * p)1617 TriaAccessor<structdim, dim, spacedim>::recursively_set_user_pointer(
1618 void *p) const
1619 {
1620 set_user_pointer(p);
1621
1622 if (this->has_children())
1623 for (unsigned int c = 0; c < this->n_children(); ++c)
1624 this->child(c)->recursively_set_user_pointer(p);
1625 }
1626
1627
1628
1629 template <int structdim, int dim, int spacedim>
1630 void
recursively_clear_user_pointer()1631 TriaAccessor<structdim, dim, spacedim>::recursively_clear_user_pointer() const
1632 {
1633 clear_user_pointer();
1634
1635 if (this->has_children())
1636 for (unsigned int c = 0; c < this->n_children(); ++c)
1637 this->child(c)->recursively_clear_user_pointer();
1638 }
1639
1640
1641
1642 template <int structdim, int dim, int spacedim>
1643 void
set_user_index(unsigned int p)1644 TriaAccessor<structdim, dim, spacedim>::set_user_index(unsigned int p) const
1645 {
1646 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1647 this->objects().user_index(this->present_index) = p;
1648 }
1649
1650
1651
1652 template <int structdim, int dim, int spacedim>
1653 void
clear_user_index()1654 TriaAccessor<structdim, dim, spacedim>::clear_user_index() const
1655 {
1656 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1657 this->objects().user_index(this->present_index) = 0;
1658 }
1659
1660
1661
1662 template <int structdim, int dim, int spacedim>
1663 unsigned int
user_index()1664 TriaAccessor<structdim, dim, spacedim>::user_index() const
1665 {
1666 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1667 return this->objects().user_index(this->present_index);
1668 }
1669
1670
1671
1672 template <int structdim, int dim, int spacedim>
1673 void
recursively_set_user_index(unsigned int p)1674 TriaAccessor<structdim, dim, spacedim>::recursively_set_user_index(
1675 unsigned int p) const
1676 {
1677 set_user_index(p);
1678
1679 if (this->has_children())
1680 for (unsigned int c = 0; c < this->n_children(); ++c)
1681 this->child(c)->recursively_set_user_index(p);
1682 }
1683
1684
1685
1686 template <int structdim, int dim, int spacedim>
1687 void
recursively_clear_user_index()1688 TriaAccessor<structdim, dim, spacedim>::recursively_clear_user_index() const
1689 {
1690 clear_user_index();
1691
1692 if (this->has_children())
1693 for (unsigned int c = 0; c < this->n_children(); ++c)
1694 this->child(c)->recursively_clear_user_index();
1695 }
1696
1697
1698
1699 template <int structdim, int dim, int spacedim>
1700 inline unsigned int
max_refinement_depth()1701 TriaAccessor<structdim, dim, spacedim>::max_refinement_depth() const
1702 {
1703 if (!this->has_children())
1704 return 0;
1705
1706 unsigned int max_depth = 1;
1707 for (unsigned int c = 0; c < n_children(); ++c)
1708 max_depth = std::max(max_depth, child(c)->max_refinement_depth() + 1);
1709 return max_depth;
1710 }
1711
1712
1713
1714 template <int structdim, int dim, int spacedim>
1715 unsigned int
number_of_children()1716 TriaAccessor<structdim, dim, spacedim>::number_of_children() const
1717 {
1718 if (!this->has_children())
1719 return 1;
1720 else
1721 {
1722 unsigned int sum = 0;
1723 for (unsigned int c = 0; c < n_children(); ++c)
1724 sum += this->child(c)->number_of_children();
1725 return sum;
1726 }
1727 }
1728
1729
1730
1731 template <int structdim, int dim, int spacedim>
1732 types::boundary_id
boundary_id()1733 TriaAccessor<structdim, dim, spacedim>::boundary_id() const
1734 {
1735 Assert(structdim < dim, ExcImpossibleInDim(dim));
1736 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1737
1738 return this->objects()
1739 .boundary_or_material_id[this->present_index]
1740 .boundary_id;
1741 }
1742
1743
1744
1745 template <int structdim, int dim, int spacedim>
1746 void
set_boundary_id(const types::boundary_id boundary_ind)1747 TriaAccessor<structdim, dim, spacedim>::set_boundary_id(
1748 const types::boundary_id boundary_ind) const
1749 {
1750 Assert(structdim < dim, ExcImpossibleInDim(dim));
1751 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1752 Assert(
1753 boundary_ind != numbers::internal_face_boundary_id,
1754 ExcMessage(
1755 "You are trying to set the boundary_id to an illegal value (numbers::internal_face_boundary_id is reserved)."));
1756 Assert(
1757 this->at_boundary(),
1758 ExcMessage(
1759 "You are trying to set the boundary_id of an internal object, which is illegal!"));
1760
1761 this->objects().boundary_or_material_id[this->present_index].boundary_id =
1762 boundary_ind;
1763 }
1764
1765
1766
1767 template <int structdim, int dim, int spacedim>
1768 void
set_boundary_id_internal(const types::boundary_id boundary_ind)1769 TriaAccessor<structdim, dim, spacedim>::set_boundary_id_internal(
1770 const types::boundary_id boundary_ind) const
1771 {
1772 Assert(structdim < dim, ExcImpossibleInDim(dim));
1773 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1774
1775 this->objects().boundary_or_material_id[this->present_index].boundary_id =
1776 boundary_ind;
1777 }
1778
1779
1780
1781 template <int structdim, int dim, int spacedim>
1782 void
set_all_boundary_ids(const types::boundary_id boundary_ind)1783 TriaAccessor<structdim, dim, spacedim>::set_all_boundary_ids(
1784 const types::boundary_id boundary_ind) const
1785 {
1786 set_boundary_id(boundary_ind);
1787
1788 switch (structdim)
1789 {
1790 case 1:
1791 // 1d objects have no sub-objects
1792 // where we have to do anything
1793 break;
1794
1795 case 2:
1796 // for boundary quads also set
1797 // boundary_id of bounding lines
1798 for (unsigned int i = 0; i < 4; ++i)
1799 this->line(i)->set_boundary_id(boundary_ind);
1800 break;
1801
1802 default:
1803 Assert(false, ExcNotImplemented());
1804 }
1805 }
1806
1807
1808
1809 template <int structdim, int dim, int spacedim>
1810 bool
at_boundary()1811 TriaAccessor<structdim, dim, spacedim>::at_boundary() const
1812 {
1813 // error checking is done
1814 // in boundary_id()
1815 return (boundary_id() != numbers::internal_face_boundary_id);
1816 }
1817
1818
1819
1820 template <int structdim, int dim, int spacedim>
1821 const Manifold<dim, spacedim> &
get_manifold()1822 TriaAccessor<structdim, dim, spacedim>::get_manifold() const
1823 {
1824 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1825 return this->tria->get_manifold(this->manifold_id());
1826 }
1827
1828
1829 template <int structdim, int dim, int spacedim>
1830 types::manifold_id
manifold_id()1831 TriaAccessor<structdim, dim, spacedim>::manifold_id() const
1832 {
1833 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1834
1835 return this->objects().manifold_id[this->present_index];
1836 }
1837
1838
1839
1840 template <int structdim, int dim, int spacedim>
1841 void
set_manifold_id(const types::manifold_id manifold_ind)1842 TriaAccessor<structdim, dim, spacedim>::set_manifold_id(
1843 const types::manifold_id manifold_ind) const
1844 {
1845 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1846
1847 this->objects().manifold_id[this->present_index] = manifold_ind;
1848 }
1849
1850
1851 template <int structdim, int dim, int spacedim>
1852 void
set_all_manifold_ids(const types::manifold_id manifold_ind)1853 TriaAccessor<structdim, dim, spacedim>::set_all_manifold_ids(
1854 const types::manifold_id manifold_ind) const
1855 {
1856 set_manifold_id(manifold_ind);
1857
1858 if (this->has_children())
1859 for (unsigned int c = 0; c < this->n_children(); ++c)
1860 this->child(c)->set_all_manifold_ids(manifold_ind);
1861
1862 switch (structdim)
1863 {
1864 case 1:
1865 if (dim == 1)
1866 {
1867 (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(0)] =
1868 manifold_ind;
1869 (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(1)] =
1870 manifold_ind;
1871 }
1872 break;
1873
1874 case 2:
1875 // for quads also set manifold_id of bounding lines
1876 for (unsigned int i = 0; i < 4; ++i)
1877 this->line(i)->set_manifold_id(manifold_ind);
1878 break;
1879 default:
1880 Assert(false, ExcNotImplemented());
1881 }
1882 }
1883
1884
1885
1886 template <int structdim, int dim, int spacedim>
1887 double
diameter()1888 TriaAccessor<structdim, dim, spacedim>::diameter() const
1889 {
1890 switch (this->reference_cell_type())
1891 {
1892 case ReferenceCell::Type::Line:
1893 return (this->vertex(1) - this->vertex(0)).norm();
1894 case ReferenceCell::Type::Tri:
1895 return std::max(std::max((this->vertex(1) - this->vertex(0)).norm(),
1896 (this->vertex(2) - this->vertex(1)).norm()),
1897 (this->vertex(2) - this->vertex(0)).norm());
1898 case ReferenceCell::Type::Quad:
1899 return std::max((this->vertex(3) - this->vertex(0)).norm(),
1900 (this->vertex(2) - this->vertex(1)).norm());
1901 case ReferenceCell::Type::Tet:
1902 return std::max(
1903 std::max(std::max((this->vertex(1) - this->vertex(0)).norm(),
1904 (this->vertex(2) - this->vertex(0)).norm()),
1905 std::max((this->vertex(2) - this->vertex(1)).norm(),
1906 (this->vertex(3) - this->vertex(0)).norm())),
1907 std::max((this->vertex(3) - this->vertex(1)).norm(),
1908 (this->vertex(3) - this->vertex(2)).norm()));
1909 case ReferenceCell::Type::Hex:
1910 return std::max(std::max((this->vertex(7) - this->vertex(0)).norm(),
1911 (this->vertex(6) - this->vertex(1)).norm()),
1912 std::max((this->vertex(2) - this->vertex(5)).norm(),
1913 (this->vertex(3) - this->vertex(4)).norm()));
1914 default:
1915 Assert(false, ExcNotImplemented());
1916 return -1e10;
1917 }
1918 }
1919
1920
1921
1922 template <int structdim, int dim, int spacedim>
1923 std::pair<Point<spacedim>, double>
enclosing_ball()1924 TriaAccessor<structdim, dim, spacedim>::enclosing_ball() const
1925 {
1926 // If the object is one dimensional,
1927 // the enclosing ball is the initial iterate
1928 // i.e., the ball's center and diameter are
1929 // the center and the diameter of the object.
1930 if (structdim == 1)
1931 return std::make_pair((this->vertex(1) + this->vertex(0)) * 0.5,
1932 (this->vertex(1) - this->vertex(0)).norm() * 0.5);
1933
1934 // The list is_initial_guess_vertex contains bool values and has
1935 // the same size as the number of vertices per object.
1936 // The entries of is_initial_guess_vertex are set true only for those
1937 // two vertices corresponding to the largest diagonal which is being used
1938 // to construct the initial ball.
1939 // We employ this mask to skip these two vertices while enlarging the ball.
1940 std::vector<bool> is_initial_guess_vertex(this->n_vertices());
1941
1942 // First let all the vertices be outside
1943 std::fill(is_initial_guess_vertex.begin(),
1944 is_initial_guess_vertex.end(),
1945 false);
1946
1947 // Get an initial guess by looking at the largest diagonal
1948 Point<spacedim> center;
1949 double radius = 0;
1950
1951 switch (structdim)
1952 {
1953 case 2:
1954 {
1955 const Point<spacedim> p30(this->vertex(3) - this->vertex(0));
1956 const Point<spacedim> p21(this->vertex(2) - this->vertex(1));
1957 if (p30.norm() > p21.norm())
1958 {
1959 center = this->vertex(0) + 0.5 * p30;
1960 radius = p30.norm() / 2.;
1961 is_initial_guess_vertex[3] = true;
1962 is_initial_guess_vertex[0] = true;
1963 }
1964 else
1965 {
1966 center = this->vertex(1) + 0.5 * p21;
1967 radius = p21.norm() / 2.;
1968 is_initial_guess_vertex[2] = true;
1969 is_initial_guess_vertex[1] = true;
1970 }
1971 break;
1972 }
1973 case 3:
1974 {
1975 const Point<spacedim> p70(this->vertex(7) - this->vertex(0));
1976 const Point<spacedim> p61(this->vertex(6) - this->vertex(1));
1977 const Point<spacedim> p25(this->vertex(2) - this->vertex(5));
1978 const Point<spacedim> p34(this->vertex(3) - this->vertex(4));
1979 const std::vector<double> diagonals = {p70.norm(),
1980 p61.norm(),
1981 p25.norm(),
1982 p34.norm()};
1983 const std::vector<double>::const_iterator it =
1984 std::max_element(diagonals.begin(), diagonals.end());
1985 if (it == diagonals.begin())
1986 {
1987 center = this->vertex(0) + 0.5 * p70;
1988 is_initial_guess_vertex[7] = true;
1989 is_initial_guess_vertex[0] = true;
1990 }
1991 else if (it == diagonals.begin() + 1)
1992 {
1993 center = this->vertex(1) + 0.5 * p61;
1994 is_initial_guess_vertex[6] = true;
1995 is_initial_guess_vertex[1] = true;
1996 }
1997 else if (it == diagonals.begin() + 2)
1998 {
1999 center = this->vertex(5) + 0.5 * p25;
2000 is_initial_guess_vertex[2] = true;
2001 is_initial_guess_vertex[5] = true;
2002 }
2003 else
2004 {
2005 center = this->vertex(4) + 0.5 * p34;
2006 is_initial_guess_vertex[3] = true;
2007 is_initial_guess_vertex[4] = true;
2008 }
2009 radius = *it * 0.5;
2010 break;
2011 }
2012 default:
2013 Assert(false, ExcNotImplemented());
2014 return std::pair<Point<spacedim>, double>();
2015 }
2016
2017 // For each vertex that is found to be geometrically outside the ball
2018 // enlarge the ball so that the new ball contains both the previous ball
2019 // and the given vertex.
2020 for (const unsigned int v : this->vertex_indices())
2021 if (!is_initial_guess_vertex[v])
2022 {
2023 const double distance = center.distance(this->vertex(v));
2024 if (distance > radius)
2025 {
2026 // we found a vertex which is outside of the ball
2027 // extend it (move center and change radius)
2028 const Point<spacedim> pCV(center - this->vertex(v));
2029 radius = (distance + radius) * 0.5;
2030 center = this->vertex(v) + pCV * (radius / distance);
2031
2032 // Now the new ball constructed in this block
2033 // encloses the vertex (v) that was found to be geometrically
2034 // outside the old ball.
2035 }
2036 }
2037 #ifdef DEBUG
2038 bool all_vertices_within_ball = true;
2039
2040 // Set all_vertices_within_ball false if any of the vertices of the object
2041 // are geometrically outside the ball
2042 for (const unsigned int v : this->vertex_indices())
2043 if (center.distance(this->vertex(v)) >
2044 radius + 100. * std::numeric_limits<double>::epsilon())
2045 {
2046 all_vertices_within_ball = false;
2047 break;
2048 }
2049 // If all the vertices are not within the ball throw error
2050 Assert(all_vertices_within_ball, ExcInternalError());
2051 #endif
2052 return std::make_pair(center, radius);
2053 }
2054
2055
2056 template <int structdim, int dim, int spacedim>
2057 double
minimum_vertex_distance()2058 TriaAccessor<structdim, dim, spacedim>::minimum_vertex_distance() const
2059 {
2060 switch (structdim)
2061 {
2062 case 1:
2063 return (this->vertex(1) - this->vertex(0)).norm();
2064 case 2:
2065 case 3:
2066 {
2067 double min = std::numeric_limits<double>::max();
2068 for (const unsigned int i : this->vertex_indices())
2069 for (unsigned int j = i + 1; j < this->n_vertices(); ++j)
2070 min = std::min(min,
2071 (this->vertex(i) - this->vertex(j)) *
2072 (this->vertex(i) - this->vertex(j)));
2073 return std::sqrt(min);
2074 }
2075 default:
2076 Assert(false, ExcNotImplemented());
2077 return -1e10;
2078 }
2079 }
2080
2081
2082 template <int structdim, int dim, int spacedim>
2083 bool
is_translation_of(const TriaIterator<TriaAccessor<structdim,dim,spacedim>> & o)2084 TriaAccessor<structdim, dim, spacedim>::is_translation_of(
2085 const TriaIterator<TriaAccessor<structdim, dim, spacedim>> &o) const
2086 {
2087 // go through the vertices and check... The
2088 // cell is a translation of the previous
2089 // one in case the distance between the
2090 // individual vertices in the two cell is
2091 // the same for all the vertices. So do the
2092 // check by first getting the distance on
2093 // the first vertex, and then checking
2094 // whether all others have the same down to
2095 // rounding errors (we have to be careful
2096 // here because the calculation of the
2097 // displacement between one cell and the
2098 // next can already result in the loss of
2099 // one or two digits), so we choose 1e-12
2100 // times the distance between the zeroth
2101 // vertices here.
2102 bool is_translation = true;
2103 const Tensor<1, spacedim> dist = o->vertex(0) - this->vertex(0);
2104 const double tol_square = 1e-24 * dist.norm_square();
2105 for (unsigned int i = 1; i < this->n_vertices(); ++i)
2106 {
2107 const Tensor<1, spacedim> dist_new =
2108 (o->vertex(i) - this->vertex(i)) - dist;
2109 if (dist_new.norm_square() > tol_square)
2110 {
2111 is_translation = false;
2112 break;
2113 }
2114 }
2115 return is_translation;
2116 }
2117
2118
2119
2120 template <int structdim, int dim, int spacedim>
2121 unsigned int
n_vertices()2122 TriaAccessor<structdim, dim, spacedim>::n_vertices() const
2123 {
2124 return this->reference_cell_info().n_vertices();
2125 }
2126
2127
2128
2129 template <int structdim, int dim, int spacedim>
2130 unsigned int
n_lines()2131 TriaAccessor<structdim, dim, spacedim>::n_lines() const
2132 {
2133 return this->reference_cell_info().n_lines();
2134 }
2135
2136
2137
2138 template <int structdim, int dim, int spacedim>
2139 unsigned int
n_faces()2140 TriaAccessor<structdim, dim, spacedim>::n_faces() const
2141 {
2142 AssertDimension(structdim, dim);
2143
2144 return this->reference_cell_info().n_faces();
2145 }
2146
2147
2148
2149 template <int structdim, int dim, int spacedim>
2150 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
vertex_indices()2151 TriaAccessor<structdim, dim, spacedim>::vertex_indices() const
2152 {
2153 return {0U, n_vertices()};
2154 }
2155
2156
2157
2158 template <int structdim, int dim, int spacedim>
2159 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
line_indices()2160 TriaAccessor<structdim, dim, spacedim>::line_indices() const
2161 {
2162 return {0U, n_lines()};
2163 }
2164
2165
2166
2167 template <int structdim, int dim, int spacedim>
2168 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
face_indices()2169 TriaAccessor<structdim, dim, spacedim>::face_indices() const
2170 {
2171 return {0U, n_faces()};
2172 }
2173
2174
2175
2176 template <int structdim, int dim, int spacedim>
2177 inline const ReferenceCell::internal::Info::Base &
reference_cell_info()2178 TriaAccessor<structdim, dim, spacedim>::reference_cell_info() const
2179 {
2180 if (structdim == 0)
2181 return ReferenceCell::internal::Info::get_cell(ReferenceCell::Type::Vertex);
2182 else if (structdim == 1)
2183 return ReferenceCell::internal::Info::get_cell(ReferenceCell::Type::Line);
2184 else
2185 return ReferenceCell::internal::Info::get_cell(this->reference_cell_type());
2186 }
2187
2188
2189 /*----------------- Functions: TriaAccessor<0,dim,spacedim> -----------------*/
2190
2191 template <int dim, int spacedim>
TriaAccessor(const Triangulation<dim,spacedim> * tria,const unsigned int vertex_index)2192 inline TriaAccessor<0, dim, spacedim>::TriaAccessor(
2193 const Triangulation<dim, spacedim> *tria,
2194 const unsigned int vertex_index)
2195 : tria(tria)
2196 , global_vertex_index(vertex_index)
2197 {}
2198
2199
2200
2201 template <int dim, int spacedim>
TriaAccessor(const Triangulation<dim,spacedim> * tria,const int,const int index,const AccessorData *)2202 inline TriaAccessor<0, dim, spacedim>::TriaAccessor(
2203 const Triangulation<dim, spacedim> *tria,
2204 const int /*level*/,
2205 const int index,
2206 const AccessorData *)
2207 : tria(tria)
2208 , global_vertex_index(index)
2209 {}
2210
2211
2212
2213 template <int dim, int spacedim>
2214 template <int structdim2, int dim2, int spacedim2>
TriaAccessor(const TriaAccessor<structdim2,dim2,spacedim2> &)2215 inline TriaAccessor<0, dim, spacedim>::TriaAccessor(
2216 const TriaAccessor<structdim2, dim2, spacedim2> &)
2217 : tria(nullptr)
2218 , global_vertex_index(numbers::invalid_unsigned_int)
2219 {
2220 Assert(false, ExcImpossibleInDim(0));
2221 }
2222
2223
2224
2225 template <int dim, int spacedim>
2226 template <int structdim2, int dim2, int spacedim2>
TriaAccessor(const InvalidAccessor<structdim2,dim2,spacedim2> &)2227 inline TriaAccessor<0, dim, spacedim>::TriaAccessor(
2228 const InvalidAccessor<structdim2, dim2, spacedim2> &)
2229 : tria(nullptr)
2230 , global_vertex_index(numbers::invalid_unsigned_int)
2231 {
2232 Assert(false, ExcImpossibleInDim(0));
2233 }
2234
2235
2236
2237 template <int dim, int spacedim>
2238 inline void
copy_from(const TriaAccessor & t)2239 TriaAccessor<0, dim, spacedim>::copy_from(const TriaAccessor &t)
2240 {
2241 tria = t.tria;
2242 global_vertex_index = t.global_vertex_index;
2243 }
2244
2245
2246
2247 template <int dim, int spacedim>
2248 inline bool
2249 TriaAccessor<0, dim, spacedim>::
2250 operator<(const TriaAccessor<0, dim, spacedim> &other) const
2251 {
2252 Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators());
2253
2254 return (global_vertex_index < other.global_vertex_index);
2255 }
2256
2257
2258
2259 template <int dim, int spacedim>
2260 inline IteratorState::IteratorStates
state()2261 TriaAccessor<0, dim, spacedim>::state() const
2262 {
2263 if (global_vertex_index != numbers::invalid_unsigned_int)
2264 return IteratorState::valid;
2265 else
2266 return IteratorState::past_the_end;
2267 }
2268
2269
2270
2271 template <int dim, int spacedim>
2272 inline int
level()2273 TriaAccessor<0, dim, spacedim>::level()
2274 {
2275 return 0;
2276 }
2277
2278
2279
2280 template <int dim, int spacedim>
2281 inline int
index()2282 TriaAccessor<0, dim, spacedim>::index() const
2283 {
2284 return global_vertex_index;
2285 }
2286
2287
2288
2289 template <int dim, int spacedim>
2290 inline const Triangulation<dim, spacedim> &
get_triangulation()2291 TriaAccessor<0, dim, spacedim>::get_triangulation() const
2292 {
2293 return *tria;
2294 }
2295
2296
2297
2298 template <int dim, int spacedim>
2299 inline void
2300 TriaAccessor<0, dim, spacedim>::operator++()
2301 {
2302 ++global_vertex_index;
2303 if (global_vertex_index >= tria->n_vertices())
2304 global_vertex_index = numbers::invalid_unsigned_int;
2305 }
2306
2307
2308
2309 template <int dim, int spacedim>
2310 inline void
2311 TriaAccessor<0, dim, spacedim>::operator--()
2312 {
2313 if (global_vertex_index != numbers::invalid_unsigned_int)
2314 {
2315 if (global_vertex_index != 0)
2316 --global_vertex_index;
2317 else
2318 global_vertex_index = numbers::invalid_unsigned_int;
2319 }
2320 }
2321
2322
2323
2324 template <int dim, int spacedim>
2325 inline bool
2326 TriaAccessor<0, dim, spacedim>::operator==(const TriaAccessor &t) const
2327 {
2328 const bool result =
2329 ((tria == t.tria) && (global_vertex_index == t.global_vertex_index));
2330
2331 return result;
2332 }
2333
2334
2335
2336 template <int dim, int spacedim>
2337 inline bool
2338 TriaAccessor<0, dim, spacedim>::operator!=(const TriaAccessor &t) const
2339 {
2340 return !(*this == t);
2341 }
2342
2343
2344
2345 template <int dim, int spacedim>
2346 inline unsigned int
vertex_index(const unsigned int)2347 TriaAccessor<0, dim, spacedim>::vertex_index(const unsigned int) const
2348 {
2349 return global_vertex_index;
2350 }
2351
2352
2353
2354 template <int dim, int spacedim>
2355 inline Point<spacedim> &
vertex(const unsigned int)2356 TriaAccessor<0, dim, spacedim>::vertex(const unsigned int) const
2357 {
2358 return const_cast<Point<spacedim> &>(
2359 this->tria->vertices[global_vertex_index]);
2360 }
2361
2362
2363
2364 template <int dim, int spacedim>
2365 inline typename dealii::internal::TriangulationImplementation::
2366 Iterators<dim, spacedim>::line_iterator
line(const unsigned int)2367 TriaAccessor<0, dim, spacedim>::line(const unsigned int)
2368 {
2369 return typename dealii::internal::TriangulationImplementation::
2370 Iterators<dim, spacedim>::line_iterator();
2371 }
2372
2373
2374
2375 template <int dim, int spacedim>
2376 inline unsigned int
line_index(const unsigned int)2377 TriaAccessor<0, dim, spacedim>::line_index(const unsigned int)
2378 {
2379 Assert(false, ExcImpossibleInDim(0));
2380 return numbers::invalid_unsigned_int;
2381 }
2382
2383
2384
2385 template <int dim, int spacedim>
2386 inline typename dealii::internal::TriangulationImplementation::
2387 Iterators<dim, spacedim>::quad_iterator
quad(const unsigned int)2388 TriaAccessor<0, dim, spacedim>::quad(const unsigned int)
2389 {
2390 return typename dealii::internal::TriangulationImplementation::
2391 Iterators<dim, spacedim>::quad_iterator();
2392 }
2393
2394
2395
2396 template <int dim, int spacedim>
2397 inline unsigned int
quad_index(const unsigned int)2398 TriaAccessor<0, dim, spacedim>::quad_index(const unsigned int)
2399 {
2400 Assert(false, ExcImpossibleInDim(0));
2401 return numbers::invalid_unsigned_int;
2402 }
2403
2404
2405
2406 template <int dim, int spacedim>
2407 inline double
diameter()2408 TriaAccessor<0, dim, spacedim>::diameter() const
2409 {
2410 return 0.;
2411 }
2412
2413
2414
2415 template <int dim, int spacedim>
2416 inline double
extent_in_direction(const unsigned int)2417 TriaAccessor<0, dim, spacedim>::extent_in_direction(const unsigned int) const
2418 {
2419 return 0.;
2420 }
2421
2422
2423
2424 template <int dim, int spacedim>
2425 inline Point<spacedim>
center(const bool,const bool)2426 TriaAccessor<0, dim, spacedim>::center(const bool, const bool) const
2427 {
2428 return this->tria->vertices[global_vertex_index];
2429 }
2430
2431
2432
2433 template <int dim, int spacedim>
2434 inline double
measure()2435 TriaAccessor<0, dim, spacedim>::measure() const
2436 {
2437 return 0.;
2438 }
2439
2440
2441
2442 template <int dim, int spacedim>
2443 inline bool
face_orientation(const unsigned int)2444 TriaAccessor<0, dim, spacedim>::face_orientation(const unsigned int /*face*/)
2445 {
2446 return false;
2447 }
2448
2449
2450
2451 template <int dim, int spacedim>
2452 inline bool
face_flip(const unsigned int)2453 TriaAccessor<0, dim, spacedim>::face_flip(const unsigned int /*face*/)
2454 {
2455 return false;
2456 }
2457
2458
2459
2460 template <int dim, int spacedim>
2461 inline bool
face_rotation(const unsigned int)2462 TriaAccessor<0, dim, spacedim>::face_rotation(const unsigned int /*face*/)
2463 {
2464 return false;
2465 }
2466
2467
2468
2469 template <int dim, int spacedim>
2470 inline bool
line_orientation(const unsigned int)2471 TriaAccessor<0, dim, spacedim>::line_orientation(const unsigned int /*line*/)
2472 {
2473 return false;
2474 }
2475
2476
2477
2478 template <int dim, int spacedim>
2479 inline bool
has_children()2480 TriaAccessor<0, dim, spacedim>::has_children()
2481 {
2482 return false;
2483 }
2484
2485
2486
2487 template <int dim, int spacedim>
2488 inline unsigned int
n_children()2489 TriaAccessor<0, dim, spacedim>::n_children()
2490 {
2491 return 0;
2492 }
2493
2494
2495
2496 template <int dim, int spacedim>
2497 inline unsigned int
number_of_children()2498 TriaAccessor<0, dim, spacedim>::number_of_children()
2499 {
2500 return 0;
2501 }
2502
2503
2504
2505 template <int dim, int spacedim>
2506 inline unsigned int
max_refinement_depth()2507 TriaAccessor<0, dim, spacedim>::max_refinement_depth()
2508 {
2509 return 0;
2510 }
2511
2512
2513
2514 template <int dim, int spacedim>
2515 inline unsigned int
child_iterator_to_index(const TriaIterator<TriaAccessor<0,dim,spacedim>> &)2516 TriaAccessor<0, dim, spacedim>::child_iterator_to_index(
2517 const TriaIterator<TriaAccessor<0, dim, spacedim>> &)
2518 {
2519 return numbers::invalid_unsigned_int;
2520 }
2521
2522
2523
2524 template <int dim, int spacedim>
2525 inline TriaIterator<TriaAccessor<0, dim, spacedim>>
child(const unsigned int)2526 TriaAccessor<0, dim, spacedim>::child(const unsigned int)
2527 {
2528 return TriaIterator<TriaAccessor<0, dim, spacedim>>();
2529 }
2530
2531
2532
2533 template <int dim, int spacedim>
2534 inline TriaIterator<TriaAccessor<0, dim, spacedim>>
isotropic_child(const unsigned int)2535 TriaAccessor<0, dim, spacedim>::isotropic_child(const unsigned int)
2536 {
2537 return TriaIterator<TriaAccessor<0, dim, spacedim>>();
2538 }
2539
2540
2541
2542 template <int dim, int spacedim>
2543 inline RefinementCase<0>
refinement_case()2544 TriaAccessor<0, dim, spacedim>::refinement_case()
2545 {
2546 return {RefinementPossibilities<0>::no_refinement};
2547 }
2548
2549
2550
2551 template <int dim, int spacedim>
2552 inline int
child_index(const unsigned int)2553 TriaAccessor<0, dim, spacedim>::child_index(const unsigned int)
2554 {
2555 return -1;
2556 }
2557
2558
2559
2560 template <int dim, int spacedim>
2561 inline int
isotropic_child_index(const unsigned int)2562 TriaAccessor<0, dim, spacedim>::isotropic_child_index(const unsigned int)
2563 {
2564 return -1;
2565 }
2566
2567
2568
2569 template <int dim, int spacedim>
2570 inline bool
used()2571 TriaAccessor<0, dim, spacedim>::used() const
2572 {
2573 return tria->vertex_used(global_vertex_index);
2574 }
2575
2576
2577
2578 /*------------------- Functions: TriaAccessor<0,1,spacedim> -----------------*/
2579
2580 template <int spacedim>
TriaAccessor(const Triangulation<1,spacedim> * tria,const VertexKind vertex_kind,const unsigned int vertex_index)2581 inline TriaAccessor<0, 1, spacedim>::TriaAccessor(
2582 const Triangulation<1, spacedim> *tria,
2583 const VertexKind vertex_kind,
2584 const unsigned int vertex_index)
2585 : tria(tria)
2586 , vertex_kind(vertex_kind)
2587 , global_vertex_index(vertex_index)
2588 {}
2589
2590
2591
2592 template <int spacedim>
TriaAccessor(const Triangulation<1,spacedim> * tria,const int level,const int index,const AccessorData *)2593 inline TriaAccessor<0, 1, spacedim>::TriaAccessor(
2594 const Triangulation<1, spacedim> *tria,
2595 const int level,
2596 const int index,
2597 const AccessorData *)
2598 : tria(tria)
2599 , vertex_kind(interior_vertex)
2600 , global_vertex_index(numbers::invalid_unsigned_int)
2601 {
2602 // in general, calling this constructor should yield an error -- users should
2603 // instead call the one immediately above. however, if you create something
2604 // like Triangulation<1>::face_iterator() then this calls the default
2605 // constructor of the iterator which calls the accessor with argument list
2606 // (0,-2,-2,0), so in this particular case accept this call and create an
2607 // object that corresponds to the default constructed (invalid) vertex
2608 // accessor
2609 (void)level;
2610 (void)index;
2611 Assert((level == -2) && (index == -2), ExcInternalError());
2612 }
2613
2614
2615
2616 template <int spacedim>
2617 template <int structdim2, int dim2, int spacedim2>
TriaAccessor(const TriaAccessor<structdim2,dim2,spacedim2> &)2618 inline TriaAccessor<0, 1, spacedim>::TriaAccessor(
2619 const TriaAccessor<structdim2, dim2, spacedim2> &)
2620 : tria(nullptr)
2621 , vertex_kind(interior_vertex)
2622 , global_vertex_index(numbers::invalid_unsigned_int)
2623 {
2624 Assert(false, ExcImpossibleInDim(0));
2625 }
2626
2627
2628
2629 template <int spacedim>
2630 template <int structdim2, int dim2, int spacedim2>
TriaAccessor(const InvalidAccessor<structdim2,dim2,spacedim2> &)2631 inline TriaAccessor<0, 1, spacedim>::TriaAccessor(
2632 const InvalidAccessor<structdim2, dim2, spacedim2> &)
2633 : tria(nullptr)
2634 , vertex_kind(interior_vertex)
2635 , global_vertex_index(numbers::invalid_unsigned_int)
2636 {
2637 Assert(false, ExcImpossibleInDim(0));
2638 }
2639
2640
2641
2642 template <int spacedim>
2643 inline void
copy_from(const TriaAccessor & t)2644 TriaAccessor<0, 1, spacedim>::copy_from(const TriaAccessor &t)
2645 {
2646 tria = t.tria;
2647 vertex_kind = t.vertex_kind;
2648 global_vertex_index = t.global_vertex_index;
2649 }
2650
2651
2652
2653 template <int spacedim>
2654 inline bool
2655 TriaAccessor<0, 1, spacedim>::
2656 operator<(const TriaAccessor<0, 1, spacedim> &other) const
2657 {
2658 Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators());
2659
2660 return (global_vertex_index < other.global_vertex_index);
2661 }
2662
2663
2664
2665 template <int spacedim>
2666 inline IteratorState::IteratorStates
state()2667 TriaAccessor<0, 1, spacedim>::state()
2668 {
2669 return IteratorState::valid;
2670 }
2671
2672
2673 template <int spacedim>
2674 inline int
level()2675 TriaAccessor<0, 1, spacedim>::level()
2676 {
2677 return 0;
2678 }
2679
2680
2681
2682 template <int spacedim>
2683 inline int
index()2684 TriaAccessor<0, 1, spacedim>::index() const
2685 {
2686 return global_vertex_index;
2687 }
2688
2689
2690
2691 template <int spacedim>
2692 inline const Triangulation<1, spacedim> &
get_triangulation()2693 TriaAccessor<0, 1, spacedim>::get_triangulation() const
2694 {
2695 return *tria;
2696 }
2697
2698
2699
2700 template <int spacedim>
2701 inline void
2702 TriaAccessor<0, 1, spacedim>::operator++() const
2703 {
2704 Assert(false, ExcNotImplemented());
2705 }
2706
2707
2708 template <int spacedim>
2709 inline void
2710 TriaAccessor<0, 1, spacedim>::operator--() const
2711 {
2712 Assert(false, ExcNotImplemented());
2713 }
2714
2715
2716
2717 template <int spacedim>
2718 inline bool
2719 TriaAccessor<0, 1, spacedim>::operator==(const TriaAccessor &t) const
2720 {
2721 const bool result =
2722 ((tria == t.tria) && (global_vertex_index == t.global_vertex_index));
2723 // if we point to the same vertex,
2724 // make sure we know the same about
2725 // it
2726 if (result == true)
2727 Assert(vertex_kind == t.vertex_kind, ExcInternalError());
2728
2729 return result;
2730 }
2731
2732
2733
2734 template <int spacedim>
2735 inline bool
2736 TriaAccessor<0, 1, spacedim>::operator!=(const TriaAccessor &t) const
2737 {
2738 return !(*this == t);
2739 }
2740
2741
2742
2743 template <int spacedim>
2744 inline unsigned int
vertex_index(const unsigned int i)2745 TriaAccessor<0, 1, spacedim>::vertex_index(const unsigned int i) const
2746 {
2747 AssertIndexRange(i, 1);
2748 (void)i;
2749 return global_vertex_index;
2750 }
2751
2752
2753
2754 template <int spacedim>
2755 inline Point<spacedim> &
vertex(const unsigned int i)2756 TriaAccessor<0, 1, spacedim>::vertex(const unsigned int i) const
2757 {
2758 AssertIndexRange(i, 1);
2759 (void)i;
2760 return const_cast<Point<spacedim> &>(
2761 this->tria->vertices[global_vertex_index]);
2762 }
2763
2764
2765
2766 template <int spacedim>
2767 inline Point<spacedim>
center()2768 TriaAccessor<0, 1, spacedim>::center() const
2769 {
2770 return this->tria->vertices[global_vertex_index];
2771 }
2772
2773
2774
2775 template <int spacedim>
2776 inline typename dealii::internal::TriangulationImplementation::
2777 Iterators<1, spacedim>::line_iterator
line(const unsigned int)2778 TriaAccessor<0, 1, spacedim>::line(const unsigned int)
2779 {
2780 return typename dealii::internal::TriangulationImplementation::
2781 Iterators<1, spacedim>::line_iterator();
2782 }
2783
2784
2785 template <int spacedim>
2786 inline unsigned int
line_index(const unsigned int)2787 TriaAccessor<0, 1, spacedim>::line_index(const unsigned int)
2788 {
2789 Assert(false, ExcImpossibleInDim(0));
2790 return numbers::invalid_unsigned_int;
2791 }
2792
2793
2794 template <int spacedim>
2795 inline typename dealii::internal::TriangulationImplementation::
2796 Iterators<1, spacedim>::quad_iterator
quad(const unsigned int)2797 TriaAccessor<0, 1, spacedim>::quad(const unsigned int)
2798 {
2799 return typename dealii::internal::TriangulationImplementation::
2800 Iterators<1, spacedim>::quad_iterator();
2801 }
2802
2803
2804
2805 template <int spacedim>
2806 inline unsigned int
quad_index(const unsigned int)2807 TriaAccessor<0, 1, spacedim>::quad_index(const unsigned int)
2808 {
2809 Assert(false, ExcImpossibleInDim(0));
2810 return numbers::invalid_unsigned_int;
2811 }
2812
2813
2814 template <int spacedim>
2815 inline bool
at_boundary()2816 TriaAccessor<0, 1, spacedim>::at_boundary() const
2817 {
2818 return vertex_kind != interior_vertex;
2819 }
2820
2821
2822 template <int spacedim>
2823 inline types::boundary_id
boundary_id()2824 TriaAccessor<0, 1, spacedim>::boundary_id() const
2825 {
2826 switch (vertex_kind)
2827 {
2828 case left_vertex:
2829 case right_vertex:
2830 {
2831 Assert(tria->vertex_to_boundary_id_map_1d->find(
2832 this->vertex_index()) !=
2833 tria->vertex_to_boundary_id_map_1d->end(),
2834 ExcInternalError());
2835
2836 return (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()];
2837 }
2838
2839 default:
2840 return numbers::internal_face_boundary_id;
2841 }
2842 }
2843
2844
2845
2846 template <int spacedim>
2847 inline const Manifold<1, spacedim> &
get_manifold()2848 TriaAccessor<0, 1, spacedim>::get_manifold() const
2849 {
2850 return this->tria->get_manifold(this->manifold_id());
2851 }
2852
2853
2854
2855 template <int spacedim>
2856 inline types::manifold_id
manifold_id()2857 TriaAccessor<0, 1, spacedim>::manifold_id() const
2858 {
2859 if (tria->vertex_to_manifold_id_map_1d->find(this->vertex_index()) !=
2860 tria->vertex_to_manifold_id_map_1d->end())
2861 return (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()];
2862 else
2863 return numbers::flat_manifold_id;
2864 }
2865
2866
2867 template <int spacedim>
2868 inline bool
face_orientation(const unsigned int)2869 TriaAccessor<0, 1, spacedim>::face_orientation(const unsigned int /*face*/)
2870 {
2871 return false;
2872 }
2873
2874
2875
2876 template <int spacedim>
2877 inline bool
face_flip(const unsigned int)2878 TriaAccessor<0, 1, spacedim>::face_flip(const unsigned int /*face*/)
2879 {
2880 return false;
2881 }
2882
2883
2884
2885 template <int spacedim>
2886 inline bool
face_rotation(const unsigned int)2887 TriaAccessor<0, 1, spacedim>::face_rotation(const unsigned int /*face*/)
2888 {
2889 return false;
2890 }
2891
2892
2893
2894 template <int spacedim>
2895 inline bool
line_orientation(const unsigned int)2896 TriaAccessor<0, 1, spacedim>::line_orientation(const unsigned int /*line*/)
2897 {
2898 return false;
2899 }
2900
2901
2902
2903 template <int spacedim>
2904 inline bool
has_children()2905 TriaAccessor<0, 1, spacedim>::has_children()
2906 {
2907 return false;
2908 }
2909
2910
2911
2912 template <int spacedim>
2913 inline unsigned int
n_children()2914 TriaAccessor<0, 1, spacedim>::n_children()
2915 {
2916 return 0;
2917 }
2918
2919
2920
2921 template <int spacedim>
2922 inline unsigned int
number_of_children()2923 TriaAccessor<0, 1, spacedim>::number_of_children()
2924 {
2925 return 0;
2926 }
2927
2928
2929
2930 template <int spacedim>
2931 inline unsigned int
max_refinement_depth()2932 TriaAccessor<0, 1, spacedim>::max_refinement_depth()
2933 {
2934 return 0;
2935 }
2936
2937
2938
2939 template <int spacedim>
2940 inline unsigned int
child_iterator_to_index(const TriaIterator<TriaAccessor<0,1,spacedim>> &)2941 TriaAccessor<0, 1, spacedim>::child_iterator_to_index(
2942 const TriaIterator<TriaAccessor<0, 1, spacedim>> &)
2943 {
2944 return numbers::invalid_unsigned_int;
2945 }
2946
2947
2948
2949 template <int spacedim>
2950 inline TriaIterator<TriaAccessor<0, 1, spacedim>>
child(const unsigned int)2951 TriaAccessor<0, 1, spacedim>::child(const unsigned int)
2952 {
2953 return TriaIterator<TriaAccessor<0, 1, spacedim>>();
2954 }
2955
2956
2957 template <int spacedim>
2958 inline TriaIterator<TriaAccessor<0, 1, spacedim>>
isotropic_child(const unsigned int)2959 TriaAccessor<0, 1, spacedim>::isotropic_child(const unsigned int)
2960 {
2961 return TriaIterator<TriaAccessor<0, 1, spacedim>>();
2962 }
2963
2964
2965 template <int spacedim>
2966 inline RefinementCase<0>
refinement_case()2967 TriaAccessor<0, 1, spacedim>::refinement_case()
2968 {
2969 return {RefinementPossibilities<0>::no_refinement};
2970 }
2971
2972 template <int spacedim>
2973 inline int
child_index(const unsigned int)2974 TriaAccessor<0, 1, spacedim>::child_index(const unsigned int)
2975 {
2976 return -1;
2977 }
2978
2979
2980 template <int spacedim>
2981 inline int
isotropic_child_index(const unsigned int)2982 TriaAccessor<0, 1, spacedim>::isotropic_child_index(const unsigned int)
2983 {
2984 return -1;
2985 }
2986
2987
2988
2989 template <int spacedim>
2990 inline void
set_boundary_id(const types::boundary_id b)2991 TriaAccessor<0, 1, spacedim>::set_boundary_id(const types::boundary_id b)
2992 {
2993 Assert(tria->vertex_to_boundary_id_map_1d->find(this->vertex_index()) !=
2994 tria->vertex_to_boundary_id_map_1d->end(),
2995 ExcInternalError());
2996
2997 (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()] = b;
2998 }
2999
3000
3001
3002 template <int spacedim>
3003 inline void
set_manifold_id(const types::manifold_id b)3004 TriaAccessor<0, 1, spacedim>::set_manifold_id(const types::manifold_id b)
3005 {
3006 (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()] = b;
3007 }
3008
3009
3010
3011 template <int spacedim>
3012 inline void
set_all_boundary_ids(const types::boundary_id b)3013 TriaAccessor<0, 1, spacedim>::set_all_boundary_ids(const types::boundary_id b)
3014 {
3015 set_boundary_id(b);
3016 }
3017
3018
3019
3020 template <int spacedim>
3021 inline void
set_all_manifold_ids(const types::manifold_id b)3022 TriaAccessor<0, 1, spacedim>::set_all_manifold_ids(const types::manifold_id b)
3023 {
3024 set_manifold_id(b);
3025 }
3026
3027
3028
3029 template <int spacedim>
3030 inline bool
used()3031 TriaAccessor<0, 1, spacedim>::used() const
3032 {
3033 return tria->vertex_used(global_vertex_index);
3034 }
3035
3036
3037
3038 template <int spacedim>
3039 inline ReferenceCell::Type
reference_cell_type()3040 TriaAccessor<0, 1, spacedim>::reference_cell_type() const
3041 {
3042 return ReferenceCell::Type::Vertex;
3043 }
3044
3045
3046
3047 template <int spacedim>
3048 unsigned int
n_vertices()3049 TriaAccessor<0, 1, spacedim>::n_vertices() const
3050 {
3051 return 1;
3052 }
3053
3054
3055
3056 template <int spacedim>
3057 unsigned int
n_lines()3058 TriaAccessor<0, 1, spacedim>::n_lines() const
3059 {
3060 return 0;
3061 }
3062
3063
3064
3065 template <int spacedim>
3066 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
vertex_indices()3067 TriaAccessor<0, 1, spacedim>::vertex_indices() const
3068 {
3069 return {0U, n_vertices()};
3070 }
3071
3072
3073
3074 template <int spacedim>
3075 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
line_indices()3076 TriaAccessor<0, 1, spacedim>::line_indices() const
3077 {
3078 return {0U, n_lines()};
3079 }
3080
3081 /*------------------ Functions: CellAccessor<dim,spacedim> ------------------*/
3082
3083
3084 template <int dim, int spacedim>
CellAccessor(const Triangulation<dim,spacedim> * parent,const int level,const int index,const AccessorData * local_data)3085 inline CellAccessor<dim, spacedim>::CellAccessor(
3086 const Triangulation<dim, spacedim> *parent,
3087 const int level,
3088 const int index,
3089 const AccessorData * local_data)
3090 : TriaAccessor<dim, dim, spacedim>(parent, level, index, local_data)
3091 {}
3092
3093
3094
3095 template <int dim, int spacedim>
CellAccessor(const TriaAccessor<dim,dim,spacedim> & cell_accessor)3096 inline CellAccessor<dim, spacedim>::CellAccessor(
3097 const TriaAccessor<dim, dim, spacedim> &cell_accessor)
3098 : TriaAccessor<dim, dim, spacedim>(
3099 static_cast<const TriaAccessor<dim, dim, spacedim> &>(cell_accessor))
3100 {}
3101
3102
3103
3104 namespace internal
3105 {
3106 namespace CellAccessorImplementation
3107 {
3108 template <int spacedim>
3109 inline dealii::TriaIterator<dealii::TriaAccessor<0, 1, spacedim>>
get_face(const dealii::CellAccessor<1,spacedim> & cell,const unsigned int i)3110 get_face(const dealii::CellAccessor<1, spacedim> &cell,
3111 const unsigned int i)
3112 {
3113 dealii::TriaAccessor<0, 1, spacedim> a(
3114 &cell.get_triangulation(),
3115 ((i == 0) && cell.at_boundary(0) ?
3116 dealii::TriaAccessor<0, 1, spacedim>::left_vertex :
3117 ((i == 1) && cell.at_boundary(1) ?
3118 dealii::TriaAccessor<0, 1, spacedim>::right_vertex :
3119 dealii::TriaAccessor<0, 1, spacedim>::interior_vertex)),
3120 cell.vertex_index(i));
3121 return dealii::TriaIterator<dealii::TriaAccessor<0, 1, spacedim>>(a);
3122 }
3123
3124
3125 template <int spacedim>
3126 inline dealii::TriaIterator<dealii::TriaAccessor<1, 2, spacedim>>
get_face(const dealii::CellAccessor<2,spacedim> & cell,const unsigned int i)3127 get_face(const dealii::CellAccessor<2, spacedim> &cell,
3128 const unsigned int i)
3129 {
3130 return cell.line(i);
3131 }
3132
3133
3134 template <int spacedim>
3135 inline dealii::TriaIterator<dealii::TriaAccessor<2, 3, spacedim>>
get_face(const dealii::CellAccessor<3,spacedim> & cell,const unsigned int i)3136 get_face(const dealii::CellAccessor<3, spacedim> &cell,
3137 const unsigned int i)
3138 {
3139 return cell.quad(i);
3140 }
3141 } // namespace CellAccessorImplementation
3142 } // namespace internal
3143
3144
3145
3146 template <int dim, int spacedim>
3147 inline TriaIterator<CellAccessor<dim, spacedim>>
child(const unsigned int i)3148 CellAccessor<dim, spacedim>::child(const unsigned int i) const
3149 {
3150 TriaIterator<CellAccessor<dim, spacedim>> q(this->tria,
3151 this->present_level + 1,
3152 this->child_index(i));
3153
3154 Assert((q.state() == IteratorState::past_the_end) || q->used(),
3155 ExcInternalError());
3156
3157 return q;
3158 }
3159
3160
3161
3162 template <int dim, int spacedim>
3163 inline TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
face(const unsigned int i)3164 CellAccessor<dim, spacedim>::face(const unsigned int i) const
3165 {
3166 return dealii::internal::CellAccessorImplementation::get_face(*this, i);
3167 }
3168
3169
3170
3171 template <int dim, int spacedim>
3172 inline unsigned int
face_iterator_to_index(const TriaIterator<TriaAccessor<dim-1,dim,spacedim>> & face)3173 CellAccessor<dim, spacedim>::face_iterator_to_index(
3174 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> &face) const
3175 {
3176 for (const unsigned int face_n : this->face_indices())
3177 if (this->face(face_n) == face)
3178 return face_n;
3179
3180 Assert(false,
3181 ExcMessage("The given face is not a face of the current cell."));
3182 return numbers::invalid_unsigned_int;
3183 }
3184
3185
3186
3187 template <int dim, int spacedim>
3188 inline boost::container::small_vector<
3189 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
3190 GeometryInfo<dim>::faces_per_cell>
face_iterators()3191 CellAccessor<dim, spacedim>::face_iterators() const
3192 {
3193 boost::container::small_vector<
3194 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
3195 GeometryInfo<dim>::faces_per_cell>
3196 face_iterators(this->n_faces());
3197
3198 for (unsigned int i : this->face_indices())
3199 face_iterators[i] =
3200 dealii::internal::CellAccessorImplementation::get_face(*this, i);
3201
3202 return face_iterators;
3203 }
3204
3205
3206
3207 template <int dim, int spacedim>
3208 inline unsigned int
face_index(const unsigned int i)3209 CellAccessor<dim, spacedim>::face_index(const unsigned int i) const
3210 {
3211 switch (dim)
3212 {
3213 case 1:
3214 {
3215 return this->vertex_index(i);
3216 }
3217
3218 case 2:
3219 return this->line_index(i);
3220
3221 case 3:
3222 return this->quad_index(i);
3223
3224 default:
3225 return numbers::invalid_unsigned_int;
3226 }
3227 }
3228
3229
3230
3231 template <int dim, int spacedim>
3232 inline int
neighbor_index(const unsigned int face_no)3233 CellAccessor<dim, spacedim>::neighbor_index(const unsigned int face_no) const
3234 {
3235 AssertIndexRange(face_no, this->n_faces());
3236 return this->tria->levels[this->present_level]
3237 ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell +
3238 face_no]
3239 .second;
3240 }
3241
3242
3243
3244 template <int dim, int spacedim>
3245 inline int
neighbor_level(const unsigned int face_no)3246 CellAccessor<dim, spacedim>::neighbor_level(const unsigned int face_no) const
3247 {
3248 AssertIndexRange(face_no, this->n_faces());
3249 return this->tria->levels[this->present_level]
3250 ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell +
3251 face_no]
3252 .first;
3253 }
3254
3255
3256
3257 template <int dim, int spacedim>
3258 inline RefinementCase<dim>
refine_flag_set()3259 CellAccessor<dim, spacedim>::refine_flag_set() const
3260 {
3261 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
3262 // cells flagged for refinement must be active
3263 // (the @p set_refine_flag function checks this,
3264 // but activity may change when refinement is
3265 // executed and for some reason the refine
3266 // flag is not cleared).
3267 Assert(this->is_active() || !this->tria->levels[this->present_level]
3268 ->refine_flags[this->present_index],
3269 ExcRefineCellNotActive());
3270 return RefinementCase<dim>(
3271 this->tria->levels[this->present_level]->refine_flags[this->present_index]);
3272 }
3273
3274
3275
3276 template <int dim, int spacedim>
3277 inline void
set_refine_flag(const RefinementCase<dim> refinement_case)3278 CellAccessor<dim, spacedim>::set_refine_flag(
3279 const RefinementCase<dim> refinement_case) const
3280 {
3281 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
3282 Assert(!coarsen_flag_set(), ExcCellFlaggedForCoarsening());
3283
3284 this->tria->levels[this->present_level]->refine_flags[this->present_index] =
3285 refinement_case;
3286 }
3287
3288
3289
3290 template <int dim, int spacedim>
3291 inline void
clear_refine_flag()3292 CellAccessor<dim, spacedim>::clear_refine_flag() const
3293 {
3294 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
3295 this->tria->levels[this->present_level]->refine_flags[this->present_index] =
3296 RefinementCase<dim>::no_refinement;
3297 }
3298
3299
3300
3301 template <int dim, int spacedim>
3302 inline bool
flag_for_face_refinement(const unsigned int face_no,const RefinementCase<dim-1> & face_refinement_case)3303 CellAccessor<dim, spacedim>::flag_for_face_refinement(
3304 const unsigned int face_no,
3305 const RefinementCase<dim - 1> &face_refinement_case) const
3306 {
3307 Assert(dim > 1, ExcImpossibleInDim(dim));
3308 AssertIndexRange(face_no, this->n_faces());
3309 AssertIndexRange(face_refinement_case,
3310 RefinementCase<dim>::isotropic_refinement + 1);
3311
3312 // the new refinement case is a combination
3313 // of the minimum required one for the given
3314 // face refinement and the already existing
3315 // flagged refinement case
3316 RefinementCase<dim> old_ref_case = refine_flag_set();
3317 RefinementCase<dim> new_ref_case =
3318 (old_ref_case |
3319 GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
3320 face_refinement_case,
3321 face_no,
3322 this->face_orientation(face_no),
3323 this->face_flip(face_no),
3324 this->face_rotation(face_no)));
3325 set_refine_flag(new_ref_case);
3326 // return, whether we had to change the
3327 // refinement flag
3328 return new_ref_case != old_ref_case;
3329 }
3330
3331
3332
3333 template <int dim, int spacedim>
3334 inline bool
flag_for_line_refinement(const unsigned int line_no)3335 CellAccessor<dim, spacedim>::flag_for_line_refinement(
3336 const unsigned int line_no) const
3337 {
3338 Assert(dim > 1, ExcImpossibleInDim(dim));
3339 AssertIndexRange(line_no, this->n_lines());
3340
3341 // the new refinement case is a combination
3342 // of the minimum required one for the given
3343 // line refinement and the already existing
3344 // flagged refinement case
3345 RefinementCase<dim>
3346 old_ref_case = refine_flag_set(),
3347 new_ref_case =
3348 old_ref_case |
3349 GeometryInfo<dim>::min_cell_refinement_case_for_line_refinement(line_no);
3350 set_refine_flag(new_ref_case);
3351 // return, whether we had to change the
3352 // refinement flag
3353 return new_ref_case != old_ref_case;
3354 }
3355
3356
3357
3358 template <>
3359 inline dealii::internal::SubfaceCase<1>
subface_case(const unsigned int)3360 CellAccessor<1>::subface_case(const unsigned int) const
3361 {
3362 return dealii::internal::SubfaceCase<1>::case_none;
3363 }
3364
3365 template <>
3366 inline dealii::internal::SubfaceCase<1>
subface_case(const unsigned int)3367 CellAccessor<1, 2>::subface_case(const unsigned int) const
3368 {
3369 return dealii::internal::SubfaceCase<1>::case_none;
3370 }
3371
3372
3373 template <>
3374 inline dealii::internal::SubfaceCase<1>
subface_case(const unsigned int)3375 CellAccessor<1, 3>::subface_case(const unsigned int) const
3376 {
3377 return dealii::internal::SubfaceCase<1>::case_none;
3378 }
3379
3380
3381 template <>
3382 inline dealii::internal::SubfaceCase<2>
subface_case(const unsigned int face_no)3383 CellAccessor<2>::subface_case(const unsigned int face_no) const
3384 {
3385 Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive());
3386 AssertIndexRange(face_no, this->n_faces());
3387 return ((face(face_no)->has_children()) ?
3388 dealii::internal::SubfaceCase<2>::case_x :
3389 dealii::internal::SubfaceCase<2>::case_none);
3390 }
3391
3392 template <>
3393 inline dealii::internal::SubfaceCase<2>
subface_case(const unsigned int face_no)3394 CellAccessor<2, 3>::subface_case(const unsigned int face_no) const
3395 {
3396 Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive());
3397 AssertIndexRange(face_no, this->n_faces());
3398 return ((face(face_no)->has_children()) ?
3399 dealii::internal::SubfaceCase<2>::case_x :
3400 dealii::internal::SubfaceCase<2>::case_none);
3401 }
3402
3403
3404 template <>
3405 inline dealii::internal::SubfaceCase<3>
subface_case(const unsigned int face_no)3406 CellAccessor<3>::subface_case(const unsigned int face_no) const
3407 {
3408 Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive());
3409 AssertIndexRange(face_no, this->n_faces());
3410 switch (static_cast<std::uint8_t>(face(face_no)->refinement_case()))
3411 {
3412 case RefinementCase<3>::no_refinement:
3413 return dealii::internal::SubfaceCase<3>::case_none;
3414 case RefinementCase<3>::cut_x:
3415 if (face(face_no)->child(0)->has_children())
3416 {
3417 Assert(face(face_no)->child(0)->refinement_case() ==
3418 RefinementCase<2>::cut_y,
3419 ExcInternalError());
3420 if (face(face_no)->child(1)->has_children())
3421 {
3422 Assert(face(face_no)->child(1)->refinement_case() ==
3423 RefinementCase<2>::cut_y,
3424 ExcInternalError());
3425 return dealii::internal::SubfaceCase<3>::case_x1y2y;
3426 }
3427 else
3428 return dealii::internal::SubfaceCase<3>::case_x1y;
3429 }
3430 else
3431 {
3432 if (face(face_no)->child(1)->has_children())
3433 {
3434 Assert(face(face_no)->child(1)->refinement_case() ==
3435 RefinementCase<2>::cut_y,
3436 ExcInternalError());
3437 return dealii::internal::SubfaceCase<3>::case_x2y;
3438 }
3439 else
3440 return dealii::internal::SubfaceCase<3>::case_x;
3441 }
3442 case RefinementCase<3>::cut_y:
3443 if (face(face_no)->child(0)->has_children())
3444 {
3445 Assert(face(face_no)->child(0)->refinement_case() ==
3446 RefinementCase<2>::cut_x,
3447 ExcInternalError());
3448 if (face(face_no)->child(1)->has_children())
3449 {
3450 Assert(face(face_no)->child(1)->refinement_case() ==
3451 RefinementCase<2>::cut_x,
3452 ExcInternalError());
3453 return dealii::internal::SubfaceCase<3>::case_y1x2x;
3454 }
3455 else
3456 return dealii::internal::SubfaceCase<3>::case_y1x;
3457 }
3458 else
3459 {
3460 if (face(face_no)->child(1)->has_children())
3461 {
3462 Assert(face(face_no)->child(1)->refinement_case() ==
3463 RefinementCase<2>::cut_x,
3464 ExcInternalError());
3465 return dealii::internal::SubfaceCase<3>::case_y2x;
3466 }
3467 else
3468 return dealii::internal::SubfaceCase<3>::case_y;
3469 }
3470 case RefinementCase<3>::cut_xy:
3471 return dealii::internal::SubfaceCase<3>::case_xy;
3472 default:
3473 Assert(false, ExcInternalError());
3474 }
3475 // we should never get here
3476 return dealii::internal::SubfaceCase<3>::case_none;
3477 }
3478
3479
3480
3481 template <int dim, int spacedim>
3482 inline bool
coarsen_flag_set()3483 CellAccessor<dim, spacedim>::coarsen_flag_set() const
3484 {
3485 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
3486 // cells flagged for coarsening must be active
3487 // (the @p set_refine_flag function checks this,
3488 // but activity may change when refinement is
3489 // executed and for some reason the refine
3490 // flag is not cleared).
3491 Assert(this->is_active() || !this->tria->levels[this->present_level]
3492 ->coarsen_flags[this->present_index],
3493 ExcRefineCellNotActive());
3494 return this->tria->levels[this->present_level]
3495 ->coarsen_flags[this->present_index];
3496 }
3497
3498
3499
3500 template <int dim, int spacedim>
3501 inline void
set_coarsen_flag()3502 CellAccessor<dim, spacedim>::set_coarsen_flag() const
3503 {
3504 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
3505 Assert(!refine_flag_set(), ExcCellFlaggedForRefinement());
3506
3507 this->tria->levels[this->present_level]->coarsen_flags[this->present_index] =
3508 true;
3509 }
3510
3511
3512
3513 template <int dim, int spacedim>
3514 inline void
clear_coarsen_flag()3515 CellAccessor<dim, spacedim>::clear_coarsen_flag() const
3516 {
3517 Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
3518 this->tria->levels[this->present_level]->coarsen_flags[this->present_index] =
3519 false;
3520 }
3521
3522
3523
3524 template <int dim, int spacedim>
3525 inline TriaIterator<CellAccessor<dim, spacedim>>
neighbor(const unsigned int face_no)3526 CellAccessor<dim, spacedim>::neighbor(const unsigned int face_no) const
3527 {
3528 TriaIterator<CellAccessor<dim, spacedim>> q(this->tria,
3529 neighbor_level(face_no),
3530 neighbor_index(face_no));
3531
3532 Assert((q.state() == IteratorState::past_the_end) || q->used(),
3533 ExcInternalError());
3534
3535 return q;
3536 }
3537
3538
3539
3540 template <int dim, int spacedim>
3541 inline bool
active()3542 CellAccessor<dim, spacedim>::active() const
3543 {
3544 return !this->has_children();
3545 }
3546
3547
3548
3549 template <int dim, int spacedim>
3550 inline bool
is_active()3551 CellAccessor<dim, spacedim>::is_active() const
3552 {
3553 return !this->has_children();
3554 }
3555
3556
3557
3558 template <int dim, int spacedim>
3559 inline bool
is_locally_owned()3560 CellAccessor<dim, spacedim>::is_locally_owned() const
3561 {
3562 Assert(this->is_active(),
3563 ExcMessage("is_locally_owned() can only be called on active cells!"));
3564 #ifndef DEAL_II_WITH_MPI
3565 return true;
3566 #else
3567
3568 // Serial triangulations report invalid_subdomain_id as their locally owned
3569 // subdomain, so the first condition checks whether we have a serial
3570 // triangulation, in which case all cells are locally owned. The second
3571 // condition compares the subdomain id in the parallel case.
3572 return (this->tria->locally_owned_subdomain() ==
3573 numbers::invalid_subdomain_id ||
3574 this->subdomain_id() == this->tria->locally_owned_subdomain());
3575
3576 #endif
3577 }
3578
3579
3580 template <int dim, int spacedim>
3581 inline bool
is_locally_owned_on_level()3582 CellAccessor<dim, spacedim>::is_locally_owned_on_level() const
3583 {
3584 #ifndef DEAL_II_WITH_MPI
3585 return true;
3586 #else
3587
3588 // Serial triangulations report invalid_subdomain_id as their locally owned
3589 // subdomain, so the first condition checks whether we have a serial
3590 // triangulation, in which case all cells are locally owned. The second
3591 // condition compares the subdomain id in the parallel case.
3592 return (this->tria->locally_owned_subdomain() ==
3593 numbers::invalid_subdomain_id ||
3594 this->level_subdomain_id() == this->tria->locally_owned_subdomain());
3595
3596 #endif
3597 }
3598
3599
3600 template <int dim, int spacedim>
3601 inline bool
is_ghost()3602 CellAccessor<dim, spacedim>::is_ghost() const
3603 {
3604 Assert(this->is_active(),
3605 ExcMessage("is_ghost() can only be called on active cells!"));
3606 if (this->has_children())
3607 return false;
3608
3609 #ifndef DEAL_II_WITH_MPI
3610 return false;
3611 #else
3612
3613 // Serial triangulations report invalid_subdomain_id as their locally owned
3614 // subdomain, so the first condition rules out that case as all cells to a
3615 // serial triangulation are locally owned and none is ghosted. The second
3616 // and third conditions check whether the cell's subdomain is not the
3617 // locally owned one and not artificial.
3618 return (this->tria->locally_owned_subdomain() !=
3619 numbers::invalid_subdomain_id &&
3620 this->subdomain_id() != this->tria->locally_owned_subdomain() &&
3621 this->subdomain_id() != numbers::artificial_subdomain_id);
3622
3623 #endif
3624 }
3625
3626
3627
3628 template <int dim, int spacedim>
3629 inline bool
is_artificial()3630 CellAccessor<dim, spacedim>::is_artificial() const
3631 {
3632 Assert(this->is_active(),
3633 ExcMessage("is_artificial() can only be called on active cells!"));
3634 #ifndef DEAL_II_WITH_MPI
3635 return false;
3636 #else
3637
3638 // Serial triangulations report invalid_subdomain_id as their locally owned
3639 // subdomain, so the first condition rules out that case as all cells to a
3640 // serial triangulation are locally owned and none is artificial.
3641 return (this->tria->locally_owned_subdomain() !=
3642 numbers::invalid_subdomain_id &&
3643 this->subdomain_id() == numbers::artificial_subdomain_id);
3644
3645 #endif
3646 }
3647
3648
3649
3650 template <int dim, int spacedim>
3651 inline types::subdomain_id
subdomain_id()3652 CellAccessor<dim, spacedim>::subdomain_id() const
3653 {
3654 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
3655 Assert(this->is_active(),
3656 ExcMessage("subdomain_id() can only be called on active cells!"));
3657 return this->tria->levels[this->present_level]
3658 ->subdomain_ids[this->present_index];
3659 }
3660
3661
3662
3663 template <int dim, int spacedim>
3664 inline unsigned int
neighbor_face_no(const unsigned int neighbor)3665 CellAccessor<dim, spacedim>::neighbor_face_no(const unsigned int neighbor) const
3666 {
3667 const unsigned int n2 = neighbor_of_neighbor_internal(neighbor);
3668 if (n2 != numbers::invalid_unsigned_int)
3669 // return this value as the
3670 // neighbor is not coarser
3671 return n2;
3672 else
3673 // the neighbor is coarser
3674 return neighbor_of_coarser_neighbor(neighbor).first;
3675 }
3676
3677
3678
3679 template <int dim, int spacedim>
3680 inline bool
is_level_cell()3681 CellAccessor<dim, spacedim>::is_level_cell()
3682 {
3683 return false;
3684 }
3685
3686
3687 DEAL_II_NAMESPACE_CLOSE
3688
3689 #endif
3690