M = r * r.T

[1,0] instead of [0,1].

There are two versions of the Ogden material formulation, one with mu / alpha and one with 2 mu / alpha^2. Change to the second version due automatic consistency with linear elasticity (mu is always equal to the initial shear modulus)

Implemented material models should only contain the essential code. A multiplicative split of the deformation gradient should be applied via a decorator.

• add decorator
• simplify material models
• check tests

because the gradient fails in rare cases. Improve it by adding a full perturbation matrix (and not only two diagonal entries).

the areastretch should be areastretch = sum1(sqrt(diag(r.T @ Cs @ r)) ** q * w) (without q-root)

replace inv(C) with cof(C) = det(C) * inv(C) because da = cof(F) * dA. In an isochoric-volumetric splitted setting this doesn't matter as det(C) = 1 but this should be changed.

A given fun should be augmented by the (u/p)-framework.

Although a bit of an overkill; implement a linear-elastic strain energy function which can be used in MaterialHyperelastic. It takes the symmetric part of the displacement gradient as strain variable, on which a quadratic potential is described.

Code:

from matadi import MaterialHyperelastic
from matadi.math import eye, sym, trace

def linear_elastic(F, mu, lmbda):
    strain = sym(F - eye(3))
    return mu * trace(strain @ strain) + lmbda / 2 * trace(strain) ** 2

mat = MaterialHyperelastic(fun=linear_elastic, mu=1.0, lmbda=2.0)

TODO:

• add eye and sym to math

It is based on strains, not deformation tensors

In the meantime, use MaterialTensor() instead.

this enables future support for microsphere models with state variables

as micro-sphere model template

Lab

Current behavior

Currently the material is analyzed as is. That means the volume change is introduced through the material. This is also called the nearly-incompressible framework.

New Feature

Add a new Lab feature to run incompressible loadcases. Here the incompressibility is enforced by a constraint eqation (stretch_1 * stretch_2 * stretch_3 = 1). It doesn't matter if the material contains a volumetric part or not because is is deactivated by the constraint.

FElupe>5.3.1 calls materials with

# allow (unused) state variables to be passed as the c-
W = umat.function([deformation_gradient, statevars])
P, statevars_new = umat.gradient([deformation_gradient, statevars]) # also returns the updated state variables
A = umat.hessian([deformation_gradient, statevars])

This is just a small change for matADi with a wrapper for the function, gradient and hessian methods for hyperelastic materials.

class MaterialHyperelastic:
    def __init__(self, fun, **kwargs):
        F = Variable("F", 3, 3)
        self.x = [F]
        self.fun = fun
        self.kwargs = kwargs
        self.W = Material(self.x, self._fun_wrapper, kwargs=self.kwargs)
        self.gradient_vector_product = self.W.gradient_vector_product
        self.hessian_vector_product = self.W.hessian_vector_product

    def _fun_wrapper(self, x, **kwargs):
        return self.fun(x[0], **kwargs)

    def function(self, x, *args, **kwargs):
        return self.W.function(x[:1], *args, **kwargs)

    def gradient(self, x, *args, **kwargs):
        return [*self.W.gradient(x[:1], *args, **kwargs), None]

    def hessian(self, x, *args, **kwargs):
        return self.W.hessian(x[:1], *args, **kwargs)

implementation is unstable

introduce a new class ThreeFieldVariationPlaneStrain for plane strain hyperelastic materials

Hyperelastic models may be unstable even if the force-stretch relation seems ok (see .

Stability

For each state of homogenous deformation (uniaxial, biaxial and planar loading)
a) the elasticity matrix has to be evaluated
b) converted to shape of (6,6) (Voigt-notation) in order to invert the elasticity matrix
c) and the sign of the resulting stretch in direction 1 as linear solution of a unit load in direction 1 is checked (1 = stable and 0 = unstable)

Additional checks:

• volume ratio J = det(F)
• slope of force per undeformed area vs. stretch

This is really important!

From a user side, the unstable regions should simply be plotted with dashed lines (vs. solid lines in stable regions).

Create a composite material out of a list of Material with same input arguments and results will be summed up.

Matadi's dot-product is taken over from casadi where the first argument is transposed. However in matadi the dot-product should be

A = dot(B, C) should be A(i,j) = B(i,k) C(k,j)

and not

