Unify diameter implementations

diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h
index b2166d14241fd511fa3e709076dd34ee41ff9307..ca2b06c927b7fc4fd1f06d470eaef007b875cc5a 100644 (file)
--- a/include/deal.II/grid/tria_accessor.h
+++ b/include/deal.II/grid/tria_accessor.h
@@ -53,6 +53,9 @@ namespace parallel
 
 template <int dim, int spacedim>
 class Manifold;
+
+template <int dim, int spacedim>
+class Mapping;
 #endif
 
 namespace internal
@@ -124,7 +127,6 @@ namespace internal
     {
       using type = int;
     };
-
   } // namespace TriaAccessorImplementation
 } // namespace internal
 template <int structdim, int dim, int spacedim>
@@ -3717,6 +3719,14 @@ public:
   CellId
   id() const;
 
+  using TriaAccessor<dim, dim, spacedim>::diameter;
+
+  /**
+   * The same as TriaAccessor::diameter() but also taking a Mapping class.
+   */
+  double
+  diameter(const Mapping<dim, spacedim> &mapping) const;
+
   /**
    * @}
    */

diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h
index 07624ba4a4749bbf2e4a3b8647816b7f1b7c0bfc..773a2180cd0bcbf3ca611ff8d8383bce19ae886e 100644 (file)
--- 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 <int dim, int spacedim>
+  inline double
+  diameter(
+    const boost::container::small_vector<Point<spacedim>,
+                                         GeometryInfo<dim>::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 <int structdim, int dim, int spacedim>
@@ -1885,62 +1958,26 @@ template <int structdim, int dim, int spacedim>
 double
 TriaAccessor<structdim, dim, spacedim>::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<Point<spacedim>,
+                                 GeometryInfo<structdim>::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<structdim, spacedim>(vertices);
+}
+
+
+
+template <int dim, int spacedim>
+double
+CellAccessor<dim, spacedim>::diameter(
+  const Mapping<dim, spacedim> &mapping) const
+{
+  return internal::diameter<dim, spacedim>(
+    mapping.get_vertices(typename Triangulation<dim, spacedim>::cell_iterator(
+      this->tria, this->level(), this->index())));
 }
 
 

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 9ad875bb2b264f0b4e13c999fe24ab6b0822c822..2900ae3e34079ef02ab7ec66724df2d6decde508 100644 (file)
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -3099,49 +3099,7 @@ namespace GridTools
     diameter(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
              const Mapping<dim, spacedim> &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