adtzlr / matadi

Material Definition with Automatic Differentiation
GNU General Public License v3.0
21 stars 2 forks source link

Add statevars materialtensor #87

Closed adtzlr closed 2 years ago

adtzlr commented 2 years ago

fixes #86

codecov-commenter commented 2 years ago

Codecov Report

Merging #87 (e3949e3) into main (9f70d97) will increase coverage by 0.26%. The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##             main      #87      +/-   ##
==========================================
+ Coverage   99.60%   99.86%   +0.26%     
==========================================
  Files          20       20              
  Lines         765      768       +3     
==========================================
+ Hits          762      767       +5     
+ Misses          3        1       -2     
Impacted Files Coverage Δ
matadi/_material.py 99.16% <100.00%> (+0.02%) :arrow_up:
matadi/math.py 100.00% <0.00%> (+8.69%) :arrow_up:

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 9f70d97...e3949e3. Read the comment docs.

adtzlr commented 2 years ago

Example code:

from matadi import MaterialTensor, Variable
from matadi.math import det, dev, inv
import numpy as np

defgrad = np.random.rand(3, 3, 5, 100) - 0.5
for a in range(3):
    defgrad[a, a] += 1.0

pressure = np.random.rand(5, 100)
statevars = np.random.rand(5, 16, 5, 100)

F = Variable("F", 3, 3)
p = Variable("p")
z = Variable("z", 5, 16)

def nh(x, C10=0.5, bulk=5000):
    """Neo-Hookean material formulation with first Piola-Kirchhoff stress
    and volumetric constraint equation (u/p formulation)."""

    F, p, z = x

    J = det(F)
    C = F.T @ F

    S = 2 * C10 * dev(J ** (-2 / 3) * C) @ inv(C) + p * J * inv(C)
    constraint = (J - 1) - p / bulk

    return F @ S, constraint, z, z[:3, :]

# (u/p) framework
# ===============
# 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)

NH = MaterialTensor(x=[F, p, z], fun=nh, statevars=1, triu=True)

W = NH.function([defgrad, pressure, statevars])
P = NH.gradient([defgrad, pressure, statevars])