Add math::FourthOrderTensor #22190

xuchenhan-tri · 2024-11-15T02:06:10Z

Make the abstraction for 3x3x3x3 fourth-order tensors and move it from fem::internal namespace to math namespace to prepare it to be used for MPM.

This change is

xuchenhan-tri

+(release notes: none) +@amcastro-tri for feature review please.

Reviewable status: LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

amcastro-tri

I like this idea a lot @xuchenhan-tri.
Right now however it looks more like a convenience struct and not a first class citizen that allows you to reuse common and easy to get wrong tensor operations in both FEM and MPM.
I gave you below a few suggestions. I'm thinking they'll prove useful to write your MPM code without having to rethink how to transalte math to indexes all over again. PTAL

Reviewed 10 of 18 files at r1, all commit messages.
Reviewable status: 9 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

math/fourth_order_tensor.h line 12 at r1 (raw file):

/* 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

I don't think this detail should show in the class's public API.
Consider moving to codcument MatrixType.

math/fourth_order_tensor.h line 46 at r1 (raw file):

   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,

I love this function in that it abstract the internal layout of FourthOrderTensor.

math/fourth_order_tensor.h line 82 at r1 (raw file):

   the class documentation.
   @pre 0 <= i, j < 9. */
  const T& operator()(int i, int j) const {

If needed, I'd document this differently, I'd just say this accesses entry Aᵢ₀ₖ₀.
Maybe unicode is not good for this, I meant to write, this accesses Aᵢⱼₖₗ for the special case j=l=0.

Ditto below.

multibody/fem/linear_corotated_model.cc line 89 at r1 (raw file):

  auto& local_dPdF = (*dPdF);
  /* Add in μ * δₐᵢδⱼᵦ. */
  local_dPdF.set_data(mu_ * Eigen::Matrix<T, 9, 9>::Identity());

ditto, this could be placed in a factory?

multibody/fem/linear_constitutive_model.cc line 28 at r1 (raw file):

        ∂εᵢⱼ/∂Fₖₗ = 0.5 * δᵢₖ δⱼₗ  + 0.5 * δᵢₗ δₖⱼ.
    Plugging in, we get:
        ∂Pᵢⱼ/∂Fₖₗ = μ * (δᵢₖδⱼₗ + δᵢₗ δⱼₖ) +  λ * δₖₗ * δᵢⱼ.

maybe FourthOrderTesnsor::Identity()?

Code quote:

(δᵢₖδⱼₗ + δᵢₗ δⱼₖ)

multibody/fem/linear_constitutive_model.cc line 28 at r1 (raw file):

        ∂εᵢⱼ/∂Fₖₗ = 0.5 * δᵢₖ δⱼₗ  + 0.5 * δᵢₗ δₖⱼ.
    Plugging in, we get:
        ∂Pᵢⱼ/∂Fₖₗ = μ * (δᵢₖδⱼₗ + δᵢₗ δⱼₖ) +  λ * δₖₗ * δᵢⱼ.

Maybe a factory ForuthOrderTensor::ProductOfSecondOrderIdentities()? (or better name)

Code quote:

δₖₗ * δᵢⱼ

multibody/fem/linear_constitutive_model.cc line 33 at r1 (raw file):

  /* First term. */
  dPdF_.set_data(mu_ * Eigen::Matrix<T, 9, 9>::Identity());
  for (int k = 0; k < 3; ++k) {

I belive this code could live in a nice factory inside FourthOrderTensor so that you can reuse it.

multibody/fem/corotated_model.cc line 58 at r1 (raw file):

  /* The contribution from derivatives of Jm1. */
  local_dPdF.mutable_data().noalias() =
      lambda_ * flat_JFinvT * flat_JFinvT.transpose();

it'd be nice we could express this as a tensor product. JFinvT is represented as a Matrix3, you could have a function TensorProduct(const Matrix3& A, const Matrix3& B, FourthOrderTensor*) that abstract the layout of the 4th order tensor. It'd avoid having to wrestle with these maps very time you need this operation. I imagine you'll need it in MPM.

multibody/fem/corotated_model.cc line 60 at r1 (raw file):

      lambda_ * flat_JFinvT * flat_JFinvT.transpose();
  /* The contribution from derivatives of F. */
  local_dPdF.mutable_data().diagonal().array() += 2.0 * mu_;

ditto here, you could have a factory FourthOrderTensor::MakeDiagonal(const T&) to abstract layout.

xuchenhan-tri

Thanks for the review! Comments addressed, PTAL.

Reviewable status: 9 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @xuchenhan-tri)

math/fourth_order_tensor.h line 12 at r1 (raw file):