Add math::FourthOrderTensor

math/BUILD.bazel

@@ -23,6 +23,7 @@ drake_cc_package_library(
         ":discrete_lyapunov_equation",
         ":eigen_sparse_triplet",
         ":evenly_distributed_pts_on_sphere",
+        ":fourth_order_tensor",
         ":geometric_transform",
         ":gradient",
         ":gray_code",
@@ -134,6 +135,16 @@ drake_cc_library(
     ],
 )
 
+drake_cc_library(
+    name = "fourth_order_tensor",
+    srcs = ["fourth_order_tensor.cc"],
+    hdrs = ["fourth_order_tensor.h"],
+    deps = [
+        "//common:default_scalars",
+        "//common:essential",
+    ],
+)
+
 # TODO(jwnimmer-tri) Improved name for this library, "pose_representations"?
 drake_cc_library(
     name = "geometric_transform",
@@ -433,6 +444,14 @@ drake_cc_googletest(
     ],
 )
 
+drake_cc_googletest(
+    name = "fourth_order_tensor_test",
+    deps = [
+        ":fourth_order_tensor",
+        "//common/test_utilities:eigen_matrix_compare",
+    ],
+)
+
 drake_cc_googletest(
     name = "jacobian_test",
     deps = [

math/fourth_order_tensor.cc

@@ -0,0 +1,33 @@
+#include "drake/math/fourth_order_tensor.h"
+
+#include "drake/common/default_scalars.h"
+
+namespace drake {
+namespace math {
+namespace internal {
+
+template <typename T>
+FourthOrderTensor<T>::FourthOrderTensor(const MatrixType& data) : data_(data) {}
+
+template <typename T>
+FourthOrderTensor<T>::FourthOrderTensor() = default;
+
+template <typename T>
+void FourthOrderTensor<T>::ContractWithVectors(
+    const Eigen::Ref<const Vector3<T>>& u,
+    const Eigen::Ref<const Vector3<T>>& v, EigenPtr<Matrix3<T>> B) const {
+  B->setZero();
+  for (int l = 0; l < 3; ++l) {
+    for (int j = 0; j < 3; ++j) {
+      *B += data_.template block<3, 3>(3 * j, 3 * l) * u(j) * v(l);
+    }
+  }
+}
+
+}  // namespace internal
+}  // namespace math
+}  // namespace drake
+
+DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS(
+    class ::drake::math::internal::FourthOrderTensor);
+template class drake::math::internal::FourthOrderTensor<float>;

math/fourth_order_tensor.h

@@ -0,0 +1,108 @@
+#pragma once
+
+#include "drake/common/drake_throw.h"
+#include "drake/common/eigen_types.h"
+
+namespace drake {
+namespace math {
+namespace internal {
+
+/* This class provides functionalities related to 4th-order tensors of
+ dimension 3*3*3*3. The tensor is represented using a a 9*9 matrix that is
+ organized as following
+
+                  l = 1       l = 2       l = 3
+              -------------------------------------
+              |           |           |           |
+    j = 1     |   Aᵢ₁ₖ₁   |   Aᵢ₁ₖ₂   |   Aᵢ₁ₖ₃   |
+              |           |           |           |
+              -------------------------------------
+              |           |           |           |
+    j = 2     |   Aᵢ₂ₖ₁   |   Aᵢ₂ₖ₂   |   Aᵢ₂ₖ₃   |
+              |           |           |           |
+              -------------------------------------
+              |           |           |           |
+    j = 3     |   Aᵢ₃ₖ₁   |   Aᵢ₃ₖ₂   |   Aᵢ₃ₖ₃   |
+              |           |           |           |
+              -------------------------------------
+ Namely the ik-th entry in the jl-th block corresponds to the value Aᵢⱼₖₗ.
+ @tparam float, double, AutoDiffXd. */
+template <typename T>
+class FourthOrderTensor {
+ public:
+  DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(FourthOrderTensor)
+  using MatrixType = Eigen::Matrix<T, 9, 9>;
+
+  /* Constructs a 4th-order tensor represented by the given matrix using the
+   convention layed out in the class documentation. */
+  explicit FourthOrderTensor(const MatrixType& data);
+
+  /* Constructs a zero 4th-order tensor. */
+  FourthOrderTensor();
+
+  /* Performs contraction between this 4th-order tensor A and two vectors u and
+   v and outputs 2nd order tensor B. In Einstein notation, the contraction being
+   done is Bᵢₖ = uⱼ Aᵢⱼₖₗ vₗ. */
+  void ContractWithVectors(const Eigen::Ref<const Vector3<T>>& u,
+                           const Eigen::Ref<const Vector3<T>>& v,
+                           EigenPtr<Matrix3<T>> B) const;
+
+  /* Returns this 4th-order tensor encoded as a matrix according to the class
+   documentation. */
+  const MatrixType& data() const { return data_; }
+
+  /* Returns this 4th-order tensor encoded as a mutable matrix according to the
+   class documentation. */
+  MatrixType& mutable_data() { return data_; }
+
+  /* Returns the value of the 4th-order tensor at the given indices.
+   @pre 0 <= i, j, k, l < 3. */
+  const T& operator()(int i, int j, int k, int l) const {
+    DRAKE_ASSERT(0 <= i && i < 3);
+    DRAKE_ASSERT(0 <= j && j < 3);
+    DRAKE_ASSERT(0 <= k && k < 3);
+    DRAKE_ASSERT(0 <= l && l < 3);
+    return data_(3 * j + i, 3 * l + k);
+  }
+
+  /* Returns the value of the 4th-order tensor at the given indices.
+   @pre 0 <= i, j, k, l < 3. */
+  T& operator()(int i, int j, int k, int l) {
+    DRAKE_ASSERT(0 <= i && i < 3);
+    DRAKE_ASSERT(0 <= j && j < 3);
+    DRAKE_ASSERT(0 <= k && k < 3);
+    DRAKE_ASSERT(0 <= l && l < 3);
+    return data_(3 * j + i, 3 * l + k);
+  }
+
+  /* Returns the value of the 4th-order tensor at the given indices,
+   interpreted as indices into the 9x9 matrix, using the convention layed out in
+   the class documentation.
+   @pre 0 <= i, j < 9. */
+  const T& operator()(int i, int j) const {
+    DRAKE_ASSERT(0 <= i && i < 9);
+    DRAKE_ASSERT(0 <= j && j < 9);
+    return data_(i, j);
+  }
+
+  /* Returns the value of the 4th-order tensor at the given indices,
+   interpreted as indices into the 9x9 matrix, using the convention layed out in
+   the class documentation.
+   @pre 0 <= i, j < 9. */
+  T& operator()(int i, int j) {
+    DRAKE_ASSERT(0 <= i && i < 9);
+    DRAKE_ASSERT(0 <= j && j < 9);
+    return data_(i, j);
+  }
+
+  /* Sets the data of this 4th-order tensor to the given matrix, using the
+   convention layed out in the class documentation. */
+  void set_data(const MatrixType& data) { data_ = data; }
+
+ private:
+  MatrixType data_{MatrixType::Zero()};
+};
+
+}  // namespace internal
+}  // namespace math
+}  // namespace drake

math/test/fourth_order_tensor_test.cc

@@ -0,0 +1,86 @@
+#include "drake/math/fourth_order_tensor.h"
+
+#include "drake/common/eigen_types.h"
+#include "drake/common/test_utilities/eigen_matrix_compare.h"
+
+using Eigen::Matrix3d;
+using Eigen::Vector3d;
+
+namespace drake {
+namespace math {
+namespace internal {
+namespace {
+
+Eigen::Matrix<double, 9, 9> MakeArbitraryMatrix() {
+  Eigen::Matrix<double, 9, 9> data;
+  for (int i = 0; i < 9; ++i) {
+    for (int j = 0; j < 9; ++j) {
+      data(i, j) = i + j;
+    }
+  }
+  return data;
+}
+
+GTEST_TEST(FourthOrderTensorTest, DefaultConstructor) {
+  FourthOrderTensor<double> tensor;
+  EXPECT_TRUE(tensor.data().isZero());
+  const Vector3d u(1.0, 2.0, 3.0);
+  const Vector3d v(4.0, 5.0, 6.0);
+  Matrix3d B;
+
+  tensor.ContractWithVectors(u, v, &B);
+  EXPECT_TRUE(B.isZero());
+}
+
+GTEST_TEST(FourthOrderTensorTest, ConstructWithData) {
+  const Eigen::Matrix<double, 9, 9> data = MakeArbitraryMatrix();
+  FourthOrderTensor<double> tensor(data);
+  /* Getter and mutable getter. */
+  EXPECT_EQ(tensor.data(), data);
+  EXPECT_EQ(tensor.data(), tensor.mutable_data());
+  tensor.mutable_data() = data.transpose();
+  EXPECT_EQ(tensor.data(), data.transpose());
+  /* Settor */
+  tensor.set_data(data);
+  EXPECT_EQ(tensor.data(), data);
+  /* Operator with four indices. */
+  EXPECT_EQ(tensor(0, 0, 0, 0), data(0, 0));
+  EXPECT_EQ(tensor(1, 1, 1, 1), data(4, 4));
+  /* Operator with two indices. */
+  EXPECT_EQ(tensor(0, 0), data(0, 0));
+  EXPECT_EQ(tensor(3, 3), data(3, 3));
+}
+
+GTEST_TEST(FourthOrderTensorTest, ContractWithVectors) {
+  FourthOrderTensor<double> tensor(MakeArbitraryMatrix());
+
+  /* If any vector input is zero, the contraction is zero. */
+  Vector3d u = Vector3d::Zero();
+  Vector3d v(4.0, 5.0, 6.0);
+  Matrix3d B;
+  tensor.ContractWithVectors(u, v, &B);
+  EXPECT_TRUE(B.isZero());
+  tensor.ContractWithVectors(v, u, &B);
+  EXPECT_TRUE(B.isZero());
+
+  /* If the 9x9 representation of the tensor has a repeating pattern in the
+   blocks, the contraction is a multiple of that block. */
+  Matrix3d block;
+  block << 1, 2, 3, 4, 5, 6, 7, 8, 9;
+  Eigen::Matrix<double, 9, 9> data;
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      data.block<3, 3>(3 * i, 3 * j) = block;
+    }
+  }
+  tensor = FourthOrderTensor<double>(data);
+  u << 1.0, 2.0, 3.0;
+  v << 4.0, 5.0, 6.0;
+  tensor.ContractWithVectors(u, v, &B);
+  EXPECT_TRUE(CompareMatrices(B, block * (u * v.transpose()).sum()));
+}
+
+}  // namespace
+}  // namespace internal
+}  // namespace math
+}  // namespace drake