A = dot(B, C) should be A(i,j) = B(k, i) C(k,j).

as a follow-up of #84

because MaterialTensorGeneral also uses a fun argument.

A Function should serve as a base class for any Material but without calculating gradient and hessian.

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

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 kinnala/scikit-fem#839 (comment)

• allow a list as output of a user defined function
• add triu argument for gradient evaluation (like hessian-method for mixed-field hyperelastic materials)
  allow state variables (same as default variables but for which no gradient is evaluated) --> as seperate issue/PR

With the above enhancements it will be possible to code a u/p - Formulation:

# W(u, p) = w(u) + p (J - 1) - p^2 / (2 K)
# dWdE = dwdE + p J inv(C) (= S)
# dWdp = (J - 1) - p / K (= constraint)

If the hyperelastic materials would be based on C = F.T @ F and C would be stored as a vector in Voigt-notation, this gives a speed-up of 2 in stress evaluation and of 6 in elasticity evaluation! However, the code gets much more complicated.

change it to keyword arguments

f1 = fiber(F, E, angle=angle, k=k, axis=axis, compression=compression)

this is used for the non-affine tube part

and be consistent with MaterialTensor. Also, make compress consistent for function output with gradient and hessian outputs.

• add triu argument for hessian output
• add statevars argument
• improve compress
• update tests

I2 is not trace(C @ C), it's (trace(C)**2 - trace(C @ C)) / 2 (the second invariant from the characteristic polynomial) for invariant based material formulations.

• Mooney-Rivlin
• Third-Order-Deformation

...so they can be used with SolidBodyNearlyIncompressible of FElupe.

The driving force of the 1d micro-sphere model has to be integrated w.r.t. the stretch, in order to obtain the 1d strain energy function. The Padé approximation of the inverse Langevin function will be integrated w.r.t. the stretch. Here is the result of WolframAlpha

integrate (3*N - λ^2)/(N-λ^2) dλ

which gives:

λ + 2 sqrt(N) tanh^(-1)(λ/sqrt(N))

This is transformed into a Python function:

def langevin(stretch, mu, N):
    """Langevin model (Padé approximation) given by the free energy 
    of a single chain as a function of the stretch."""

    return mu * (stretch + 2 * sqrt(N) * atanh(stretch / sqrt(N)))

Originally posted by @adtzlr in #28 (comment)

These two classes should have the same syntax as the 3d hyperelastic material class. All material formulations should be compatible.

The constraint equation needs to be divided by d2W/dJdJ in order to be consistent with the derivation of a potential.

see https://numpy.org/doc/stable/reference/generated/numpy.zeros_like.html

with support for

• state variables
• current and initial deformation gradient

Introduce a tensor valued material definition by a stress-like function in terms of a tensor valued kinematic variable. This is more or less a copy and paste thing; replace gradient -> jacobian and remove hessian because these only operate on scalar valued functions.

Implement the affine, non-affine, stretch and area-stretch micro-sphere models for rubber elasticity [1].

Tasks

• Add numeric integration scheme on the surface of a sphere [2]
• Add affine stretch model
• Add affine tube model
• Add non-affine stretch model
• Add non-affine tube model
• Add one-dimensional micro-sphere model (langevin and gauss model)
• Add links to README
• Add tests

References

[1] C. Miehe, “A micro-macro approach to rubber-like materials. Part I: the non-affine micro-sphere model of rubber elasticity,” Journal of the Mechanics and Physics of Solids, vol. 52, no. 11. Elsevier BV, pp. 2617–2660, Nov. 2004.

[2] P. Bažant and B. H. Oh, “Efficient Numerical Integration on the Surface of a Sphere,” ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, vol. 66, no. 1. Wiley, pp. 37–49, 1986.

as a general material based on the first Piola Kirchhoff stress.

like Ogden & Roxburgh (doi: 10.1098/rspa.1999.0431)

as described in

A. V. Shutov, R. Landgraf, and J. Ihlemann, An explicit solution for implicit time stepping in multiplicative finite strain viscoelasticity, Computer Methods in Applied Mechanics and Engineering, vol. 265. Elsevier BV, pp. 213–225, Oct. 2013. doi: 10.1016/j.cma.2013.07.004.

see also the version of the paper on arXiv: https://arxiv.org/abs/1304.3380