From: Peter Munch Date: Thu, 3 Dec 2020 06:36:38 +0000 (+0100) Subject: Unify diameter implementations X-Git-Tag: v9.3.0-rc1~799^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11301%2Fhead;p=dealii.git Unify diameter implementations --- diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index b2166d1424..ca2b06c927 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -53,6 +53,9 @@ namespace parallel template class Manifold; + +template +class Mapping; #endif namespace internal @@ -124,7 +127,6 @@ namespace internal { using type = int; }; - } // namespace TriaAccessorImplementation } // namespace internal template @@ -3717,6 +3719,14 @@ public: CellId id() const; + using TriaAccessor::diameter; + + /** + * The same as TriaAccessor::diameter() but also taking a Mapping class. + */ + double + diameter(const Mapping &mapping) const; + /** * @} */ diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index 07624ba4a4..773a2180cd 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -38,6 +38,79 @@ DEAL_II_ENABLE_EXTRA_DIAGNOSTICS DEAL_II_NAMESPACE_OPEN +namespace internal +{ + /** + * 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) + { + switch (ReferenceCell::n_vertices_to_type(dim, vertices.size())) + { + case ReferenceCell::Type::Line: + // Return the distance between the two vertices + return (vertices[1] - vertices[0]).norm(); + + case ReferenceCell::Type::Tri: + // 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()}); + + case ReferenceCell::Type::Quad: + // Return the longer one of the two diagonals of the quadrilateral + return std::max({(vertices[3] - vertices[0]).norm(), + (vertices[2] - vertices[1]).norm()}); + + case ReferenceCell::Type::Tet: + // 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()}); + + case ReferenceCell::Type::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()}); + + case ReferenceCell::Type::Hex: + // 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()}); + default: + Assert(false, ExcNotImplemented()); + return -1e10; + } + } +} // namespace internal + + /*--------------------- Functions: TriaAccessorBase -------------------------*/ template @@ -1885,62 +1958,26 @@ template double TriaAccessor::diameter() const { - switch (this->reference_cell_type()) - { - case ReferenceCell::Type::Line: - // Return the distance between the two vertices - return (this->vertex(1) - this->vertex(0)).norm(); + boost::container::small_vector, + GeometryInfo::vertices_per_cell> + vertices(this->n_vertices()); - case ReferenceCell::Type::Tri: - // Return the longest of the three edges - return std::max({(this->vertex(1) - this->vertex(0)).norm(), - (this->vertex(2) - this->vertex(1)).norm(), - (this->vertex(2) - this->vertex(0)).norm()}); - - case ReferenceCell::Type::Quad: - // Return the longer one of the two diagonals of the quadrilateral - return std::max({(this->vertex(3) - this->vertex(0)).norm(), - (this->vertex(2) - this->vertex(1)).norm()}); - - case ReferenceCell::Type::Tet: - // Return the longest of the six edges of the tetrahedron - return std::max({(this->vertex(1) - this->vertex(0)).norm(), - (this->vertex(2) - this->vertex(0)).norm(), - (this->vertex(2) - this->vertex(1)).norm(), - (this->vertex(3) - this->vertex(0)).norm(), - (this->vertex(3) - this->vertex(1)).norm(), - (this->vertex(3) - this->vertex(2)).norm()}); - - case ReferenceCell::Type::Wedge: - // Return ... - return std::max({// the longest of the 2*3=6 diagonals of the three - // quadrilateral sides of the wedge or ... - (this->vertex(4) - this->vertex(0)).norm(), - (this->vertex(3) - this->vertex(1)).norm(), - (this->vertex(5) - this->vertex(1)).norm(), - (this->vertex(4) - this->vertex(2)).norm(), - (this->vertex(5) - this->vertex(0)).norm(), - (this->vertex(3) - this->vertex(2)).norm(), - // the longest of the 3*2=6 edges of the two triangular - // faces of the wedge - (this->vertex(1) - this->vertex(0)).norm(), - (this->vertex(2) - this->vertex(1)).norm(), - (this->vertex(2) - this->vertex(0)).norm(), - (this->vertex(4) - this->vertex(3)).norm(), - (this->vertex(5) - this->vertex(4)).norm(), - (this->vertex(5) - this->vertex(3)).norm()}); - - case ReferenceCell::Type::Hex: - // Return the longest of the four diagonals of the hexahedron - return std::max({(this->vertex(7) - this->vertex(0)).norm(), - (this->vertex(6) - this->vertex(1)).norm(), - (this->vertex(2) - this->vertex(5)).norm(), - (this->vertex(3) - this->vertex(4)).norm()}); + for (unsigned int v = 0; v < vertices.size(); ++v) + vertices[v] = this->vertex(v); - default: - Assert(false, ExcNotImplemented()); - return -1e10; - } + return internal::diameter(vertices); +} + + + +template +double +CellAccessor::diameter( + const Mapping &mapping) const +{ + return internal::diameter( + mapping.get_vertices(typename Triangulation::cell_iterator( + this->tria, this->level(), this->index()))); } diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 9ad875bb2b..2900ae3e34 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -3099,49 +3099,7 @@ namespace GridTools diameter(const typename Triangulation::cell_iterator &cell, const Mapping &mapping) { - // see also TriaAccessor::diameter() - - const auto vertices = mapping.get_vertices(cell); - switch (cell->reference_cell_type()) - { - case ReferenceCell::Type::Line: - return (vertices[1] - vertices[0]).norm(); - case ReferenceCell::Type::Tri: - return std::max({(vertices[1] - vertices[0]).norm(), - (vertices[2] - vertices[1]).norm(), - (vertices[2] - vertices[0]).norm()}); - case ReferenceCell::Type::Quad: - return std::max({(vertices[3] - vertices[0]).norm(), - (vertices[2] - vertices[1]).norm()}); - case ReferenceCell::Type::Tet: - 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()}); - case ReferenceCell::Type::Wedge: - return std::max({(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(), - (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()}); - case ReferenceCell::Type::Hex: - return std::max({(vertices[7] - vertices[0]).norm(), - (vertices[6] - vertices[1]).norm(), - (vertices[2] - vertices[5]).norm(), - (vertices[3] - vertices[4]).norm()}); - default: - Assert(false, ExcNotImplemented()); - return -1e10; - } + return cell->diameter(mapping); } } // namespace internal