Add math::FourthOrderTensor

[xuchenhan-tri]

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):

Previously, amcastro-tri (Alejandro Castro) wrote…

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

Done.

---

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

Previously, amcastro-tri (Alejandro Castro) wrote…

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

Done.

---

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

Previously, amcastro-tri (Alejandro Castro) wrote…

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.

That's not what this function does.
I removed this sugar instead of trying to explain what it does.

---

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

Previously, amcastro-tri (Alejandro Castro) wrote…

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.

Done

---

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

Previously, amcastro-tri (Alejandro Castro) wrote…

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

This is not the diagonal of the tensor, but rather a diagonal along one of the 2d planes. There are a few combinations like T_{iijj}, T_{ijij}, T_{ijji}. The one I'm interested in is T_{ijij}. I named it major diagonal, borrowing the name from major symmetry.

---

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

Previously, amcastro-tri (Alejandro Castro) wrote…

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

Done

---

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

Previously, amcastro-tri (Alejandro Castro) wrote…

maybe FourthOrderTesnsor::Identity()?

Done

---

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

Previously, amcastro-tri (Alejandro Castro) wrote…

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

Done

---

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

Previously, amcastro-tri (Alejandro Castro) wrote…

ditto, this could be placed in a factory?

Done

[amcastro-tri]

This is great Xuchen, thanks! Just one question remaining on SetAsOuterProduct()

Reviewed 3 of 18 files at r1, 7 of 8 files at r2.
Reviewable status: 6 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)

---

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

Previously, xuchenhan-tri wrote…

This is not the diagonal of the tensor, but rather a diagonal along one of the 2d planes. There are a few combinations like T_{iijj}, T_{ijij}, T_{ijji}. The one I'm interested in is T_{ijij}. I named it major diagonal, borrowing the name from major symmetry.

Sounds good!. I figured a better name was in order. I like that.

---

