adtzlr / matadi Goto Github PK
View Code? Open in Web Editor NEWMaterial Definition with Automatic Differentiation
License: GNU General Public License v3.0
Material Definition with Automatic Differentiation
License: GNU General Public License v3.0
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)
and be consistent with MaterialTensor
. Also, make compress
consistent for function output with gradient and hessian outputs.
triu
argument for hessian outputstatevars
argumentcompress
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.
In the meantime, use MaterialTensor()
instead.
Hyperelastic models may be unstable even if the force-stretch relation seems ok (see .
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:
J = det(F)
This is really important!
From a user side, the unstable regions should simply be plotted with dashed lines (vs. solid lines in stable regions).
Implement the affine, non-affine, stretch and area-stretch micro-sphere models for rubber elasticity [1].
[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.
the areastretch
should be areastretch = sum1(sqrt(diag(r.T @ Cs @ r)) ** q * w)
(without q-root)
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)
because MaterialTensorGeneral
also uses a fun
argument.
A Function
should serve as a base class for any Material
but without calculating gradient and hessian.
It is based on strains, not deformation tensors
Create a composite material out of a list of Material
with same input arguments and results will be summed up.
A given fun
should be augmented by the (u/p)-framework.
Implemented material models should only contain the essential code. A multiplicative split of the deformation gradient should be applied via a decorator.
with support for
like Ogden & Roxburgh (doi: 10.1098/rspa.1999.0431)
implementation is unstable
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)
.
...so they can be used with SolidBodyNearlyIncompressible
of FElupe.
as a follow-up of #84
The constraint equation needs to be divided by d2W/dJdJ
in order to be consistent with the derivation of a potential.
this is used for the non-affine tube part
introduce a new class ThreeFieldVariationPlaneStrain
for plane strain hyperelastic materials
[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)
These two classes should have the same syntax as the 3d hyperelastic material class. All material formulations should be compatible.
because the gradient fails in rare cases. Improve it by adding a full perturbation matrix (and not only two diagonal entries).
this enables future support for microsphere models with state variables
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.
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.
M = r * r.T
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.
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)
as a general material based on the first Piola Kirchhoff stress.
triu
argument for gradient evaluation (like hessian-method for mixed-field hyperelastic materials)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)
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)
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
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.
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:
eye
and sym
to math
as micro-sphere model template
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.