From 0bce841203ab2cd9051f4d77a4a21a081bc1ec44 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 17 Dec 2024 18:00:05 -0700 Subject: [PATCH] Just include a file directly. --- include/deal.II/grid/tria_accessor.h | 4011 +++++++++++++++- .../deal.II/grid/tria_accessor.templates.h | 4024 +---------------- source/grid/tria_accessor.cc | 1 - 3 files changed, 4011 insertions(+), 4025 deletions(-) diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index d69ebf5c44..39586e1489 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -22,14 +22,22 @@ #include #include #include +#include #include #include +#include +#include +#include +#include #include #include +#include #include +#include +#include #include @@ -4338,9 +4346,4008 @@ TriaAccessor<3, 3, 3>::set_all_manifold_ids(const types::manifold_id) const; #endif // DOXYGEN + + +namespace internal +{ + namespace TriaAccessorImplementation + { + /** + * Compute the diameter for a given set of vertices. The vertices are + * interpreted, depending on their count, as the vertices of a particular + * reference cell. + */ + template + inline double + diameter( + const boost::container::small_vector, + GeometryInfo::vertices_per_cell> + vertices) + { + const ReferenceCell reference_cell = + ReferenceCell::n_vertices_to_type(dim, vertices.size()); + + if (reference_cell == ReferenceCells::Line) + // Return the distance between the two vertices + return (vertices[1] - vertices[0]).norm(); + else if (reference_cell == ReferenceCells::Triangle) + // Return the longest of the three edges + return std::max({(vertices[1] - vertices[0]).norm(), + (vertices[2] - vertices[1]).norm(), + (vertices[2] - vertices[0]).norm()}); + else if (reference_cell == ReferenceCells::Quadrilateral) + // Return the longer one of the two diagonals of the quadrilateral + return std::max({(vertices[3] - vertices[0]).norm(), + (vertices[2] - vertices[1]).norm()}); + else if (reference_cell == ReferenceCells::Tetrahedron) + // Return the longest of the six edges of the tetrahedron + return std::max({(vertices[1] - vertices[0]).norm(), + (vertices[2] - vertices[0]).norm(), + (vertices[2] - vertices[1]).norm(), + (vertices[3] - vertices[0]).norm(), + (vertices[3] - vertices[1]).norm(), + (vertices[3] - vertices[2]).norm()}); + else if (reference_cell == ReferenceCells::Pyramid) + // Return ... + return std::max({// the longest diagonal of the quadrilateral base + // of the pyramid or ... + (vertices[3] - vertices[0]).norm(), + (vertices[2] - vertices[1]).norm(), + // the longest edge connected with the apex of the + // pyramid + (vertices[4] - vertices[0]).norm(), + (vertices[4] - vertices[1]).norm(), + (vertices[4] - vertices[2]).norm(), + (vertices[4] - vertices[3]).norm()}); + else if (reference_cell == ReferenceCells::Wedge) + // Return ... + return std::max({// the longest of the 2*3=6 diagonals of the three + // quadrilateral sides of the wedge or ... + (vertices[4] - vertices[0]).norm(), + (vertices[3] - vertices[1]).norm(), + (vertices[5] - vertices[1]).norm(), + (vertices[4] - vertices[2]).norm(), + (vertices[5] - vertices[0]).norm(), + (vertices[3] - vertices[2]).norm(), + // the longest of the 3*2=6 edges of the two + // triangular faces of the wedge + (vertices[1] - vertices[0]).norm(), + (vertices[2] - vertices[1]).norm(), + (vertices[2] - vertices[0]).norm(), + (vertices[4] - vertices[3]).norm(), + (vertices[5] - vertices[4]).norm(), + (vertices[5] - vertices[3]).norm()}); + else if (reference_cell == ReferenceCells::Hexahedron) + // Return the longest of the four diagonals of the hexahedron + return std::max({(vertices[7] - vertices[0]).norm(), + (vertices[6] - vertices[1]).norm(), + (vertices[2] - vertices[5]).norm(), + (vertices[3] - vertices[4]).norm()}); + + DEAL_II_NOT_IMPLEMENTED(); + return -1e10; + } + } // namespace TriaAccessorImplementation +} // namespace internal + + +/*--------------------- Functions: TriaAccessorBase -------------------------*/ + +template +inline TriaAccessorBase::TriaAccessorBase( + const Triangulation *tria, + const int level, + const int index, + const AccessorData *) + : present_level((structdim == dim) ? level : 0) + , present_index(index) + , tria(tria) +{ + // non-cells have no level, so a 0 + // should have been passed, or a -1 + // for an end-iterator, or -2 for + // an invalid (default constructed) + // iterator + if (structdim != dim) + { + Assert((level == 0) || (level == -1) || (level == -2), + ExcInternalError()); + } +} + + +template +inline TriaAccessorBase::TriaAccessorBase( + const TriaAccessorBase &a) + : present_level(a.present_level) + , present_index(a.present_index) + , tria(a.tria) +{} + + +template +inline void +TriaAccessorBase::copy_from( + const TriaAccessorBase &a) +{ + present_level = a.present_level; + present_index = a.present_index; + tria = a.tria; + + if (structdim != dim) + { + Assert((present_level == 0) || (present_level == -1) || + (present_level == -2), + ExcInternalError()); + } +} + + + +template +inline TriaAccessorBase & +TriaAccessorBase::operator=( + const TriaAccessorBase &a) +{ + present_level = a.present_level; + present_index = a.present_index; + tria = a.tria; + + if (structdim != dim) + { + Assert((present_level == 0) || (present_level == -1) || + (present_level == -2), + ExcInternalError()); + } + return *this; +} + + + +template +inline bool +TriaAccessorBase::operator==( + const TriaAccessorBase &a) const +{ + Assert(tria == a.tria || tria == nullptr || a.tria == nullptr, + TriaAccessorExceptions::ExcCantCompareIterators()); + return ((tria == a.tria) && (present_level == a.present_level) && + (present_index == a.present_index)); +} + + + +template +inline bool +TriaAccessorBase::operator!=( + const TriaAccessorBase &a) const +{ + Assert(tria == a.tria || tria == nullptr || a.tria == nullptr, + TriaAccessorExceptions::ExcCantCompareIterators()); + return ((tria != a.tria) || (present_level != a.present_level) || + (present_index != a.present_index)); +} + + + +template +inline bool +TriaAccessorBase::operator<( + const TriaAccessorBase &other) const +{ + Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators()); + + if (present_level != other.present_level) + return (present_level < other.present_level); + + return (present_index < other.present_index); +} + + + +template +inline int +TriaAccessorBase::level() const +{ + // This is always zero or invalid + // if the object is not a cell + return present_level; +} + + + +template +inline int +TriaAccessorBase::index() const +{ + return present_index; +} + + + +template +inline IteratorState::IteratorStates +TriaAccessorBase::state() const +{ + if ((present_level >= 0) && (present_index >= 0)) + return IteratorState::valid; + else if (present_index == -1) + return IteratorState::past_the_end; + else + return IteratorState::invalid; +} + + + +template +inline const Triangulation & +TriaAccessorBase::get_triangulation() const +{ + return *tria; +} + + + +template +inline void +TriaAccessorBase::operator++() +{ + // this iterator is used for + // objects without level + ++this->present_index; + + if (structdim != dim) + { + // is index still in the range of + // the vector? (note that we don't + // have to set the level, since + // dim!=1 and the object therefore + // has no level) + if (this->present_index >= static_cast(objects().n_objects())) + this->present_index = -1; + } + else + { + while (this->present_index >= + static_cast( + this->tria->levels[this->present_level]->cells.n_objects())) + { + // no -> go one level up until we find + // one with more than zero cells + ++this->present_level; + this->present_index = 0; + // highest level reached? + if (this->present_level >= + static_cast(this->tria->levels.size())) + { + // return with past the end pointer + this->present_level = this->present_index = -1; + return; + } + } + } +} + + +template +inline void +TriaAccessorBase::operator--() +{ + // same as operator++ + --this->present_index; + + if (structdim != dim) + { + if (this->present_index < 0) + this->present_index = -1; + } + else + { + while (this->present_index < 0) + { + // no -> go one level down + --this->present_level; + // lowest level reached? + if (this->present_level == -1) + { + // return with past the end pointer + this->present_level = this->present_index = -1; + return; + } + // else + this->present_index = + this->tria->levels[this->present_level]->cells.n_objects() - 1; + } + } +} + + + +template +inline dealii::internal::TriangulationImplementation::TriaObjects & +TriaAccessorBase::objects() const +{ + if (structdim == dim) + return this->tria->levels[this->present_level]->cells; + + if (structdim == 1 && dim > 1) + return this->tria->faces->lines; + + if (structdim == 2 && dim > 2) + return this->tria->faces->quads; + + DEAL_II_ASSERT_UNREACHABLE(); + + return this->tria->levels[this->present_level]->cells; +} + + + +/*---------------------- Functions: InvalidAccessor -------------------------*/ + +template +InvalidAccessor::InvalidAccessor(const void *, + const int, + const int, + const AccessorData *) +{ + Assert(false, + ExcMessage("You are attempting an invalid conversion between " + "iterator/accessor types. The constructor you call " + "only exists to make certain template constructs " + "easier to write as dimension independent code but " + "the conversion is not valid in the current context.")); +} + + + +template +InvalidAccessor::InvalidAccessor( + const InvalidAccessor &) +{ + Assert(false, + ExcMessage("You are attempting an invalid conversion between " + "iterator/accessor types. The constructor you call " + "only exists to make certain template constructs " + "easier to write as dimension independent code but " + "the conversion is not valid in the current context.")); +} + + + +template +void +InvalidAccessor::copy_from(const InvalidAccessor &) +{ + // nothing to do here. we could + // throw an exception but we can't + // get here without first creating + // an object which would have + // already thrown +} + + + +template +bool +InvalidAccessor::operator==( + const InvalidAccessor &) const +{ + // nothing to do here. we could + // throw an exception but we can't + // get here without first creating + // an object which would have + // already thrown + return false; +} + + + +template +bool +InvalidAccessor::operator!=( + const InvalidAccessor &) const +{ + // nothing to do here. we could + // throw an exception but we can't + // get here without first creating + // an object which would have + // already thrown + return true; +} + + + +template +bool +InvalidAccessor::used() const +{ + // nothing to do here. we could + // throw an exception but we can't + // get here without first creating + // an object which would have + // already thrown + return false; +} + + + +template +bool +InvalidAccessor::has_children() const +{ + // nothing to do here. we could + // throw an exception but we can't + // get here without first creating + // an object which would have + // already thrown + return false; +} + + + +template +void +InvalidAccessor::operator++() const +{} + + + +template +void +InvalidAccessor::operator--() const +{} + + + +template +types::manifold_id +InvalidAccessor::manifold_id() const +{ + return numbers::flat_manifold_id; +} + + + +template +unsigned int +InvalidAccessor::user_index() const +{ + return numbers::invalid_unsigned_int; +} + + + +template +void +InvalidAccessor::set_user_index( + const unsigned int) const +{ + Assert(false, + ExcMessage("You are trying to set the user index of an " + "invalid object.")); +} + + + +template +void +InvalidAccessor::set_manifold_id( + const types::manifold_id) const +{ + Assert(false, + ExcMessage("You are trying to set the manifold id of an " + "invalid object.")); +} + + + +template +inline Point & +InvalidAccessor::vertex(const unsigned int) const +{ + // nothing to do here. we could throw an exception but we can't get here + // without first creating an object which would have already thrown + static Point invalid_vertex; + return invalid_vertex; +} + + +template +inline void * +InvalidAccessor::line(const unsigned int) const +{ + // nothing to do here. we could throw an exception but we can't get here + // without first creating an object which would have already thrown + return nullptr; +} + + + +template +inline void * +InvalidAccessor::quad(const unsigned int) const +{ + // nothing to do here. we could throw an exception but we can't get here + // without first creating an object which would have already thrown + return nullptr; +} + + +/*------------------------ Functions: TriaAccessor ---------------------------*/ + + +namespace internal +{ + namespace TriaAccessorImplementation + { + // make sure that if in the following we + // write TriaAccessor + // we mean the *class* + // dealii::TriaAccessor, not the + // enclosing namespace + // dealii::internal::TriaAccessor + using dealii::TriaAccessor; + + /** + * A class with the same purpose as the similarly named class of the + * Triangulation class. See there for more information. + */ + struct Implementation + { + /** + * Implementation of the function of some name in the parent class. + */ + template + inline static unsigned int + line_index(const TriaAccessor<1, dim, spacedim> &, const unsigned int) + { + Assert(false, + ExcMessage("You can't ask for the index of a line bounding " + "a one-dimensional cell because it is not " + "bounded by lines.")); + return numbers::invalid_unsigned_int; + } + + + template + inline static unsigned int + line_index(const TriaAccessor<2, dim, spacedim> &accessor, + const unsigned int i) + { + constexpr unsigned int max_faces_per_cell = 4; + return accessor.objects() + .cells[accessor.present_index * max_faces_per_cell + i]; + } + + + inline static unsigned int + line_index(const TriaAccessor<3, 3, 3> &accessor, const unsigned int i) + { + const auto [face_index, line_index] = + accessor.reference_cell().standard_line_to_face_and_line_index(i); + const auto line_within_face_index = + accessor.reference_cell().standard_to_real_face_line( + line_index, + face_index, + accessor.combined_face_orientation(face_index)); + + return accessor.quad(face_index)->line_index(line_within_face_index); + } + + + + /** + * Implementation of the function of some name in the parent class. + */ + template + inline static unsigned int + quad_index(const TriaAccessor &, + const unsigned int) + { + Assert(false, + ExcMessage("You can't ask for the index of a quad bounding " + "a one- or two-dimensional cell because it is not " + "bounded by quads.")); + return numbers::invalid_unsigned_int; + } + + + inline static unsigned int + quad_index(const TriaAccessor<3, 3, 3> &accessor, const unsigned int i) + { + constexpr unsigned int max_faces_per_cell = 6; + return accessor.tria->levels[accessor.present_level] + ->cells.cells[accessor.present_index * max_faces_per_cell + i]; + } + + + + /** + * Implementation of the function of some name in the parent class. + */ + template + inline static void + set_combined_face_orientation( + const TriaAccessor &accessor, + const unsigned int face, + const unsigned char combined_orientation) + { + Assert(structdim == dim, + ExcMessage("This function can only be used on objects that are " + "cells and not on objects which bound cells.")); + AssertIndexRange(face, accessor.n_faces()); + + if (dim == 1) + Assert(combined_orientation == + ReferenceCell::default_combined_face_orientation(), + ExcMessage("In 1d, faces do not have an orientation, so the " + "only valid value is the default.")); + else if (dim == 2) + Assert(combined_orientation == + ReferenceCell::default_combined_face_orientation() || + combined_orientation == + ReferenceCell::reversed_combined_line_orientation(), + ExcMessage( + "In 2d, the only valid values of the combined orientation " + "are the standard orientation or the reversed line " + "orientation.")); + + // face_orientations is not set up in 1d + if (dim != 1) + accessor.tria->levels[accessor.present_level] + ->face_orientations.set_combined_orientation( + accessor.present_index * GeometryInfo::faces_per_cell + face, + combined_orientation); + } + + /** + * Implementation of the function of some name in the parent class. + */ + template + inline static void + set_line_orientation(const TriaAccessor<1, dim, spacedim> &, + const unsigned int, + const bool) + { + DEAL_II_ASSERT_UNREACHABLE(); + } + + + template + inline static void + set_line_orientation(const TriaAccessor<2, 2, spacedim> &, + const unsigned int, + const bool) + { + // quads in 2d have no + // non-standard orientation + DEAL_II_ASSERT_UNREACHABLE(); + } + + + template + inline static void + set_line_orientation(const TriaAccessor<2, 3, spacedim> &accessor, + const unsigned int line, + const bool value) + { + Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed()); + AssertIndexRange(line, accessor.n_lines()); + Assert(accessor.present_index * GeometryInfo<3>::lines_per_face + line < + accessor.tria->faces->quads_line_orientations.size(), + ExcInternalError()); + + // quads as part of 3d hexes can have non-standard orientation + accessor.tria->faces->quads_line_orientations + [accessor.present_index * GeometryInfo<3>::lines_per_face + line] = + value; + } + + + inline static void + set_line_orientation(const TriaAccessor<3, 3, 3> &, + const unsigned int, + const bool) + { + // it seems like we don't need this + // one + DEAL_II_NOT_IMPLEMENTED(); + } + + + /** + * Implementation of the function of same name in the enclosing class. + */ + template + inline static unsigned int + vertex_index(const TriaAccessor<1, dim, spacedim> &accessor, + const unsigned int corner) + { + constexpr unsigned int max_faces_per_cell = 2; + return accessor.objects() + .cells[accessor.present_index * max_faces_per_cell + corner]; + } + + + template + inline static unsigned int + vertex_index(const TriaAccessor<2, dim, spacedim> &accessor, + const unsigned int corner) + { + const auto [line_index, vertex_index] = + accessor.reference_cell().standard_vertex_to_face_and_vertex_index( + corner); + const auto vertex_within_line_index = + accessor.reference_cell().standard_to_real_face_vertex( + vertex_index, + line_index, + accessor.line_orientation(line_index) == true ? + ReferenceCell::default_combined_face_orientation() : + ReferenceCell::reversed_combined_line_orientation()); + + return accessor.line(line_index) + ->vertex_index(vertex_within_line_index); + } + + + + inline static unsigned int + vertex_index(const TriaAccessor<3, 3, 3> &accessor, + const unsigned int corner) + { + const auto [face_index, vertex_index] = + accessor.reference_cell().standard_vertex_to_face_and_vertex_index( + corner); + const auto vertex_within_face_index = + accessor.reference_cell().standard_to_real_face_vertex( + vertex_index, + face_index, + accessor.combined_face_orientation(face_index)); + + return accessor.quad(face_index) + ->vertex_index(vertex_within_face_index); + } + + + + template + static std::array + get_line_indices_of_cell(const TriaAccessor<1, dim, spacedim> &) + { + DEAL_II_ASSERT_UNREACHABLE(); + return {}; + } + + + + template + static std::array + get_line_indices_of_cell(const TriaAccessor<2, dim, spacedim> &cell) + { + // For 2d cells the access cell->line_orientation() is already + // efficient + std::array line_indices = {}; + for (const unsigned int line : cell.line_indices()) + line_indices[line] = cell.line_index(line); + return line_indices; + } + + /** + * A helper function to provide faster access to cell->line_index() in + * 3d + */ + template + static std::array + get_line_indices_of_cell( + const TriaAccessor &cell) + { + std::array line_indices = {}; + + // For hexahedra, the classical access via quads -> lines is too + // inefficient. Unroll this code here to allow the compiler to inline + // the necessary functions. + const auto ref_cell = cell.reference_cell(); + if (ref_cell == ReferenceCells::Hexahedron) + { + for (unsigned int f = 4; f < 6; ++f) + { + const unsigned char orientation = + cell.get_triangulation() + .levels[cell.level()] + ->face_orientations.get_combined_orientation( + cell.index() * GeometryInfo<3>::faces_per_cell + f); + + // It might seem superfluous to spell out the four indices + // that get later consumed by a for loop over these four + // elements; however, for the compiler it is easier to inline + // the statement of standard_to_real_face_line() when next to + // each other, as opposed to be interleaved with a + // line_index() call. + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, f, orientation), + ref_cell.standard_to_real_face_line(1, f, orientation), + ref_cell.standard_to_real_face_line(2, f, orientation), + ref_cell.standard_to_real_face_line(3, f, orientation)}}; + const auto quad = cell.quad(f); + for (unsigned int l = 0; l < 4; ++l) + line_indices[4 * (f - 4) + l] = + quad->line_index(my_indices[l]); + } + for (unsigned int f = 0; f < 2; ++f) + { + const unsigned char orientation = + cell.get_triangulation() + .levels[cell.level()] + ->face_orientations.get_combined_orientation( + cell.index() * GeometryInfo<3>::faces_per_cell + f); + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, f, orientation), + ref_cell.standard_to_real_face_line(1, f, orientation)}}; + const auto quad = cell.quad(f); + line_indices[8 + f] = quad->line_index(my_indices[0]); + line_indices[10 + f] = quad->line_index(my_indices[1]); + } + } + else if (ref_cell == ReferenceCells::Tetrahedron) + { + std::array orientations{ + {cell.combined_face_orientation(0), + cell.combined_face_orientation(1), + cell.combined_face_orientation(2)}}; + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, 0, orientations[0]), + ref_cell.standard_to_real_face_line(1, 0, orientations[0]), + ref_cell.standard_to_real_face_line(2, 0, orientations[0]), + ref_cell.standard_to_real_face_line(1, 1, orientations[1]), + ref_cell.standard_to_real_face_line(2, 1, orientations[1]), + ref_cell.standard_to_real_face_line(1, 2, orientations[2])}}; + line_indices[0] = cell.quad(0)->line_index(my_indices[0]); + line_indices[1] = cell.quad(0)->line_index(my_indices[1]); + line_indices[2] = cell.quad(0)->line_index(my_indices[2]); + line_indices[3] = cell.quad(1)->line_index(my_indices[3]); + line_indices[4] = cell.quad(1)->line_index(my_indices[4]); + line_indices[5] = cell.quad(2)->line_index(my_indices[5]); + } + else + // For other shapes (wedges, pyramids), we do not currently + // implement an optimized function. + for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l) + line_indices[l] = cell.line_index(l); + + return line_indices; + } + + + + /** + * A helper function to provide faster access to + * cell->line_orientation(), 1d specialization + */ + template + static std::array + get_line_orientations_of_cell(const TriaAccessor<1, dim, spacedim> &) + { + DEAL_II_ASSERT_UNREACHABLE(); + return {}; + } + + + + /** + * A helper function to provide faster access to + * cell->line_orientation(), 2d specialization + */ + template + static std::array + get_line_orientations_of_cell(const TriaAccessor<2, dim, spacedim> &cell) + { + // For 2d cells the access cell->line_orientation() is already + // efficient + std::array line_orientations = {}; + for (const unsigned int line : cell.line_indices()) + line_orientations[line] = + cell.line_orientation(line) == true ? + ReferenceCell::default_combined_face_orientation() : + ReferenceCell::reversed_combined_line_orientation(); + return line_orientations; + } + + + + /** + * A helper function to provide faster access to + * cell->line_orientation(), 3d specialization + */ + template + static std::array + get_line_orientations_of_cell(const TriaAccessor<3, dim, spacedim> &cell) + { + std::array line_orientations = {}; + + // For hexahedra, the classical access via quads -> lines is too + // inefficient. Unroll this code here to allow the compiler to inline + // the necessary functions. + const auto ref_cell = cell.reference_cell(); + if (ref_cell == ReferenceCells::Hexahedron) + { + for (unsigned int f = 4; f < 6; ++f) + { + const unsigned char orientation = + cell.get_triangulation() + .levels[cell.level()] + ->face_orientations.get_combined_orientation( + cell.index() * GeometryInfo<3>::faces_per_cell + f); + + // It might seem superfluous to spell out the four indices and + // orientations that get later consumed by a for loop over + // these four elements; however, for the compiler it is easier + // to inline the statement of standard_to_real_face_line() + // when next to each other, as opposed to be interleaved with + // a line_index() call. + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, f, orientation), + ref_cell.standard_to_real_face_line(1, f, orientation), + ref_cell.standard_to_real_face_line(2, f, orientation), + ref_cell.standard_to_real_face_line(3, f, orientation)}}; + const auto quad = cell.quad(f); + const std::array my_orientations{ + {ref_cell.standard_vs_true_line_orientation( + 0, f, orientation, quad->line_orientation(my_indices[0])), + ref_cell.standard_vs_true_line_orientation( + 1, f, orientation, quad->line_orientation(my_indices[1])), + ref_cell.standard_vs_true_line_orientation( + 2, f, orientation, quad->line_orientation(my_indices[2])), + ref_cell.standard_vs_true_line_orientation( + 3, + f, + orientation, + quad->line_orientation(my_indices[3]))}}; + for (unsigned int l = 0; l < 4; ++l) + line_orientations[4 * (f - 4) + l] = my_orientations[l]; + } + for (unsigned int f = 0; f < 2; ++f) + { + const unsigned char orientation = + cell.get_triangulation() + .levels[cell.level()] + ->face_orientations.get_combined_orientation( + cell.index() * GeometryInfo<3>::faces_per_cell + f); + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, f, orientation), + ref_cell.standard_to_real_face_line(1, f, orientation)}}; + const auto quad = cell.quad(f); + const std::array my_orientations{ + {ref_cell.standard_vs_true_line_orientation( + 0, f, orientation, quad->line_orientation(my_indices[0])), + ref_cell.standard_vs_true_line_orientation( + 1, + f, + orientation, + quad->line_orientation(my_indices[1]))}}; + line_orientations[8 + f] = my_orientations[0]; + line_orientations[10 + f] = my_orientations[1]; + } + } + else if (ref_cell == ReferenceCells::Tetrahedron) + { + std::array orientations{ + {cell.combined_face_orientation(0), + cell.combined_face_orientation(1), + cell.combined_face_orientation(2)}}; + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, 0, orientations[0]), + ref_cell.standard_to_real_face_line(1, 0, orientations[0]), + ref_cell.standard_to_real_face_line(2, 0, orientations[0]), + ref_cell.standard_to_real_face_line(1, 1, orientations[1]), + ref_cell.standard_to_real_face_line(2, 1, orientations[1]), + ref_cell.standard_to_real_face_line(1, 2, orientations[2])}}; + line_orientations[0] = ref_cell.standard_vs_true_line_orientation( + 0, + 0, + orientations[0], + cell.quad(0)->line_orientation(my_indices[0])); + line_orientations[1] = ref_cell.standard_vs_true_line_orientation( + 1, + 0, + orientations[0], + cell.quad(0)->line_orientation(my_indices[1])); + line_orientations[2] = ref_cell.standard_vs_true_line_orientation( + 2, + 0, + orientations[0], + cell.quad(0)->line_orientation(my_indices[2])); + line_orientations[3] = ref_cell.standard_vs_true_line_orientation( + 1, + 1, + orientations[1], + cell.quad(1)->line_orientation(my_indices[3])); + line_orientations[4] = ref_cell.standard_vs_true_line_orientation( + 2, + 1, + orientations[1], + cell.quad(1)->line_orientation(my_indices[4])); + line_orientations[5] = ref_cell.standard_vs_true_line_orientation( + 1, + 2, + orientations[2], + cell.quad(2)->line_orientation(my_indices[5])); + } + else + // For other shapes (wedges, pyramids), we do not currently + // implement an optimized function + for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l) + line_orientations[l] = cell.line_orientation(l); + + return line_orientations; + } + }; + } // namespace TriaAccessorImplementation +} // namespace internal + + + +template +inline TriaAccessor::TriaAccessor( + const Triangulation *parent, + const int level, + const int index, + const AccessorData *local_data) + : TriaAccessorBase(parent, level, index, local_data) +{} + + + +template +inline bool +TriaAccessor::used() const +{ + Assert(this->state() == IteratorState::valid, + TriaAccessorExceptions::ExcDereferenceInvalidObject( + *this)); + return this->objects().used[this->present_index]; +} + + + +template +inline TriaIterator> +TriaAccessor::vertex_iterator( + const unsigned int i) const +{ + return TriaIterator>(this->tria, + 0, + vertex_index(i)); +} + + + +template +inline ReferenceCell +TriaAccessor::reference_cell() const +{ + if (structdim == 0) + return ReferenceCells::Vertex; + else if (structdim == 1) + return ReferenceCells::Line; + else if (structdim == dim) + return this->tria->levels[this->present_level] + ->reference_cell[this->present_index]; + else + return this->tria->faces->get_quad_type(this->present_index); +} + + + +template +inline unsigned int +TriaAccessor::vertex_index( + const unsigned int corner) const +{ + AssertIndexRange(corner, this->n_vertices()); + + if (structdim == dim) + { + constexpr unsigned int max_vertices_per_cell = 1 << dim; + const std::size_t my_index = + static_cast(this->present_index) * max_vertices_per_cell; + AssertIndexRange(my_index + corner, + this->tria->levels[this->present_level] + ->cell_vertex_indices_cache.size()); + const unsigned int vertex_index = + this->tria->levels[this->present_level] + ->cell_vertex_indices_cache[my_index + corner]; + Assert(vertex_index != numbers::invalid_unsigned_int, ExcInternalError()); + return vertex_index; + } + else + return dealii::internal::TriaAccessorImplementation::Implementation:: + vertex_index(*this, corner); +} + + + +template +inline Point & +TriaAccessor::vertex(const unsigned int i) const +{ + return const_cast &>(this->tria->vertices[vertex_index(i)]); +} + + + +template +inline typename dealii::internal::TriangulationImplementation:: + Iterators::line_iterator + TriaAccessor::line(const unsigned int i) const +{ + // checks happen in line_index + return typename dealii::internal::TriangulationImplementation:: + Iterators::line_iterator(this->tria, 0, line_index(i)); +} + + + +template +inline unsigned int +TriaAccessor::line_index(const unsigned int i) const +{ + AssertIndexRange(i, this->n_lines()); + + return dealii::internal::TriaAccessorImplementation::Implementation:: + line_index(*this, i); +} + + + +template +inline typename dealii::internal::TriangulationImplementation:: + Iterators::quad_iterator + TriaAccessor::quad(const unsigned int i) const +{ + // checks happen in quad_index + return typename dealii::internal::TriangulationImplementation:: + Iterators::quad_iterator(this->tria, 0, quad_index(i)); +} + + + +template +inline unsigned int +TriaAccessor::quad_index(const unsigned int i) const +{ + return dealii::internal::TriaAccessorImplementation::Implementation:: + quad_index(*this, i); +} + + + +template +inline unsigned char +TriaAccessor::combined_face_orientation( + const unsigned int face) const +{ + Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); + AssertIndexRange(face, n_faces()); + Assert(structdim == dim, + ExcMessage("This function can only be used on objects " + "that are cells, but not on faces or edges " + "that bound cells.")); + // work around a bogus GCC-9 warning which considers face unused except in 3d + (void)face; + + if constexpr (structdim == 1) + return ReferenceCell::default_combined_face_orientation(); + else if constexpr (structdim == 2) + return this->line_orientation(face) == true ? + ReferenceCell::default_combined_face_orientation() : + ReferenceCell::reversed_combined_line_orientation(); + else + return this->tria->levels[this->present_level] + ->face_orientations.get_combined_orientation( + this->present_index * GeometryInfo::faces_per_cell + face); +} + + + +template +inline bool +TriaAccessor::face_orientation( + const unsigned int face) const +{ + Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); + AssertIndexRange(face, n_faces()); + Assert(structdim == dim, + ExcMessage("This function can only be used on objects " + "that are cells, but not on faces or edges " + "that bound cells.")); + // work around a bogus GCC-9 warning which considers face unused in 1d + (void)face; + + if constexpr (structdim == 1) + // in 1d 'faces' are vertices and those are always consistently oriented + return true; + else if constexpr (structdim == 2) + return this->line_orientation(face); + else + return this->tria->levels[this->present_level] + ->face_orientations.get_orientation( + this->present_index * GeometryInfo::faces_per_cell + face); +} + + + +template +inline bool +TriaAccessor::face_flip(const unsigned int face) const +{ + Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); + Assert(structdim == dim, + ExcMessage("This function can only be used on objects " + "that are cells, but not on faces or edges " + "that bound cells.")); + AssertIndexRange(face, n_faces()); + // work around a bogus GCC-9 warning which considers face unused except in 3d + (void)face; + + if constexpr (structdim == 3) + return this->tria->levels[this->present_level]->face_orientations.get_flip( + this->present_index * GeometryInfo<3>::faces_per_cell + face); + else + // In 1d and 2d, face_flip is always false as faces can only be + // 'flipped' in 3d. + return false; +} + + +template +inline bool +TriaAccessor::face_rotation( + const unsigned int face) const +{ + Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); + Assert(structdim == dim, + ExcMessage("This function can only be used on objects " + "that are cells, but not on faces or edges " + "that bound cells.")); + AssertIndexRange(face, n_faces()); + // work around a bogus GCC-9 warning which considers face unused except in 3d + (void)face; + + if constexpr (structdim == 3) + return this->tria->levels[this->present_level] + ->face_orientations.get_rotation( + this->present_index * GeometryInfo<3>::faces_per_cell + face); + else + // In 1d and 2d, face_rotation is always false as faces can only be + // 'rotated' in 3d. + return false; +} + + + +template +inline bool +TriaAccessor::line_orientation( + const unsigned int line) const +{ + Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); + AssertIndexRange(line, this->n_lines()); + // work around a bogus GCC-9 warning which considers line unused in 1d + (void)line; + + if constexpr (structdim == 1) + return true; + else if constexpr (structdim == 2 && dim == 2) + // lines in 2d are faces + { + // if all elements are quads (or if we have a very special consistently + // oriented triangular mesh) then we do not store this array + if (this->tria->levels[this->present_level] + ->face_orientations.n_objects() == 0) + { + return true; + } + else + { + return this->tria->levels[this->present_level] + ->face_orientations.get_orientation( + this->present_index * GeometryInfo::faces_per_cell + + line); + } + } + else if constexpr (structdim == 2 && dim == 3) + { + // line orientations in 3d are stored in their own array + Assert(this->present_index * GeometryInfo<3>::lines_per_face + line < + this->tria->faces->quads_line_orientations.size(), + ExcInternalError()); + return this->tria->faces->quads_line_orientations + [this->present_index * GeometryInfo<3>::lines_per_face + line]; + } + else if constexpr (structdim == 3 && dim == 3) + { + const auto reference_cell = this->reference_cell(); + // First pick a face on which this line is a part of, and the + // index of the line within. + const auto [face_index, line_index] = + reference_cell.standard_line_to_face_and_line_index(line); + const auto line_within_face_index = + reference_cell.standard_to_real_face_line( + line_index, face_index, this->combined_face_orientation(face_index)); + + // Then query how that line is oriented within that face: + return reference_cell.standard_vs_true_line_orientation( + line_index, + face_index, + this->combined_face_orientation(face_index), + this->quad(face_index)->line_orientation(line_within_face_index)); + } + else + { + DEAL_II_ASSERT_UNREACHABLE(); + return false; + } +} + + + +template +inline void +TriaAccessor::set_line_orientation( + const unsigned int line, + const bool value) const +{ + Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); + AssertIndexRange(line, this->n_lines()); + + dealii::internal::TriaAccessorImplementation::Implementation:: + set_line_orientation(*this, line, value); +} + + + +template +inline void +TriaAccessor::set_combined_face_orientation( + const unsigned int face, + const unsigned char combined_orientation) const +{ + Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); + AssertIndexRange(face, this->n_faces()); + + dealii::internal::TriaAccessorImplementation::Implementation:: + set_combined_face_orientation(*this, face, combined_orientation); +} + + + +template +void +TriaAccessor::set_used_flag() const +{ + Assert(this->state() == IteratorState::valid, + TriaAccessorExceptions::ExcDereferenceInvalidObject( + *this)); + this->objects().used[this->present_index] = true; +} + + + +template +void +TriaAccessor::clear_used_flag() const +{ + Assert(this->state() == IteratorState::valid, + TriaAccessorExceptions::ExcDereferenceInvalidObject( + *this)); + this->objects().used[this->present_index] = false; +} + + +template +int +TriaAccessor::child_index(const unsigned int i) const +{ + Assert(has_children(), TriaAccessorExceptions::ExcCellHasNoChildren()); + AssertIndexRange(i, n_children()); + + // each set of two children are stored + // consecutively, so we only have to find + // the location of the set of children + const unsigned int n_sets_of_two = + GeometryInfo::max_children_per_cell / 2; + return this->objects().children[n_sets_of_two * this->present_index + i / 2] + + i % 2; +} + + + +template +int +TriaAccessor::isotropic_child_index( + const unsigned int i) const +{ + AssertIndexRange(i, GeometryInfo::max_children_per_cell); + + switch (structdim) + { + case 1: + return child_index(i); + case 2: + { + const RefinementCase<2> this_refinement_case( + static_cast(refinement_case())); + + Assert(this_refinement_case != RefinementCase<2>::no_refinement, + TriaAccessorExceptions::ExcCellHasNoChildren()); + + if (this_refinement_case == RefinementCase<2>::cut_xy) + return child_index(i); + else if ((this_refinement_case == RefinementCase<2>::cut_x) && + (child(i % 2)->refinement_case() == + RefinementCase<2>::cut_y)) + return child(i % 2)->child_index(i / 2); + else if ((this_refinement_case == RefinementCase<2>::cut_y) && + (child(i / 2)->refinement_case() == + RefinementCase<2>::cut_x)) + return child(i / 2)->child_index(i % 2); + else + Assert( + false, + ExcMessage( + "This cell has no grandchildren equivalent to isotropic refinement")); + break; + } + + case 3: + DEAL_II_NOT_IMPLEMENTED(); + } + return -1; +} + + + +template +RefinementCase +TriaAccessor::refinement_case() const +{ + Assert(this->state() == IteratorState::valid, + TriaAccessorExceptions::ExcDereferenceInvalidObject( + *this)); + + switch (structdim) + { + case 1: + return (RefinementCase( + this->objects().children[this->present_index] != -1 ? + // cast the branches + // here first to uchar + // and then (above) to + // RefinementCase + // so that the + // conversion is valid + // even for the case + // structdim>1 (for + // which this part of + // the code is dead + // anyway) + static_cast(RefinementCase<1>::cut_x) : + static_cast(RefinementCase<1>::no_refinement))); + + default: + Assert(static_cast(this->present_index) < + this->objects().refinement_cases.size(), + ExcIndexRange(this->present_index, + 0, + this->objects().refinement_cases.size())); + + return (static_cast>( + this->objects().refinement_cases[this->present_index])); + } +} + + + +template +inline TriaIterator> +TriaAccessor::child(const unsigned int i) const + +{ + // checking of 'i' happens in child_index + const TriaIterator> q( + this->tria, (dim == structdim ? this->level() + 1 : 0), child_index(i)); + + Assert((q.state() == IteratorState::past_the_end) || q->used(), + ExcInternalError()); + + return q; +} + + + +template +inline unsigned int +TriaAccessor::child_iterator_to_index( + const TriaIterator> &child) const +{ + const auto n_children = this->n_children(); + for (unsigned int child_n = 0; child_n < n_children; ++child_n) + if (this->child(child_n) == child) + return child_n; + + Assert(false, + ExcMessage("The given child is not a child of the current object.")); + return numbers::invalid_unsigned_int; +} + + + +template +inline TriaIterator> +TriaAccessor::isotropic_child( + const unsigned int i) const +{ + // checking of 'i' happens in child() or + // child_index() called below + switch (structdim) + { + case 1: + // no anisotropic refinement in 1d + return child(i); + + case 2: + { + const RefinementCase<2> this_refinement_case( + static_cast(refinement_case())); + + Assert(this_refinement_case != RefinementCase<2>::no_refinement, + TriaAccessorExceptions::ExcCellHasNoChildren()); + + if (this_refinement_case == RefinementCase<2>::cut_xy) + return child(i); + else if ((this_refinement_case == RefinementCase<2>::cut_x) && + (child(i % 2)->refinement_case() == + RefinementCase<2>::cut_y)) + return child(i % 2)->child(i / 2); + else if ((this_refinement_case == RefinementCase<2>::cut_y) && + (child(i / 2)->refinement_case() == + RefinementCase<2>::cut_x)) + return child(i / 2)->child(i % 2); + else + Assert( + false, + ExcMessage( + "This cell has no grandchildren equivalent to isotropic refinement")); + break; + } + + default: + DEAL_II_NOT_IMPLEMENTED(); + } + // we don't get here but have to return + // something... + return child(0); +} + + + +template +inline bool +TriaAccessor::has_children() const +{ + Assert(this->state() == IteratorState::valid, + TriaAccessorExceptions::ExcDereferenceInvalidObject( + *this)); + + // each set of two children are stored + // consecutively, so we only have to find + // the location of the set of children + const unsigned int n_sets_of_two = + GeometryInfo::max_children_per_cell / 2; + return (this->objects().children[n_sets_of_two * this->present_index] != -1); +} + + + +template +inline unsigned int +TriaAccessor::n_children() const +{ + Assert(this->state() == IteratorState::valid, + TriaAccessorExceptions::ExcDereferenceInvalidObject( + *this)); + if (reference_cell() == ReferenceCells::Tetrahedron) + return GeometryInfo::n_children( + RefinementCase::isotropic_refinement); + else + return GeometryInfo::n_children(refinement_case()); +} + + + +template +inline void +TriaAccessor::set_refinement_case( + const RefinementCase &refinement_case) const +{ + Assert(this->state() == IteratorState::valid, + TriaAccessorExceptions::ExcDereferenceInvalidObject( + *this)); + Assert(static_cast(this->present_index) < + this->objects().refinement_cases.size(), + ExcIndexRange(this->present_index, + 0, + this->objects().refinement_cases.size())); + + this->objects().refinement_cases[this->present_index] = refinement_case; +} + + +template +inline void +TriaAccessor::clear_refinement_case() const +{ + Assert(this->state() == IteratorState::valid, + TriaAccessorExceptions::ExcDereferenceInvalidObject( + *this)); + Assert(static_cast(this->present_index) < + this->objects().refinement_cases.size(), + ExcIndexRange(this->present_index, + 0, + this->objects().refinement_cases.size())); + + this->objects().refinement_cases[this->present_index] = + RefinementCase::no_refinement; +} + + + +template +void +TriaAccessor::set_children(const unsigned int i, + const int index) const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + Assert(i % 2 == 0, TriaAccessorExceptions::ExcSetOnlyEvenChildren(i)); + + // each set of two children are stored + // consecutively, so we only have to find + // the location of the set of children + const unsigned int n_sets_of_two = + GeometryInfo::max_children_per_cell / 2; + + Assert( + // clearing the child index for a cell + (index == -1) || + // if setting the child index for the i'th child (with i==0), + // then the index must be a non-negative number + (i == 0 && !this->has_children() && (index >= 0)) || + // if setting the child index for the i'th child (with i>0), + // then the previously stored index must be the invalid + // index + (i > 0 && this->has_children() && (index >= 0) && + this->objects().children[n_sets_of_two * this->present_index + i / 2] == + -1), + TriaAccessorExceptions::ExcCantSetChildren(index)); + + this->objects().children[n_sets_of_two * this->present_index + i / 2] = index; +} + + + +template +void +TriaAccessor::clear_children() const +{ + // each set of two children are stored + // consecutively, so we only have to find + // the location of the set of children + const unsigned int n_sets_of_two = + GeometryInfo::max_children_per_cell / 2; + + for (unsigned int i = 0; i < n_sets_of_two; ++i) + set_children(2 * i, -1); +} + + + +template +inline bool +TriaAccessor::user_flag_set() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + return this->objects().user_flags[this->present_index]; +} + + + +template +inline void +TriaAccessor::set_user_flag() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + this->objects().user_flags[this->present_index] = true; +} + + + +template +inline void +TriaAccessor::clear_user_flag() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + this->objects().user_flags[this->present_index] = false; +} + + + +template +void +TriaAccessor::recursively_set_user_flag() const +{ + set_user_flag(); + + if (this->has_children()) + for (unsigned int c = 0; c < this->n_children(); ++c) + this->child(c)->recursively_set_user_flag(); +} + + + +template +void +TriaAccessor::recursively_clear_user_flag() const +{ + clear_user_flag(); + + if (this->has_children()) + for (unsigned int c = 0; c < this->n_children(); ++c) + this->child(c)->recursively_clear_user_flag(); +} + + + +template +void +TriaAccessor::clear_user_data() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + this->objects().clear_user_data(this->present_index); +} + + + +template +void +TriaAccessor::set_user_pointer(void *p) const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + this->objects().user_pointer(this->present_index) = p; +} + + + +template +void +TriaAccessor::clear_user_pointer() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + this->objects().user_pointer(this->present_index) = nullptr; +} + + + +template +void * +TriaAccessor::user_pointer() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + return this->objects().user_pointer(this->present_index); +} + + + +template +void +TriaAccessor::recursively_set_user_pointer( + void *p) const +{ + set_user_pointer(p); + + if (this->has_children()) + for (unsigned int c = 0; c < this->n_children(); ++c) + this->child(c)->recursively_set_user_pointer(p); +} + + + +template +void +TriaAccessor::recursively_clear_user_pointer() const +{ + clear_user_pointer(); + + if (this->has_children()) + for (unsigned int c = 0; c < this->n_children(); ++c) + this->child(c)->recursively_clear_user_pointer(); +} + + + +template +void +TriaAccessor::set_user_index( + const unsigned int p) const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + this->objects().user_index(this->present_index) = p; +} + + + +template +void +TriaAccessor::clear_user_index() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + this->objects().user_index(this->present_index) = 0; +} + + + +template +unsigned int +TriaAccessor::user_index() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + return this->objects().user_index(this->present_index); +} + + + +template +void +TriaAccessor::recursively_set_user_index( + const unsigned int p) const +{ + set_user_index(p); + + if (this->has_children()) + for (unsigned int c = 0; c < this->n_children(); ++c) + this->child(c)->recursively_set_user_index(p); +} + + + +template +void +TriaAccessor::recursively_clear_user_index() const +{ + clear_user_index(); + + if (this->has_children()) + for (unsigned int c = 0; c < this->n_children(); ++c) + this->child(c)->recursively_clear_user_index(); +} + + + +template +inline unsigned int +TriaAccessor::max_refinement_depth() const +{ + if (!this->has_children()) + return 0; + + unsigned int max_depth = 1; + for (unsigned int c = 0; c < n_children(); ++c) + max_depth = std::max(max_depth, child(c)->max_refinement_depth() + 1); + return max_depth; +} + + + +template +unsigned int +TriaAccessor::n_active_descendants() const +{ + if (!this->has_children()) + return 1; + else + { + unsigned int sum = 0; + for (unsigned int c = 0; c < n_children(); ++c) + sum += this->child(c)->n_active_descendants(); + return sum; + } +} + + + +template +types::boundary_id +TriaAccessor::boundary_id() const +{ + Assert(structdim < dim, ExcImpossibleInDim(dim)); + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + + return this->objects() + .boundary_or_material_id[this->present_index] + .boundary_id; +} + + + +template +void +TriaAccessor::set_boundary_id( + const types::boundary_id boundary_ind) const +{ + Assert(structdim < dim, ExcImpossibleInDim(dim)); + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + Assert(boundary_ind != numbers::internal_face_boundary_id, + ExcMessage("You are trying to set the boundary_id to an invalid " + "value (numbers::internal_face_boundary_id is reserved).")); + Assert(this->at_boundary(), + ExcMessage("You are trying to set the boundary_id of an " + "internal object, which is not allowed!")); + + this->objects().boundary_or_material_id[this->present_index].boundary_id = + boundary_ind; +} + + + +template +void +TriaAccessor::set_boundary_id_internal( + const types::boundary_id boundary_ind) const +{ + Assert(structdim < dim, ExcImpossibleInDim(dim)); + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + + this->objects().boundary_or_material_id[this->present_index].boundary_id = + boundary_ind; +} + + + +template +void +TriaAccessor::set_all_boundary_ids( + const types::boundary_id boundary_ind) const +{ + set_boundary_id(boundary_ind); + + switch (structdim) + { + case 1: + // 1d objects have no sub-objects + // where we have to do anything + break; + + case 2: + // for boundary quads also set + // boundary_id of bounding lines + for (unsigned int i = 0; i < this->n_lines(); ++i) + this->line(i)->set_boundary_id(boundary_ind); + break; + + default: + DEAL_II_NOT_IMPLEMENTED(); + } +} + + + +template +bool +TriaAccessor::at_boundary() const +{ + // error checking is done + // in boundary_id() + return (boundary_id() != numbers::internal_face_boundary_id); +} + + + +template +const Manifold & +TriaAccessor::get_manifold() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + return this->tria->get_manifold(this->manifold_id()); +} + + +template +types::manifold_id +TriaAccessor::manifold_id() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + + return this->objects().manifold_id[this->present_index]; +} + + + +template +void +TriaAccessor::set_manifold_id( + const types::manifold_id manifold_ind) const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + + this->objects().manifold_id[this->present_index] = manifold_ind; +} + + +template +void +TriaAccessor::set_all_manifold_ids( + const types::manifold_id manifold_ind) const +{ + set_manifold_id(manifold_ind); + + if (this->has_children()) + for (unsigned int c = 0; c < this->n_children(); ++c) + this->child(c)->set_all_manifold_ids(manifold_ind); + + switch (structdim) + { + case 1: + if (dim == 1) + { + (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(0)] = + manifold_ind; + (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(1)] = + manifold_ind; + } + break; + + case 2: + // for quads/simplices also set manifold_id of bounding lines + for (unsigned int i = 0; i < this->n_lines(); ++i) + this->line(i)->set_manifold_id(manifold_ind); + break; + default: + DEAL_II_NOT_IMPLEMENTED(); + } +} + + + +template +double +TriaAccessor::diameter() const +{ + boost::container::small_vector, + GeometryInfo::vertices_per_cell> + vertices(this->n_vertices()); + + for (unsigned int v = 0; v < vertices.size(); ++v) + vertices[v] = this->vertex(v); + + return internal::TriaAccessorImplementation::diameter( + vertices); +} + + + +template +double +CellAccessor::diameter( + const Mapping &mapping) const +{ + return internal::TriaAccessorImplementation::diameter( + mapping.get_vertices(typename Triangulation::cell_iterator( + this->tria, this->level(), this->index()))); +} + + + +template +std::pair, double> +TriaAccessor::enclosing_ball() const +{ + // If the object is one dimensional, + // the enclosing ball is the initial iterate + // i.e., the ball's center and diameter are + // the center and the diameter of the object. + if (structdim == 1) + return std::make_pair((this->vertex(1) + this->vertex(0)) * 0.5, + (this->vertex(1) - this->vertex(0)).norm() * 0.5); + + // The list is_initial_guess_vertex contains bool values and has + // the same size as the number of vertices per object. + // The entries of is_initial_guess_vertex are set true only for those + // two vertices corresponding to the largest diagonal which is being used + // to construct the initial ball. + // We employ this mask to skip these two vertices while enlarging the ball. + std::vector is_initial_guess_vertex(this->n_vertices()); + + // First let all the vertices be outside + std::fill(is_initial_guess_vertex.begin(), + is_initial_guess_vertex.end(), + false); + + // Get an initial guess by looking at the largest diagonal + Point center; + double radius = 0; + + switch (structdim) + { + case 2: + { + const Point p30(this->vertex(3) - this->vertex(0)); + const Point p21(this->vertex(2) - this->vertex(1)); + if (p30.norm() > p21.norm()) + { + center = this->vertex(0) + 0.5 * p30; + radius = p30.norm() / 2.; + is_initial_guess_vertex[3] = true; + is_initial_guess_vertex[0] = true; + } + else + { + center = this->vertex(1) + 0.5 * p21; + radius = p21.norm() / 2.; + is_initial_guess_vertex[2] = true; + is_initial_guess_vertex[1] = true; + } + break; + } + case 3: + { + const Point p70(this->vertex(7) - this->vertex(0)); + const Point p61(this->vertex(6) - this->vertex(1)); + const Point p25(this->vertex(2) - this->vertex(5)); + const Point p34(this->vertex(3) - this->vertex(4)); + const std::vector diagonals = {p70.norm(), + p61.norm(), + p25.norm(), + p34.norm()}; + const std::vector::const_iterator it = + std::max_element(diagonals.begin(), diagonals.end()); + if (it == diagonals.begin()) + { + center = this->vertex(0) + 0.5 * p70; + is_initial_guess_vertex[7] = true; + is_initial_guess_vertex[0] = true; + } + else if (it == diagonals.begin() + 1) + { + center = this->vertex(1) + 0.5 * p61; + is_initial_guess_vertex[6] = true; + is_initial_guess_vertex[1] = true; + } + else if (it == diagonals.begin() + 2) + { + center = this->vertex(5) + 0.5 * p25; + is_initial_guess_vertex[2] = true; + is_initial_guess_vertex[5] = true; + } + else + { + center = this->vertex(4) + 0.5 * p34; + is_initial_guess_vertex[3] = true; + is_initial_guess_vertex[4] = true; + } + radius = *it * 0.5; + break; + } + default: + DEAL_II_NOT_IMPLEMENTED(); + return std::pair, double>(); + } + + // For each vertex that is found to be geometrically outside the ball + // enlarge the ball so that the new ball contains both the previous ball + // and the given vertex. + for (const unsigned int v : this->vertex_indices()) + if (!is_initial_guess_vertex[v]) + { + const double distance = center.distance(this->vertex(v)); + if (distance > radius) + { + // we found a vertex which is outside of the ball + // extend it (move center and change radius) + const Point pCV(center - this->vertex(v)); + radius = (distance + radius) * 0.5; + center = this->vertex(v) + pCV * (radius / distance); + + // Now the new ball constructed in this block + // encloses the vertex (v) that was found to be geometrically + // outside the old ball. + } + } +#ifdef DEBUG + bool all_vertices_within_ball = true; + + // Set all_vertices_within_ball false if any of the vertices of the object + // are geometrically outside the ball + for (const unsigned int v : this->vertex_indices()) + if (center.distance(this->vertex(v)) > + radius + 100. * std::numeric_limits::epsilon()) + { + all_vertices_within_ball = false; + break; + } + // If all the vertices are not within the ball throw error + Assert(all_vertices_within_ball, ExcInternalError()); +#endif + return std::make_pair(center, radius); +} + + +template +double +TriaAccessor::minimum_vertex_distance() const +{ + switch (structdim) + { + case 1: + return (this->vertex(1) - this->vertex(0)).norm(); + case 2: + case 3: + { + double min = std::numeric_limits::max(); + for (const unsigned int i : this->vertex_indices()) + for (unsigned int j = i + 1; j < this->n_vertices(); ++j) + min = std::min(min, + (this->vertex(i) - this->vertex(j)) * + (this->vertex(i) - this->vertex(j))); + return std::sqrt(min); + } + default: + DEAL_II_NOT_IMPLEMENTED(); + return -1e10; + } +} + + +template +bool +TriaAccessor::is_translation_of( + const TriaIterator> &o) const +{ + // go through the vertices and check... The + // cell is a translation of the previous + // one in case the distance between the + // individual vertices in the two cell is + // the same for all the vertices. So do the + // check by first getting the distance on + // the first vertex, and then checking + // whether all others have the same down to + // rounding errors (we have to be careful + // here because the calculation of the + // displacement between one cell and the + // next can already result in the loss of + // one or two digits), so we choose 1e-12 + // times the distance between the zeroth + // vertices here. + bool is_translation = true; + const Tensor<1, spacedim> dist = o->vertex(0) - this->vertex(0); + const double tol_square = 1e-24 * dist.norm_square(); + for (unsigned int i = 1; i < this->n_vertices(); ++i) + { + const Tensor<1, spacedim> dist_new = + (o->vertex(i) - this->vertex(i)) - dist; + if (dist_new.norm_square() > tol_square) + { + is_translation = false; + break; + } + } + return is_translation; +} + + + +template +unsigned int +TriaAccessor::n_vertices() const +{ + return this->reference_cell().n_vertices(); +} + + + +template +unsigned int +TriaAccessor::n_lines() const +{ + return this->reference_cell().n_lines(); +} + + + +template +unsigned int +TriaAccessor::n_faces() const +{ + Assert(structdim == dim, + ExcMessage("This function can only be used on objects " + "that are cells, but not on faces or edges " + "that bound cells.")); + + return this->reference_cell().n_faces(); +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor::vertex_indices() const +{ + return {0U, n_vertices()}; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor::line_indices() const +{ + return {0U, n_lines()}; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor::face_indices() const +{ + return {0U, n_faces()}; +} + + + +/*----------------- Functions: TriaAccessor<0,dim,spacedim> -----------------*/ + +template +inline TriaAccessor<0, dim, spacedim>::TriaAccessor( + const Triangulation *tria, + const unsigned int vertex_index) + : tria(tria) + , global_vertex_index(vertex_index) +{} + + + +template +inline TriaAccessor<0, dim, spacedim>::TriaAccessor( + const Triangulation *tria, + const int /*level*/, + const int index, + const AccessorData *) + : tria(tria) + , global_vertex_index(index) +{} + + + +template +template +inline TriaAccessor<0, dim, spacedim>::TriaAccessor( + const TriaAccessor &) + : tria(nullptr) + , global_vertex_index(numbers::invalid_unsigned_int) +{ + Assert(false, ExcImpossibleInDim(0)); +} + + + +template +template +inline TriaAccessor<0, dim, spacedim>::TriaAccessor( + const InvalidAccessor &) + : tria(nullptr) + , global_vertex_index(numbers::invalid_unsigned_int) +{ + Assert(false, ExcImpossibleInDim(0)); +} + + + +template +inline void +TriaAccessor<0, dim, spacedim>::copy_from(const TriaAccessor &t) +{ + tria = t.tria; + global_vertex_index = t.global_vertex_index; +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::operator<( + const TriaAccessor<0, dim, spacedim> &other) const +{ + Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators()); + + return (global_vertex_index < other.global_vertex_index); +} + + + +template +inline IteratorState::IteratorStates +TriaAccessor<0, dim, spacedim>::state() const +{ + if (global_vertex_index != numbers::invalid_unsigned_int) + return IteratorState::valid; + else + return IteratorState::past_the_end; +} + + + +template +inline int +TriaAccessor<0, dim, spacedim>::level() +{ + return 0; +} + + + +template +inline int +TriaAccessor<0, dim, spacedim>::index() const +{ + return global_vertex_index; +} + + + +template +inline const Triangulation & +TriaAccessor<0, dim, spacedim>::get_triangulation() const +{ + return *tria; +} + + + +template +inline void +TriaAccessor<0, dim, spacedim>::operator++() +{ + ++global_vertex_index; + if (global_vertex_index >= tria->n_vertices()) + global_vertex_index = numbers::invalid_unsigned_int; +} + + + +template +inline void +TriaAccessor<0, dim, spacedim>::operator--() +{ + if (global_vertex_index != numbers::invalid_unsigned_int) + { + if (global_vertex_index != 0) + --global_vertex_index; + else + global_vertex_index = numbers::invalid_unsigned_int; + } +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::operator==(const TriaAccessor &t) const +{ + const bool result = + ((tria == t.tria) && (global_vertex_index == t.global_vertex_index)); + + return result; +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::operator!=(const TriaAccessor &t) const +{ + return !(*this == t); +} + + + +template +inline unsigned int +TriaAccessor<0, dim, spacedim>::vertex_index(const unsigned int) const +{ + return global_vertex_index; +} + + + +template +inline Point & +TriaAccessor<0, dim, spacedim>::vertex(const unsigned int) const +{ + return const_cast &>( + this->tria->vertices[global_vertex_index]); +} + + + +template +inline typename dealii::internal::TriangulationImplementation:: + Iterators::line_iterator + TriaAccessor<0, dim, spacedim>::line(const unsigned int) +{ + return typename dealii::internal::TriangulationImplementation:: + Iterators::line_iterator(); +} + + + +template +inline unsigned int +TriaAccessor<0, dim, spacedim>::line_index(const unsigned int) +{ + Assert(false, ExcImpossibleInDim(0)); + return numbers::invalid_unsigned_int; +} + + + +template +inline typename dealii::internal::TriangulationImplementation:: + Iterators::quad_iterator + TriaAccessor<0, dim, spacedim>::quad(const unsigned int) +{ + return typename dealii::internal::TriangulationImplementation:: + Iterators::quad_iterator(); +} + + + +template +inline unsigned int +TriaAccessor<0, dim, spacedim>::quad_index(const unsigned int) +{ + Assert(false, ExcImpossibleInDim(0)); + return numbers::invalid_unsigned_int; +} + + + +template +inline double +TriaAccessor<0, dim, spacedim>::diameter() const +{ + return 0.; +} + + + +template +inline double +TriaAccessor<0, dim, spacedim>::extent_in_direction(const unsigned int) const +{ + return 0.; +} + + + +template +inline Point +TriaAccessor<0, dim, spacedim>::center(const bool, const bool) const +{ + return this->tria->vertices[global_vertex_index]; +} + + + +template +inline double +TriaAccessor<0, dim, spacedim>::measure() const +{ + return 0.; +} + + + +template +inline unsigned char +TriaAccessor<0, dim, spacedim>::combined_face_orientation( + const unsigned int /*face*/) +{ + return 0; +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::face_orientation(const unsigned int /*face*/) +{ + return false; +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::face_flip(const unsigned int /*face*/) +{ + return false; +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::face_rotation(const unsigned int /*face*/) +{ + return false; +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::line_orientation(const unsigned int /*line*/) +{ + return false; +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::has_children() +{ + return false; +} + + + +template +inline unsigned int +TriaAccessor<0, dim, spacedim>::n_children() +{ + return 0; +} + + + +template +inline unsigned int +TriaAccessor<0, dim, spacedim>::n_active_descendants() +{ + return 0; +} + + + +template +inline unsigned int +TriaAccessor<0, dim, spacedim>::max_refinement_depth() +{ + return 0; +} + + + +template +inline unsigned int +TriaAccessor<0, dim, spacedim>::child_iterator_to_index( + const TriaIterator> &) +{ + return numbers::invalid_unsigned_int; +} + + + +template +inline TriaIterator> +TriaAccessor<0, dim, spacedim>::child(const unsigned int) +{ + return TriaIterator>(); +} + + + +template +inline TriaIterator> +TriaAccessor<0, dim, spacedim>::isotropic_child(const unsigned int) +{ + return TriaIterator>(); +} + + + +template +inline RefinementCase<0> +TriaAccessor<0, dim, spacedim>::refinement_case() +{ + return {RefinementPossibilities<0>::no_refinement}; +} + + + +template +inline int +TriaAccessor<0, dim, spacedim>::child_index(const unsigned int) +{ + return -1; +} + + + +template +inline int +TriaAccessor<0, dim, spacedim>::isotropic_child_index(const unsigned int) +{ + return -1; +} + + + +template +inline bool +TriaAccessor<0, dim, spacedim>::used() const +{ + return tria->vertex_used(global_vertex_index); +} + + + +/*------------------- Functions: TriaAccessor<0,1,spacedim> -----------------*/ + +template +inline TriaAccessor<0, 1, spacedim>::TriaAccessor( + const Triangulation<1, spacedim> *tria, + const VertexKind vertex_kind, + const unsigned int vertex_index) + : tria(tria) + , vertex_kind(vertex_kind) + , global_vertex_index(vertex_index) +{} + + + +template +inline TriaAccessor<0, 1, spacedim>::TriaAccessor( + const Triangulation<1, spacedim> *tria, + const int level, + const int index, + const AccessorData *) + : tria(tria) + , vertex_kind(interior_vertex) + , global_vertex_index(numbers::invalid_unsigned_int) +{ + // in general, calling this constructor should yield an error -- users should + // instead call the one immediately above. however, if you create something + // like Triangulation<1>::face_iterator() then this calls the default + // constructor of the iterator which calls the accessor with argument list + // (0,-2,-2,0), so in this particular case accept this call and create an + // object that corresponds to the default constructed (invalid) vertex + // accessor + (void)level; + (void)index; + Assert((level == -2) && (index == -2), + ExcMessage( + "This constructor can not be called for face iterators in 1d, " + "except to default-construct iterator objects.")); +} + + + +template +template +inline TriaAccessor<0, 1, spacedim>::TriaAccessor( + const TriaAccessor &) + : tria(nullptr) + , vertex_kind(interior_vertex) + , global_vertex_index(numbers::invalid_unsigned_int) +{ + Assert(false, ExcImpossibleInDim(0)); +} + + + +template +template +inline TriaAccessor<0, 1, spacedim>::TriaAccessor( + const InvalidAccessor &) + : tria(nullptr) + , vertex_kind(interior_vertex) + , global_vertex_index(numbers::invalid_unsigned_int) +{ + Assert(false, ExcImpossibleInDim(0)); +} + + + +template +inline void +TriaAccessor<0, 1, spacedim>::copy_from(const TriaAccessor &t) +{ + tria = t.tria; + vertex_kind = t.vertex_kind; + global_vertex_index = t.global_vertex_index; +} + + + +template +inline void +TriaAccessor<0, 1, spacedim>::copy_from( + const TriaAccessorBase<0, 1, spacedim> &) +{ + // We cannot convert from TriaAccessorBase to + // TriaAccessor<0,1,spacedim> because the latter is not derived from + // the former. We should never get here. + DEAL_II_ASSERT_UNREACHABLE(); +} + + + +template +inline bool +TriaAccessor<0, 1, spacedim>::operator<( + const TriaAccessor<0, 1, spacedim> &other) const +{ + Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators()); + + return (global_vertex_index < other.global_vertex_index); +} + + + +template +inline IteratorState::IteratorStates +TriaAccessor<0, 1, spacedim>::state() +{ + return IteratorState::valid; +} + + +template +inline int +TriaAccessor<0, 1, spacedim>::level() +{ + return 0; +} + + + +template +inline int +TriaAccessor<0, 1, spacedim>::index() const +{ + return global_vertex_index; +} + + + +template +inline const Triangulation<1, spacedim> & +TriaAccessor<0, 1, spacedim>::get_triangulation() const +{ + return *tria; +} + + + +template +inline void +TriaAccessor<0, 1, spacedim>::operator++() const +{ + DEAL_II_NOT_IMPLEMENTED(); +} + + +template +inline void +TriaAccessor<0, 1, spacedim>::operator--() const +{ + DEAL_II_NOT_IMPLEMENTED(); +} + + + +template +inline bool +TriaAccessor<0, 1, spacedim>::operator==(const TriaAccessor &t) const +{ + const bool result = + ((tria == t.tria) && (global_vertex_index == t.global_vertex_index)); + // if we point to the same vertex, + // make sure we know the same about + // it + if (result == true) + Assert(vertex_kind == t.vertex_kind, ExcInternalError()); + + return result; +} + + + +template +inline bool +TriaAccessor<0, 1, spacedim>::operator!=(const TriaAccessor &t) const +{ + return !(*this == t); +} + + + +template +inline unsigned int +TriaAccessor<0, 1, spacedim>::vertex_index(const unsigned int i) const +{ + AssertIndexRange(i, 1); + (void)i; + return global_vertex_index; +} + + + +template +inline Point & +TriaAccessor<0, 1, spacedim>::vertex(const unsigned int i) const +{ + AssertIndexRange(i, 1); + (void)i; + return const_cast &>( + this->tria->vertices[global_vertex_index]); +} + + + +template +inline Point +TriaAccessor<0, 1, spacedim>::center() const +{ + return this->tria->vertices[global_vertex_index]; +} + + + +template +inline typename dealii::internal::TriangulationImplementation:: + Iterators<1, spacedim>::line_iterator + TriaAccessor<0, 1, spacedim>::line(const unsigned int) +{ + return {}; +} + + +template +inline unsigned int +TriaAccessor<0, 1, spacedim>::line_index(const unsigned int) +{ + Assert(false, ExcImpossibleInDim(0)); + return numbers::invalid_unsigned_int; +} + + +template +inline typename dealii::internal::TriangulationImplementation:: + Iterators<1, spacedim>::quad_iterator + TriaAccessor<0, 1, spacedim>::quad(const unsigned int) +{ + return {}; +} + + + +template +inline unsigned int +TriaAccessor<0, 1, spacedim>::quad_index(const unsigned int) +{ + Assert(false, ExcImpossibleInDim(0)); + return numbers::invalid_unsigned_int; +} + + +template +inline bool +TriaAccessor<0, 1, spacedim>::at_boundary() const +{ + return vertex_kind != interior_vertex; +} + + +template +inline types::boundary_id +TriaAccessor<0, 1, spacedim>::boundary_id() const +{ + switch (vertex_kind) + { + case left_vertex: + case right_vertex: + { + Assert(tria->vertex_to_boundary_id_map_1d->find( + this->vertex_index()) != + tria->vertex_to_boundary_id_map_1d->end(), + ExcInternalError()); + + return (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()]; + } + + default: + return numbers::internal_face_boundary_id; + } +} + + + +template +inline const Manifold<1, spacedim> & +TriaAccessor<0, 1, spacedim>::get_manifold() const +{ + return this->tria->get_manifold(this->manifold_id()); +} + + + +template +inline types::manifold_id +TriaAccessor<0, 1, spacedim>::manifold_id() const +{ + if (tria->vertex_to_manifold_id_map_1d->find(this->vertex_index()) != + tria->vertex_to_manifold_id_map_1d->end()) + return (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()]; + else + return numbers::flat_manifold_id; +} + + +template +inline unsigned char +TriaAccessor<0, 1, spacedim>::combined_face_orientation( + const unsigned int /*face*/) +{ + return 0; +} + + +template +inline bool +TriaAccessor<0, 1, spacedim>::face_orientation(const unsigned int /*face*/) +{ + return false; +} + + + +template +inline bool +TriaAccessor<0, 1, spacedim>::face_flip(const unsigned int /*face*/) +{ + return false; +} + + + +template +inline bool +TriaAccessor<0, 1, spacedim>::face_rotation(const unsigned int /*face*/) +{ + return false; +} + + + +template +inline bool +TriaAccessor<0, 1, spacedim>::line_orientation(const unsigned int /*line*/) +{ + return false; +} + + + +template +inline bool +TriaAccessor<0, 1, spacedim>::has_children() +{ + return false; +} + + + +template +inline unsigned int +TriaAccessor<0, 1, spacedim>::n_children() +{ + return 0; +} + + + +template +inline unsigned int +TriaAccessor<0, 1, spacedim>::n_active_descendants() +{ + return 0; +} + + + +template +inline unsigned int +TriaAccessor<0, 1, spacedim>::max_refinement_depth() +{ + return 0; +} + + + +template +inline unsigned int +TriaAccessor<0, 1, spacedim>::child_iterator_to_index( + const TriaIterator> &) +{ + return numbers::invalid_unsigned_int; +} + + + +template +inline TriaIterator> +TriaAccessor<0, 1, spacedim>::child(const unsigned int) +{ + return TriaIterator>(); +} + + +template +inline TriaIterator> +TriaAccessor<0, 1, spacedim>::isotropic_child(const unsigned int) +{ + return TriaIterator>(); +} + + +template +inline RefinementCase<0> +TriaAccessor<0, 1, spacedim>::refinement_case() +{ + return {RefinementPossibilities<0>::no_refinement}; +} + +template +inline int +TriaAccessor<0, 1, spacedim>::child_index(const unsigned int) +{ + return -1; +} + + +template +inline int +TriaAccessor<0, 1, spacedim>::isotropic_child_index(const unsigned int) +{ + return -1; +} + + + +template +inline void +TriaAccessor<0, 1, spacedim>::set_boundary_id(const types::boundary_id b) const +{ + Assert(tria->vertex_to_boundary_id_map_1d->find(this->vertex_index()) != + tria->vertex_to_boundary_id_map_1d->end(), + ExcMessage("You can't set the boundary_id of a face of a cell that is " + "not actually at the boundary.")); + + (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()] = b; +} + + + +template +inline void +TriaAccessor<0, 1, spacedim>::set_manifold_id(const types::manifold_id b) +{ + (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()] = b; +} + + + +template +inline void +TriaAccessor<0, 1, spacedim>::set_all_boundary_ids( + const types::boundary_id b) const +{ + set_boundary_id(b); +} + + + +template +inline void +TriaAccessor<0, 1, spacedim>::set_all_manifold_ids(const types::manifold_id b) +{ + set_manifold_id(b); +} + + + +template +inline bool +TriaAccessor<0, 1, spacedim>::used() const +{ + return tria->vertex_used(global_vertex_index); +} + + + +template +inline ReferenceCell +TriaAccessor<0, 1, spacedim>::reference_cell() const +{ + return ReferenceCells::Vertex; +} + + + +template +unsigned int +TriaAccessor<0, 1, spacedim>::n_vertices() const +{ + return 1; +} + + + +template +unsigned int +TriaAccessor<0, 1, spacedim>::n_lines() const +{ + return 0; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor<0, 1, spacedim>::vertex_indices() const +{ + return {0U, n_vertices()}; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor<0, 1, spacedim>::line_indices() const +{ + return {0U, n_lines()}; +} + +/*------------------ Functions: CellAccessor ------------------*/ + + +template +inline CellAccessor::CellAccessor( + const Triangulation *parent, + const int level, + const int index, + const AccessorData *local_data) + : TriaAccessor(parent, level, index, local_data) +{} + + + +template +inline CellAccessor::CellAccessor( + const TriaAccessor &cell_accessor) + : TriaAccessor( + static_cast &>(cell_accessor)) +{} + + + +template +inline TriaIterator> +CellAccessor::child(const unsigned int i) const +{ + TriaIterator> q(this->tria, + this->present_level + 1, + this->child_index(i)); + + Assert((q.state() == IteratorState::past_the_end) || q->used(), + ExcInternalError()); + + return q; +} + + + +template +inline boost::container::small_vector>, + GeometryInfo::max_children_per_cell> +CellAccessor::child_iterators() const +{ + boost::container::small_vector>, + GeometryInfo::max_children_per_cell> + child_iterators(this->n_children()); + + for (unsigned int i = 0; i < this->n_children(); ++i) + child_iterators[i] = this->child(i); + + return child_iterators; +} + + + +template +inline TriaIterator> +CellAccessor::face(const unsigned int i) const +{ + AssertIndexRange(i, this->n_faces()); + if constexpr (dim == 1) + { + TriaAccessor<0, 1, spacedim> a( + &this->get_triangulation(), + ((i == 0) && at_boundary(0) ? + TriaAccessor<0, 1, spacedim>::left_vertex : + ((i == 1) && at_boundary(1) ? + TriaAccessor<0, 1, spacedim>::right_vertex : + TriaAccessor<0, 1, spacedim>::interior_vertex)), + internal::TriaAccessorImplementation::Implementation::vertex_index( + *this, i)); + return TriaIterator>(a); + } + else if constexpr (dim == 2) + return this->line(i); + else if constexpr (dim == 3) + return this->quad(i); + else + { + Assert(false, ExcNotImplemented()); + return {}; + } +} + + + +template +inline unsigned int +CellAccessor::face_iterator_to_index( + const TriaIterator> &face) const +{ + for (const unsigned int face_n : this->face_indices()) + if (this->face(face_n) == face) + return face_n; + + Assert(false, + ExcMessage("The given face is not a face of the current cell.")); + return numbers::invalid_unsigned_int; +} + + + +template +inline boost::container::small_vector< + TriaIterator>, + GeometryInfo::faces_per_cell> +CellAccessor::face_iterators() const +{ + boost::container::small_vector< + TriaIterator>, + GeometryInfo::faces_per_cell> + face_iterators(this->n_faces()); + + for (const unsigned int i : this->face_indices()) + face_iterators[i] = this->face(i); + + return face_iterators; +} + + + +template +inline unsigned int +CellAccessor::face_index(const unsigned int i) const +{ + switch (dim) + { + case 1: + { + return this->vertex_index(i); + } + + case 2: + return this->line_index(i); + + case 3: + return this->quad_index(i); + + default: + return numbers::invalid_unsigned_int; + } +} + + + +template +inline int +CellAccessor::neighbor_index(const unsigned int face_no) const +{ + AssertIndexRange(face_no, this->n_faces()); + return this->tria->levels[this->present_level] + ->neighbors[this->present_index * GeometryInfo::faces_per_cell + + face_no] + .second; +} + + + +template +inline int +CellAccessor::neighbor_level(const unsigned int face_no) const +{ + AssertIndexRange(face_no, this->n_faces()); + return this->tria->levels[this->present_level] + ->neighbors[this->present_index * GeometryInfo::faces_per_cell + + face_no] + .first; +} + + + +template +inline RefinementCase +CellAccessor::refine_flag_set() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + // cells flagged for refinement must be active + // (the @p set_refine_flag function checks this, + // but activity may change when refinement is + // executed and for some reason the refine + // flag is not cleared). + Assert(this->is_active() || !this->tria->levels[this->present_level] + ->refine_flags[this->present_index], + ExcRefineCellNotActive()); + return RefinementCase( + this->tria->levels[this->present_level]->refine_flags[this->present_index]); +} + + + +template +inline void +CellAccessor::set_refine_flag( + const RefinementCase refinement_case) const +{ + Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); + Assert(!coarsen_flag_set(), ExcCellFlaggedForCoarsening()); + + this->tria->levels[this->present_level]->refine_flags[this->present_index] = + refinement_case; +} + + + +template +inline void +CellAccessor::clear_refine_flag() const +{ + Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); + this->tria->levels[this->present_level]->refine_flags[this->present_index] = + RefinementCase::no_refinement; +} + + +template +inline std::uint8_t +CellAccessor::refine_choice() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + if (this->tria->levels[this->present_level]->refine_choice.size() == 0) + return 0U; + return this->tria->levels[this->present_level] + ->refine_choice[this->present_index]; +} + + +template +inline void +CellAccessor::set_refine_choice( + const std::uint8_t refinement_choice) const +{ + Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); + if (this->tria->levels[this->present_level]->refine_choice.size() != 0) + this->tria->levels[this->present_level] + ->refine_choice[this->present_index] = refinement_choice; +} + + +template +inline void +CellAccessor::clear_refine_choice() const +{ + Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); + if (this->tria->levels[this->present_level]->refine_choice.size() != 0) + this->tria->levels[this->present_level] + ->refine_choice[this->present_index] = + static_cast(IsotropicRefinementChoice::isotropic_refinement); +} + + +template +inline bool +CellAccessor::flag_for_face_refinement( + const unsigned int face_no, + const RefinementCase &face_refinement_case) const +{ + Assert(dim > 1, ExcImpossibleInDim(dim)); + AssertIndexRange(face_no, this->n_faces()); + AssertIndexRange(face_refinement_case, + RefinementCase::isotropic_refinement + 1); + + // the new refinement case is a combination + // of the minimum required one for the given + // face refinement and the already existing + // flagged refinement case + RefinementCase old_ref_case = refine_flag_set(); + RefinementCase new_ref_case = + (old_ref_case | + GeometryInfo::min_cell_refinement_case_for_face_refinement( + face_refinement_case, + face_no, + this->face_orientation(face_no), + this->face_flip(face_no), + this->face_rotation(face_no))); + set_refine_flag(new_ref_case); + // return, whether we had to change the + // refinement flag + return new_ref_case != old_ref_case; +} + + + +template +inline bool +CellAccessor::flag_for_line_refinement( + const unsigned int line_no) const +{ + Assert(dim > 1, ExcImpossibleInDim(dim)); + AssertIndexRange(line_no, this->n_lines()); + + // the new refinement case is a combination + // of the minimum required one for the given + // line refinement and the already existing + // flagged refinement case + RefinementCase + old_ref_case = refine_flag_set(), + new_ref_case = + old_ref_case | + GeometryInfo::min_cell_refinement_case_for_line_refinement(line_no); + set_refine_flag(new_ref_case); + // return, whether we had to change the + // refinement flag + return new_ref_case != old_ref_case; +} + + + +template <> +inline dealii::internal::SubfaceCase<1> +CellAccessor<1>::subface_case(const unsigned int) const +{ + return dealii::internal::SubfaceCase<1>::case_none; +} + +template <> +inline dealii::internal::SubfaceCase<1> +CellAccessor<1, 2>::subface_case(const unsigned int) const +{ + return dealii::internal::SubfaceCase<1>::case_none; +} + + +template <> +inline dealii::internal::SubfaceCase<1> +CellAccessor<1, 3>::subface_case(const unsigned int) const +{ + return dealii::internal::SubfaceCase<1>::case_none; +} + + +template <> +inline dealii::internal::SubfaceCase<2> +CellAccessor<2>::subface_case(const unsigned int face_no) const +{ + Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); + AssertIndexRange(face_no, this->n_faces()); + return ((face(face_no)->has_children()) ? + dealii::internal::SubfaceCase<2>::case_x : + dealii::internal::SubfaceCase<2>::case_none); +} + +template <> +inline dealii::internal::SubfaceCase<2> +CellAccessor<2, 3>::subface_case(const unsigned int face_no) const +{ + Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); + AssertIndexRange(face_no, this->n_faces()); + return ((face(face_no)->has_children()) ? + dealii::internal::SubfaceCase<2>::case_x : + dealii::internal::SubfaceCase<2>::case_none); +} + + +template <> +inline dealii::internal::SubfaceCase<3> +CellAccessor<3>::subface_case(const unsigned int face_no) const +{ + Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); + AssertIndexRange(face_no, this->n_faces()); + switch (static_cast(face(face_no)->refinement_case())) + { + case RefinementCase<3>::no_refinement: + return dealii::internal::SubfaceCase<3>::case_none; + case RefinementCase<3>::cut_x: + if (face(face_no)->child(0)->has_children()) + { + Assert(face(face_no)->child(0)->refinement_case() == + RefinementCase<2>::cut_y, + ExcInternalError()); + if (face(face_no)->child(1)->has_children()) + { + Assert(face(face_no)->child(1)->refinement_case() == + RefinementCase<2>::cut_y, + ExcInternalError()); + return dealii::internal::SubfaceCase<3>::case_x1y2y; + } + else + return dealii::internal::SubfaceCase<3>::case_x1y; + } + else + { + if (face(face_no)->child(1)->has_children()) + { + Assert(face(face_no)->child(1)->refinement_case() == + RefinementCase<2>::cut_y, + ExcInternalError()); + return dealii::internal::SubfaceCase<3>::case_x2y; + } + else + return dealii::internal::SubfaceCase<3>::case_x; + } + case RefinementCase<3>::cut_y: + if (face(face_no)->child(0)->has_children()) + { + Assert(face(face_no)->child(0)->refinement_case() == + RefinementCase<2>::cut_x, + ExcInternalError()); + if (face(face_no)->child(1)->has_children()) + { + Assert(face(face_no)->child(1)->refinement_case() == + RefinementCase<2>::cut_x, + ExcInternalError()); + return dealii::internal::SubfaceCase<3>::case_y1x2x; + } + else + return dealii::internal::SubfaceCase<3>::case_y1x; + } + else + { + if (face(face_no)->child(1)->has_children()) + { + Assert(face(face_no)->child(1)->refinement_case() == + RefinementCase<2>::cut_x, + ExcInternalError()); + return dealii::internal::SubfaceCase<3>::case_y2x; + } + else + return dealii::internal::SubfaceCase<3>::case_y; + } + case RefinementCase<3>::cut_xy: + return dealii::internal::SubfaceCase<3>::case_xy; + default: + DEAL_II_ASSERT_UNREACHABLE(); + } + // we should never get here + return dealii::internal::SubfaceCase<3>::case_none; +} + + + +template +inline bool +CellAccessor::coarsen_flag_set() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + // cells flagged for coarsening must be active + // (the @p set_refine_flag function checks this, + // but activity may change when refinement is + // executed and for some reason the refine + // flag is not cleared). + Assert(this->is_active() || !this->tria->levels[this->present_level] + ->coarsen_flags[this->present_index], + ExcRefineCellNotActive()); + return this->tria->levels[this->present_level] + ->coarsen_flags[this->present_index]; +} + + + +template +inline void +CellAccessor::set_coarsen_flag() const +{ + Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); + Assert(!refine_flag_set(), ExcCellFlaggedForRefinement()); + + this->tria->levels[this->present_level]->coarsen_flags[this->present_index] = + true; +} + + + +template +inline void +CellAccessor::clear_coarsen_flag() const +{ + Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); + this->tria->levels[this->present_level]->coarsen_flags[this->present_index] = + false; +} + + + +template +inline TriaIterator> +CellAccessor::neighbor(const unsigned int face_no) const +{ + TriaIterator> q(this->tria, + neighbor_level(face_no), + neighbor_index(face_no)); + + Assert((q.state() == IteratorState::past_the_end) || q->used(), + ExcInternalError()); + + return q; +} + + + +template +inline bool +CellAccessor::is_active() const +{ + return !this->has_children(); +} + + + +template +inline bool +CellAccessor::is_locally_owned() const +{ + Assert(this->is_active(), + ExcMessage("is_locally_owned() can only be called on active cells!")); +#ifndef DEAL_II_WITH_MPI + return true; +#else + + // Serial triangulations report invalid_subdomain_id as their locally owned + // subdomain, so the first condition checks whether we have a serial + // triangulation, in which case all cells are locally owned. The second + // condition compares the subdomain id in the parallel case. + const types::subdomain_id locally_owned_subdomain = + this->tria->locally_owned_subdomain(); + return (locally_owned_subdomain == numbers::invalid_subdomain_id || + this->subdomain_id() == locally_owned_subdomain); + +#endif +} + + +template +inline bool +CellAccessor::is_locally_owned_on_level() const +{ +#ifndef DEAL_II_WITH_MPI + return true; +#else + + // Serial triangulations report invalid_subdomain_id as their locally owned + // subdomain, so the first condition checks whether we have a serial + // triangulation, in which case all cells are locally owned. The second + // condition compares the subdomain id in the parallel case. + const types::subdomain_id locally_owned_subdomain = + this->tria->locally_owned_subdomain(); + return (locally_owned_subdomain == numbers::invalid_subdomain_id || + this->level_subdomain_id() == locally_owned_subdomain); + +#endif +} + + +template +inline bool +CellAccessor::is_ghost() const +{ + Assert(this->is_active(), + ExcMessage("is_ghost() can only be called on active cells!")); + if (this->has_children()) + return false; + +#ifndef DEAL_II_WITH_MPI + return false; +#else + + // Serial triangulations report invalid_subdomain_id as their locally owned + // subdomain, so the first condition rules out that case as all cells to a + // serial triangulation are locally owned and none is ghosted. The second + // and third conditions check whether the cell's subdomain is not the + // locally owned one and not artificial. + const types::subdomain_id locally_owned_subdomain = + this->tria->locally_owned_subdomain(); + const types::subdomain_id subdomain_id = this->subdomain_id(); + return (locally_owned_subdomain != numbers::invalid_subdomain_id && + subdomain_id != locally_owned_subdomain && + subdomain_id != numbers::artificial_subdomain_id); + +#endif +} + + +template +inline bool +CellAccessor::is_ghost_on_level() const +{ +#ifndef DEAL_II_WITH_MPI + return false; +#else + + // Serial triangulations report invalid_subdomain_id as their locally owned + // subdomain, so the first condition checks whether we have a serial + // triangulation, in which case all cells are locally owned. The second + // condition compares the subdomain id in the parallel case. + const types::subdomain_id locally_owned_subdomain = + this->tria->locally_owned_subdomain(); + const types::subdomain_id subdomain_id = this->level_subdomain_id(); + return (locally_owned_subdomain != numbers::invalid_subdomain_id && + subdomain_id != locally_owned_subdomain && + subdomain_id != numbers::artificial_subdomain_id); + +#endif +} + + + +template +inline bool +CellAccessor::is_artificial() const +{ + Assert(this->is_active(), + ExcMessage("is_artificial() can only be called on active cells!")); +#ifndef DEAL_II_WITH_MPI + return false; +#else + + // Serial triangulations report invalid_subdomain_id as their locally owned + // subdomain, so the first condition rules out that case as all cells to a + // serial triangulation are locally owned and none is artificial. + return (this->tria->locally_owned_subdomain() != + numbers::invalid_subdomain_id && + this->subdomain_id() == numbers::artificial_subdomain_id); + +#endif +} + + + +template +inline bool +CellAccessor::is_artificial_on_level() const +{ +#ifndef DEAL_II_WITH_MPI + return false; +#else + return (this->tria->locally_owned_subdomain() != + numbers::invalid_subdomain_id && + this->level_subdomain_id() == numbers::artificial_subdomain_id); +#endif +} + + + +template +inline types::subdomain_id +CellAccessor::subdomain_id() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + Assert(this->is_active(), + ExcMessage("subdomain_id() can only be called on active cells!")); + return this->tria->levels[this->present_level] + ->subdomain_ids[this->present_index]; +} + + + +template +inline types::subdomain_id +CellAccessor::level_subdomain_id() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + return this->tria->levels[this->present_level] + ->level_subdomain_ids[this->present_index]; +} + + + +template +inline unsigned int +CellAccessor::neighbor_face_no(const unsigned int neighbor) const +{ + const unsigned int n2 = neighbor_of_neighbor_internal(neighbor); + if (n2 != numbers::invalid_unsigned_int) + // return this value as the + // neighbor is not coarser + return n2; + else + // the neighbor is coarser + return neighbor_of_coarser_neighbor(neighbor).first; +} + + + +template +inline bool +CellAccessor::is_level_cell() +{ + return false; +} + + + +template +inline unsigned int +CellAccessor::active_cell_index() const +{ + Assert(this->is_active(), TriaAccessorExceptions::ExcCellNotActive()); + return this->tria->levels[this->present_level] + ->active_cell_indices[this->present_index]; +} + + + +template +inline types::global_cell_index +CellAccessor::global_active_cell_index() const +{ + Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); + Assert(this->is_active(), + ExcMessage( + "global_active_cell_index() can only be called on active cells!")); + + return this->tria->levels[this->present_level] + ->global_active_cell_indices[this->present_index]; +} + + + +template +inline types::global_cell_index +CellAccessor::global_level_cell_index() const +{ + return this->tria->levels[this->present_level] + ->global_level_cell_indices[this->present_index]; +} + DEAL_II_NAMESPACE_CLOSE -// include more templates in debug and optimized mode -#include #endif diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index 686e9802eb..6ab9563bd4 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -15,4027 +15,7 @@ #ifndef dealii_tria_accessor_templates_h #define dealii_tria_accessor_templates_h - -#include - -#include -#include - -#include -#include -#include -#include -#include -#include - -#include - -#include -#include - -DEAL_II_NAMESPACE_OPEN - - -namespace internal -{ - namespace TriaAccessorImplementation - { - /** - * Compute the diameter for a given set of vertices. The vertices are - * interpreted, depending on their count, as the vertices of a particular - * reference cell. - */ - template - inline double - diameter( - const boost::container::small_vector, - GeometryInfo::vertices_per_cell> - vertices) - { - const ReferenceCell reference_cell = - ReferenceCell::n_vertices_to_type(dim, vertices.size()); - - if (reference_cell == ReferenceCells::Line) - // Return the distance between the two vertices - return (vertices[1] - vertices[0]).norm(); - else if (reference_cell == ReferenceCells::Triangle) - // Return the longest of the three edges - return std::max({(vertices[1] - vertices[0]).norm(), - (vertices[2] - vertices[1]).norm(), - (vertices[2] - vertices[0]).norm()}); - else if (reference_cell == ReferenceCells::Quadrilateral) - // Return the longer one of the two diagonals of the quadrilateral - return std::max({(vertices[3] - vertices[0]).norm(), - (vertices[2] - vertices[1]).norm()}); - else if (reference_cell == ReferenceCells::Tetrahedron) - // Return the longest of the six edges of the tetrahedron - return std::max({(vertices[1] - vertices[0]).norm(), - (vertices[2] - vertices[0]).norm(), - (vertices[2] - vertices[1]).norm(), - (vertices[3] - vertices[0]).norm(), - (vertices[3] - vertices[1]).norm(), - (vertices[3] - vertices[2]).norm()}); - else if (reference_cell == ReferenceCells::Pyramid) - // Return ... - return std::max({// the longest diagonal of the quadrilateral base - // of the pyramid or ... - (vertices[3] - vertices[0]).norm(), - (vertices[2] - vertices[1]).norm(), - // the longest edge connected with the apex of the - // pyramid - (vertices[4] - vertices[0]).norm(), - (vertices[4] - vertices[1]).norm(), - (vertices[4] - vertices[2]).norm(), - (vertices[4] - vertices[3]).norm()}); - else if (reference_cell == ReferenceCells::Wedge) - // Return ... - return std::max({// the longest of the 2*3=6 diagonals of the three - // quadrilateral sides of the wedge or ... - (vertices[4] - vertices[0]).norm(), - (vertices[3] - vertices[1]).norm(), - (vertices[5] - vertices[1]).norm(), - (vertices[4] - vertices[2]).norm(), - (vertices[5] - vertices[0]).norm(), - (vertices[3] - vertices[2]).norm(), - // the longest of the 3*2=6 edges of the two - // triangular faces of the wedge - (vertices[1] - vertices[0]).norm(), - (vertices[2] - vertices[1]).norm(), - (vertices[2] - vertices[0]).norm(), - (vertices[4] - vertices[3]).norm(), - (vertices[5] - vertices[4]).norm(), - (vertices[5] - vertices[3]).norm()}); - else if (reference_cell == ReferenceCells::Hexahedron) - // Return the longest of the four diagonals of the hexahedron - return std::max({(vertices[7] - vertices[0]).norm(), - (vertices[6] - vertices[1]).norm(), - (vertices[2] - vertices[5]).norm(), - (vertices[3] - vertices[4]).norm()}); - - DEAL_II_NOT_IMPLEMENTED(); - return -1e10; - } - } // namespace TriaAccessorImplementation -} // namespace internal - - -/*--------------------- Functions: TriaAccessorBase -------------------------*/ - -template -inline TriaAccessorBase::TriaAccessorBase( - const Triangulation *tria, - const int level, - const int index, - const AccessorData *) - : present_level((structdim == dim) ? level : 0) - , present_index(index) - , tria(tria) -{ - // non-cells have no level, so a 0 - // should have been passed, or a -1 - // for an end-iterator, or -2 for - // an invalid (default constructed) - // iterator - if (structdim != dim) - { - Assert((level == 0) || (level == -1) || (level == -2), - ExcInternalError()); - } -} - - -template -inline TriaAccessorBase::TriaAccessorBase( - const TriaAccessorBase &a) - : present_level(a.present_level) - , present_index(a.present_index) - , tria(a.tria) -{} - - -template -inline void -TriaAccessorBase::copy_from( - const TriaAccessorBase &a) -{ - present_level = a.present_level; - present_index = a.present_index; - tria = a.tria; - - if (structdim != dim) - { - Assert((present_level == 0) || (present_level == -1) || - (present_level == -2), - ExcInternalError()); - } -} - - - -template -inline TriaAccessorBase & -TriaAccessorBase::operator=( - const TriaAccessorBase &a) -{ - present_level = a.present_level; - present_index = a.present_index; - tria = a.tria; - - if (structdim != dim) - { - Assert((present_level == 0) || (present_level == -1) || - (present_level == -2), - ExcInternalError()); - } - return *this; -} - - - -template -inline bool -TriaAccessorBase::operator==( - const TriaAccessorBase &a) const -{ - Assert(tria == a.tria || tria == nullptr || a.tria == nullptr, - TriaAccessorExceptions::ExcCantCompareIterators()); - return ((tria == a.tria) && (present_level == a.present_level) && - (present_index == a.present_index)); -} - - - -template -inline bool -TriaAccessorBase::operator!=( - const TriaAccessorBase &a) const -{ - Assert(tria == a.tria || tria == nullptr || a.tria == nullptr, - TriaAccessorExceptions::ExcCantCompareIterators()); - return ((tria != a.tria) || (present_level != a.present_level) || - (present_index != a.present_index)); -} - - - -template -inline bool -TriaAccessorBase::operator<( - const TriaAccessorBase &other) const -{ - Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators()); - - if (present_level != other.present_level) - return (present_level < other.present_level); - - return (present_index < other.present_index); -} - - - -template -inline int -TriaAccessorBase::level() const -{ - // This is always zero or invalid - // if the object is not a cell - return present_level; -} - - - -template -inline int -TriaAccessorBase::index() const -{ - return present_index; -} - - - -template -inline IteratorState::IteratorStates -TriaAccessorBase::state() const -{ - if ((present_level >= 0) && (present_index >= 0)) - return IteratorState::valid; - else if (present_index == -1) - return IteratorState::past_the_end; - else - return IteratorState::invalid; -} - - - -template -inline const Triangulation & -TriaAccessorBase::get_triangulation() const -{ - return *tria; -} - - - -template -inline void -TriaAccessorBase::operator++() -{ - // this iterator is used for - // objects without level - ++this->present_index; - - if (structdim != dim) - { - // is index still in the range of - // the vector? (note that we don't - // have to set the level, since - // dim!=1 and the object therefore - // has no level) - if (this->present_index >= static_cast(objects().n_objects())) - this->present_index = -1; - } - else - { - while (this->present_index >= - static_cast( - this->tria->levels[this->present_level]->cells.n_objects())) - { - // no -> go one level up until we find - // one with more than zero cells - ++this->present_level; - this->present_index = 0; - // highest level reached? - if (this->present_level >= - static_cast(this->tria->levels.size())) - { - // return with past the end pointer - this->present_level = this->present_index = -1; - return; - } - } - } -} - - -template -inline void -TriaAccessorBase::operator--() -{ - // same as operator++ - --this->present_index; - - if (structdim != dim) - { - if (this->present_index < 0) - this->present_index = -1; - } - else - { - while (this->present_index < 0) - { - // no -> go one level down - --this->present_level; - // lowest level reached? - if (this->present_level == -1) - { - // return with past the end pointer - this->present_level = this->present_index = -1; - return; - } - // else - this->present_index = - this->tria->levels[this->present_level]->cells.n_objects() - 1; - } - } -} - - - -template -inline dealii::internal::TriangulationImplementation::TriaObjects & -TriaAccessorBase::objects() const -{ - if (structdim == dim) - return this->tria->levels[this->present_level]->cells; - - if (structdim == 1 && dim > 1) - return this->tria->faces->lines; - - if (structdim == 2 && dim > 2) - return this->tria->faces->quads; - - DEAL_II_ASSERT_UNREACHABLE(); - - return this->tria->levels[this->present_level]->cells; -} - - - -/*---------------------- Functions: InvalidAccessor -------------------------*/ - -template -InvalidAccessor::InvalidAccessor(const void *, - const int, - const int, - const AccessorData *) -{ - Assert(false, - ExcMessage("You are attempting an invalid conversion between " - "iterator/accessor types. The constructor you call " - "only exists to make certain template constructs " - "easier to write as dimension independent code but " - "the conversion is not valid in the current context.")); -} - - - -template -InvalidAccessor::InvalidAccessor( - const InvalidAccessor &) -{ - Assert(false, - ExcMessage("You are attempting an invalid conversion between " - "iterator/accessor types. The constructor you call " - "only exists to make certain template constructs " - "easier to write as dimension independent code but " - "the conversion is not valid in the current context.")); -} - - - -template -void -InvalidAccessor::copy_from(const InvalidAccessor &) -{ - // nothing to do here. we could - // throw an exception but we can't - // get here without first creating - // an object which would have - // already thrown -} - - - -template -bool -InvalidAccessor::operator==( - const InvalidAccessor &) const -{ - // nothing to do here. we could - // throw an exception but we can't - // get here without first creating - // an object which would have - // already thrown - return false; -} - - - -template -bool -InvalidAccessor::operator!=( - const InvalidAccessor &) const -{ - // nothing to do here. we could - // throw an exception but we can't - // get here without first creating - // an object which would have - // already thrown - return true; -} - - - -template -bool -InvalidAccessor::used() const -{ - // nothing to do here. we could - // throw an exception but we can't - // get here without first creating - // an object which would have - // already thrown - return false; -} - - - -template -bool -InvalidAccessor::has_children() const -{ - // nothing to do here. we could - // throw an exception but we can't - // get here without first creating - // an object which would have - // already thrown - return false; -} - - - -template -void -InvalidAccessor::operator++() const -{} - - - -template -void -InvalidAccessor::operator--() const -{} - - - -template -types::manifold_id -InvalidAccessor::manifold_id() const -{ - return numbers::flat_manifold_id; -} - - - -template -unsigned int -InvalidAccessor::user_index() const -{ - return numbers::invalid_unsigned_int; -} - - - -template -void -InvalidAccessor::set_user_index( - const unsigned int) const -{ - Assert(false, - ExcMessage("You are trying to set the user index of an " - "invalid object.")); -} - - - -template -void -InvalidAccessor::set_manifold_id( - const types::manifold_id) const -{ - Assert(false, - ExcMessage("You are trying to set the manifold id of an " - "invalid object.")); -} - - - -template -inline Point & -InvalidAccessor::vertex(const unsigned int) const -{ - // nothing to do here. we could throw an exception but we can't get here - // without first creating an object which would have already thrown - static Point invalid_vertex; - return invalid_vertex; -} - - -template -inline void * -InvalidAccessor::line(const unsigned int) const -{ - // nothing to do here. we could throw an exception but we can't get here - // without first creating an object which would have already thrown - return nullptr; -} - - - -template -inline void * -InvalidAccessor::quad(const unsigned int) const -{ - // nothing to do here. we could throw an exception but we can't get here - // without first creating an object which would have already thrown - return nullptr; -} - - -/*------------------------ Functions: TriaAccessor ---------------------------*/ - - -namespace internal -{ - namespace TriaAccessorImplementation - { - // make sure that if in the following we - // write TriaAccessor - // we mean the *class* - // dealii::TriaAccessor, not the - // enclosing namespace - // dealii::internal::TriaAccessor - using dealii::TriaAccessor; - - /** - * A class with the same purpose as the similarly named class of the - * Triangulation class. See there for more information. - */ - struct Implementation - { - /** - * Implementation of the function of some name in the parent class. - */ - template - inline static unsigned int - line_index(const TriaAccessor<1, dim, spacedim> &, const unsigned int) - { - Assert(false, - ExcMessage("You can't ask for the index of a line bounding " - "a one-dimensional cell because it is not " - "bounded by lines.")); - return numbers::invalid_unsigned_int; - } - - - template - inline static unsigned int - line_index(const TriaAccessor<2, dim, spacedim> &accessor, - const unsigned int i) - { - constexpr unsigned int max_faces_per_cell = 4; - return accessor.objects() - .cells[accessor.present_index * max_faces_per_cell + i]; - } - - - inline static unsigned int - line_index(const TriaAccessor<3, 3, 3> &accessor, const unsigned int i) - { - const auto [face_index, line_index] = - accessor.reference_cell().standard_line_to_face_and_line_index(i); - const auto line_within_face_index = - accessor.reference_cell().standard_to_real_face_line( - line_index, - face_index, - accessor.combined_face_orientation(face_index)); - - return accessor.quad(face_index)->line_index(line_within_face_index); - } - - - - /** - * Implementation of the function of some name in the parent class. - */ - template - inline static unsigned int - quad_index(const TriaAccessor &, - const unsigned int) - { - Assert(false, - ExcMessage("You can't ask for the index of a quad bounding " - "a one- or two-dimensional cell because it is not " - "bounded by quads.")); - return numbers::invalid_unsigned_int; - } - - - inline static unsigned int - quad_index(const TriaAccessor<3, 3, 3> &accessor, const unsigned int i) - { - constexpr unsigned int max_faces_per_cell = 6; - return accessor.tria->levels[accessor.present_level] - ->cells.cells[accessor.present_index * max_faces_per_cell + i]; - } - - - - /** - * Implementation of the function of some name in the parent class. - */ - template - inline static void - set_combined_face_orientation( - const TriaAccessor &accessor, - const unsigned int face, - const unsigned char combined_orientation) - { - Assert(structdim == dim, - ExcMessage("This function can only be used on objects that are " - "cells and not on objects which bound cells.")); - AssertIndexRange(face, accessor.n_faces()); - - if (dim == 1) - Assert(combined_orientation == - ReferenceCell::default_combined_face_orientation(), - ExcMessage("In 1d, faces do not have an orientation, so the " - "only valid value is the default.")); - else if (dim == 2) - Assert(combined_orientation == - ReferenceCell::default_combined_face_orientation() || - combined_orientation == - ReferenceCell::reversed_combined_line_orientation(), - ExcMessage( - "In 2d, the only valid values of the combined orientation " - "are the standard orientation or the reversed line " - "orientation.")); - - // face_orientations is not set up in 1d - if (dim != 1) - accessor.tria->levels[accessor.present_level] - ->face_orientations.set_combined_orientation( - accessor.present_index * GeometryInfo::faces_per_cell + face, - combined_orientation); - } - - /** - * Implementation of the function of some name in the parent class. - */ - template - inline static void - set_line_orientation(const TriaAccessor<1, dim, spacedim> &, - const unsigned int, - const bool) - { - DEAL_II_ASSERT_UNREACHABLE(); - } - - - template - inline static void - set_line_orientation(const TriaAccessor<2, 2, spacedim> &, - const unsigned int, - const bool) - { - // quads in 2d have no - // non-standard orientation - DEAL_II_ASSERT_UNREACHABLE(); - } - - - template - inline static void - set_line_orientation(const TriaAccessor<2, 3, spacedim> &accessor, - const unsigned int line, - const bool value) - { - Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(line, accessor.n_lines()); - Assert(accessor.present_index * GeometryInfo<3>::lines_per_face + line < - accessor.tria->faces->quads_line_orientations.size(), - ExcInternalError()); - - // quads as part of 3d hexes can have non-standard orientation - accessor.tria->faces->quads_line_orientations - [accessor.present_index * GeometryInfo<3>::lines_per_face + line] = - value; - } - - - inline static void - set_line_orientation(const TriaAccessor<3, 3, 3> &, - const unsigned int, - const bool) - { - // it seems like we don't need this - // one - DEAL_II_NOT_IMPLEMENTED(); - } - - - /** - * Implementation of the function of same name in the enclosing class. - */ - template - inline static unsigned int - vertex_index(const TriaAccessor<1, dim, spacedim> &accessor, - const unsigned int corner) - { - constexpr unsigned int max_faces_per_cell = 2; - return accessor.objects() - .cells[accessor.present_index * max_faces_per_cell + corner]; - } - - - template - inline static unsigned int - vertex_index(const TriaAccessor<2, dim, spacedim> &accessor, - const unsigned int corner) - { - const auto [line_index, vertex_index] = - accessor.reference_cell().standard_vertex_to_face_and_vertex_index( - corner); - const auto vertex_within_line_index = - accessor.reference_cell().standard_to_real_face_vertex( - vertex_index, - line_index, - accessor.line_orientation(line_index) == true ? - ReferenceCell::default_combined_face_orientation() : - ReferenceCell::reversed_combined_line_orientation()); - - return accessor.line(line_index) - ->vertex_index(vertex_within_line_index); - } - - - - inline static unsigned int - vertex_index(const TriaAccessor<3, 3, 3> &accessor, - const unsigned int corner) - { - const auto [face_index, vertex_index] = - accessor.reference_cell().standard_vertex_to_face_and_vertex_index( - corner); - const auto vertex_within_face_index = - accessor.reference_cell().standard_to_real_face_vertex( - vertex_index, - face_index, - accessor.combined_face_orientation(face_index)); - - return accessor.quad(face_index) - ->vertex_index(vertex_within_face_index); - } - - - - template - static std::array - get_line_indices_of_cell(const TriaAccessor<1, dim, spacedim> &) - { - DEAL_II_ASSERT_UNREACHABLE(); - return {}; - } - - - - template - static std::array - get_line_indices_of_cell(const TriaAccessor<2, dim, spacedim> &cell) - { - // For 2d cells the access cell->line_orientation() is already - // efficient - std::array line_indices = {}; - for (const unsigned int line : cell.line_indices()) - line_indices[line] = cell.line_index(line); - return line_indices; - } - - /** - * A helper function to provide faster access to cell->line_index() in - * 3d - */ - template - static std::array - get_line_indices_of_cell( - const TriaAccessor &cell) - { - std::array line_indices = {}; - - // For hexahedra, the classical access via quads -> lines is too - // inefficient. Unroll this code here to allow the compiler to inline - // the necessary functions. - const auto ref_cell = cell.reference_cell(); - if (ref_cell == ReferenceCells::Hexahedron) - { - for (unsigned int f = 4; f < 6; ++f) - { - const unsigned char orientation = - cell.get_triangulation() - .levels[cell.level()] - ->face_orientations.get_combined_orientation( - cell.index() * GeometryInfo<3>::faces_per_cell + f); - - // It might seem superfluous to spell out the four indices - // that get later consumed by a for loop over these four - // elements; however, for the compiler it is easier to inline - // the statement of standard_to_real_face_line() when next to - // each other, as opposed to be interleaved with a - // line_index() call. - const std::array my_indices{ - {ref_cell.standard_to_real_face_line(0, f, orientation), - ref_cell.standard_to_real_face_line(1, f, orientation), - ref_cell.standard_to_real_face_line(2, f, orientation), - ref_cell.standard_to_real_face_line(3, f, orientation)}}; - const auto quad = cell.quad(f); - for (unsigned int l = 0; l < 4; ++l) - line_indices[4 * (f - 4) + l] = - quad->line_index(my_indices[l]); - } - for (unsigned int f = 0; f < 2; ++f) - { - const unsigned char orientation = - cell.get_triangulation() - .levels[cell.level()] - ->face_orientations.get_combined_orientation( - cell.index() * GeometryInfo<3>::faces_per_cell + f); - const std::array my_indices{ - {ref_cell.standard_to_real_face_line(0, f, orientation), - ref_cell.standard_to_real_face_line(1, f, orientation)}}; - const auto quad = cell.quad(f); - line_indices[8 + f] = quad->line_index(my_indices[0]); - line_indices[10 + f] = quad->line_index(my_indices[1]); - } - } - else if (ref_cell == ReferenceCells::Tetrahedron) - { - std::array orientations{ - {cell.combined_face_orientation(0), - cell.combined_face_orientation(1), - cell.combined_face_orientation(2)}}; - const std::array my_indices{ - {ref_cell.standard_to_real_face_line(0, 0, orientations[0]), - ref_cell.standard_to_real_face_line(1, 0, orientations[0]), - ref_cell.standard_to_real_face_line(2, 0, orientations[0]), - ref_cell.standard_to_real_face_line(1, 1, orientations[1]), - ref_cell.standard_to_real_face_line(2, 1, orientations[1]), - ref_cell.standard_to_real_face_line(1, 2, orientations[2])}}; - line_indices[0] = cell.quad(0)->line_index(my_indices[0]); - line_indices[1] = cell.quad(0)->line_index(my_indices[1]); - line_indices[2] = cell.quad(0)->line_index(my_indices[2]); - line_indices[3] = cell.quad(1)->line_index(my_indices[3]); - line_indices[4] = cell.quad(1)->line_index(my_indices[4]); - line_indices[5] = cell.quad(2)->line_index(my_indices[5]); - } - else - // For other shapes (wedges, pyramids), we do not currently - // implement an optimized function. - for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l) - line_indices[l] = cell.line_index(l); - - return line_indices; - } - - - - /** - * A helper function to provide faster access to - * cell->line_orientation(), 1d specialization - */ - template - static std::array - get_line_orientations_of_cell(const TriaAccessor<1, dim, spacedim> &) - { - DEAL_II_ASSERT_UNREACHABLE(); - return {}; - } - - - - /** - * A helper function to provide faster access to - * cell->line_orientation(), 2d specialization - */ - template - static std::array - get_line_orientations_of_cell(const TriaAccessor<2, dim, spacedim> &cell) - { - // For 2d cells the access cell->line_orientation() is already - // efficient - std::array line_orientations = {}; - for (const unsigned int line : cell.line_indices()) - line_orientations[line] = - cell.line_orientation(line) == true ? - ReferenceCell::default_combined_face_orientation() : - ReferenceCell::reversed_combined_line_orientation(); - return line_orientations; - } - - - - /** - * A helper function to provide faster access to - * cell->line_orientation(), 3d specialization - */ - template - static std::array - get_line_orientations_of_cell(const TriaAccessor<3, dim, spacedim> &cell) - { - std::array line_orientations = {}; - - // For hexahedra, the classical access via quads -> lines is too - // inefficient. Unroll this code here to allow the compiler to inline - // the necessary functions. - const auto ref_cell = cell.reference_cell(); - if (ref_cell == ReferenceCells::Hexahedron) - { - for (unsigned int f = 4; f < 6; ++f) - { - const unsigned char orientation = - cell.get_triangulation() - .levels[cell.level()] - ->face_orientations.get_combined_orientation( - cell.index() * GeometryInfo<3>::faces_per_cell + f); - - // It might seem superfluous to spell out the four indices and - // orientations that get later consumed by a for loop over - // these four elements; however, for the compiler it is easier - // to inline the statement of standard_to_real_face_line() - // when next to each other, as opposed to be interleaved with - // a line_index() call. - const std::array my_indices{ - {ref_cell.standard_to_real_face_line(0, f, orientation), - ref_cell.standard_to_real_face_line(1, f, orientation), - ref_cell.standard_to_real_face_line(2, f, orientation), - ref_cell.standard_to_real_face_line(3, f, orientation)}}; - const auto quad = cell.quad(f); - const std::array my_orientations{ - {ref_cell.standard_vs_true_line_orientation( - 0, f, orientation, quad->line_orientation(my_indices[0])), - ref_cell.standard_vs_true_line_orientation( - 1, f, orientation, quad->line_orientation(my_indices[1])), - ref_cell.standard_vs_true_line_orientation( - 2, f, orientation, quad->line_orientation(my_indices[2])), - ref_cell.standard_vs_true_line_orientation( - 3, - f, - orientation, - quad->line_orientation(my_indices[3]))}}; - for (unsigned int l = 0; l < 4; ++l) - line_orientations[4 * (f - 4) + l] = my_orientations[l]; - } - for (unsigned int f = 0; f < 2; ++f) - { - const unsigned char orientation = - cell.get_triangulation() - .levels[cell.level()] - ->face_orientations.get_combined_orientation( - cell.index() * GeometryInfo<3>::faces_per_cell + f); - const std::array my_indices{ - {ref_cell.standard_to_real_face_line(0, f, orientation), - ref_cell.standard_to_real_face_line(1, f, orientation)}}; - const auto quad = cell.quad(f); - const std::array my_orientations{ - {ref_cell.standard_vs_true_line_orientation( - 0, f, orientation, quad->line_orientation(my_indices[0])), - ref_cell.standard_vs_true_line_orientation( - 1, - f, - orientation, - quad->line_orientation(my_indices[1]))}}; - line_orientations[8 + f] = my_orientations[0]; - line_orientations[10 + f] = my_orientations[1]; - } - } - else if (ref_cell == ReferenceCells::Tetrahedron) - { - std::array orientations{ - {cell.combined_face_orientation(0), - cell.combined_face_orientation(1), - cell.combined_face_orientation(2)}}; - const std::array my_indices{ - {ref_cell.standard_to_real_face_line(0, 0, orientations[0]), - ref_cell.standard_to_real_face_line(1, 0, orientations[0]), - ref_cell.standard_to_real_face_line(2, 0, orientations[0]), - ref_cell.standard_to_real_face_line(1, 1, orientations[1]), - ref_cell.standard_to_real_face_line(2, 1, orientations[1]), - ref_cell.standard_to_real_face_line(1, 2, orientations[2])}}; - line_orientations[0] = ref_cell.standard_vs_true_line_orientation( - 0, - 0, - orientations[0], - cell.quad(0)->line_orientation(my_indices[0])); - line_orientations[1] = ref_cell.standard_vs_true_line_orientation( - 1, - 0, - orientations[0], - cell.quad(0)->line_orientation(my_indices[1])); - line_orientations[2] = ref_cell.standard_vs_true_line_orientation( - 2, - 0, - orientations[0], - cell.quad(0)->line_orientation(my_indices[2])); - line_orientations[3] = ref_cell.standard_vs_true_line_orientation( - 1, - 1, - orientations[1], - cell.quad(1)->line_orientation(my_indices[3])); - line_orientations[4] = ref_cell.standard_vs_true_line_orientation( - 2, - 1, - orientations[1], - cell.quad(1)->line_orientation(my_indices[4])); - line_orientations[5] = ref_cell.standard_vs_true_line_orientation( - 1, - 2, - orientations[2], - cell.quad(2)->line_orientation(my_indices[5])); - } - else - // For other shapes (wedges, pyramids), we do not currently - // implement an optimized function - for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l) - line_orientations[l] = cell.line_orientation(l); - - return line_orientations; - } - }; - } // namespace TriaAccessorImplementation -} // namespace internal - - - -template -inline TriaAccessor::TriaAccessor( - const Triangulation *parent, - const int level, - const int index, - const AccessorData *local_data) - : TriaAccessorBase(parent, level, index, local_data) -{} - - - -template -inline bool -TriaAccessor::used() const -{ - Assert(this->state() == IteratorState::valid, - TriaAccessorExceptions::ExcDereferenceInvalidObject( - *this)); - return this->objects().used[this->present_index]; -} - - - -template -inline TriaIterator> -TriaAccessor::vertex_iterator( - const unsigned int i) const -{ - return TriaIterator>(this->tria, - 0, - vertex_index(i)); -} - - - -template -inline ReferenceCell -TriaAccessor::reference_cell() const -{ - if (structdim == 0) - return ReferenceCells::Vertex; - else if (structdim == 1) - return ReferenceCells::Line; - else if (structdim == dim) - return this->tria->levels[this->present_level] - ->reference_cell[this->present_index]; - else - return this->tria->faces->get_quad_type(this->present_index); -} - - - -template -inline unsigned int -TriaAccessor::vertex_index( - const unsigned int corner) const -{ - AssertIndexRange(corner, this->n_vertices()); - - if (structdim == dim) - { - constexpr unsigned int max_vertices_per_cell = 1 << dim; - const std::size_t my_index = - static_cast(this->present_index) * max_vertices_per_cell; - AssertIndexRange(my_index + corner, - this->tria->levels[this->present_level] - ->cell_vertex_indices_cache.size()); - const unsigned int vertex_index = - this->tria->levels[this->present_level] - ->cell_vertex_indices_cache[my_index + corner]; - Assert(vertex_index != numbers::invalid_unsigned_int, ExcInternalError()); - return vertex_index; - } - else - return dealii::internal::TriaAccessorImplementation::Implementation:: - vertex_index(*this, corner); -} - - - -template -inline Point & -TriaAccessor::vertex(const unsigned int i) const -{ - return const_cast &>(this->tria->vertices[vertex_index(i)]); -} - - - -template -inline typename dealii::internal::TriangulationImplementation:: - Iterators::line_iterator - TriaAccessor::line(const unsigned int i) const -{ - // checks happen in line_index - return typename dealii::internal::TriangulationImplementation:: - Iterators::line_iterator(this->tria, 0, line_index(i)); -} - - - -template -inline unsigned int -TriaAccessor::line_index(const unsigned int i) const -{ - AssertIndexRange(i, this->n_lines()); - - return dealii::internal::TriaAccessorImplementation::Implementation:: - line_index(*this, i); -} - - - -template -inline typename dealii::internal::TriangulationImplementation:: - Iterators::quad_iterator - TriaAccessor::quad(const unsigned int i) const -{ - // checks happen in quad_index - return typename dealii::internal::TriangulationImplementation:: - Iterators::quad_iterator(this->tria, 0, quad_index(i)); -} - - - -template -inline unsigned int -TriaAccessor::quad_index(const unsigned int i) const -{ - return dealii::internal::TriaAccessorImplementation::Implementation:: - quad_index(*this, i); -} - - - -template -inline unsigned char -TriaAccessor::combined_face_orientation( - const unsigned int face) const -{ - Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(face, n_faces()); - Assert(structdim == dim, - ExcMessage("This function can only be used on objects " - "that are cells, but not on faces or edges " - "that bound cells.")); - // work around a bogus GCC-9 warning which considers face unused except in 3d - (void)face; - - if constexpr (structdim == 1) - return ReferenceCell::default_combined_face_orientation(); - else if constexpr (structdim == 2) - return this->line_orientation(face) == true ? - ReferenceCell::default_combined_face_orientation() : - ReferenceCell::reversed_combined_line_orientation(); - else - return this->tria->levels[this->present_level] - ->face_orientations.get_combined_orientation( - this->present_index * GeometryInfo::faces_per_cell + face); -} - - - -template -inline bool -TriaAccessor::face_orientation( - const unsigned int face) const -{ - Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(face, n_faces()); - Assert(structdim == dim, - ExcMessage("This function can only be used on objects " - "that are cells, but not on faces or edges " - "that bound cells.")); - // work around a bogus GCC-9 warning which considers face unused in 1d - (void)face; - - if constexpr (structdim == 1) - // in 1d 'faces' are vertices and those are always consistently oriented - return true; - else if constexpr (structdim == 2) - return this->line_orientation(face); - else - return this->tria->levels[this->present_level] - ->face_orientations.get_orientation( - this->present_index * GeometryInfo::faces_per_cell + face); -} - - - -template -inline bool -TriaAccessor::face_flip(const unsigned int face) const -{ - Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - Assert(structdim == dim, - ExcMessage("This function can only be used on objects " - "that are cells, but not on faces or edges " - "that bound cells.")); - AssertIndexRange(face, n_faces()); - // work around a bogus GCC-9 warning which considers face unused except in 3d - (void)face; - - if constexpr (structdim == 3) - return this->tria->levels[this->present_level]->face_orientations.get_flip( - this->present_index * GeometryInfo<3>::faces_per_cell + face); - else - // In 1d and 2d, face_flip is always false as faces can only be - // 'flipped' in 3d. - return false; -} - - -template -inline bool -TriaAccessor::face_rotation( - const unsigned int face) const -{ - Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - Assert(structdim == dim, - ExcMessage("This function can only be used on objects " - "that are cells, but not on faces or edges " - "that bound cells.")); - AssertIndexRange(face, n_faces()); - // work around a bogus GCC-9 warning which considers face unused except in 3d - (void)face; - - if constexpr (structdim == 3) - return this->tria->levels[this->present_level] - ->face_orientations.get_rotation( - this->present_index * GeometryInfo<3>::faces_per_cell + face); - else - // In 1d and 2d, face_rotation is always false as faces can only be - // 'rotated' in 3d. - return false; -} - - - -template -inline bool -TriaAccessor::line_orientation( - const unsigned int line) const -{ - Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(line, this->n_lines()); - // work around a bogus GCC-9 warning which considers line unused in 1d - (void)line; - - if constexpr (structdim == 1) - return true; - else if constexpr (structdim == 2 && dim == 2) - // lines in 2d are faces - { - // if all elements are quads (or if we have a very special consistently - // oriented triangular mesh) then we do not store this array - if (this->tria->levels[this->present_level] - ->face_orientations.n_objects() == 0) - { - return true; - } - else - { - return this->tria->levels[this->present_level] - ->face_orientations.get_orientation( - this->present_index * GeometryInfo::faces_per_cell + - line); - } - } - else if constexpr (structdim == 2 && dim == 3) - { - // line orientations in 3d are stored in their own array - Assert(this->present_index * GeometryInfo<3>::lines_per_face + line < - this->tria->faces->quads_line_orientations.size(), - ExcInternalError()); - return this->tria->faces->quads_line_orientations - [this->present_index * GeometryInfo<3>::lines_per_face + line]; - } - else if constexpr (structdim == 3 && dim == 3) - { - const auto reference_cell = this->reference_cell(); - // First pick a face on which this line is a part of, and the - // index of the line within. - const auto [face_index, line_index] = - reference_cell.standard_line_to_face_and_line_index(line); - const auto line_within_face_index = - reference_cell.standard_to_real_face_line( - line_index, face_index, this->combined_face_orientation(face_index)); - - // Then query how that line is oriented within that face: - return reference_cell.standard_vs_true_line_orientation( - line_index, - face_index, - this->combined_face_orientation(face_index), - this->quad(face_index)->line_orientation(line_within_face_index)); - } - else - { - DEAL_II_ASSERT_UNREACHABLE(); - return false; - } -} - - - -template -inline void -TriaAccessor::set_line_orientation( - const unsigned int line, - const bool value) const -{ - Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(line, this->n_lines()); - - dealii::internal::TriaAccessorImplementation::Implementation:: - set_line_orientation(*this, line, value); -} - - - -template -inline void -TriaAccessor::set_combined_face_orientation( - const unsigned int face, - const unsigned char combined_orientation) const -{ - Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(face, this->n_faces()); - - dealii::internal::TriaAccessorImplementation::Implementation:: - set_combined_face_orientation(*this, face, combined_orientation); -} - - - -template -void -TriaAccessor::set_used_flag() const -{ - Assert(this->state() == IteratorState::valid, - TriaAccessorExceptions::ExcDereferenceInvalidObject( - *this)); - this->objects().used[this->present_index] = true; -} - - - -template -void -TriaAccessor::clear_used_flag() const -{ - Assert(this->state() == IteratorState::valid, - TriaAccessorExceptions::ExcDereferenceInvalidObject( - *this)); - this->objects().used[this->present_index] = false; -} - - -template -int -TriaAccessor::child_index(const unsigned int i) const -{ - Assert(has_children(), TriaAccessorExceptions::ExcCellHasNoChildren()); - AssertIndexRange(i, n_children()); - - // each set of two children are stored - // consecutively, so we only have to find - // the location of the set of children - const unsigned int n_sets_of_two = - GeometryInfo::max_children_per_cell / 2; - return this->objects().children[n_sets_of_two * this->present_index + i / 2] + - i % 2; -} - - - -template -int -TriaAccessor::isotropic_child_index( - const unsigned int i) const -{ - AssertIndexRange(i, GeometryInfo::max_children_per_cell); - - switch (structdim) - { - case 1: - return child_index(i); - case 2: - { - const RefinementCase<2> this_refinement_case( - static_cast(refinement_case())); - - Assert(this_refinement_case != RefinementCase<2>::no_refinement, - TriaAccessorExceptions::ExcCellHasNoChildren()); - - if (this_refinement_case == RefinementCase<2>::cut_xy) - return child_index(i); - else if ((this_refinement_case == RefinementCase<2>::cut_x) && - (child(i % 2)->refinement_case() == - RefinementCase<2>::cut_y)) - return child(i % 2)->child_index(i / 2); - else if ((this_refinement_case == RefinementCase<2>::cut_y) && - (child(i / 2)->refinement_case() == - RefinementCase<2>::cut_x)) - return child(i / 2)->child_index(i % 2); - else - Assert( - false, - ExcMessage( - "This cell has no grandchildren equivalent to isotropic refinement")); - break; - } - - case 3: - DEAL_II_NOT_IMPLEMENTED(); - } - return -1; -} - - - -template -RefinementCase -TriaAccessor::refinement_case() const -{ - Assert(this->state() == IteratorState::valid, - TriaAccessorExceptions::ExcDereferenceInvalidObject( - *this)); - - switch (structdim) - { - case 1: - return (RefinementCase( - this->objects().children[this->present_index] != -1 ? - // cast the branches - // here first to uchar - // and then (above) to - // RefinementCase - // so that the - // conversion is valid - // even for the case - // structdim>1 (for - // which this part of - // the code is dead - // anyway) - static_cast(RefinementCase<1>::cut_x) : - static_cast(RefinementCase<1>::no_refinement))); - - default: - Assert(static_cast(this->present_index) < - this->objects().refinement_cases.size(), - ExcIndexRange(this->present_index, - 0, - this->objects().refinement_cases.size())); - - return (static_cast>( - this->objects().refinement_cases[this->present_index])); - } -} - - - -template -inline TriaIterator> -TriaAccessor::child(const unsigned int i) const - -{ - // checking of 'i' happens in child_index - const TriaIterator> q( - this->tria, (dim == structdim ? this->level() + 1 : 0), child_index(i)); - - Assert((q.state() == IteratorState::past_the_end) || q->used(), - ExcInternalError()); - - return q; -} - - - -template -inline unsigned int -TriaAccessor::child_iterator_to_index( - const TriaIterator> &child) const -{ - const auto n_children = this->n_children(); - for (unsigned int child_n = 0; child_n < n_children; ++child_n) - if (this->child(child_n) == child) - return child_n; - - Assert(false, - ExcMessage("The given child is not a child of the current object.")); - return numbers::invalid_unsigned_int; -} - - - -template -inline TriaIterator> -TriaAccessor::isotropic_child( - const unsigned int i) const -{ - // checking of 'i' happens in child() or - // child_index() called below - switch (structdim) - { - case 1: - // no anisotropic refinement in 1d - return child(i); - - case 2: - { - const RefinementCase<2> this_refinement_case( - static_cast(refinement_case())); - - Assert(this_refinement_case != RefinementCase<2>::no_refinement, - TriaAccessorExceptions::ExcCellHasNoChildren()); - - if (this_refinement_case == RefinementCase<2>::cut_xy) - return child(i); - else if ((this_refinement_case == RefinementCase<2>::cut_x) && - (child(i % 2)->refinement_case() == - RefinementCase<2>::cut_y)) - return child(i % 2)->child(i / 2); - else if ((this_refinement_case == RefinementCase<2>::cut_y) && - (child(i / 2)->refinement_case() == - RefinementCase<2>::cut_x)) - return child(i / 2)->child(i % 2); - else - Assert( - false, - ExcMessage( - "This cell has no grandchildren equivalent to isotropic refinement")); - break; - } - - default: - DEAL_II_NOT_IMPLEMENTED(); - } - // we don't get here but have to return - // something... - return child(0); -} - - - -template -inline bool -TriaAccessor::has_children() const -{ - Assert(this->state() == IteratorState::valid, - TriaAccessorExceptions::ExcDereferenceInvalidObject( - *this)); - - // each set of two children are stored - // consecutively, so we only have to find - // the location of the set of children - const unsigned int n_sets_of_two = - GeometryInfo::max_children_per_cell / 2; - return (this->objects().children[n_sets_of_two * this->present_index] != -1); -} - - - -template -inline unsigned int -TriaAccessor::n_children() const -{ - Assert(this->state() == IteratorState::valid, - TriaAccessorExceptions::ExcDereferenceInvalidObject( - *this)); - if (reference_cell() == ReferenceCells::Tetrahedron) - return GeometryInfo::n_children( - RefinementCase::isotropic_refinement); - else - return GeometryInfo::n_children(refinement_case()); -} - - - -template -inline void -TriaAccessor::set_refinement_case( - const RefinementCase &refinement_case) const -{ - Assert(this->state() == IteratorState::valid, - TriaAccessorExceptions::ExcDereferenceInvalidObject( - *this)); - Assert(static_cast(this->present_index) < - this->objects().refinement_cases.size(), - ExcIndexRange(this->present_index, - 0, - this->objects().refinement_cases.size())); - - this->objects().refinement_cases[this->present_index] = refinement_case; -} - - -template -inline void -TriaAccessor::clear_refinement_case() const -{ - Assert(this->state() == IteratorState::valid, - TriaAccessorExceptions::ExcDereferenceInvalidObject( - *this)); - Assert(static_cast(this->present_index) < - this->objects().refinement_cases.size(), - ExcIndexRange(this->present_index, - 0, - this->objects().refinement_cases.size())); - - this->objects().refinement_cases[this->present_index] = - RefinementCase::no_refinement; -} - - - -template -void -TriaAccessor::set_children(const unsigned int i, - const int index) const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - Assert(i % 2 == 0, TriaAccessorExceptions::ExcSetOnlyEvenChildren(i)); - - // each set of two children are stored - // consecutively, so we only have to find - // the location of the set of children - const unsigned int n_sets_of_two = - GeometryInfo::max_children_per_cell / 2; - - Assert( - // clearing the child index for a cell - (index == -1) || - // if setting the child index for the i'th child (with i==0), - // then the index must be a non-negative number - (i == 0 && !this->has_children() && (index >= 0)) || - // if setting the child index for the i'th child (with i>0), - // then the previously stored index must be the invalid - // index - (i > 0 && this->has_children() && (index >= 0) && - this->objects().children[n_sets_of_two * this->present_index + i / 2] == - -1), - TriaAccessorExceptions::ExcCantSetChildren(index)); - - this->objects().children[n_sets_of_two * this->present_index + i / 2] = index; -} - - - -template -void -TriaAccessor::clear_children() const -{ - // each set of two children are stored - // consecutively, so we only have to find - // the location of the set of children - const unsigned int n_sets_of_two = - GeometryInfo::max_children_per_cell / 2; - - for (unsigned int i = 0; i < n_sets_of_two; ++i) - set_children(2 * i, -1); -} - - - -template -inline bool -TriaAccessor::user_flag_set() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - return this->objects().user_flags[this->present_index]; -} - - - -template -inline void -TriaAccessor::set_user_flag() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - this->objects().user_flags[this->present_index] = true; -} - - - -template -inline void -TriaAccessor::clear_user_flag() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - this->objects().user_flags[this->present_index] = false; -} - - - -template -void -TriaAccessor::recursively_set_user_flag() const -{ - set_user_flag(); - - if (this->has_children()) - for (unsigned int c = 0; c < this->n_children(); ++c) - this->child(c)->recursively_set_user_flag(); -} - - - -template -void -TriaAccessor::recursively_clear_user_flag() const -{ - clear_user_flag(); - - if (this->has_children()) - for (unsigned int c = 0; c < this->n_children(); ++c) - this->child(c)->recursively_clear_user_flag(); -} - - - -template -void -TriaAccessor::clear_user_data() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - this->objects().clear_user_data(this->present_index); -} - - - -template -void -TriaAccessor::set_user_pointer(void *p) const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - this->objects().user_pointer(this->present_index) = p; -} - - - -template -void -TriaAccessor::clear_user_pointer() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - this->objects().user_pointer(this->present_index) = nullptr; -} - - - -template -void * -TriaAccessor::user_pointer() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - return this->objects().user_pointer(this->present_index); -} - - - -template -void -TriaAccessor::recursively_set_user_pointer( - void *p) const -{ - set_user_pointer(p); - - if (this->has_children()) - for (unsigned int c = 0; c < this->n_children(); ++c) - this->child(c)->recursively_set_user_pointer(p); -} - - - -template -void -TriaAccessor::recursively_clear_user_pointer() const -{ - clear_user_pointer(); - - if (this->has_children()) - for (unsigned int c = 0; c < this->n_children(); ++c) - this->child(c)->recursively_clear_user_pointer(); -} - - - -template -void -TriaAccessor::set_user_index( - const unsigned int p) const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - this->objects().user_index(this->present_index) = p; -} - - - -template -void -TriaAccessor::clear_user_index() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - this->objects().user_index(this->present_index) = 0; -} - - - -template -unsigned int -TriaAccessor::user_index() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - return this->objects().user_index(this->present_index); -} - - - -template -void -TriaAccessor::recursively_set_user_index( - const unsigned int p) const -{ - set_user_index(p); - - if (this->has_children()) - for (unsigned int c = 0; c < this->n_children(); ++c) - this->child(c)->recursively_set_user_index(p); -} - - - -template -void -TriaAccessor::recursively_clear_user_index() const -{ - clear_user_index(); - - if (this->has_children()) - for (unsigned int c = 0; c < this->n_children(); ++c) - this->child(c)->recursively_clear_user_index(); -} - - - -template -inline unsigned int -TriaAccessor::max_refinement_depth() const -{ - if (!this->has_children()) - return 0; - - unsigned int max_depth = 1; - for (unsigned int c = 0; c < n_children(); ++c) - max_depth = std::max(max_depth, child(c)->max_refinement_depth() + 1); - return max_depth; -} - - - -template -unsigned int -TriaAccessor::n_active_descendants() const -{ - if (!this->has_children()) - return 1; - else - { - unsigned int sum = 0; - for (unsigned int c = 0; c < n_children(); ++c) - sum += this->child(c)->n_active_descendants(); - return sum; - } -} - - - -template -types::boundary_id -TriaAccessor::boundary_id() const -{ - Assert(structdim < dim, ExcImpossibleInDim(dim)); - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - - return this->objects() - .boundary_or_material_id[this->present_index] - .boundary_id; -} - - - -template -void -TriaAccessor::set_boundary_id( - const types::boundary_id boundary_ind) const -{ - Assert(structdim < dim, ExcImpossibleInDim(dim)); - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - Assert(boundary_ind != numbers::internal_face_boundary_id, - ExcMessage("You are trying to set the boundary_id to an invalid " - "value (numbers::internal_face_boundary_id is reserved).")); - Assert(this->at_boundary(), - ExcMessage("You are trying to set the boundary_id of an " - "internal object, which is not allowed!")); - - this->objects().boundary_or_material_id[this->present_index].boundary_id = - boundary_ind; -} - - - -template -void -TriaAccessor::set_boundary_id_internal( - const types::boundary_id boundary_ind) const -{ - Assert(structdim < dim, ExcImpossibleInDim(dim)); - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - - this->objects().boundary_or_material_id[this->present_index].boundary_id = - boundary_ind; -} - - - -template -void -TriaAccessor::set_all_boundary_ids( - const types::boundary_id boundary_ind) const -{ - set_boundary_id(boundary_ind); - - switch (structdim) - { - case 1: - // 1d objects have no sub-objects - // where we have to do anything - break; - - case 2: - // for boundary quads also set - // boundary_id of bounding lines - for (unsigned int i = 0; i < this->n_lines(); ++i) - this->line(i)->set_boundary_id(boundary_ind); - break; - - default: - DEAL_II_NOT_IMPLEMENTED(); - } -} - - - -template -bool -TriaAccessor::at_boundary() const -{ - // error checking is done - // in boundary_id() - return (boundary_id() != numbers::internal_face_boundary_id); -} - - - -template -const Manifold & -TriaAccessor::get_manifold() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - return this->tria->get_manifold(this->manifold_id()); -} - - -template -types::manifold_id -TriaAccessor::manifold_id() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - - return this->objects().manifold_id[this->present_index]; -} - - - -template -void -TriaAccessor::set_manifold_id( - const types::manifold_id manifold_ind) const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - - this->objects().manifold_id[this->present_index] = manifold_ind; -} - - -template -void -TriaAccessor::set_all_manifold_ids( - const types::manifold_id manifold_ind) const -{ - set_manifold_id(manifold_ind); - - if (this->has_children()) - for (unsigned int c = 0; c < this->n_children(); ++c) - this->child(c)->set_all_manifold_ids(manifold_ind); - - switch (structdim) - { - case 1: - if (dim == 1) - { - (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(0)] = - manifold_ind; - (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(1)] = - manifold_ind; - } - break; - - case 2: - // for quads/simplices also set manifold_id of bounding lines - for (unsigned int i = 0; i < this->n_lines(); ++i) - this->line(i)->set_manifold_id(manifold_ind); - break; - default: - DEAL_II_NOT_IMPLEMENTED(); - } -} - - - -template -double -TriaAccessor::diameter() const -{ - boost::container::small_vector, - GeometryInfo::vertices_per_cell> - vertices(this->n_vertices()); - - for (unsigned int v = 0; v < vertices.size(); ++v) - vertices[v] = this->vertex(v); - - return internal::TriaAccessorImplementation::diameter( - vertices); -} - - - -template -double -CellAccessor::diameter( - const Mapping &mapping) const -{ - return internal::TriaAccessorImplementation::diameter( - mapping.get_vertices(typename Triangulation::cell_iterator( - this->tria, this->level(), this->index()))); -} - - - -template -std::pair, double> -TriaAccessor::enclosing_ball() const -{ - // If the object is one dimensional, - // the enclosing ball is the initial iterate - // i.e., the ball's center and diameter are - // the center and the diameter of the object. - if (structdim == 1) - return std::make_pair((this->vertex(1) + this->vertex(0)) * 0.5, - (this->vertex(1) - this->vertex(0)).norm() * 0.5); - - // The list is_initial_guess_vertex contains bool values and has - // the same size as the number of vertices per object. - // The entries of is_initial_guess_vertex are set true only for those - // two vertices corresponding to the largest diagonal which is being used - // to construct the initial ball. - // We employ this mask to skip these two vertices while enlarging the ball. - std::vector is_initial_guess_vertex(this->n_vertices()); - - // First let all the vertices be outside - std::fill(is_initial_guess_vertex.begin(), - is_initial_guess_vertex.end(), - false); - - // Get an initial guess by looking at the largest diagonal - Point center; - double radius = 0; - - switch (structdim) - { - case 2: - { - const Point p30(this->vertex(3) - this->vertex(0)); - const Point p21(this->vertex(2) - this->vertex(1)); - if (p30.norm() > p21.norm()) - { - center = this->vertex(0) + 0.5 * p30; - radius = p30.norm() / 2.; - is_initial_guess_vertex[3] = true; - is_initial_guess_vertex[0] = true; - } - else - { - center = this->vertex(1) + 0.5 * p21; - radius = p21.norm() / 2.; - is_initial_guess_vertex[2] = true; - is_initial_guess_vertex[1] = true; - } - break; - } - case 3: - { - const Point p70(this->vertex(7) - this->vertex(0)); - const Point p61(this->vertex(6) - this->vertex(1)); - const Point p25(this->vertex(2) - this->vertex(5)); - const Point p34(this->vertex(3) - this->vertex(4)); - const std::vector diagonals = {p70.norm(), - p61.norm(), - p25.norm(), - p34.norm()}; - const std::vector::const_iterator it = - std::max_element(diagonals.begin(), diagonals.end()); - if (it == diagonals.begin()) - { - center = this->vertex(0) + 0.5 * p70; - is_initial_guess_vertex[7] = true; - is_initial_guess_vertex[0] = true; - } - else if (it == diagonals.begin() + 1) - { - center = this->vertex(1) + 0.5 * p61; - is_initial_guess_vertex[6] = true; - is_initial_guess_vertex[1] = true; - } - else if (it == diagonals.begin() + 2) - { - center = this->vertex(5) + 0.5 * p25; - is_initial_guess_vertex[2] = true; - is_initial_guess_vertex[5] = true; - } - else - { - center = this->vertex(4) + 0.5 * p34; - is_initial_guess_vertex[3] = true; - is_initial_guess_vertex[4] = true; - } - radius = *it * 0.5; - break; - } - default: - DEAL_II_NOT_IMPLEMENTED(); - return std::pair, double>(); - } - - // For each vertex that is found to be geometrically outside the ball - // enlarge the ball so that the new ball contains both the previous ball - // and the given vertex. - for (const unsigned int v : this->vertex_indices()) - if (!is_initial_guess_vertex[v]) - { - const double distance = center.distance(this->vertex(v)); - if (distance > radius) - { - // we found a vertex which is outside of the ball - // extend it (move center and change radius) - const Point pCV(center - this->vertex(v)); - radius = (distance + radius) * 0.5; - center = this->vertex(v) + pCV * (radius / distance); - - // Now the new ball constructed in this block - // encloses the vertex (v) that was found to be geometrically - // outside the old ball. - } - } -#ifdef DEBUG - bool all_vertices_within_ball = true; - - // Set all_vertices_within_ball false if any of the vertices of the object - // are geometrically outside the ball - for (const unsigned int v : this->vertex_indices()) - if (center.distance(this->vertex(v)) > - radius + 100. * std::numeric_limits::epsilon()) - { - all_vertices_within_ball = false; - break; - } - // If all the vertices are not within the ball throw error - Assert(all_vertices_within_ball, ExcInternalError()); -#endif - return std::make_pair(center, radius); -} - - -template -double -TriaAccessor::minimum_vertex_distance() const -{ - switch (structdim) - { - case 1: - return (this->vertex(1) - this->vertex(0)).norm(); - case 2: - case 3: - { - double min = std::numeric_limits::max(); - for (const unsigned int i : this->vertex_indices()) - for (unsigned int j = i + 1; j < this->n_vertices(); ++j) - min = std::min(min, - (this->vertex(i) - this->vertex(j)) * - (this->vertex(i) - this->vertex(j))); - return std::sqrt(min); - } - default: - DEAL_II_NOT_IMPLEMENTED(); - return -1e10; - } -} - - -template -bool -TriaAccessor::is_translation_of( - const TriaIterator> &o) const -{ - // go through the vertices and check... The - // cell is a translation of the previous - // one in case the distance between the - // individual vertices in the two cell is - // the same for all the vertices. So do the - // check by first getting the distance on - // the first vertex, and then checking - // whether all others have the same down to - // rounding errors (we have to be careful - // here because the calculation of the - // displacement between one cell and the - // next can already result in the loss of - // one or two digits), so we choose 1e-12 - // times the distance between the zeroth - // vertices here. - bool is_translation = true; - const Tensor<1, spacedim> dist = o->vertex(0) - this->vertex(0); - const double tol_square = 1e-24 * dist.norm_square(); - for (unsigned int i = 1; i < this->n_vertices(); ++i) - { - const Tensor<1, spacedim> dist_new = - (o->vertex(i) - this->vertex(i)) - dist; - if (dist_new.norm_square() > tol_square) - { - is_translation = false; - break; - } - } - return is_translation; -} - - - -template -unsigned int -TriaAccessor::n_vertices() const -{ - return this->reference_cell().n_vertices(); -} - - - -template -unsigned int -TriaAccessor::n_lines() const -{ - return this->reference_cell().n_lines(); -} - - - -template -unsigned int -TriaAccessor::n_faces() const -{ - Assert(structdim == dim, - ExcMessage("This function can only be used on objects " - "that are cells, but not on faces or edges " - "that bound cells.")); - - return this->reference_cell().n_faces(); -} - - - -template -std_cxx20::ranges::iota_view -TriaAccessor::vertex_indices() const -{ - return {0U, n_vertices()}; -} - - - -template -std_cxx20::ranges::iota_view -TriaAccessor::line_indices() const -{ - return {0U, n_lines()}; -} - - - -template -std_cxx20::ranges::iota_view -TriaAccessor::face_indices() const -{ - return {0U, n_faces()}; -} - - - -/*----------------- Functions: TriaAccessor<0,dim,spacedim> -----------------*/ - -template -inline TriaAccessor<0, dim, spacedim>::TriaAccessor( - const Triangulation *tria, - const unsigned int vertex_index) - : tria(tria) - , global_vertex_index(vertex_index) -{} - - - -template -inline TriaAccessor<0, dim, spacedim>::TriaAccessor( - const Triangulation *tria, - const int /*level*/, - const int index, - const AccessorData *) - : tria(tria) - , global_vertex_index(index) -{} - - - -template -template -inline TriaAccessor<0, dim, spacedim>::TriaAccessor( - const TriaAccessor &) - : tria(nullptr) - , global_vertex_index(numbers::invalid_unsigned_int) -{ - Assert(false, ExcImpossibleInDim(0)); -} - - - -template -template -inline TriaAccessor<0, dim, spacedim>::TriaAccessor( - const InvalidAccessor &) - : tria(nullptr) - , global_vertex_index(numbers::invalid_unsigned_int) -{ - Assert(false, ExcImpossibleInDim(0)); -} - - - -template -inline void -TriaAccessor<0, dim, spacedim>::copy_from(const TriaAccessor &t) -{ - tria = t.tria; - global_vertex_index = t.global_vertex_index; -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::operator<( - const TriaAccessor<0, dim, spacedim> &other) const -{ - Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators()); - - return (global_vertex_index < other.global_vertex_index); -} - - - -template -inline IteratorState::IteratorStates -TriaAccessor<0, dim, spacedim>::state() const -{ - if (global_vertex_index != numbers::invalid_unsigned_int) - return IteratorState::valid; - else - return IteratorState::past_the_end; -} - - - -template -inline int -TriaAccessor<0, dim, spacedim>::level() -{ - return 0; -} - - - -template -inline int -TriaAccessor<0, dim, spacedim>::index() const -{ - return global_vertex_index; -} - - - -template -inline const Triangulation & -TriaAccessor<0, dim, spacedim>::get_triangulation() const -{ - return *tria; -} - - - -template -inline void -TriaAccessor<0, dim, spacedim>::operator++() -{ - ++global_vertex_index; - if (global_vertex_index >= tria->n_vertices()) - global_vertex_index = numbers::invalid_unsigned_int; -} - - - -template -inline void -TriaAccessor<0, dim, spacedim>::operator--() -{ - if (global_vertex_index != numbers::invalid_unsigned_int) - { - if (global_vertex_index != 0) - --global_vertex_index; - else - global_vertex_index = numbers::invalid_unsigned_int; - } -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::operator==(const TriaAccessor &t) const -{ - const bool result = - ((tria == t.tria) && (global_vertex_index == t.global_vertex_index)); - - return result; -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::operator!=(const TriaAccessor &t) const -{ - return !(*this == t); -} - - - -template -inline unsigned int -TriaAccessor<0, dim, spacedim>::vertex_index(const unsigned int) const -{ - return global_vertex_index; -} - - - -template -inline Point & -TriaAccessor<0, dim, spacedim>::vertex(const unsigned int) const -{ - return const_cast &>( - this->tria->vertices[global_vertex_index]); -} - - - -template -inline typename dealii::internal::TriangulationImplementation:: - Iterators::line_iterator - TriaAccessor<0, dim, spacedim>::line(const unsigned int) -{ - return typename dealii::internal::TriangulationImplementation:: - Iterators::line_iterator(); -} - - - -template -inline unsigned int -TriaAccessor<0, dim, spacedim>::line_index(const unsigned int) -{ - Assert(false, ExcImpossibleInDim(0)); - return numbers::invalid_unsigned_int; -} - - - -template -inline typename dealii::internal::TriangulationImplementation:: - Iterators::quad_iterator - TriaAccessor<0, dim, spacedim>::quad(const unsigned int) -{ - return typename dealii::internal::TriangulationImplementation:: - Iterators::quad_iterator(); -} - - - -template -inline unsigned int -TriaAccessor<0, dim, spacedim>::quad_index(const unsigned int) -{ - Assert(false, ExcImpossibleInDim(0)); - return numbers::invalid_unsigned_int; -} - - - -template -inline double -TriaAccessor<0, dim, spacedim>::diameter() const -{ - return 0.; -} - - - -template -inline double -TriaAccessor<0, dim, spacedim>::extent_in_direction(const unsigned int) const -{ - return 0.; -} - - - -template -inline Point -TriaAccessor<0, dim, spacedim>::center(const bool, const bool) const -{ - return this->tria->vertices[global_vertex_index]; -} - - - -template -inline double -TriaAccessor<0, dim, spacedim>::measure() const -{ - return 0.; -} - - - -template -inline unsigned char -TriaAccessor<0, dim, spacedim>::combined_face_orientation( - const unsigned int /*face*/) -{ - return 0; -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::face_orientation(const unsigned int /*face*/) -{ - return false; -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::face_flip(const unsigned int /*face*/) -{ - return false; -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::face_rotation(const unsigned int /*face*/) -{ - return false; -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::line_orientation(const unsigned int /*line*/) -{ - return false; -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::has_children() -{ - return false; -} - - - -template -inline unsigned int -TriaAccessor<0, dim, spacedim>::n_children() -{ - return 0; -} - - - -template -inline unsigned int -TriaAccessor<0, dim, spacedim>::n_active_descendants() -{ - return 0; -} - - - -template -inline unsigned int -TriaAccessor<0, dim, spacedim>::max_refinement_depth() -{ - return 0; -} - - - -template -inline unsigned int -TriaAccessor<0, dim, spacedim>::child_iterator_to_index( - const TriaIterator> &) -{ - return numbers::invalid_unsigned_int; -} - - - -template -inline TriaIterator> -TriaAccessor<0, dim, spacedim>::child(const unsigned int) -{ - return TriaIterator>(); -} - - - -template -inline TriaIterator> -TriaAccessor<0, dim, spacedim>::isotropic_child(const unsigned int) -{ - return TriaIterator>(); -} - - - -template -inline RefinementCase<0> -TriaAccessor<0, dim, spacedim>::refinement_case() -{ - return {RefinementPossibilities<0>::no_refinement}; -} - - - -template -inline int -TriaAccessor<0, dim, spacedim>::child_index(const unsigned int) -{ - return -1; -} - - - -template -inline int -TriaAccessor<0, dim, spacedim>::isotropic_child_index(const unsigned int) -{ - return -1; -} - - - -template -inline bool -TriaAccessor<0, dim, spacedim>::used() const -{ - return tria->vertex_used(global_vertex_index); -} - - - -/*------------------- Functions: TriaAccessor<0,1,spacedim> -----------------*/ - -template -inline TriaAccessor<0, 1, spacedim>::TriaAccessor( - const Triangulation<1, spacedim> *tria, - const VertexKind vertex_kind, - const unsigned int vertex_index) - : tria(tria) - , vertex_kind(vertex_kind) - , global_vertex_index(vertex_index) -{} - - - -template -inline TriaAccessor<0, 1, spacedim>::TriaAccessor( - const Triangulation<1, spacedim> *tria, - const int level, - const int index, - const AccessorData *) - : tria(tria) - , vertex_kind(interior_vertex) - , global_vertex_index(numbers::invalid_unsigned_int) -{ - // in general, calling this constructor should yield an error -- users should - // instead call the one immediately above. however, if you create something - // like Triangulation<1>::face_iterator() then this calls the default - // constructor of the iterator which calls the accessor with argument list - // (0,-2,-2,0), so in this particular case accept this call and create an - // object that corresponds to the default constructed (invalid) vertex - // accessor - (void)level; - (void)index; - Assert((level == -2) && (index == -2), - ExcMessage( - "This constructor can not be called for face iterators in 1d, " - "except to default-construct iterator objects.")); -} - - - -template -template -inline TriaAccessor<0, 1, spacedim>::TriaAccessor( - const TriaAccessor &) - : tria(nullptr) - , vertex_kind(interior_vertex) - , global_vertex_index(numbers::invalid_unsigned_int) -{ - Assert(false, ExcImpossibleInDim(0)); -} - - - -template -template -inline TriaAccessor<0, 1, spacedim>::TriaAccessor( - const InvalidAccessor &) - : tria(nullptr) - , vertex_kind(interior_vertex) - , global_vertex_index(numbers::invalid_unsigned_int) -{ - Assert(false, ExcImpossibleInDim(0)); -} - - - -template -inline void -TriaAccessor<0, 1, spacedim>::copy_from(const TriaAccessor &t) -{ - tria = t.tria; - vertex_kind = t.vertex_kind; - global_vertex_index = t.global_vertex_index; -} - - - -template -inline void -TriaAccessor<0, 1, spacedim>::copy_from( - const TriaAccessorBase<0, 1, spacedim> &) -{ - // We cannot convert from TriaAccessorBase to - // TriaAccessor<0,1,spacedim> because the latter is not derived from - // the former. We should never get here. - DEAL_II_ASSERT_UNREACHABLE(); -} - - - -template -inline bool -TriaAccessor<0, 1, spacedim>::operator<( - const TriaAccessor<0, 1, spacedim> &other) const -{ - Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators()); - - return (global_vertex_index < other.global_vertex_index); -} - - - -template -inline IteratorState::IteratorStates -TriaAccessor<0, 1, spacedim>::state() -{ - return IteratorState::valid; -} - - -template -inline int -TriaAccessor<0, 1, spacedim>::level() -{ - return 0; -} - - - -template -inline int -TriaAccessor<0, 1, spacedim>::index() const -{ - return global_vertex_index; -} - - - -template -inline const Triangulation<1, spacedim> & -TriaAccessor<0, 1, spacedim>::get_triangulation() const -{ - return *tria; -} - - - -template -inline void -TriaAccessor<0, 1, spacedim>::operator++() const -{ - DEAL_II_NOT_IMPLEMENTED(); -} - - -template -inline void -TriaAccessor<0, 1, spacedim>::operator--() const -{ - DEAL_II_NOT_IMPLEMENTED(); -} - - - -template -inline bool -TriaAccessor<0, 1, spacedim>::operator==(const TriaAccessor &t) const -{ - const bool result = - ((tria == t.tria) && (global_vertex_index == t.global_vertex_index)); - // if we point to the same vertex, - // make sure we know the same about - // it - if (result == true) - Assert(vertex_kind == t.vertex_kind, ExcInternalError()); - - return result; -} - - - -template -inline bool -TriaAccessor<0, 1, spacedim>::operator!=(const TriaAccessor &t) const -{ - return !(*this == t); -} - - - -template -inline unsigned int -TriaAccessor<0, 1, spacedim>::vertex_index(const unsigned int i) const -{ - AssertIndexRange(i, 1); - (void)i; - return global_vertex_index; -} - - - -template -inline Point & -TriaAccessor<0, 1, spacedim>::vertex(const unsigned int i) const -{ - AssertIndexRange(i, 1); - (void)i; - return const_cast &>( - this->tria->vertices[global_vertex_index]); -} - - - -template -inline Point -TriaAccessor<0, 1, spacedim>::center() const -{ - return this->tria->vertices[global_vertex_index]; -} - - - -template -inline typename dealii::internal::TriangulationImplementation:: - Iterators<1, spacedim>::line_iterator - TriaAccessor<0, 1, spacedim>::line(const unsigned int) -{ - return {}; -} - - -template -inline unsigned int -TriaAccessor<0, 1, spacedim>::line_index(const unsigned int) -{ - Assert(false, ExcImpossibleInDim(0)); - return numbers::invalid_unsigned_int; -} - - -template -inline typename dealii::internal::TriangulationImplementation:: - Iterators<1, spacedim>::quad_iterator - TriaAccessor<0, 1, spacedim>::quad(const unsigned int) -{ - return {}; -} - - - -template -inline unsigned int -TriaAccessor<0, 1, spacedim>::quad_index(const unsigned int) -{ - Assert(false, ExcImpossibleInDim(0)); - return numbers::invalid_unsigned_int; -} - - -template -inline bool -TriaAccessor<0, 1, spacedim>::at_boundary() const -{ - return vertex_kind != interior_vertex; -} - - -template -inline types::boundary_id -TriaAccessor<0, 1, spacedim>::boundary_id() const -{ - switch (vertex_kind) - { - case left_vertex: - case right_vertex: - { - Assert(tria->vertex_to_boundary_id_map_1d->find( - this->vertex_index()) != - tria->vertex_to_boundary_id_map_1d->end(), - ExcInternalError()); - - return (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()]; - } - - default: - return numbers::internal_face_boundary_id; - } -} - - - -template -inline const Manifold<1, spacedim> & -TriaAccessor<0, 1, spacedim>::get_manifold() const -{ - return this->tria->get_manifold(this->manifold_id()); -} - - - -template -inline types::manifold_id -TriaAccessor<0, 1, spacedim>::manifold_id() const -{ - if (tria->vertex_to_manifold_id_map_1d->find(this->vertex_index()) != - tria->vertex_to_manifold_id_map_1d->end()) - return (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()]; - else - return numbers::flat_manifold_id; -} - - -template -inline unsigned char -TriaAccessor<0, 1, spacedim>::combined_face_orientation( - const unsigned int /*face*/) -{ - return 0; -} - - -template -inline bool -TriaAccessor<0, 1, spacedim>::face_orientation(const unsigned int /*face*/) -{ - return false; -} - - - -template -inline bool -TriaAccessor<0, 1, spacedim>::face_flip(const unsigned int /*face*/) -{ - return false; -} - - - -template -inline bool -TriaAccessor<0, 1, spacedim>::face_rotation(const unsigned int /*face*/) -{ - return false; -} - - - -template -inline bool -TriaAccessor<0, 1, spacedim>::line_orientation(const unsigned int /*line*/) -{ - return false; -} - - - -template -inline bool -TriaAccessor<0, 1, spacedim>::has_children() -{ - return false; -} - - - -template -inline unsigned int -TriaAccessor<0, 1, spacedim>::n_children() -{ - return 0; -} - - - -template -inline unsigned int -TriaAccessor<0, 1, spacedim>::n_active_descendants() -{ - return 0; -} - - - -template -inline unsigned int -TriaAccessor<0, 1, spacedim>::max_refinement_depth() -{ - return 0; -} - - - -template -inline unsigned int -TriaAccessor<0, 1, spacedim>::child_iterator_to_index( - const TriaIterator> &) -{ - return numbers::invalid_unsigned_int; -} - - - -template -inline TriaIterator> -TriaAccessor<0, 1, spacedim>::child(const unsigned int) -{ - return TriaIterator>(); -} - - -template -inline TriaIterator> -TriaAccessor<0, 1, spacedim>::isotropic_child(const unsigned int) -{ - return TriaIterator>(); -} - - -template -inline RefinementCase<0> -TriaAccessor<0, 1, spacedim>::refinement_case() -{ - return {RefinementPossibilities<0>::no_refinement}; -} - -template -inline int -TriaAccessor<0, 1, spacedim>::child_index(const unsigned int) -{ - return -1; -} - - -template -inline int -TriaAccessor<0, 1, spacedim>::isotropic_child_index(const unsigned int) -{ - return -1; -} - - - -template -inline void -TriaAccessor<0, 1, spacedim>::set_boundary_id(const types::boundary_id b) const -{ - Assert(tria->vertex_to_boundary_id_map_1d->find(this->vertex_index()) != - tria->vertex_to_boundary_id_map_1d->end(), - ExcMessage("You can't set the boundary_id of a face of a cell that is " - "not actually at the boundary.")); - - (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()] = b; -} - - - -template -inline void -TriaAccessor<0, 1, spacedim>::set_manifold_id(const types::manifold_id b) -{ - (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()] = b; -} - - - -template -inline void -TriaAccessor<0, 1, spacedim>::set_all_boundary_ids( - const types::boundary_id b) const -{ - set_boundary_id(b); -} - - - -template -inline void -TriaAccessor<0, 1, spacedim>::set_all_manifold_ids(const types::manifold_id b) -{ - set_manifold_id(b); -} - - - -template -inline bool -TriaAccessor<0, 1, spacedim>::used() const -{ - return tria->vertex_used(global_vertex_index); -} - - - -template -inline ReferenceCell -TriaAccessor<0, 1, spacedim>::reference_cell() const -{ - return ReferenceCells::Vertex; -} - - - -template -unsigned int -TriaAccessor<0, 1, spacedim>::n_vertices() const -{ - return 1; -} - - - -template -unsigned int -TriaAccessor<0, 1, spacedim>::n_lines() const -{ - return 0; -} - - - -template -std_cxx20::ranges::iota_view -TriaAccessor<0, 1, spacedim>::vertex_indices() const -{ - return {0U, n_vertices()}; -} - - - -template -std_cxx20::ranges::iota_view -TriaAccessor<0, 1, spacedim>::line_indices() const -{ - return {0U, n_lines()}; -} - -/*------------------ Functions: CellAccessor ------------------*/ - - -template -inline CellAccessor::CellAccessor( - const Triangulation *parent, - const int level, - const int index, - const AccessorData *local_data) - : TriaAccessor(parent, level, index, local_data) -{} - - - -template -inline CellAccessor::CellAccessor( - const TriaAccessor &cell_accessor) - : TriaAccessor( - static_cast &>(cell_accessor)) -{} - - - -template -inline TriaIterator> -CellAccessor::child(const unsigned int i) const -{ - TriaIterator> q(this->tria, - this->present_level + 1, - this->child_index(i)); - - Assert((q.state() == IteratorState::past_the_end) || q->used(), - ExcInternalError()); - - return q; -} - - - -template -inline boost::container::small_vector>, - GeometryInfo::max_children_per_cell> -CellAccessor::child_iterators() const -{ - boost::container::small_vector>, - GeometryInfo::max_children_per_cell> - child_iterators(this->n_children()); - - for (unsigned int i = 0; i < this->n_children(); ++i) - child_iterators[i] = this->child(i); - - return child_iterators; -} - - - -template -inline TriaIterator> -CellAccessor::face(const unsigned int i) const -{ - AssertIndexRange(i, this->n_faces()); - if constexpr (dim == 1) - { - TriaAccessor<0, 1, spacedim> a( - &this->get_triangulation(), - ((i == 0) && at_boundary(0) ? - TriaAccessor<0, 1, spacedim>::left_vertex : - ((i == 1) && at_boundary(1) ? - TriaAccessor<0, 1, spacedim>::right_vertex : - TriaAccessor<0, 1, spacedim>::interior_vertex)), - internal::TriaAccessorImplementation::Implementation::vertex_index( - *this, i)); - return TriaIterator>(a); - } - else if constexpr (dim == 2) - return this->line(i); - else if constexpr (dim == 3) - return this->quad(i); - else - { - Assert(false, ExcNotImplemented()); - return {}; - } -} - - - -template -inline unsigned int -CellAccessor::face_iterator_to_index( - const TriaIterator> &face) const -{ - for (const unsigned int face_n : this->face_indices()) - if (this->face(face_n) == face) - return face_n; - - Assert(false, - ExcMessage("The given face is not a face of the current cell.")); - return numbers::invalid_unsigned_int; -} - - - -template -inline boost::container::small_vector< - TriaIterator>, - GeometryInfo::faces_per_cell> -CellAccessor::face_iterators() const -{ - boost::container::small_vector< - TriaIterator>, - GeometryInfo::faces_per_cell> - face_iterators(this->n_faces()); - - for (const unsigned int i : this->face_indices()) - face_iterators[i] = this->face(i); - - return face_iterators; -} - - - -template -inline unsigned int -CellAccessor::face_index(const unsigned int i) const -{ - switch (dim) - { - case 1: - { - return this->vertex_index(i); - } - - case 2: - return this->line_index(i); - - case 3: - return this->quad_index(i); - - default: - return numbers::invalid_unsigned_int; - } -} - - - -template -inline int -CellAccessor::neighbor_index(const unsigned int face_no) const -{ - AssertIndexRange(face_no, this->n_faces()); - return this->tria->levels[this->present_level] - ->neighbors[this->present_index * GeometryInfo::faces_per_cell + - face_no] - .second; -} - - - -template -inline int -CellAccessor::neighbor_level(const unsigned int face_no) const -{ - AssertIndexRange(face_no, this->n_faces()); - return this->tria->levels[this->present_level] - ->neighbors[this->present_index * GeometryInfo::faces_per_cell + - face_no] - .first; -} - - - -template -inline RefinementCase -CellAccessor::refine_flag_set() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - // cells flagged for refinement must be active - // (the @p set_refine_flag function checks this, - // but activity may change when refinement is - // executed and for some reason the refine - // flag is not cleared). - Assert(this->is_active() || !this->tria->levels[this->present_level] - ->refine_flags[this->present_index], - ExcRefineCellNotActive()); - return RefinementCase( - this->tria->levels[this->present_level]->refine_flags[this->present_index]); -} - - - -template -inline void -CellAccessor::set_refine_flag( - const RefinementCase refinement_case) const -{ - Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); - Assert(!coarsen_flag_set(), ExcCellFlaggedForCoarsening()); - - this->tria->levels[this->present_level]->refine_flags[this->present_index] = - refinement_case; -} - - - -template -inline void -CellAccessor::clear_refine_flag() const -{ - Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); - this->tria->levels[this->present_level]->refine_flags[this->present_index] = - RefinementCase::no_refinement; -} - - -template -inline std::uint8_t -CellAccessor::refine_choice() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - if (this->tria->levels[this->present_level]->refine_choice.size() == 0) - return 0U; - return this->tria->levels[this->present_level] - ->refine_choice[this->present_index]; -} - - -template -inline void -CellAccessor::set_refine_choice( - const std::uint8_t refinement_choice) const -{ - Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); - if (this->tria->levels[this->present_level]->refine_choice.size() != 0) - this->tria->levels[this->present_level] - ->refine_choice[this->present_index] = refinement_choice; -} - - -template -inline void -CellAccessor::clear_refine_choice() const -{ - Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); - if (this->tria->levels[this->present_level]->refine_choice.size() != 0) - this->tria->levels[this->present_level] - ->refine_choice[this->present_index] = - static_cast(IsotropicRefinementChoice::isotropic_refinement); -} - - -template -inline bool -CellAccessor::flag_for_face_refinement( - const unsigned int face_no, - const RefinementCase &face_refinement_case) const -{ - Assert(dim > 1, ExcImpossibleInDim(dim)); - AssertIndexRange(face_no, this->n_faces()); - AssertIndexRange(face_refinement_case, - RefinementCase::isotropic_refinement + 1); - - // the new refinement case is a combination - // of the minimum required one for the given - // face refinement and the already existing - // flagged refinement case - RefinementCase old_ref_case = refine_flag_set(); - RefinementCase new_ref_case = - (old_ref_case | - GeometryInfo::min_cell_refinement_case_for_face_refinement( - face_refinement_case, - face_no, - this->face_orientation(face_no), - this->face_flip(face_no), - this->face_rotation(face_no))); - set_refine_flag(new_ref_case); - // return, whether we had to change the - // refinement flag - return new_ref_case != old_ref_case; -} - - - -template -inline bool -CellAccessor::flag_for_line_refinement( - const unsigned int line_no) const -{ - Assert(dim > 1, ExcImpossibleInDim(dim)); - AssertIndexRange(line_no, this->n_lines()); - - // the new refinement case is a combination - // of the minimum required one for the given - // line refinement and the already existing - // flagged refinement case - RefinementCase - old_ref_case = refine_flag_set(), - new_ref_case = - old_ref_case | - GeometryInfo::min_cell_refinement_case_for_line_refinement(line_no); - set_refine_flag(new_ref_case); - // return, whether we had to change the - // refinement flag - return new_ref_case != old_ref_case; -} - - - -template <> -inline dealii::internal::SubfaceCase<1> -CellAccessor<1>::subface_case(const unsigned int) const -{ - return dealii::internal::SubfaceCase<1>::case_none; -} - -template <> -inline dealii::internal::SubfaceCase<1> -CellAccessor<1, 2>::subface_case(const unsigned int) const -{ - return dealii::internal::SubfaceCase<1>::case_none; -} - - -template <> -inline dealii::internal::SubfaceCase<1> -CellAccessor<1, 3>::subface_case(const unsigned int) const -{ - return dealii::internal::SubfaceCase<1>::case_none; -} - - -template <> -inline dealii::internal::SubfaceCase<2> -CellAccessor<2>::subface_case(const unsigned int face_no) const -{ - Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); - AssertIndexRange(face_no, this->n_faces()); - return ((face(face_no)->has_children()) ? - dealii::internal::SubfaceCase<2>::case_x : - dealii::internal::SubfaceCase<2>::case_none); -} - -template <> -inline dealii::internal::SubfaceCase<2> -CellAccessor<2, 3>::subface_case(const unsigned int face_no) const -{ - Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); - AssertIndexRange(face_no, this->n_faces()); - return ((face(face_no)->has_children()) ? - dealii::internal::SubfaceCase<2>::case_x : - dealii::internal::SubfaceCase<2>::case_none); -} - - -template <> -inline dealii::internal::SubfaceCase<3> -CellAccessor<3>::subface_case(const unsigned int face_no) const -{ - Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); - AssertIndexRange(face_no, this->n_faces()); - switch (static_cast(face(face_no)->refinement_case())) - { - case RefinementCase<3>::no_refinement: - return dealii::internal::SubfaceCase<3>::case_none; - case RefinementCase<3>::cut_x: - if (face(face_no)->child(0)->has_children()) - { - Assert(face(face_no)->child(0)->refinement_case() == - RefinementCase<2>::cut_y, - ExcInternalError()); - if (face(face_no)->child(1)->has_children()) - { - Assert(face(face_no)->child(1)->refinement_case() == - RefinementCase<2>::cut_y, - ExcInternalError()); - return dealii::internal::SubfaceCase<3>::case_x1y2y; - } - else - return dealii::internal::SubfaceCase<3>::case_x1y; - } - else - { - if (face(face_no)->child(1)->has_children()) - { - Assert(face(face_no)->child(1)->refinement_case() == - RefinementCase<2>::cut_y, - ExcInternalError()); - return dealii::internal::SubfaceCase<3>::case_x2y; - } - else - return dealii::internal::SubfaceCase<3>::case_x; - } - case RefinementCase<3>::cut_y: - if (face(face_no)->child(0)->has_children()) - { - Assert(face(face_no)->child(0)->refinement_case() == - RefinementCase<2>::cut_x, - ExcInternalError()); - if (face(face_no)->child(1)->has_children()) - { - Assert(face(face_no)->child(1)->refinement_case() == - RefinementCase<2>::cut_x, - ExcInternalError()); - return dealii::internal::SubfaceCase<3>::case_y1x2x; - } - else - return dealii::internal::SubfaceCase<3>::case_y1x; - } - else - { - if (face(face_no)->child(1)->has_children()) - { - Assert(face(face_no)->child(1)->refinement_case() == - RefinementCase<2>::cut_x, - ExcInternalError()); - return dealii::internal::SubfaceCase<3>::case_y2x; - } - else - return dealii::internal::SubfaceCase<3>::case_y; - } - case RefinementCase<3>::cut_xy: - return dealii::internal::SubfaceCase<3>::case_xy; - default: - DEAL_II_ASSERT_UNREACHABLE(); - } - // we should never get here - return dealii::internal::SubfaceCase<3>::case_none; -} - - - -template -inline bool -CellAccessor::coarsen_flag_set() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - // cells flagged for coarsening must be active - // (the @p set_refine_flag function checks this, - // but activity may change when refinement is - // executed and for some reason the refine - // flag is not cleared). - Assert(this->is_active() || !this->tria->levels[this->present_level] - ->coarsen_flags[this->present_index], - ExcRefineCellNotActive()); - return this->tria->levels[this->present_level] - ->coarsen_flags[this->present_index]; -} - - - -template -inline void -CellAccessor::set_coarsen_flag() const -{ - Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); - Assert(!refine_flag_set(), ExcCellFlaggedForRefinement()); - - this->tria->levels[this->present_level]->coarsen_flags[this->present_index] = - true; -} - - - -template -inline void -CellAccessor::clear_coarsen_flag() const -{ - Assert(this->used() && this->is_active(), ExcRefineCellNotActive()); - this->tria->levels[this->present_level]->coarsen_flags[this->present_index] = - false; -} - - - -template -inline TriaIterator> -CellAccessor::neighbor(const unsigned int face_no) const -{ - TriaIterator> q(this->tria, - neighbor_level(face_no), - neighbor_index(face_no)); - - Assert((q.state() == IteratorState::past_the_end) || q->used(), - ExcInternalError()); - - return q; -} - - - -template -inline bool -CellAccessor::is_active() const -{ - return !this->has_children(); -} - - - -template -inline bool -CellAccessor::is_locally_owned() const -{ - Assert(this->is_active(), - ExcMessage("is_locally_owned() can only be called on active cells!")); -#ifndef DEAL_II_WITH_MPI - return true; -#else - - // Serial triangulations report invalid_subdomain_id as their locally owned - // subdomain, so the first condition checks whether we have a serial - // triangulation, in which case all cells are locally owned. The second - // condition compares the subdomain id in the parallel case. - const types::subdomain_id locally_owned_subdomain = - this->tria->locally_owned_subdomain(); - return (locally_owned_subdomain == numbers::invalid_subdomain_id || - this->subdomain_id() == locally_owned_subdomain); - -#endif -} - - -template -inline bool -CellAccessor::is_locally_owned_on_level() const -{ -#ifndef DEAL_II_WITH_MPI - return true; -#else - - // Serial triangulations report invalid_subdomain_id as their locally owned - // subdomain, so the first condition checks whether we have a serial - // triangulation, in which case all cells are locally owned. The second - // condition compares the subdomain id in the parallel case. - const types::subdomain_id locally_owned_subdomain = - this->tria->locally_owned_subdomain(); - return (locally_owned_subdomain == numbers::invalid_subdomain_id || - this->level_subdomain_id() == locally_owned_subdomain); - -#endif -} - - -template -inline bool -CellAccessor::is_ghost() const -{ - Assert(this->is_active(), - ExcMessage("is_ghost() can only be called on active cells!")); - if (this->has_children()) - return false; - -#ifndef DEAL_II_WITH_MPI - return false; -#else - - // Serial triangulations report invalid_subdomain_id as their locally owned - // subdomain, so the first condition rules out that case as all cells to a - // serial triangulation are locally owned and none is ghosted. The second - // and third conditions check whether the cell's subdomain is not the - // locally owned one and not artificial. - const types::subdomain_id locally_owned_subdomain = - this->tria->locally_owned_subdomain(); - const types::subdomain_id subdomain_id = this->subdomain_id(); - return (locally_owned_subdomain != numbers::invalid_subdomain_id && - subdomain_id != locally_owned_subdomain && - subdomain_id != numbers::artificial_subdomain_id); - -#endif -} - - -template -inline bool -CellAccessor::is_ghost_on_level() const -{ -#ifndef DEAL_II_WITH_MPI - return false; -#else - - // Serial triangulations report invalid_subdomain_id as their locally owned - // subdomain, so the first condition checks whether we have a serial - // triangulation, in which case all cells are locally owned. The second - // condition compares the subdomain id in the parallel case. - const types::subdomain_id locally_owned_subdomain = - this->tria->locally_owned_subdomain(); - const types::subdomain_id subdomain_id = this->level_subdomain_id(); - return (locally_owned_subdomain != numbers::invalid_subdomain_id && - subdomain_id != locally_owned_subdomain && - subdomain_id != numbers::artificial_subdomain_id); - -#endif -} - - - -template -inline bool -CellAccessor::is_artificial() const -{ - Assert(this->is_active(), - ExcMessage("is_artificial() can only be called on active cells!")); -#ifndef DEAL_II_WITH_MPI - return false; -#else - - // Serial triangulations report invalid_subdomain_id as their locally owned - // subdomain, so the first condition rules out that case as all cells to a - // serial triangulation are locally owned and none is artificial. - return (this->tria->locally_owned_subdomain() != - numbers::invalid_subdomain_id && - this->subdomain_id() == numbers::artificial_subdomain_id); - -#endif -} - - - -template -inline bool -CellAccessor::is_artificial_on_level() const -{ -#ifndef DEAL_II_WITH_MPI - return false; -#else - return (this->tria->locally_owned_subdomain() != - numbers::invalid_subdomain_id && - this->level_subdomain_id() == numbers::artificial_subdomain_id); -#endif -} - - - -template -inline types::subdomain_id -CellAccessor::subdomain_id() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - Assert(this->is_active(), - ExcMessage("subdomain_id() can only be called on active cells!")); - return this->tria->levels[this->present_level] - ->subdomain_ids[this->present_index]; -} - - - -template -inline types::subdomain_id -CellAccessor::level_subdomain_id() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - return this->tria->levels[this->present_level] - ->level_subdomain_ids[this->present_index]; -} - - - -template -inline unsigned int -CellAccessor::neighbor_face_no(const unsigned int neighbor) const -{ - const unsigned int n2 = neighbor_of_neighbor_internal(neighbor); - if (n2 != numbers::invalid_unsigned_int) - // return this value as the - // neighbor is not coarser - return n2; - else - // the neighbor is coarser - return neighbor_of_coarser_neighbor(neighbor).first; -} - - - -template -inline bool -CellAccessor::is_level_cell() -{ - return false; -} - - - -template -inline unsigned int -CellAccessor::active_cell_index() const -{ - Assert(this->is_active(), TriaAccessorExceptions::ExcCellNotActive()); - return this->tria->levels[this->present_level] - ->active_cell_indices[this->present_index]; -} - - - -template -inline types::global_cell_index -CellAccessor::global_active_cell_index() const -{ - Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - Assert(this->is_active(), - ExcMessage( - "global_active_cell_index() can only be called on active cells!")); - - return this->tria->levels[this->present_level] - ->global_active_cell_indices[this->present_index]; -} - - - -template -inline types::global_cell_index -CellAccessor::global_level_cell_index() const -{ - return this->tria->levels[this->present_level] - ->global_level_cell_indices[this->present_index]; -} - - -DEAL_II_NAMESPACE_CLOSE +#warning \ + "The use of this header file is deprecated. Just include ." #endif diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 157ba26555..7a2bf6f9c4 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -24,7 +24,6 @@ #include #include #include -#include #include #include #include -- 2.39.5