[geometry] Add precomputed values to VolumeMesh and MeshFieldLinear (#22097)

Author: joemasterjohn

bindings/pydrake/geometry/geometry_py_meshes.cc

@@ -122,6 +122,10 @@ void DoScalarDependentDefinitions(py::module m, T) {
             cls_doc.tetrahedra.doc)
         .def("num_elements", &Class::num_elements, cls_doc.num_elements.doc)
         .def("num_vertices", &Class::num_vertices, cls_doc.num_vertices.doc)
+        .def("inward_normal", &Class::inward_normal, py::arg("e"), py::arg("f"),
+            cls_doc.inward_normal.doc)
+        .def("edge_vector", &Class::edge_vector, py::arg("e"), py::arg("a"),
+            py::arg("b"), cls_doc.edge_vector.doc)
         .def("CalcTetrahedronVolume", &Class::CalcTetrahedronVolume,
             py::arg("e"), cls_doc.CalcTetrahedronVolume.doc)
         .def("CalcVolume", &Class::CalcVolume, cls_doc.CalcVolume.doc)
@@ -177,7 +181,8 @@ void DoScalarIndependentDefinitions(py::module m) {
         .def(py::init<int, int, int, int>(), py::arg("v0"), py::arg("v1"),
             py::arg("v2"), py::arg("v3"), cls_doc.ctor.doc_4args)
         // TODO(SeanCurtis-TRI): Bind constructor that takes array of ints.
-        .def("vertex", &Class::vertex, py::arg("i"), cls_doc.vertex.doc);
+        .def("vertex", &Class::vertex, py::arg("i"), cls_doc.vertex.doc)
+        .def("num_vertices", &Class::num_vertices, cls_doc.num_vertices.doc);
     DefCopyAndDeepCopy(&cls);
   }
 }

bindings/pydrake/geometry/test/meshes_test.py

@@ -1,5 +1,6 @@
 import pydrake.geometry as mut
 import pydrake.geometry._testing as mut_testing
+from pydrake.common.test_utilities import numpy_compare
 
 import copy
 import unittest
@@ -141,6 +142,8 @@ def test_volume_mesh(self):
         self.assertIsInstance(dut.vertex(v=0), np.ndarray)
         self.assertIsInstance(dut.num_elements(), int)
         self.assertIsInstance(dut.num_vertices(), int)
+        self.assertIsInstance(dut.inward_normal(e=0, f=0), np.ndarray)
+        self.assertIsInstance(dut.edge_vector(e=0, a=0, b=1), np.ndarray)
 
         # Sanity check some calculations
         self.assertAlmostEqual(dut.CalcTetrahedronVolume(e=1), 1/6.0,
@@ -154,9 +157,19 @@ def test_volume_mesh(self):
         self.assertIsInstance(dut.tetrahedra()[0], mut.VolumeElement)
         self.assertEqual(len(dut.vertices()), 5)
 
+        numpy_compare.assert_float_equal(
+            dut.inward_normal(e=0, f=3), (-1., 0., 0.))
+        numpy_compare.assert_float_equal(
+            dut.inward_normal(e=1, f=3), (1., 0., 0.))
+        numpy_compare.assert_float_equal(
+            dut.edge_vector(e=0, a=1, b=3), (-1., 0., 0.))
+        numpy_compare.assert_float_equal(
+            dut.edge_vector(e=1, a=1, b=3), (1., 0., 0.))
+
         # Now check the VolumeElement bindings.
         tetrahedron0 = dut.element(e=0)
         self.assertEqual(tetrahedron0.vertex(i=0), 2)
+        self.assertEqual(tetrahedron0.num_vertices(), 4)
 
     def test_convert_volume_to_surface_mesh(self):
         # Use the volume mesh from `test_volume_mesh()`.

geometry/proximity/mesh_field_linear.h

@@ -198,6 +198,7 @@ class MeshFieldLinear {
                      static_cast<int>(values_at_Mo_.size()));
       }
     }
+    CalcMinAndMaxValues();
   }
 
   /** (Advanced) Constructor variant which receives the pre-computed,
@@ -224,6 +225,7 @@ class MeshFieldLinear {
     DRAKE_DEMAND(static_cast<int>(gradients_.size()) == mesh_->num_elements());
 
     CalcValueAtMeshOriginForAllElements();
+    CalcMinAndMaxValues();
   }
 
   /** @returns true iff the gradient field could not be computed, and the mesh
@@ -238,10 +240,40 @@ class MeshFieldLinear {
    @pre v ∈ [0, this->mesh().num_vertices()).
    */
   const T& EvaluateAtVertex(int v) const {
-    DRAKE_ASSERT(v >= 0 && v < mesh_->num_vertices());
+    DRAKE_DEMAND(v >= 0 && v < mesh_->num_vertices());
     return values_[v];
   }
 
