Implement gradient-vector- and hessian-vector-product

[adtzlr]

Enable evaluations of hyperelastic forms without explicit gradients and hessians evaluated by AD; only jacobian - vector products (jtimes) are generated.

Material definition with Automatic Differentiation:

import casadi as ca

class NeoHooke:

    def __init__(self, mu, bulk):

        F  = ca.SX.sym("F", 3, 3)
        dF = ca.SX.sym("dF", 3, 3)
        DF = ca.SX.sym("DF", 3, 3)

        C = ca.transpose(F) @ F
        J = ca.det(F)

        W = mu * (J**(-2/3) * ca.trace(C) - 3) + bulk * (J - 1)**2/2

        dW  = ca.jtimes( W, F, dF)
        DdW = ca.jtimes(dW, F, DF)

        self._gvp = ca.Function("g", [F, dF], [dW])
        self._hvp = ca.Function("h", [F, dF, DF], [DdW])

    def ravel(self, T):
        return T.reshape(T.shape[0], -1, order="F")

    def gradient_vector_product(self, F, dF):
        n = F.shape[-2] * F.shape[-1]
        gvp = self._gvp.map(n, "thread", 4)
        return np.array(
                gvp(self.ravel(F), self.ravel(dF))
            ).reshape(F.shape[-2:], order="F")

    def hessian_vector_product(self, F, dF, DF):
        n = F.shape[-2] * F.shape[-1]
        hvp = self._hvp.map(n, "thread", 4)
        return np.array(
                hvp(self.ravel(F), self.ravel(dF), self.ravel(DF))
            ).reshape(F.shape[-2:], order="F")

And here is the defintion of (bi-) linear forms with the above material class in scikit-fem:

```python umat = NeoHooke(mu=1.0, bulk=2.0) def deformation_gradient(w): dudX = grad(w["displacement"]) F = dudX + identity(dudX) return F @LinearForm(nthreads=4) def L(v, w): F = deformation_gradient(w) return umat.gradient_vector_product(F, grad(v)) @BilinearForm(nthreads=4) def a(u, v, w): F = deformation_gradient(w) return umat.hessian_vector_product(F, grad(v), grad(u)) ```

Originally posted by @adtzlr in https://github.com/kinnala/scikit-fem/discussions/839#discussioncomment-1887176

[adtzlr]

With the merged PR the following code is possible:

from matadi import MaterialHyperelastic
from matadi.models import neo_hooke

mat = MaterialHyperelastic(neo_hooke, C10=0.5)

# init some random deformation gradients
import numpy as np

defgrad = np.random.rand(3, 3, 5, 100) - 0.5
dF = np.random.rand(3, 3, 5, 100) - 0.5

for a in range(3):
    defgrad[a, a] += 1.0

dWdF = mat.gradient([defgrad])
dW = mat.gradient_vector_product([defgrad], [dF])

dW_check = np.einsum("ij...,ij...->...", dWdF[0], dF)

assert np.allclose(dW_check, dW[0])