Implement TriaAccessor::measure() and TriaAccessor::diameter() for simplex

diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h
index 1ee14383773ba5ac6f0c5f1f3f5d7f08c9e332bf..a8e940da23f7f4d77a6e15ad2353a8917f4b3386 100644 (file)
--- a/include/deal.II/grid/tria_accessor.templates.h
+++ b/include/deal.II/grid/tria_accessor.templates.h
@@ -1885,14 +1885,26 @@ template <int structdim, int dim, int spacedim>
 double
 TriaAccessor<structdim, dim, spacedim>::diameter() const
 {
-  switch (structdim)
+  switch (this->reference_cell_type())
     {
-      case 1:
+      case ReferenceCell::Type::Line:
         return (this->vertex(1) - this->vertex(0)).norm();
-      case 2:
+      case ReferenceCell::Type::Tri:
+        return std::max(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 std::max((this->vertex(3) - this->vertex(0)).norm(),
                         (this->vertex(2) - this->vertex(1)).norm());
-      case 3:
+      case ReferenceCell::Type::Tet:
+        return std::max(
+          std::max(std::max((this->vertex(1) - this->vertex(0)).norm(),
+                            (this->vertex(2) - this->vertex(0)).norm()),
+                   std::max((this->vertex(2) - this->vertex(1)).norm(),
+                            (this->vertex(3) - this->vertex(0)).norm())),
+          std::max((this->vertex(3) - this->vertex(1)).norm(),
+                   (this->vertex(3) - this->vertex(2)).norm()));
+      case ReferenceCell::Type::Hex:
         return std::max(std::max((this->vertex(7) - this->vertex(0)).norm(),
                                  (this->vertex(6) - this->vertex(1)).norm()),
                         std::max((this->vertex(2) - this->vertex(5)).norm(),

diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc
index 126ea9ec8a2d9adc6ee13fa4faf6b4f0595b7ce0..1bdc089ebe3c99877f19dc48614ca0ace848596d 100644 (file)
--- a/source/grid/grid_tools_nontemplates.cc
+++ b/source/grid/grid_tools_nontemplates.cc
@@ -45,6 +45,20 @@ namespace GridTools
   cell_measure<2>(const std::vector<Point<2>> &        all_vertices,
                   const ArrayView<const unsigned int> &vertex_indices)
   {
+    if (vertex_indices.size() == 3) // triangle
+      {
+        const double x[3] = {all_vertices[vertex_indices[0]](0),
+                             all_vertices[vertex_indices[1]](0),
+                             all_vertices[vertex_indices[2]](0)};
+
+        const double y[3] = {all_vertices[vertex_indices[0]](1),
+                             all_vertices[vertex_indices[1]](1),
+                             all_vertices[vertex_indices[2]](1)};
+
+        return 0.5 * std::abs((x[0] - x[2]) * (y[1] - y[0]) -
+                              (x[0] - x[1]) * (y[2] - y[0]));
+      }
+
     AssertDimension(vertex_indices.size(), GeometryInfo<2>::vertices_per_cell);
 
     /*
@@ -105,6 +119,16 @@ namespace GridTools
   cell_measure<3>(const std::vector<Point<3>> &        all_vertices,
                   const ArrayView<const unsigned int> &vertex_indices)
   {
+    if (vertex_indices.size() == 4) // tetrahedron
+      {
+        const auto &a = all_vertices[vertex_indices[0]];
+        const auto &b = all_vertices[vertex_indices[1]];
+        const auto &c = all_vertices[vertex_indices[2]];
+        const auto &d = all_vertices[vertex_indices[3]];
+
+        return (1.0 / 6.0) * std::abs((a - d) * cross_product_3d(b - d, c - d));
+      }
+
     AssertDimension(vertex_indices.size(), GeometryInfo<3>::vertices_per_cell);
     // note that this is the
     // cell_measure based on the new

diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index a693d61ccb6a24db315838d90e02796d58623e74..d14bfcdb6e8882b99e8e3fb81ba7538efcea7e6b 100644 (file)
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -1269,7 +1269,8 @@ namespace
       vertex_indices[i] = accessor.vertex_index(i);
 
     return GridTools::cell_measure<2>(
-      accessor.get_triangulation().get_vertices(), vertex_indices);
+      accessor.get_triangulation().get_vertices(),
+      ArrayView<unsigned int>(vertex_indices, accessor.n_vertices()));
   }
 
 
@@ -1281,7 +1282,8 @@ namespace
       vertex_indices[i] = accessor.vertex_index(i);
 
     return GridTools::cell_measure<3>(
-      accessor.get_triangulation().get_vertices(), vertex_indices);
+      accessor.get_triangulation().get_vertices(),
+      ArrayView<unsigned int>(vertex_indices, accessor.n_vertices()));
   }
 
 

diff --git a/tests/simplex/cell_measure_01.cc b/tests/simplex/cell_measure_01.cc
new file mode 100644 (file)
index 0000000..08a57fd
--- /dev/null
+++ b/tests/simplex/cell_measure_01.cc
@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test TriaAccessor::measure() and TriaAccessor::diameter().
+
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+template <int dim, int spacedim>
+void
+process(const std::vector<Point<spacedim>> &vertices,
+        const std::vector<CellData<dim>> &  cells)
+{
+  Triangulation<dim, spacedim> tria;
+  tria.create_triangulation(vertices, cells, SubCellData());
+
+  deallog << "dim=" << dim << " spacedim=" << spacedim << ":" << std::endl;
+  for (const auto &cell : tria.active_cell_iterators())
+    {
+      deallog << "measure:  " << cell->measure() << std::endl;
+      deallog << "diameter: " << cell->diameter() << std::endl;
+    }
+  deallog << std::endl;
+}
+
+template <int dim>
+void
+test()
+{
+  Assert(false, ExcNotImplemented());
+}
+
+template <>
+void
+test<2>()
+{
+  const int dim      = 2;
+  const int spacedim = 2;
+
+  std::vector<Point<spacedim>> vertices;
+  vertices.emplace_back(0, 0);
+  vertices.emplace_back(1, 0);
+  vertices.emplace_back(0, 1);
+
+  std::vector<CellData<dim>> cells;
+  CellData<dim>              cell;
+  cell.vertices = {0, 1, 2};
+  cells.push_back(cell);
+
+  process(vertices, cells);
+}
+
+template <>
+void
+test<3>()
+{
+  const int dim      = 3;
+  const int spacedim = 3;
+
+  std::vector<Point<spacedim>> vertices;
+  vertices.emplace_back(0, 0, 0);
+  vertices.emplace_back(1, 0, 0);
+  vertices.emplace_back(0, 1, 0);
+  vertices.emplace_back(0, 0, 1);
+
+  std::vector<CellData<dim>> cells;
+  CellData<dim>              cell;
+  cell.vertices = {0, 1, 2, 3};
+  cells.push_back(cell);
+
+  process(vertices, cells);
+}
+
+int
+main()
+{
+  initlog();
+
+  test<2>();
+  test<3>();
+}

diff --git a/tests/simplex/cell_measure_01.with_simplex_support=on.output b/tests/simplex/cell_measure_01.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..a812711
--- /dev/null
+++ b/tests/simplex/cell_measure_01.with_simplex_support=on.output
@@ -0,0 +1,9 @@
+
+DEAL::dim=2 spacedim=2:
+DEAL::measure:  0.500000
+DEAL::diameter: 1.41421
+DEAL::
+DEAL::dim=3 spacedim=3:
+DEAL::measure:  0.166667
+DEAL::diameter: 1.41421
+DEAL::