+  /** (Advanced) Evaluates the minimum field value on an element.
+   @param e The index of the element.
+   @pre e ∈ [0, this->mesh().num_elements()).
+   */
+  const T& EvaluateMin(int e) const {
+    DRAKE_DEMAND(e >= 0 && e < mesh_->num_elements());
+    return min_values_[e];
+  }
+
+  /** (Advanced) Evaluates the maximum field value on an element.
+   @param e The index of the element.
+   @pre e ∈ [0, this->mesh().num_elements()).
+   */
+  const T& EvaluateMax(int e) const {
+    DRAKE_DEMAND(e >= 0 && e < mesh_->num_elements());
+    return max_values_[e];
+  }
+
+  /** (Advanced) Evaluates the linear function associated with element e at the
+   mesh origin Mo.
+   @param e The index of the element.
+   @pre e ∈ [0, this->mesh().num_elements()).
+   @pre The field has valid gradients.
+  */
+  const T& EvaluateAtMo(int e) const {
+    DRAKE_DEMAND(e >= 0 && e < mesh_->num_elements());
+    DRAKE_DEMAND(e < ssize(values_at_Mo_));
+    return values_at_Mo_[e];
+  }
+
   /** Evaluates the field value at a location on an element.
 
    The return type depends on both the field's scalar type `T` and the
@@ -361,8 +393,14 @@ class MeshFieldLinear {
     return new_mesh_field;
   }
 
+  /** The mesh M to which this field refers. */
   const MeshType& mesh() const { return *mesh_; }
+  /** The field value at each vertex. */
   const std::vector<T>& values() const { return values_; }
+  /** The minimum field value on each element. */
+  const std::vector<T>& min_values() const { return min_values_; }
+  /** The maximum field value on each element. */
+  const std::vector<T>& max_values() const { return max_values_; }
 
   // TODO(#12173): Consider NaN==NaN to be true in equality tests.
   /** Checks to see whether the given MeshFieldLinear object is equal via deep
@@ -416,6 +454,25 @@ class MeshFieldLinear {
     }
   }
 
+  void CalcMinAndMaxValues() {
+    min_values_.clear();
+    max_values_.clear();
+    min_values_.reserve(this->mesh().num_elements());
+    max_values_.reserve(this->mesh().num_elements());
+    for (int e = 0; e < this->mesh().num_elements(); ++e) {
+      T min, max;
+      min = max = values_[this->mesh().element(e).vertex(0)];
+
+      for (int i = 1; i < mesh().element(e).num_vertices(); ++i) {
+        min = std::min(min, values_[this->mesh().element(e).vertex(i)]);
+        max = std::max(max, values_[this->mesh().element(e).vertex(i)]);
+      }
+
+      min_values_.emplace_back(min);
+      max_values_.emplace_back(max);
+    }
+  }
+
   std::optional<Vector3<T>> MaybeCalcGradientVector(int e) const {
     // In the case of the PolygonSurfaceMesh, where kVertexPerElement is marked
     // as "indeterminate" (aka -1), we'll simply use the first three vertices.
@@ -454,6 +511,12 @@ class MeshFieldLinear {
   // The field values are indexed in the same way as vertices, i.e.,
   // values_[i] is the field value for the mesh vertices_[i].
   std::vector<T> values_;
+  // Stores the minimum value of the field for each element, i.e.,
+  // min_values_[i] is the minimum field value on the domain of elements_[i].
+  std::vector<T> min_values_;
+  // Stores the maximum value of the field for each element, i.e.,
+  // min_values_[i] is the minimum field value on the domain of elements_[i].
+  std::vector<T> max_values_;
   // The gradients are indexed in the same way as elements, i.e.,
   // gradients_[i] is the gradient vector on elements_[i]. The elements could
   // be tetrahedra for VolumeMesh or triangles for TriangleSurfaceMesh.

geometry/proximity/test/mesh_field_linear_test.cc

@@ -68,6 +68,23 @@ GTEST_TEST(MeshFieldLinearTest, EvaluateAtVertex) {
   EXPECT_EQ(mesh_field->EvaluateAtVertex(3), 3);
 }
 
+GTEST_TEST(MeshFieldLinearTest, EvaluateMinAndMaxAndMo) {
+  auto mesh = GenerateMesh<double>();
+  constexpr double kEps = std::numeric_limits<double>::epsilon();
+
+  // Arbitrary values such that min and max are unique on each element.
+  std::vector<double> e_values = {1., -2., 2., 3.};
+  auto mesh_field =
+      std::make_unique<MeshFieldLinear<double, TriangleSurfaceMesh<double>>>(
+          std::move(e_values), mesh.get());
+  EXPECT_EQ(mesh_field->EvaluateMin(0), -2);
+  EXPECT_EQ(mesh_field->EvaluateMin(1), 1);
+  EXPECT_EQ(mesh_field->EvaluateMax(0), 2);
+  EXPECT_EQ(mesh_field->EvaluateMax(1), 3);
+  EXPECT_NEAR(mesh_field->EvaluateAtMo(0), 1, kEps);
+  EXPECT_NEAR(mesh_field->EvaluateAtMo(1), 1, kEps);
+}
+
 // Tests CloneAndSetMesh(). We use `double` and TriangleSurfaceMesh<double> as
 // representative arguments for type parameters.
 GTEST_TEST(MeshFieldLinearTest, TestDoCloneWithMesh) {