RobotLocomotion · jwnimmer-tri · Oct 6, 2025 · Oct 6, 2025
diff --git a/bindings/generated_docstrings/geometry_proximity.h b/bindings/generated_docstrings/geometry_proximity.h
@@ -485,7 +485,10 @@ Parameter ``p_MQ``:
 
 Raises:
     RuntimeError if the field does not have gradients defined *and*
-    the MeshType doesn't support Barycentric coordinates.)""";
+    the MeshType doesn't support Barycentric coordinates.
+
+Template parameter ``C``:
+    must be either ``double`` or ``AutoDiffXd``.)""";
         } EvaluateCartesian;
         // Symbol: drake::geometry::MeshFieldLinear::EvaluateGradient
         struct /* EvaluateGradient */ {
@@ -871,7 +874,10 @@ elements (three or more sides).)""";
 R"""(See TriangleSurfaceMesh::CalcBaryCentric(). This implementation is
 provided to maintain compatibility with MeshFieldLinear. However, it
 only throws. PolygonSurfaceMesh does not support barycentric
-coordinates.)""";
+coordinates.
+
+Template parameter ``C``:
+    must be either ``double`` or ``AutoDiffXd``.)""";
         } CalcBarycentric;
         // Symbol: drake::geometry::PolygonSurfaceMesh::CalcBoundingBox
         struct /* CalcBoundingBox */ {
@@ -1292,7 +1298,7 @@ Parameter ``t``:
     The index of a triangle.
 
 Returns ``b_Q``:
-    ' The barycentric coordinates of Q' (projection of Q onto `t`'s
+    ' The barycentric coordinates of Q' (projection of Q onto ``t`'s
     plane) relative to triangle t.
 
 Note:
@@ -1301,7 +1307,10 @@ Returns ``b_Q``:
     negative.
 
 Precondition:
-    t ∈ {0, 1, 2,..., num_triangles()-1}.)""";
+    t ∈ {0, 1, 2,..., num_triangles()-1}.
+
+Template parameter ``C``:
+    must be either `double`` or ``AutoDiffXd``.)""";
         } CalcBarycentric;
         // Symbol: drake::geometry::TriangleSurfaceMesh::CalcBoundingBox
         struct /* CalcBoundingBox */ {
@@ -1616,7 +1625,10 @@ Parameter ``e``:
 Note:
     If p_MQ is outside the tetrahedral element, the barycentric
     coordinates (b₀, b₁, b₂, b₃) still satisfy b₀ + b₁ + b₂ + b₃ = 1;
-    however, some bᵢ will be negative.)""";
+    however, some bᵢ will be negative.
+
+Template parameter ``C``:
+    must be either ``double`` or ``AutoDiffXd``.)""";
         } CalcBarycentric;
         // Symbol: drake::geometry::VolumeMesh::CalcGradientVectorOfLinearField
         struct /* CalcGradientVectorOfLinearField */ {

diff --git a/geometry/proximity/mesh_field_linear.h b/geometry/proximity/mesh_field_linear.h
@@ -330,10 +330,12 @@ class MeshFieldLinear {
                coordinates. M is the frame of the mesh.
    @throws std::exception if the field does not have gradients defined _and_ the
            MeshType doesn't support Barycentric coordinates.
-   */
+   @tparam C must be either `double` or `AutoDiffXd`. */
   template <typename C>
   promoted_numerical_t<C, T> EvaluateCartesian(int e,
-                                               const Vector3<C>& p_MQ) const {
+                                               const Vector3<C>& p_MQ) const
+    requires scalar_predicate<C>::is_bool
+  {
     if (is_gradient_field_degenerate_) {
       throw std::runtime_error("Gradient field is degenerate.");
     }

diff --git a/geometry/proximity/polygon_surface_mesh.cc b/geometry/proximity/polygon_surface_mesh.cc
@@ -52,6 +52,18 @@ void PolygonSurfaceMesh<T>::ReverseFaceWinding() {
   }
 }
 
+template <typename T>
+template <typename C>
+typename PolygonSurfaceMesh<T>::template Barycentric<promoted_numerical_t<T, C>>
+PolygonSurfaceMesh<T>::CalcBarycentric(const Vector3<C>& p_MQ, int p) const
+  requires scalar_predicate<C>::is_bool
+{
+  unused(p_MQ, p);
+  throw std::runtime_error(
+      "PolygonSurfaceMesh::CalcBarycentric(): PolygonSurfaceMesh does not have "
+      "barycentric coordinates.");
+}
+
 template <typename T>
 std::pair<Vector3<T>, Vector3<T>> PolygonSurfaceMesh<T>::CalcBoundingBox()
     const {
@@ -243,5 +255,8 @@ void PolygonSurfaceMesh<T>::SetAllPositions(
 DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS(
     class PolygonSurfaceMesh);
 
+DRAKE_DEFINE_FUNCTION_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS(
+    (&PolygonSurfaceMesh<T>::template CalcBarycentric<U>));
+
 }  // namespace geometry
 }  // namespace drake
diff --git a/geometry/proximity/polygon_surface_mesh.h b/geometry/proximity/polygon_surface_mesh.h
@@ -233,15 +233,12 @@ class PolygonSurfaceMesh {
 
   /** See TriangleSurfaceMesh::CalcBaryCentric(). This implementation is
    provided to maintain compatibility with MeshFieldLinear. However, it only
-   throws. %PolygonSurfaceMesh does not support barycentric coordinates. */
+   throws. %PolygonSurfaceMesh does not support barycentric coordinates.
+   @tparam C must be either `double` or `AutoDiffXd`. */
   template <typename C>
   Barycentric<promoted_numerical_t<T, C>> CalcBarycentric(
-      const Vector3<C>& p_MQ, int p) const {
-    unused(p_MQ, p);
-    throw std::runtime_error(
-        "PolygonSurfaceMesh::CalcBarycentric(): PolygonSurfaceMesh does not "
-        "have barycentric coordinates.");
-  }
+      const Vector3<C>& p_MQ, int p) const
+    requires scalar_predicate<C>::is_bool;
 
   // TODO(DamrongGuoy): Consider using an oriented bounding box in obb.h.
   //  Currently we have a problem that PolygonSurfaceMesh and its vertices are

diff --git a/geometry/proximity/triangle_surface_mesh.cc b/geometry/proximity/triangle_surface_mesh.cc
@@ -15,6 +15,69 @@ void TriangleSurfaceMesh<T>::ReverseFaceWinding() {
   }
 }
 
+template <typename T>
+template <typename C>
+typename TriangleSurfaceMesh<T>::template Barycentric<
+    promoted_numerical_t<T, C>>
+TriangleSurfaceMesh<T>::CalcBarycentric(const Vector3<C>& p_MQ, int t) const
+  requires scalar_predicate<C>::is_bool
+{
+  const Vector3<T>& v0 = vertex(element(t).vertex(0));
+  const Vector3<T>& v1 = vertex(element(t).vertex(1));
+  const Vector3<T>& v2 = vertex(element(t).vertex(2));
+  // Translate the triangle to the origin to simplify calculations;
+  // barycentric coordinates stay the same.
+  //     u⃗i = v⃗i - v0
+  //     p_MR = p_MQ - v0
+  //
+  // Consider R' on the spanning plane through the origin, u1, u2:
+  //     R' = b₀*u0 + b₁*u1 + b₂*u2
+  //        = 0 + b₁*u1 + b₂*u2
+  //        = b₁*u1 + b₂*u2
+  //
+  // Solve for b₁, b₂ that give R' "closest" to R in the least square sense:
+  //
+  //      |      ||b1|
+  //      |u⃗1  u⃗2||b2| ~ R'
+  //      |      |
+  //
+  // return Barycentric (1-b₁-b₂, b₁, b₂)
+  //
+  using ReturnType = promoted_numerical_t<T, C>;
+  Eigen::Matrix<ReturnType, 3, 2> A;
+  A.col(0) << v1 - v0;
+  A.col(1) << v2 - v0;
+  Vector2<ReturnType> solution = A.colPivHouseholderQr().solve(p_MQ - v0);
+
+  const ReturnType& b1 = solution(0);
+  const ReturnType& b2 = solution(1);
+  const ReturnType b0 = T(1.) - b1 - b2;
+  return {b0, b1, b2};
+}
+// TODO(DamrongGuoy): Investigate alternative calculation suggested by
+//  Alejandro Castro:
+// 1. Starting with the same ui and p_MR.
+// 2. Calculate the unit normal vector n to the spanning plane S through
+//    the origin, u1, and u2.
+//        n = u1.cross(u2).normalize().
+// 3. Project p_MR to p_MR' on the plane S,
+//        p_MR' = p_MR - (p_MR.dot(n))*n
+//
+// Now we have p_MR' = b₀*u⃗0 + b₁*u⃗1 + b₂*u⃗2 by barycentric coordinates.
+//                   =   0   + b₁*u1 + b₂*u2
+//
+// 5. Solve for b₁ and b₂.
+//        (b₁*u1 + b₂*u2).dot(u1) = p_MR'.dot(u1)
+//        (b₁*u1 + b₂*u2).dot(u2) = p_MR'.dot(u2)
+//    Therefore, the 2x2 system:
+//        |u1.dot(u1)  u2.dot(u1)||b1| = |p_MR'.dot(u1)|
+//        |u1.dot(u2)  u2.dot(u2)||b2|   |p_MR'.dot(u2)|
+//
+// 6. return Barycentric(1-b₁-b₂, b₁, b₂)
+//
+// Optimization: save n, and the inverse of matrix |uᵢ.dot(uⱼ)| for later.
+//
+
 template <typename T>
 void TriangleSurfaceMesh<T>::SetAllPositions(
     const Eigen::Ref<const VectorX<T>>& p_MVs) {
@@ -35,5 +98,8 @@ void TriangleSurfaceMesh<T>::SetAllPositions(
 DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS(
     class TriangleSurfaceMesh);
 
+DRAKE_DEFINE_FUNCTION_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS(
+    (&TriangleSurfaceMesh<T>::template CalcBarycentric<U>));
+
 }  // namespace geometry
 }  // namespace drake
diff --git a/geometry/proximity/triangle_surface_mesh.h b/geometry/proximity/triangle_surface_mesh.h
@@ -287,65 +287,11 @@ class TriangleSurfaceMesh {
           (b₀, b₁, b₂) still satisfy b₀ + b₁ + b₂ = 1; however, some bᵢ will be
           negative.
    @pre t ∈ {0, 1, 2,..., num_triangles()-1}.
-   */
+   @tparam C must be either `double` or `AutoDiffXd`. */
   template <typename C>
   Barycentric<promoted_numerical_t<T, C>> CalcBarycentric(
-      const Vector3<C>& p_MQ, int t) const {
-    const Vector3<T>& v0 = vertex(element(t).vertex(0));
-    const Vector3<T>& v1 = vertex(element(t).vertex(1));
-    const Vector3<T>& v2 = vertex(element(t).vertex(2));
-    // Translate the triangle to the origin to simplify calculations;
-    // barycentric coordinates stay the same.
-    //     u⃗i = v⃗i - v0
-    //     p_MR = p_MQ - v0
-    //
-    // Consider R' on the spanning plane through the origin, u1, u2:
-    //     R' = b₀*u0 + b₁*u1 + b₂*u2
-    //        = 0 + b₁*u1 + b₂*u2
-    //        = b₁*u1 + b₂*u2
-    //
-    // Solve for b₁, b₂ that give R' "closest" to R in the least square sense:
-    //
-    //      |      ||b1|
-    //      |u⃗1  u⃗2||b2| ~ R'
-    //      |      |
-    //
-    // return Barycentric (1-b₁-b₂, b₁, b₂)
-    //
-    using ReturnType = promoted_numerical_t<T, C>;
-    Eigen::Matrix<ReturnType, 3, 2> A;
-    A.col(0) << v1 - v0;
-    A.col(1) << v2 - v0;
-    Vector2<ReturnType> solution = A.colPivHouseholderQr().solve(p_MQ - v0);
-
-    const ReturnType& b1 = solution(0);
-    const ReturnType& b2 = solution(1);
-    const ReturnType b0 = T(1.) - b1 - b2;
-    return {b0, b1, b2};
-  }
-  // TODO(DamrongGuoy): Investigate alternative calculation suggested by
-  //  Alejandro Castro:
-  // 1. Starting with the same ui and p_MR.
-  // 2. Calculate the unit normal vector n to the spanning plane S through
-  //    the origin, u1, and u2.
-  //        n = u1.cross(u2).normalize().
-  // 3. Project p_MR to p_MR' on the plane S,
-  //        p_MR' = p_MR - (p_MR.dot(n))*n
-  //
-  // Now we have p_MR' = b₀*u⃗0 + b₁*u⃗1 + b₂*u⃗2 by barycentric coordinates.
-  //                   =   0   + b₁*u1 + b₂*u2
-  //
-  // 5. Solve for b₁ and b₂.
-  //        (b₁*u1 + b₂*u2).dot(u1) = p_MR'.dot(u1)
-  //        (b₁*u1 + b₂*u2).dot(u2) = p_MR'.dot(u2)
-  //    Therefore, the 2x2 system:
-  //        |u1.dot(u1)  u2.dot(u1)||b1| = |p_MR'.dot(u1)|
-  //        |u1.dot(u2)  u2.dot(u2)||b2|   |p_MR'.dot(u2)|
-  //
-  // 6. return Barycentric(1-b₁-b₂, b₁, b₂)
-  //
-  // Optimization: save n, and the inverse of matrix |uᵢ.dot(uⱼ)| for later.
-  //
+      const Vector3<C>& p_MQ, int t) const
+    requires scalar_predicate<C>::is_bool;
 
   // TODO(DamrongGuoy): Consider using an oriented bounding box in obb.h.
   //  Currently we have a problem that TriangleSurfaceMesh and its vertices are

diff --git a/geometry/proximity/volume_mesh.cc b/geometry/proximity/volume_mesh.cc
@@ -1,6 +1,7 @@
 #include "drake/geometry/proximity/volume_mesh.h"
 
 #include "drake/common/default_scalars.h"
+#include "drake/math/linear_solve.h"
 
 namespace drake {
 namespace geometry {
@@ -15,6 +16,39 @@ VolumeMesh<T>::VolumeMesh(std::vector<VolumeElement>&& elements,
   ComputePositionDependentQuantities();
 }
 
+template <typename T>
+template <typename C>
+typename VolumeMesh<T>::template Barycentric<promoted_numerical_t<T, C>>
+VolumeMesh<T>::CalcBarycentric(const Vector3<C>& p_MQ, int e) const
+  requires scalar_predicate<C>::is_bool
+{
+  // We have two conditions to satisfy.
+  // 1. b₀ + b₁ + b₂ + b₃ = 1
+  // 2. b₀*v0 + b₁*v1 + b₂*v2 + b₃*v3 = p_M.
+  // Together they create this 4x4 linear system:
+  //
+  //      | 1  1  1  1 ||b₀|   | 1 |
+  //      | |  |  |  | ||b₁| = | | |
+  //      | v0 v1 v2 v3||b₂|   |p_M|
+  //      | |  |  |  | ||b₃|   | | |
+  //
+  // q = p_M - v0 = b₀*u0 + b₁*u1 + b₂*u2 + b₃*u3
+  //              = 0 + b₁*u1 + b₂*u2 + b₃*u3
+  using ReturnType = promoted_numerical_t<T, C>;
+  Matrix4<ReturnType> A;
+  for (int i = 0; i < 4; ++i) {
+    A.col(i) << ReturnType(1.0), vertex(element(e).vertex(i));
+  }
+  Vector4<ReturnType> b;
+  b << ReturnType(1.0), p_MQ;
+  const math::LinearSolver<Eigen::PartialPivLU, Matrix4<ReturnType>> A_lu(A);
+  const Vector4<ReturnType> b_Q = A_lu.Solve(b);
+  // TODO(DamrongGuoy): Save the inverse of the matrix instead of
+  //  calculating it on the fly. We can reduce to 3x3 system too.  See
+  //  issue #11653.
+  return b_Q;
+}
+
 template <typename T>
 void VolumeMesh<T>::TransformVertices(
     const math::RigidTransform<T>& transform) {
@@ -127,5 +161,8 @@ bool VolumeMesh<T>::Equal(const VolumeMesh<T>& mesh,
 DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS(
     class VolumeMesh);
 
+DRAKE_DEFINE_FUNCTION_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS(
+    (&VolumeMesh<T>::template CalcBarycentric<U>));
+
 }  // namespace geometry
 }  // namespace drake
diff --git a/geometry/proximity/volume_mesh.h b/geometry/proximity/volume_mesh.h
@@ -13,7 +13,6 @@
 #include "drake/common/drake_copyable.h"
 #include "drake/common/eigen_types.h"
 #include "drake/geometry/proximity/mesh_traits.h"
-#include "drake/math/linear_solve.h"
 #include "drake/math/rigid_transform.h"
 
 namespace drake {
@@ -241,36 +240,11 @@ class VolumeMesh {
    @note  If p_MQ is outside the tetrahedral element, the barycentric
           coordinates (b₀, b₁, b₂, b₃) still satisfy b₀ + b₁ + b₂ + b₃ = 1;
           however, some bᵢ will be negative.
-   */
+   @tparam C must be either `double` or `AutoDiffXd`. */
   template <typename C>
   Barycentric<promoted_numerical_t<T, C>> CalcBarycentric(
-      const Vector3<C>& p_MQ, int e) const {
-    // We have two conditions to satisfy.
-    // 1. b₀ + b₁ + b₂ + b₃ = 1
-    // 2. b₀*v0 + b₁*v1 + b₂*v2 + b₃*v3 = p_M.
-    // Together they create this 4x4 linear system:
-    //
-    //      | 1  1  1  1 ||b₀|   | 1 |
-    //      | |  |  |  | ||b₁| = | | |
-    //      | v0 v1 v2 v3||b₂|   |p_M|
-    //      | |  |  |  | ||b₃|   | | |
-    //
-    // q = p_M - v0 = b₀*u0 + b₁*u1 + b₂*u2 + b₃*u3
-    //              = 0 + b₁*u1 + b₂*u2 + b₃*u3
-    using ReturnType = promoted_numerical_t<T, C>;
-    Matrix4<ReturnType> A;
-    for (int i = 0; i < 4; ++i) {
-      A.col(i) << ReturnType(1.0), vertex(element(e).vertex(i));
-    }
-    Vector4<ReturnType> b;
-    b << ReturnType(1.0), p_MQ;
-    const math::LinearSolver<Eigen::PartialPivLU, Matrix4<ReturnType>> A_lu(A);
-    const Vector4<ReturnType> b_Q = A_lu.Solve(b);
-    // TODO(DamrongGuoy): Save the inverse of the matrix instead of
-    //  calculating it on the fly. We can reduce to 3x3 system too.  See
-    //  issue #11653.
-    return b_Q;
-  }
+      const Vector3<C>& p_MQ, int e) const
+    requires scalar_predicate<C>::is_bool;
 
   /** Checks to see whether the given VolumeMesh object is equal via deep
    comparison (up to a tolerance). NaNs are treated as not equal as per the IEEE