Extend GridTools::maximal_cell_diameter and ::minimal_cell_diameter for simplex

diff --git a/doc/news/changes/incompatibilities/20201016Munch b/doc/news/changes/incompatibilities/20201016Munch
new file mode 100644 (file)
index 0000000..eeb889e
--- /dev/null
+++ b/doc/news/changes/incompatibilities/20201016Munch
@@ -0,0 +1,4 @@
+Changed: The return type of the method Mapping::get_vertices() has been changed
+from std::array to boost::container::small_vector.  
+<br>
+(Peter Munch, 2020/10/16)

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 8b302d8ba74d26b64a78e5fb78c4bd39931d576c..aec51196304a083ad8c27e839f29e6b006bc8d77 100644 (file)
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -332,7 +332,8 @@ public:
    * information stored by the triangulation, i.e.,
    * <code>cell-@>vertex(v)</code>.
    */
-  virtual std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  virtual boost::container::small_vector<Point<spacedim>,
+                                         GeometryInfo<dim>::vertices_per_cell>
   get_vertices(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
 

diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h
index 691549b4d4a9217b8b4b18b2303fd0c5095a0957..5c2eae1d94ebb2ee4d90568441b5be3bbff2b694 100644 (file)
--- a/include/deal.II/fe/mapping_fe_field.h
+++ b/include/deal.II/fe/mapping_fe_field.h
@@ -169,7 +169,8 @@ public:
    * and constructs the position of the vertices according to the @p euler_vector
    * that was passed at construction time.
    */
-  virtual std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  virtual boost::container::small_vector<Point<spacedim>,
+                                         GeometryInfo<dim>::vertices_per_cell>
   get_vertices(const typename Triangulation<dim, spacedim>::cell_iterator &cell)
     const override;
 

diff --git a/include/deal.II/fe/mapping_q1_eulerian.h b/include/deal.II/fe/mapping_q1_eulerian.h
index 531f904a8287aba294aa6cbcc9a750a52e26425d..6a9d77df97503fdf3b66808e4448152f2c271f04 100644 (file)
--- a/include/deal.II/fe/mapping_q1_eulerian.h
+++ b/include/deal.II/fe/mapping_q1_eulerian.h
@@ -116,7 +116,8 @@ public:
    * cell but instead evaluates an externally given displacement field in
    * addition to the geometry of the cell.
    */
-  virtual std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  virtual boost::container::small_vector<Point<spacedim>,
+                                         GeometryInfo<dim>::vertices_per_cell>
   get_vertices(const typename Triangulation<dim, spacedim>::cell_iterator &cell)
     const override;
 

diff --git a/include/deal.II/fe/mapping_q_eulerian.h b/include/deal.II/fe/mapping_q_eulerian.h
index 7e5b5adab8796f9f839a09442e69bd54ac6b492a..10bd13b13f888084c9586be54c86bc880bce5423 100644 (file)
--- a/include/deal.II/fe/mapping_q_eulerian.h
+++ b/include/deal.II/fe/mapping_q_eulerian.h
@@ -120,7 +120,8 @@ public:
    * cell but instead evaluates an externally given displacement field in
    * addition to the geometry of the cell.
    */
-  virtual std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  virtual boost::container::small_vector<Point<spacedim>,
+                                         GeometryInfo<dim>::vertices_per_cell>
   get_vertices(const typename Triangulation<dim, spacedim>::cell_iterator &cell)
     const override;
 
@@ -204,7 +205,8 @@ private:
      * current cell but instead evaluates an externally given displacement
      * field in addition to the geometry of the cell.
      */
-    virtual std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+    virtual boost::container::small_vector<Point<spacedim>,
+                                           GeometryInfo<dim>::vertices_per_cell>
     get_vertices(const typename Triangulation<dim, spacedim>::cell_iterator
                    &cell) const override;
 

diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc
index 5433bf7471a06a27d7ff513c44f251dcc9d70574..3e10ef99e2cb59848bb153e9f3429da2b69f1b49 100644 (file)
--- a/source/fe/mapping.cc
+++ b/source/fe/mapping.cc
@@ -24,15 +24,17 @@ DEAL_II_NAMESPACE_OPEN
 
 
 template <int dim, int spacedim>
-std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+boost::container::small_vector<Point<spacedim>,
+                               GeometryInfo<dim>::vertices_per_cell>
 Mapping<dim, spacedim>::get_vertices(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
 {
-  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
-  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
-    {
-      vertices[i] = cell->vertex(i);
-    }
+  boost::container::small_vector<Point<spacedim>,
+                                 GeometryInfo<dim>::vertices_per_cell>
+    vertices;
+  for (const unsigned int i : cell->vertex_indices())
+    vertices.push_back(cell->vertex(i));
+
   return vertices;
 }
 

diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc
index 163719f12e6d3f35a984db0319319506f3aa4a76..50667b2f517491751988efaedb524c07c1b558f1 100644 (file)
--- a/source/fe/mapping_fe_field.cc
+++ b/source/fe/mapping_fe_field.cc
@@ -357,7 +357,8 @@ MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
 
 
 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+boost::container::small_vector<Point<spacedim>,
+                               GeometryInfo<dim>::vertices_per_cell>
 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::get_vertices(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
 {
@@ -395,7 +396,9 @@ MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::get_vertices(
   const VectorType &vector =
     uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
 
-  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
+  boost::container::small_vector<Point<spacedim>,
+                                 GeometryInfo<dim>::vertices_per_cell>
+    vertices(GeometryInfo<dim>::vertices_per_cell);
   for (unsigned int i = 0; i < dofs_per_cell; ++i)
     {
       const unsigned int comp = fe_to_real

diff --git a/source/fe/mapping_q1_eulerian.cc b/source/fe/mapping_q1_eulerian.cc
index 36ff806f398071aa24a383f182f47a880a275277..3d495c32b4776669e09555ffb8b6fc094cbb3c06 100644 (file)
--- a/source/fe/mapping_q1_eulerian.cc
+++ b/source/fe/mapping_q1_eulerian.cc
@@ -50,11 +50,14 @@ MappingQ1Eulerian<dim, VectorType, spacedim>::MappingQ1Eulerian(
 
 
 template <int dim, class VectorType, int spacedim>
-std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+boost::container::small_vector<Point<spacedim>,
+                               GeometryInfo<dim>::vertices_per_cell>
 MappingQ1Eulerian<dim, VectorType, spacedim>::get_vertices(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
 {
-  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
+  boost::container::small_vector<Point<spacedim>,
+                                 GeometryInfo<dim>::vertices_per_cell>
+    vertices(GeometryInfo<dim>::vertices_per_cell);
   // The assertions can not be in the constructor, since this would
   // require to call dof_handler.distribute_dofs(fe) *before* the mapping
   // object is constructed, which is not necessarily what we want.
@@ -105,8 +108,7 @@ std::vector<Point<spacedim>>
 MappingQ1Eulerian<dim, VectorType, spacedim>::compute_mapping_support_points(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
 {
-  const std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
-    vertices = this->get_vertices(cell);
+  const auto vertices = this->get_vertices(cell);
 
   std::vector<Point<spacedim>> a(GeometryInfo<dim>::vertices_per_cell);
   for (const unsigned int i : GeometryInfo<dim>::vertex_indices())

diff --git a/source/fe/mapping_q_eulerian.cc b/source/fe/mapping_q_eulerian.cc
index e10223f422f5b3b1405dfac13e066d3a750e7222..f44f0d796ea60b085d3eb0cbea57af679ece2938 100644 (file)
--- a/source/fe/mapping_q_eulerian.cc
+++ b/source/fe/mapping_q_eulerian.cc
@@ -121,7 +121,8 @@ MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerianGeneric::
 // .... COMPUTE MAPPING SUPPORT POINTS
 
 template <int dim, class VectorType, int spacedim>
-std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+boost::container::small_vector<Point<spacedim>,
+                               GeometryInfo<dim>::vertices_per_cell>
 MappingQEulerian<dim, VectorType, spacedim>::get_vertices(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
 {
@@ -130,11 +131,10 @@ MappingQEulerian<dim, VectorType, spacedim>::get_vertices(
     dynamic_cast<const MappingQEulerianGeneric &>(*this->qp_mapping)
       .compute_mapping_support_points(cell);
 
-  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
-    vertex_locations;
-  std::copy(a.begin(),
-            a.begin() + GeometryInfo<dim>::vertices_per_cell,
-            vertex_locations.begin());
+  boost::container::small_vector<Point<spacedim>,
+                                 GeometryInfo<dim>::vertices_per_cell>
+    vertex_locations(a.begin(),
+                     a.begin() + GeometryInfo<dim>::vertices_per_cell);
 
   return vertex_locations;
 }
@@ -142,7 +142,8 @@ MappingQEulerian<dim, VectorType, spacedim>::get_vertices(
 
 
 template <int dim, class VectorType, int spacedim>
-std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+boost::container::small_vector<Point<spacedim>,
+                               GeometryInfo<dim>::vertices_per_cell>
 MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerianGeneric::
   get_vertices(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell) const

diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index 3d878bbc04a435dd57c59a24c4e36b794b80e121..3497a2ed00e45ecf114f980aae5d5eb94ab38c2c 100644 (file)
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -2518,8 +2518,13 @@ MappingQGeneric<dim, spacedim>::transform_real_to_unit_cell(
       //   cell and that value is returned.
       // * In 3D there is no (known to the authors) exact formula, so the Newton
       //   algorithm is used.
-      const std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
-        vertices = this->get_vertices(cell);
+      const auto vertices_ = this->get_vertices(cell);
+
+      std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+        vertices;
+      for (unsigned int i = 0; i < vertices.size(); ++i)
+        vertices[i] = vertices_[i];
+
       try
         {
           switch (dim)

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 62212f88816fe1f5ab1c0bf3486bf8c434aff0fa..4b1aec7c5d702c35338783899c2c3c6b58133b73 100644 (file)
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -2988,15 +2988,29 @@ 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 (dim)
+      switch (cell->reference_cell_type())
         {
-          case 1:
+          case ReferenceCell::Type::Line:
             return (vertices[1] - vertices[0]).norm();
-          case 2:
+          case ReferenceCell::Type::Tri:
+            return std::max(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 3:
+          case ReferenceCell::Type::Tet:
+            return std::max(
+              std::max(std::max((vertices[1] - vertices[0]).norm(),
+                                (vertices[2] - vertices[0]).norm()),
+                       std::max((vertices[2] - vertices[1]).norm(),
+                                (vertices[3] - vertices[0]).norm())),
+              std::max((vertices[3] - vertices[1]).norm(),
+                       (vertices[3] - vertices[2]).norm()));
+          case ReferenceCell::Type::Hex:
             return std::max(std::max((vertices[7] - vertices[0]).norm(),
                                      (vertices[6] - vertices[1]).norm()),
                             std::max((vertices[2] - vertices[5]).norm(),

diff --git a/source/opencascade/utilities.cc b/source/opencascade/utilities.cc
index b9f11d79b94addcdc28511a4164792ab446d822a..1011d14ac68868b31fe8d7321e46cdaf63c6edcc 100644 (file)
--- a/source/opencascade/utilities.cc
+++ b/source/opencascade/utilities.cc
@@ -619,8 +619,7 @@ namespace OpenCASCADE
             visited_faces[face_index]        = false;
 
             // extract mapped vertex locations
-            std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell>
-              verts           = mapping.get_vertices(cell);
+            const auto verts  = mapping.get_vertices(cell);
             vert_to_point[v0] = verts[GeometryInfo<2>::face_to_cell_vertices(
               f, 0, true, false, false)];
             vert_to_point[v1] = verts[GeometryInfo<2>::face_to_cell_vertices(

diff --git a/tests/simplex/cell_measure_01.cc b/tests/simplex/cell_measure_01.cc
index 08a57fd5126a1025589079a7e35ec83ff3b7807f..a5112d7330075df41102b7759694d1583639904a 100644 (file)
--- a/tests/simplex/cell_measure_01.cc
+++ b/tests/simplex/cell_measure_01.cc
@@ -17,8 +17,13 @@
 
 // Test TriaAccessor::measure() and TriaAccessor::diameter().
 
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/tria.h>
 
+#include <deal.II/simplex/fe_lib.h>
+
 #include "../tests.h"
 
 template <int dim, int spacedim>
@@ -35,6 +40,13 @@ process(const std::vector<Point<spacedim>> &vertices,
       deallog << "measure:  " << cell->measure() << std::endl;
       deallog << "diameter: " << cell->diameter() << std::endl;
     }
+
+  const MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+
+  deallog << "diameter_min: " << GridTools::minimal_cell_diameter(tria, mapping)
+          << std::endl;
+  deallog << "diameter_max: " << GridTools::maximal_cell_diameter(tria, mapping)
+          << std::endl;
   deallog << std::endl;
 }
 

diff --git a/tests/simplex/cell_measure_01.with_simplex_support=on.output b/tests/simplex/cell_measure_01.with_simplex_support=on.output
index a8127115b17326558ff2ce7395c031cfa32fdb39..8fe42d75a29721b19019a6fe8df9d1afbde53c8a 100644 (file)
--- a/tests/simplex/cell_measure_01.with_simplex_support=on.output
+++ b/tests/simplex/cell_measure_01.with_simplex_support=on.output
@@ -2,8 +2,12 @@
 DEAL::dim=2 spacedim=2:
 DEAL::measure:  0.500000
 DEAL::diameter: 1.41421
+DEAL::diameter_min: 1.41421
+DEAL::diameter_max: 1.41421
 DEAL::
 DEAL::dim=3 spacedim=3:
 DEAL::measure:  0.166667
 DEAL::diameter: 1.41421
+DEAL::diameter_min: 1.41421
+DEAL::diameter_max: 1.41421
 DEAL::