multibody/fem/BUILD.bazel

@@ -83,6 +83,7 @@ drake_cc_library(
     deps = [
         "//common:essential",
         "//common:nice_type_name",
+        "//math:fourth_order_tensor",
     ],
 )
 
@@ -407,6 +408,7 @@ drake_cc_library(
     deps = [
         "//common:default_scalars",
         "//common:essential",
+        "//math:fourth_order_tensor",
     ],
 )
 

multibody/fem/constitutive_model.h

@@ -6,6 +6,7 @@
 
 #include "drake/common/eigen_types.h"
 #include "drake/common/nice_type_name.h"
+#include "drake/math/fourth_order_tensor.h"
 
 namespace drake {
 namespace multibody {
@@ -92,8 +93,8 @@ class ConstitutiveModel {
                    |           |           |           |
                    -------------------------------------
   @pre `dPdF != nullptr`. */
-  void CalcFirstPiolaStressDerivative(const Data& data,
-                                      Eigen::Matrix<T, 9, 9>* dPdF) const {
+  void CalcFirstPiolaStressDerivative(
+      const Data& data, math::internal::FourthOrderTensor<T>* dPdF) const {
     DRAKE_ASSERT(dPdF != nullptr);
     derived().CalcFirstPiolaStressDerivativeImpl(data, dPdF);
   }
@@ -123,8 +124,8 @@ class ConstitutiveModel {
                     NiceTypeName::Get(derived())));
   }
 
-  void CalcFirstPiolaStressDerivativeImpl(const Data& data,
-                                          Eigen::Matrix<T, 9, 9>* dPdF) const {
+  void CalcFirstPiolaStressDerivativeImpl(
+      const Data& data, math::internal::FourthOrderTensor<T>* dPdF) const {
     throw std::logic_error(
         fmt::format("The derived class {} must provide a shadow definition of "
                     "CalcFirstPiolaStressDerivativeImpl() to be correct.",

multibody/fem/corotated_model.cc

@@ -44,7 +44,7 @@ void CorotatedModel<T>::CalcFirstPiolaStressImpl(const Data& data,
 
 template <typename T>
 void CorotatedModel<T>::CalcFirstPiolaStressDerivativeImpl(
-    const Data& data, Eigen::Matrix<T, 9, 9>* dPdF) const {
+    const Data& data, math::internal::FourthOrderTensor<T>* dPdF) const {
   const T& Jm1 = data.Jm1();
   const Matrix3<T>& F = data.deformation_gradient();
   const Matrix3<T>& R = data.R();
@@ -54,9 +54,10 @@ void CorotatedModel<T>::CalcFirstPiolaStressDerivativeImpl(
       Eigen::Map<const Vector<T, 3 * 3>>(JFinvT.data(), 3 * 3);
   auto& local_dPdF = (*dPdF);
   /* The contribution from derivatives of Jm1. */
-  local_dPdF.noalias() = lambda_ * flat_JFinvT * flat_JFinvT.transpose();
+  local_dPdF.mutable_data().noalias() =
+      lambda_ * flat_JFinvT * flat_JFinvT.transpose();
   /* The contribution from derivatives of F. */
-  local_dPdF.diagonal().array() += 2.0 * mu_;
+  local_dPdF.mutable_data().diagonal().array() += 2.0 * mu_;
   /* The contribution from derivatives of R. */
   internal::AddScaledRotationalDerivative<T>(R, S, -2.0 * mu_, &local_dPdF);
   /* The contribution from derivatives of JFinvT. */

multibody/fem/corotated_model.h

@@ -69,8 +69,8 @@ class CorotatedModel final
 
   /* Shadows ConstitutiveModel::CalcFirstPiolaStressDerivativeImpl() as required
    by the CRTP base class. */
-  void CalcFirstPiolaStressDerivativeImpl(const Data& data,
-                                          Eigen::Matrix<T, 9, 9>* dPdF) const;
+  void CalcFirstPiolaStressDerivativeImpl(
+      const Data& data, math::internal::FourthOrderTensor<T>* dPdF) const;
 
   T E_;       // Young's modulus, N/m².
   T nu_;      // Poisson's ratio.

multibody/fem/linear_constitutive_model.cc

@@ -29,7 +29,7 @@ LinearConstitutiveModel<T>::LinearConstitutiveModel(const T& youngs_modulus,
     Keep in mind that the indices are laid out such that the ik-th entry in
     the jl-th block corresponds to the value dPᵢⱼ/dFₖₗ.  */
   /* First term. */
-  dPdF_ = mu_ * Eigen::Matrix<T, 9, 9>::Identity();
+  dPdF_.set_data(mu_ * Eigen::Matrix<T, 9, 9>::Identity());
   for (int k = 0; k < 3; ++k) {
     /* Second term. */
     for (int l = 0; l < 3; ++l) {
@@ -65,7 +65,7 @@ void LinearConstitutiveModel<T>::CalcFirstPiolaStressImpl(const Data& data,
 
 template <typename T>
 void LinearConstitutiveModel<T>::CalcFirstPiolaStressDerivativeImpl(
-    const Data&, Eigen::Matrix<T, 9, 9>* dPdF) const {
+    const Data&, math::internal::FourthOrderTensor<T>* dPdF) const {
   *dPdF = dPdF_;
 }
 

multibody/fem/linear_constitutive_model.h

@@ -71,14 +71,14 @@ class LinearConstitutiveModel final
 
   /* Shadows ConstitutiveModel::CalcFirstPiolaStressDerivativeImpl() as required
    by the CRTP base class. */
-  void CalcFirstPiolaStressDerivativeImpl(const Data& data,
-                                          Eigen::Matrix<T, 9, 9>* dPdF) const;
+  void CalcFirstPiolaStressDerivativeImpl(
+      const Data& data, math::internal::FourthOrderTensor<T>* dPdF) const;
 
   T E_;       // Young's modulus, N/m².
   T nu_;      // Poisson's ratio.
   T mu_;      // Lamé's second parameter/Shear modulus, N/m².
   T lambda_;  // Lamé's first parameter, N/m².
-  Eigen::Matrix<T, 9, 9>
+  math::internal::FourthOrderTensor<T>
       dPdF_;  // The First Piola stress derivative is constant and precomputed.
 };
 

multibody/fem/linear_corotated_model.cc

@@ -82,11 +82,11 @@ void LinearCorotatedModel<T>::CalcFirstPiolaStressImpl(const Data& data,
 */
 template <typename T>
 void LinearCorotatedModel<T>::CalcFirstPiolaStressDerivativeImpl(
-    const Data& data, Eigen::Matrix<T, 9, 9>* dPdF) const {
+    const Data& data, math::internal::FourthOrderTensor<T>* dPdF) const {
   const Matrix3<T>& R0 = data.R0();
   auto& local_dPdF = (*dPdF);
   /* Add in μ * δₐᵢδⱼᵦ. */
-  local_dPdF = mu_ * Eigen::Matrix<T, 9, 9>::Identity();
+  local_dPdF.set_data(mu_ * Eigen::Matrix<T, 9, 9>::Identity());
   for (int i = 0; i < 3; ++i) {
     for (int j = 0; j < 3; ++j) {
       for (int alpha = 0; alpha < 3; ++alpha) {