From fad03389f3ece3d566dcd5fcfa5ae84ca5d55f30 Mon Sep 17 00:00:00 2001 From: Xuchen Han Date: Thu, 14 Nov 2024 17:58:17 -0800 Subject: [PATCH] Add math::FourthOrderTensor --- math/BUILD.bazel | 19 ++++ math/fourth_order_tensor.cc | 31 +++++ math/fourth_order_tensor.h | 106 ++++++++++++++++++ math/test/fourth_order_tensor_test.cc | 76 +++++++++++++ multibody/fem/BUILD.bazel | 2 + multibody/fem/constitutive_model.h | 7 +- multibody/fem/corotated_model.cc | 7 +- multibody/fem/corotated_model.h | 4 +- multibody/fem/linear_constitutive_model.cc | 4 +- multibody/fem/linear_constitutive_model.h | 6 +- multibody/fem/linear_corotated_model.cc | 4 +- multibody/fem/linear_corotated_model.h | 4 +- multibody/fem/matrix_utilities.cc | 8 +- multibody/fem/matrix_utilities.h | 32 +----- multibody/fem/test/constitutive_model_test.cc | 2 +- .../test/constitutive_model_test_utilities.cc | 4 +- multibody/fem/test/matrix_utilities_test.cc | 6 +- multibody/fem/volumetric_element.h | 41 +------ 18 files changed, 271 insertions(+), 92 deletions(-) create mode 100644 math/fourth_order_tensor.cc create mode 100644 math/fourth_order_tensor.h create mode 100644 math/test/fourth_order_tensor_test.cc diff --git a/math/BUILD.bazel b/math/BUILD.bazel index bee2b3cbf8d1..46654f9a6697 100644 --- a/math/BUILD.bazel +++ b/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 = [ diff --git a/math/fourth_order_tensor.cc b/math/fourth_order_tensor.cc new file mode 100644 index 000000000000..6017344d5134 --- /dev/null +++ b/math/fourth_order_tensor.cc @@ -0,0 +1,31 @@ +#include "drake/math/fourth_order_tensor.h" + +#include "drake/common/default_scalars.h" + +namespace drake { +namespace math { + +template +FourthOrderTensor::FourthOrderTensor(const MatrixType& data) : data_(data) {} + +template +FourthOrderTensor::FourthOrderTensor() = default; + +template +void FourthOrderTensor::ContractWithVectors( + const Eigen::Ref>& u, + const Eigen::Ref>& v, EigenPtr> 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 math +} // namespace drake + +DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_NONSYMBOLIC_SCALARS( + class ::drake::math::FourthOrderTensor); +template class drake::math::FourthOrderTensor; diff --git a/math/fourth_order_tensor.h b/math/fourth_order_tensor.h new file mode 100644 index 000000000000..52026cfa4945 --- /dev/null +++ b/math/fourth_order_tensor.h @@ -0,0 +1,106 @@ +#pragma once + +#include "drake/common/drake_throw.h" +#include "drake/common/eigen_types.h" + +namespace drake { +namespace math { + +/** 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. + @experimental */ +template +class FourthOrderTensor { + public: + DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(FourthOrderTensor) + using MatrixType = Eigen::Matrix; + + /** 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>& u, + const Eigen::Ref>& v, + EigenPtr> B) const; + + /** Returns this fourth order tensor encoded as a matrix. */ + const MatrixType& data() const { return data_; } + + /** (Advanced) Returns this fourth order tensor encoded as a mutable matrix. + */ + 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 math +} // namespace drake diff --git a/math/test/fourth_order_tensor_test.cc b/math/test/fourth_order_tensor_test.cc new file mode 100644 index 000000000000..0b4bce02e47e --- /dev/null +++ b/math/test/fourth_order_tensor_test.cc @@ -0,0 +1,76 @@ +#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 { + +Eigen::Matrix MakeArbitraryMatrix() { + Eigen::Matrix 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 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 data = MakeArbitraryMatrix(); + FourthOrderTensor tensor(data); + EXPECT_EQ(tensor.data(), data); + tensor.set_data(data.transpose()); + EXPECT_EQ(tensor.data(), data.transpose()); + + EXPECT_EQ(tensor(0, 0, 0, 0), data(0, 0)); + EXPECT_EQ(tensor(1, 1, 1, 1), data(4, 4)); +} + +GTEST_TEST(FourthOrderTensorTest, ContractWithVectors) { + FourthOrderTensor 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 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(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 math +} // namespace drake diff --git a/multibody/fem/BUILD.bazel b/multibody/fem/BUILD.bazel index ea9a3ef84324..2c3b4608637f 100644 --- a/multibody/fem/BUILD.bazel +++ b/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", ], ) diff --git a/multibody/fem/constitutive_model.h b/multibody/fem/constitutive_model.h index 223f62f09285..ac762a551bad 100644 --- a/multibody/fem/constitutive_model.h +++ b/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 { @@ -93,7 +94,7 @@ class ConstitutiveModel { ------------------------------------- @pre `dPdF != nullptr`. */ void CalcFirstPiolaStressDerivative(const Data& data, - Eigen::Matrix* dPdF) const { + math::FourthOrderTensor* 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* dPdF) const { + void CalcFirstPiolaStressDerivativeImpl( + const Data& data, math::FourthOrderTensor* dPdF) const { throw std::logic_error( fmt::format("The derived class {} must provide a shadow definition of " "CalcFirstPiolaStressDerivativeImpl() to be correct.", diff --git a/multibody/fem/corotated_model.cc b/multibody/fem/corotated_model.cc index 3f102e82afe1..f813ec05aeee 100644 --- a/multibody/fem/corotated_model.cc +++ b/multibody/fem/corotated_model.cc @@ -44,7 +44,7 @@ void CorotatedModel::CalcFirstPiolaStressImpl(const Data& data, template void CorotatedModel::CalcFirstPiolaStressDerivativeImpl( - const Data& data, Eigen::Matrix* dPdF) const { + const Data& data, math::FourthOrderTensor* dPdF) const { const T& Jm1 = data.Jm1(); const Matrix3& F = data.deformation_gradient(); const Matrix3& R = data.R(); @@ -54,9 +54,10 @@ void CorotatedModel::CalcFirstPiolaStressDerivativeImpl( Eigen::Map>(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(R, S, -2.0 * mu_, &local_dPdF); /* The contribution from derivatives of JFinvT. */ diff --git a/multibody/fem/corotated_model.h b/multibody/fem/corotated_model.h index 9e3673fe3be7..398830905344 100644 --- a/multibody/fem/corotated_model.h +++ b/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* dPdF) const; + void CalcFirstPiolaStressDerivativeImpl( + const Data& data, math::FourthOrderTensor* dPdF) const; T E_; // Young's modulus, N/m². T nu_; // Poisson's ratio. diff --git a/multibody/fem/linear_constitutive_model.cc b/multibody/fem/linear_constitutive_model.cc index b0383816efe0..0cd04e8ecbcf 100644 --- a/multibody/fem/linear_constitutive_model.cc +++ b/multibody/fem/linear_constitutive_model.cc @@ -29,7 +29,7 @@ LinearConstitutiveModel::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::Identity(); + dPdF_.set_data(mu_ * Eigen::Matrix::Identity()); for (int k = 0; k < 3; ++k) { /* Second term. */ for (int l = 0; l < 3; ++l) { @@ -65,7 +65,7 @@ void LinearConstitutiveModel::CalcFirstPiolaStressImpl(const Data& data, template void LinearConstitutiveModel::CalcFirstPiolaStressDerivativeImpl( - const Data&, Eigen::Matrix* dPdF) const { + const Data&, math::FourthOrderTensor* dPdF) const { *dPdF = dPdF_; } diff --git a/multibody/fem/linear_constitutive_model.h b/multibody/fem/linear_constitutive_model.h index b986f85ddb89..005b326f4302 100644 --- a/multibody/fem/linear_constitutive_model.h +++ b/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* dPdF) const; + void CalcFirstPiolaStressDerivativeImpl( + const Data& data, math::FourthOrderTensor* 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 + math::FourthOrderTensor dPdF_; // The First Piola stress derivative is constant and precomputed. }; diff --git a/multibody/fem/linear_corotated_model.cc b/multibody/fem/linear_corotated_model.cc index a3fed81e495c..a67eed857a1e 100644 --- a/multibody/fem/linear_corotated_model.cc +++ b/multibody/fem/linear_corotated_model.cc @@ -82,11 +82,11 @@ void LinearCorotatedModel::CalcFirstPiolaStressImpl(const Data& data, */ template void LinearCorotatedModel::CalcFirstPiolaStressDerivativeImpl( - const Data& data, Eigen::Matrix* dPdF) const { + const Data& data, math::FourthOrderTensor* dPdF) const { const Matrix3& R0 = data.R0(); auto& local_dPdF = (*dPdF); /* Add in μ * δₐᵢδⱼᵦ. */ - local_dPdF = mu_ * Eigen::Matrix::Identity(); + local_dPdF.set_data(mu_ * Eigen::Matrix::Identity()); for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { for (int alpha = 0; alpha < 3; ++alpha) { diff --git a/multibody/fem/linear_corotated_model.h b/multibody/fem/linear_corotated_model.h index c561a73ffdbe..3c59bb90b28d 100644 --- a/multibody/fem/linear_corotated_model.h +++ b/multibody/fem/linear_corotated_model.h @@ -71,8 +71,8 @@ class LinearCorotatedModel final /* Shadows ConstitutiveModel::CalcFirstPiolaStressDerivativeImpl() as required by the CRTP base class. */ - void CalcFirstPiolaStressDerivativeImpl(const Data& data, - Eigen::Matrix* dPdF) const; + void CalcFirstPiolaStressDerivativeImpl( + const Data& data, math::FourthOrderTensor* dPdF) const; T E_; // Young's modulus, N/m². T nu_; // Poisson's ratio. diff --git a/multibody/fem/matrix_utilities.cc b/multibody/fem/matrix_utilities.cc index b0258785cc4e..7b56a4f734f4 100644 --- a/multibody/fem/matrix_utilities.cc +++ b/multibody/fem/matrix_utilities.cc @@ -50,9 +50,9 @@ void PolarDecompose(const Matrix3& F, EigenPtr> R, } template -void AddScaledRotationalDerivative( - const Matrix3& R, const Matrix3& S, const T& scale, - EigenPtr> scaled_dRdF) { +void AddScaledRotationalDerivative(const Matrix3& R, const Matrix3& S, + const T& scale, + math::FourthOrderTensor* scaled_dRdF) { /* Some notes on derivation on the derivative of the rotation matrix from polar decomposition: we start with the result from section 2 of [McAdams, 2011] about the differential of the rotation matrix, which states that δR = @@ -118,7 +118,7 @@ void CalcCofactorMatrix(const Matrix3& M, EigenPtr> cofactor) { template void AddScaledCofactorMatrixDerivative( const Matrix3& M, const T& scale, - EigenPtr> scaled_dCdM) { + math::FourthOrderTensor* scaled_dCdM) { /* See the convention for ordering the 9-by-9 derivatives at the top of the header file. */ const Matrix3 A = scale * M; diff --git a/multibody/fem/matrix_utilities.h b/multibody/fem/matrix_utilities.h index 56be4822a62c..f9afa5bdac4b 100644 --- a/multibody/fem/matrix_utilities.h +++ b/multibody/fem/matrix_utilities.h @@ -5,32 +5,13 @@ #include #include "drake/common/eigen_types.h" +#include "drake/math/fourth_order_tensor.h" namespace drake { namespace multibody { namespace fem { namespace internal { -/* Some of the following functions involve calculations about a 4th order tensor -(call it A) of dimension 3*3*3*3. We follow the following convention to flatten -the 4th order tensor into 9*9 matrices that are organized as follows: - - 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ᵢⱼₖₗ. */ - /* Calculates the condition number for the given matrix A. The condition number is calculated via k(A) = σₘₐₓ / σₘᵢₙ @@ -65,9 +46,9 @@ void PolarDecompose(const Matrix3& F, EigenPtr> R, @pre scaled_dRdF != nullptr. @tparam_nonsymbolic_scalar */ template -void AddScaledRotationalDerivative( - const Matrix3& R, const Matrix3& S, const T& scale, - EigenPtr> scaled_dRdF); +void AddScaledRotationalDerivative(const Matrix3& R, const Matrix3& S, + const T& scale, + math::FourthOrderTensor* scaled_dRdF); /* Calculates the cofactor matrix of the given input 3-by-3 matrix M. @pre cofactor != nullptr. @@ -83,9 +64,8 @@ void CalcCofactorMatrix(const Matrix3& M, EigenPtr> cofactor); @pre scaled_dCdM != nullptr. @tparam_nonsymbolic_scalar */ template -void AddScaledCofactorMatrixDerivative( - const Matrix3& M, const T& scale, - EigenPtr> scaled_dCdM); +void AddScaledCofactorMatrixDerivative(const Matrix3& M, const T& scale, + math::FourthOrderTensor* scaled_dCdM); /* Given a size 3N vector with block structure with size 3 block entries Bᵢ where i ∈ V = {0, ..., N-1} and a permutation P on V, this function builds the diff --git a/multibody/fem/test/constitutive_model_test.cc b/multibody/fem/test/constitutive_model_test.cc index 1e47d75a621c..637cde1fd190 100644 --- a/multibody/fem/test/constitutive_model_test.cc +++ b/multibody/fem/test/constitutive_model_test.cc @@ -48,7 +48,7 @@ GTEST_TEST(ConstitutiveModelTest, InvalidModel) { "CalcFirstPiolaStressImpl.. to be correct.", NiceTypeName::Get(model))); - Eigen::Matrix dPdF; + math::FourthOrderTensor dPdF; DRAKE_EXPECT_THROWS_MESSAGE( model.CalcFirstPiolaStressDerivative(data, &dPdF), fmt::format("The derived class {} must provide a shadow definition of " diff --git a/multibody/fem/test/constitutive_model_test_utilities.cc b/multibody/fem/test/constitutive_model_test_utilities.cc index df206effaf63..ad1462eb3587 100644 --- a/multibody/fem/test/constitutive_model_test_utilities.cc +++ b/multibody/fem/test/constitutive_model_test_utilities.cc @@ -160,14 +160,14 @@ void TestdPdFIsDerivativeOfP() { data.UpdateData(deformation_gradients, previous_step_deformation_gradients); Matrix3 P; model.CalcFirstPiolaStress(data, &P); - Eigen::Matrix dPdF; + math::FourthOrderTensor dPdF; model.CalcFirstPiolaStressDerivative(data, &dPdF); for (int i = 0; i < kSpaceDimension; ++i) { for (int j = 0; j < kSpaceDimension; ++j) { Matrix3d dPijdF; for (int k = 0; k < kSpaceDimension; ++k) { for (int l = 0; l < kSpaceDimension; ++l) { - dPijdF(k, l) = dPdF(3 * j + i, 3 * l + k).value(); + dPijdF(k, l) = dPdF(i, j, k, l).value(); } } EXPECT_TRUE(CompareMatrices( diff --git a/multibody/fem/test/matrix_utilities_test.cc b/multibody/fem/test/matrix_utilities_test.cc index ff49160a0ad3..8116df572ec0 100644 --- a/multibody/fem/test/matrix_utilities_test.cc +++ b/multibody/fem/test/matrix_utilities_test.cc @@ -68,8 +68,7 @@ GTEST_TEST(MatrixUtilitiesTest, AddScaledRotationalDerivative) { const Matrix3 F = MakeAutoDiffMatrix(3, 3); Matrix3 R, S; PolarDecompose(F, &R, &S); - Eigen::Matrix scaled_dRdF = - Eigen::Matrix::Zero(); + math::FourthOrderTensor scaled_dRdF; AutoDiffXd scale = 1.23; AddScaledRotationalDerivative(R, S, scale, &scaled_dRdF); for (int i = 0; i < 3; ++i) { @@ -100,8 +99,7 @@ GTEST_TEST(MatrixUtilitiesTest, AddScaledCofactorMatrixDerivative) { const Matrix3 A = MakeAutoDiffMatrix(3, 3); Matrix3 C; CalcCofactorMatrix(A, &C); - Eigen::Matrix scaled_dCdA = - Eigen::Matrix::Zero(); + math::FourthOrderTensor scaled_dCdA; AutoDiffXd scale = 1.23; AddScaledCofactorMatrixDerivative(A, scale, &scaled_dCdA); for (int i = 0; i < 3; ++i) { diff --git a/multibody/fem/volumetric_element.h b/multibody/fem/volumetric_element.h index c48563c1438f..e0784abba0d0 100644 --- a/multibody/fem/volumetric_element.h +++ b/multibody/fem/volumetric_element.h @@ -14,41 +14,6 @@ namespace multibody { namespace fem { namespace internal { -// TODO(xuchenhan-tri): Encapsulate the the memory layout of 4th order tensors -// and the contraction operation in a FourthOrderTensor class. -/* Helper function that performs a contraction between a 4th order tensor A - and two vectors u and v and returns a matrix B. In Einstein notation, the - contraction is: Bᵢₖ = uⱼ Aᵢⱼₖₗ vₗ. The 4th order tensor A of dimension - 3*3*3*3 is flattened to 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ᵢⱼₖₗ. */ -template -void PerformDoubleTensorContraction( - const Eigen::Ref>& A, - const Eigen::Ref>& u, - const Eigen::Ref>& v, EigenPtr> B) { - B->setZero(); - for (int l = 0; l < 3; ++l) { - for (int j = 0; j < 3; ++j) { - *B += A.template block<3, 3>(3 * j, 3 * l) * u(j) * v(l); - } - } -} - /* The data struct that stores per element data for VolumetricElement. See FemElement for the requirement. We define it here instead of nesting it in the traits class below due to #17109. */ @@ -74,7 +39,7 @@ struct VolumetricElementData { std::array, num_quadrature_points> P; /* The derivative of first Piola stress with respect to the deformation gradient evaluated at quadrature points. */ - std::array, num_quadrature_points> dPdF; + std::array, num_quadrature_points> dPdF; }; /* Forward declaration needed for defining the traits below. */ @@ -354,8 +319,8 @@ class VolumetricElement /* Note that the scale is negated here because the tensor contraction gives the second derivative of energy, which is the opposite of the force derivative. */ - PerformDoubleTensorContraction( - data.dPdF[q], dSdX_transpose_[q].col(a), + data.dPdF[q].ContractWithVectors( + dSdX_transpose_[q].col(a), dSdX_transpose_[q].col(b) * reference_volume_[q] * -scale, &K_ab); AccumulateMatrixBlock(K_ab, a, b, K); }