Repository: Genesis-Embodied-AI/Genesis
File: genesis/engine/materials/FEM/elastic.py (Python)

Code before:
            The deformation gradient matrix.
        actu: ti.Matrix
            The activation matrix (not used in linear elasticity).
        m_dir: ti.Matrix
            The material direction (not used in linear elasticity).
        hessian_field: ti.Matrix
            The Hessian of the energy with respect to the deformation gradient F.

        Returns
        -------
        energy: float
            The computed energy.
        gradient: ti.Matrix
            The gradient of the energy with respect to the deformation gradient F.

        Notes
        -------
        This implementation assumes small deformations and linear stress-strain relationship.
        It is adapted from the HOBAKv1 implementation for linear elasticity:
        https://github.com/theodorekim/HOBAKv1/blob/main/src/Hyperelastic/Volume/LINEAR.cpp

        """
        I = ti.Matrix.identity(dt=gs.ti_float, n=3)
        R = self.R[i_b, i_e]
        F_hat = R.transpose() @ F
        eps = 0.5 * (F_hat + F_hat.transpose()) - I
        trEps = eps.trace()
        energy = mu * eps.norm_sqr() + 0.5 * lam * trEps**2

        gradient = 2.0 * mu * R @ eps + lam * trEps * R

        # Zero out the matrix
        for i in ti.static(ti.grouped(ti.ndrange(3, 3))):
            hessian_field[i_b, i, i_e].fill(0.0)

        # Identity part
        for i, k in ti.static(ti.ndrange(3, 3)):
            hessian_field[i_b, i, i, i_e][k, k] = mu

        for i, j, alpha, beta in ti.static(ti.ndrange(3, 3, 3, 3)):
            hessian_field[i_b, j, beta, i_e][i, alpha] += mu * R[i, beta] * R[alpha, j] + lam * R[alpha, beta] * R[i, j]

        return energy, gradient

    @ti.func
    def compute_energy_gradient_linear_corotated(self, mu, lam, J, F, actu, m_dir, i_e, i_b):
        """
        Compute the energy, gradient for linear elasticity.

        Parameters
        ----------

Review comment (suggestion):
You should rather just subtract 1.0 from the diagonal terms:
```python
for i in ti.static(range(3)):
    eps[i] = eps[I] - 1.0
```

Code after:
            The deformation gradient matrix.
        actu: ti.Matrix
            The activation matrix (not used in linear elasticity).
        m_dir: ti.Matrix
            The material direction (not used in linear elasticity).
        hessian_field: ti.Matrix
            The Hessian of the energy with respect to the deformation gradient F.

        Returns
        -------
        energy: float
            The computed energy.
        gradient: ti.Matrix
            The gradient of the energy with respect to the deformation gradient F.

        Notes
        -------
        This implementation assumes small deformations and linear stress-strain relationship.
        It is adapted from the HOBAKv1 implementation for linear elasticity:
        https://github.com/theodorekim/HOBAKv1/blob/main/src/Hyperelastic/Volume/LINEAR.cpp

        """
        R = self.R[i_b, i_e]
        F_hat = R.transpose() @ F
        # E = 1/2(F_hat + F_hat.transpose()) - I
        eps = 0.5 * (F_hat + F_hat.transpose())
        for i in ti.static(range(3)):
            eps[i, i] -= 1.0
        trEps = eps.trace()
        energy = mu * eps.norm_sqr() + 0.5 * lam * trEps**2

        gradient = 2.0 * mu * R @ eps + lam * trEps * R

        # Zero out the matrix
        for i in ti.static(ti.grouped(ti.ndrange(3, 3))):
            hessian_field[i_b, i, i_e].fill(0.0)

        # Identity part
        for i, k in ti.static(ti.ndrange(3, 3)):
            hessian_field[i_b, i, i, i_e][k, k] = mu

        for i, j, alpha, beta in ti.static(ti.ndrange(3, 3, 3, 3)):
            hessian_field[i_b, j, beta, i_e][i, alpha] += mu * R[i, beta] * R[alpha, j] + lam * R[alpha, beta] * R[i, j]

        return energy, gradient

    @ti.func
    def compute_energy_gradient_linear_corotated(self, mu, lam, J, F, actu, m_dir, i_e, i_b):
        """
        Compute the energy, gradient for linear elasticity.