multibody/fem/linear_constitutive_model.cc line 32 at r2 (raw file):

    the jl-th block corresponds to the value dPᵢⱼ/dFₖₗ.  */
  /* First term. */
  dPdF_.SetAsOuterProduct(Matrix3<T>::Identity(),

in retrospect, this seems like an expensive way to compute δₖₗ * δᵢⱼ

---

math/fourth_order_tensor.cc line 33 at r2 (raw file):

  const auto M_vec = Eigen::Map<const Vector<T, 9>>(M.data(), 9);
  const auto N_vec = Eigen::Map<const Vector<T, 9>>(N.data(), 9);
  data_.noalias() = M_vec * N_vec.transpose();

double check this for me. I believe this is computing the outer product of M with the transpose of N?
Per the layout documented in the header.

That is, I believe this computes Aᵢⱼₖₗ = MᵢⱼNₗₖ instead of Aᵢⱼₖₗ = MᵢⱼNₖₗ

---

math/fourth_order_tensor.cc line 42 at r2 (raw file):

  result.data_ = half_scale * Eigen::Matrix<T, 9, 9>::Identity();
  for (int k = 0; k < 3; ++k) {
    /* Second term. */

You mean "second" in this form (δᵢₖδⱼₗ + δᵢₗδⱼₖ), but I could easily swap those.

I'd instead say: "term δᵢₗδⱼₖ."

Similarly for the "first" term.

Code quote:

 Second term.

---

math/fourth_order_tensor.cc line 43 at r2 (raw file):

  for (int k = 0; k < 3; ++k) {
    /* Second term. */
    for (int l = 0; l < 3; ++l) {

I'm curious, any reason to choose to loop over k-l insted of i-j? wondering if there is anything worth documenting here on this particular choice?

---

math/fourth_order_tensor.cc line 46 at r2 (raw file):

      const int i = l;
      const int j = k;
      result.data_(3 * j + i, 3 * l + k) += half_scale;

nit, we can avoid the addition

Suggestion:

 result.data_(3 * j + i, 3 * l + k) = half_scale;

---

math/fourth_order_tensor.cc line 54 at r2 (raw file):

template <typename T>
FourthOrderTensor<T> FourthOrderTensor<T>::MakeMajorIdentity(T scale) {
  return FourthOrderTensor<T>(scale * Eigen::Matrix<T, 9, 9>::Identity());

I believe here we'll make a temporary and pass it down to the ctor.
Consider a ctor with move semantics instead (i.e. FourthOrderTensor(MatrixType data) : data_(std::move(data)) {})

[xuchenhan-tri]

Reviewable status: 6 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.cc line 33 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

double check this for me. I believe this is computing the outer product of M with the transpose of N?
Per the layout documented in the header.

That is, I believe this computes Aᵢⱼₖₗ = MᵢⱼNₗₖ instead of Aᵢⱼₖₗ = MᵢⱼNₖₗ

The matrix lay out is Aᵢⱼₖₗ = Aₐᵦ with a = 3j + i, b = 3l + k.

We want Aₐᵦ = Aᵢⱼₖₗ = Mᵢⱼ*Nₖₗ = Mₐ * Nᵦ which is M_vec * N_vec.transpose(), where Foo_vec is the column concatenation of Foo.

This is also confirmed by the unit test.

---

math/fourth_order_tensor.cc line 42 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

You mean "second" in this form (δᵢₖδⱼₗ + δᵢₗδⱼₖ), but I could easily swap those.

I'd instead say: "term δᵢₗδⱼₖ."

Similarly for the "first" term.

Done.

---

math/fourth_order_tensor.cc line 43 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

I'm curious, any reason to choose to loop over k-l insted of i-j? wondering if there is anything worth documenting here on this particular choice?

It's better to loop over j,i because the matrix is stored column major. Done.

---

math/fourth_order_tensor.cc line 46 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

nit, we can avoid the addition

We need the addition; otherwise the entries on the diagonal are wrong.
This is also confirmed by the unit test.

---

math/fourth_order_tensor.cc line 54 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

I believe here we'll make a temporary and pass it down to the ctor.
Consider a ctor with move semantics instead (i.e. FourthOrderTensor(MatrixType data) : data_(std::move(data)) {})

I don't think the move semantic would do much here. With a fixed size matrix, the data will be copied regardless whether it's moved or copied.

---

multibody/fem/linear_constitutive_model.cc line 32 at r2 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

in retrospect, this seems like an expensive way to compute δₖₗ * δᵢⱼ

Yeah this is indeed more expensive than directly modifying the matrix. The alternatives are directly modifying the underlying matrix or coming up with another name for this identity.

On the bright, we only create constitutive model at registration time, and this is not that much more expensive.

[amcastro-tri]

Excellent!

Reviewed 1 of 18 files at r1, 1 of 8 files at r2, 1 of 1 files at r3, all commit messages.
Reviewable status: 5 unresolved discussions, 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.cc line 33 at r2 (raw file):

Previously, xuchenhan-tri wrote…

The matrix lay out is Aᵢⱼₖₗ = Aₐᵦ with a = 3j + i, b = 3l + k.

We want Aₐᵦ = Aᵢⱼₖₗ = Mᵢⱼ*Nₖₗ = Mₐ * Nᵦ which is M_vec * N_vec.transpose(), where Foo_vec is the column concatenation of Foo.

This is also confirmed by the unit test.

Ah, I just realized I flipped the indexes in my notes. Thanks for double checking.

---

math/fourth_order_tensor.cc line 43 at r2 (raw file):

Previously, xuchenhan-tri wrote…

It's better to loop over j,i because the matrix is stored column major. Done.

ah, very smart! could you add a comment stating that choice here?

---

math/fourth_order_tensor.cc line 46 at r2 (raw file):

Previously, xuchenhan-tri wrote…

We need the addition; otherwise the entries on the diagonal are wrong.
This is also confirmed by the unit test.

ah, of course, this is the symmetric one. Thanks

---

math/fourth_order_tensor.cc line 54 at r2 (raw file):

Previously, xuchenhan-tri wrote…

I don't think the move semantic would do much here. With a fixed size matrix, the data will be copied regardless whether it's moved or copied.

is that rigth?. I would've thougth the expression returned by scale*Matrix<T,9,9>::Identity() would end up in a temporary, passed by reference to the ctor, and copied into data_.

Vs with moving the expression would instantiate a temporary, which then would get immediately moved into data_, thus avoiding the copy of 81 numbers.

Marking resolved though if you think not worth doing.

---

multibody/fem/linear_constitutive_model.cc line 32 at r2 (raw file):

Previously, xuchenhan-tri wrote…

Yeah this is indeed more expensive than directly modifying the matrix. The alternatives are directly modifying the underlying matrix or coming up with another name for this identity.

On the bright, we only create constitutive model at registration time, and this is not that much more expensive.

oh, that's right, this is only linear. Probably not worth the effort. Thanks.

---

math/test/fourth_order_tensor_test.cc line 31 at r3 (raw file):

  Matrix3d B;

  tensor.ContractWithVectors(u, v, &B);

minor, this last constration test is probably not needed since you already checked all data in the tensor is indeed zero.

---

math/test/fourth_order_tensor_test.cc line 43 at r3 (raw file):

  tensor.mutable_data() = data.transpose();
  EXPECT_EQ(tensor.data(), data.transpose());
  /* Settor */

minor

Suggestion:

Setter

---

math/test/fourth_order_tensor_test.cc line 48 at r3 (raw file):

  /* Operator with four indices. */
  EXPECT_EQ(tensor(0, 0, 0, 0), data(0, 0));
  EXPECT_EQ(tensor(1, 1, 1, 1), data(4, 4));

minor, you might want to test a completely arbitrary set of indexes that is not that symmetric.

---

math/test/fourth_order_tensor_test.cc line 128 at r3 (raw file):

        for (int l = 0; l < 3; ++l) {
          double expected_entry = 0.0;
          if (i == k && j == l) expected_entry += scale;

minor, just to be crystal clear, here the += can be replace by = correct?

Code quote:

expected_entry += scale;