[Tutorial] Easy and fast assembly of bilinear forms for hyperelastic solids #839
Replies: 4 comments 12 replies
-
|
Thanks, looks great: |
Beta Was this translation helpful? Give feedback.
-
|
Would be great but please have a look at the scikit-fem specific lines of code if something could be simplified. |
Beta Was this translation helpful? Give feedback.
-
|
Alright, I modified the boundary condition and the mesh slightly so the result looks like this: |
Beta Was this translation helpful? Give feedback.
-
|
Code is here, I'll continue later: import numpy as np
import skfem as fem
from skfem.assembly import LinearForm, BilinearForm
from skfem.helpers import grad, identity, ddot, det, transpose, inv, trace, mul
mesh = fem.MeshHex.init_tensor(
np.linspace(0, 1, 40),
np.linspace(-.1, .1, 8),
np.linspace(-.1, .1, 8),
)
element = fem.ElementVectorH1(fem.ElementHex1())
basis = fem.Basis(mesh, element, intorder=1)
mu, bulk = 1.0, 2.0
def deformation_gradient(w):
dudX = grad(w["displacement"])
F = dudX + identity(dudX)
return F, inv(F)
@LinearForm(nthreads=2)
def L(v, w):
F, iF = deformation_gradient(w)
dF = grad(v)
lnJ = np.log(det(F))
return mu * ddot(F, dF) + (bulk * lnJ - mu) * ddot(transpose(iF), dF)
@BilinearForm(nthreads=2)
def a(u, v, w):
F, iF = deformation_gradient(w)
DF = grad(u)
dF = grad(v)
dFiF = mul(dF, iF)
DFiF = mul(DF, iF)
tr_DFiF_dFiF = np.einsum("ji...,ij...->...", dFiF, DFiF)
lnJ = np.log(det(F))
return mu * ddot(DF, dF) - (bulk * lnJ - mu) * tr_DFiF_dFiF \
+ bulk * trace(dFiF) * trace(DFiF)
u = basis.zeros()
du = basis.zeros()
uv = basis.interpolate(u)
dofs = {
"left": basis.get_dofs(lambda x: x[0] == 0.0),
"right": basis.get_dofs(lambda x: x[0] == 1.0),
}
I = basis.complement_dofs(dofs)
f = fem.asm(L, basis, displacement=uv)
K = fem.asm(a, basis, basis, displacement=uv)
right = 0.0
tol = 1e-10
for iteration in range(28):
c = (min(iteration, 20) + 1.) / 21.
du_D = u.copy()
x, y, z = basis.doflocs[:, dofs["right"].nodal['u^1']]
du_D[dofs["right"].nodal['u^1']] = right - du_D[dofs["right"].nodal['u^1']]
du_D[dofs["right"].nodal['u^2']] = (
y * np.cos(c * np.pi) - z * np.sin(c * np.pi) - y
- du_D[dofs["right"].nodal['u^2']]
)
du_D[dofs["right"].nodal['u^3']] = (
y * np.sin(c * np.pi) + z * np.cos(c * np.pi) - z
- du_D[dofs["right"].nodal['u^3']]
)
system = fem.condense(K, -f, x=du_D, I=I)
du = fem.solve(*system)
norm_du = np.linalg.norm(du)
u += du
uv = basis.interpolate(u)
f = L.assemble(basis, displacement=uv)
K = a.assemble(basis, displacement=uv)
print(1 + iteration, norm_du)
if norm_du < tol:
break |
Beta Was this translation helpful? Give feedback.
-
|
Looks good! A similar but a little bit more complicated example can be found in "J. Bonet and R. D. Wood, Nonlinear continuum mechanics for Finite Element Analysis. Cambridge: Cambridge University Press, 2009. " (section 10.12.5 twisting column) |
Beta Was this translation helpful? Give feedback.
-
|
I'll have to update the docstring and the code because for a compressible Neo-Hookean material the volumetric term is associated to lambda and not the bulk modulus. |
Beta Was this translation helpful? Give feedback.
-
|
Here is what I had so far: def deformation_gradient(w):
dudX = grad(w["displacement"])
F = dudX + identity(dudX)
return F, inv(F)
@LinearForm
def L(v, w):
F, iF = deformation_gradient(w)
dF = grad(v)
lnJ = np.log(det(F))
return mu * ddot(F, dF) + (bulk * lnJ - mu) * ddot(transpose(iF), dF)
@BilinearForm
def a(u, v, w):
F, iF = deformation_gradient(w)
DF = grad(u)
dF = grad(v)
dFiF = mul(dF, iF)
DFiF = mul(DF, iF)
tr_DFiF_dFiF = ddot(transpose(dFiF), DFiF)
lnJ = np.log(det(F))
return mu * ddot(DF, dF) - (bulk * lnJ - mu) * tr_DFiF_dFiF \
+ bulk * trace(dFiF) * trace(DFiF)Used the helpers instead of einsum. |
Beta Was this translation helpful? Give feedback.
-
|
Please replace def deformation_gradient(w):
dudX = grad(w["displacement"])
F = dudX + identity(dudX)
return F, inv(F)
@LinearForm
def L(v, w):
F, iF = deformation_gradient(w)
dF = grad(v)
lnJ = np.log(det(F))
return mu * ddot(F, dF) + (lmbda * lnJ - mu) * ddot(transpose(iF), dF)
@BilinearForm
def a(u, v, w):
F, iF = deformation_gradient(w)
DF = grad(u)
dF = grad(v)
dFiF = mul(dF, iF)
DFiF = mul(DF, iF)
tr_DFiF_dFiF = ddot(transpose(dFiF), DFiF)
lnJ = np.log(det(F))
return mu * ddot(DF, dF) - (lmbda * lnJ - mu) * tr_DFiF_dFiF \
+ lmbda * trace(dFiF) * trace(DFiF) |
Beta Was this translation helpful? Give feedback.
-
|
And here's the updated docstring if you'd like to use it. r"""
The strain energy density function per unit undeformed volume of the
isotropic hyperelastic Neo-Hookean material formulation is given as
.. math::
\psi(\boldsymbol{F}) = \frac{\mu}{2} ( I_1 - 3) - \mu \ln(J) +
\frac{\lambda}{2} \ln^2(J)
by the first invariant of the right Cauchy-Green deformation tensor
and the determinant of the deformation gradient (volume change).
.. math::
I_1 &= \text{tr}(\boldsymbol{F}^T \boldsymbol{F}) =
\boldsymbol{F} : \boldsymbol{F}
J &= \det(\boldsymbol{F})
The variation of the strain energy function leads to
.. math::
\delta\psi(\boldsymbol{F}) = \frac{\mu}{2} \delta I_1 - \mu \delta\ln(J)
+ \lambda \ln(J) \delta\ln(J)
with
.. math::
\delta I_1 &= 2 \boldsymbol{F} : \delta\boldsymbol{F}
\delta \ln(J) &= \boldsymbol{F}^{-T} : \delta\boldsymbol{F}
which finally results in
.. math::
\delta\psi(\boldsymbol{F}) = \mu \boldsymbol{F} : \delta\boldsymbol{F} -
(\lambda \ln(J) - \mu) \boldsymbol{F}^{-T} : \delta\boldsymbol{F}
The linearization of the variation leads to
.. math::
\Delta\delta\psi(\boldsymbol{F}) = \mu \Delta\boldsymbol{F} :
\delta\boldsymbol{F} + (\lambda \ln(J) - \mu)
\Delta\boldsymbol{F}^{-T} : \delta\boldsymbol{F}
+ \lambda \left( \boldsymbol{F}^{-T} : \delta\boldsymbol{F} \right)
\left( \boldsymbol{F}^{-T} : \Delta\boldsymbol{F} \right)
with
.. math::
\Delta\boldsymbol{F}^{-T} = -\boldsymbol{F}^{-T} \Delta\boldsymbol{F}^{T}
\boldsymbol{F}^{-T}
which may be alternatively formulated by expressions of the trace:
.. math::
\Delta\delta\psi(\boldsymbol{F}) = \mu \Delta\boldsymbol{F} :
\delta\boldsymbol{F} + (\lambda \ln(J) - \mu)
\text{tr}(\Delta\boldsymbol{F} \boldsymbol{F}^{-1}
\delta\boldsymbol{F} \boldsymbol{F}^{-1})
+ \lambda \text{tr}(\delta\boldsymbol{F} \boldsymbol{F}^{-1})
\text{tr}(\Delta\boldsymbol{F} \boldsymbol{F}^{-1})
""" |
Beta Was this translation helpful? Give feedback.
-
|
Do you really think it is a good idea to update the boundaries inside the newton iterations? In my opinion we should have two loops, one for load/prescribed displacement incrementation and another one for equilibrium (these are the newton iterations). Doing it this way the solution of every increment is valid, i.e. it is in equilibrium. If we are only interested in the max. deformed state, the final solution seems fine though. |
Beta Was this translation helpful? Give feedback.
-
|
Not good idea, no. It's just a quick hack for testing because you cannot directly go to c=1. I plan to make an animation of the twisting so I need to do two loops eventually. |
Beta Was this translation helpful? Give feedback.
-
|
This technique also enables faster evaluations of hyperelastic forms using automatic differentiation (just as an idea using casADi...). No explicit gradients and hessians are evaluated by AD; only jacobian - vector products ( Material definition with Automatic Differentiationimport 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: 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)) |
Beta Was this translation helpful? Give feedback.
-
|
I'll implement this in matADi, see progress in adtzlr/matadi#76 |
Beta Was this translation helpful? Give feedback.
-
|
Hey @kinnala: Recently, I experimented a bit with forward-mode automatic differentiation / the dual number approach, and it ended up in tensortrax. I know, Autograd/JAX is the way to go; but I decided to write everything in "calculus of variation". It always bothered me that I didn't fully understand the basic concept. Finally, I just tried
import tensortrax as tr
import tensortrax.math as tm
def neo_hooke(F):
I1 = tm.ddot(F, F)
J = tm.linalg.det(F)
lnJ = tm.log(J)
return mu / 2 * (I1 - 3) - mu * lnJ + lmbda * lnJ ** 2 / 2
@LinearForm
def L(v, w):
F = deformation_gradient(w)
return tr.gradient_vector_product(neo_hooke, ntrax=2)(F, δx=grad(v))
@BilinearForm
def a(u, v, w):
F = deformation_gradient(w)
return tr.hessian_vector_product(neo_hooke, ntrax=2)(F, δx=grad(v), Δx=grad(u))Full code:r"""Hyperelasticity.
The strain energy density function per unit undeformed volume of the
isotropic hyperelastic Neo-Hookean material formulation is given as
.. math::
\psi(\boldsymbol{F}) = \frac{\mu}{2} ( I_1 - 3) - \mu \ln(J) +
\frac{\lambda}{2} \ln^2(J)
by the first invariant of the right Cauchy-Green deformation tensor
and the determinant of the deformation gradient (volume change).
.. math::
I_1 &= \text{tr}(\boldsymbol{F}^T \boldsymbol{F}) =
\boldsymbol{F} : \boldsymbol{F}
J &= \det(\boldsymbol{F})
The variation of the strain energy function leads to
.. math::
\delta\psi(\boldsymbol{F}) = \frac{\mu}{2} \delta I_1 - \mu \delta\ln(J)
+ \lambda \ln(J) \delta\ln(J)
with
.. math::
\delta I_1 &= 2 \boldsymbol{F} : \delta\boldsymbol{F}
\delta \ln(J) &= \boldsymbol{F}^{-T} : \delta\boldsymbol{F}
which finally results in
.. math::
\delta\psi(\boldsymbol{F}) = \mu \boldsymbol{F} : \delta\boldsymbol{F} -
(\lambda \ln(J) - \mu) \boldsymbol{F}^{-T} : \delta\boldsymbol{F}
The linearization of the variation leads to
.. math::
\Delta\delta\psi(\boldsymbol{F}) = \mu \Delta\boldsymbol{F} :
\delta\boldsymbol{F} + (\lambda \ln(J) - \mu)
\Delta\boldsymbol{F}^{-T} : \delta\boldsymbol{F}
+ \lambda \left( \boldsymbol{F}^{-T} : \delta\boldsymbol{F} \right)
\left( \boldsymbol{F}^{-T} : \Delta\boldsymbol{F} \right)
with
.. math::
\Delta\boldsymbol{F}^{-T} = -\boldsymbol{F}^{-T} \Delta\boldsymbol{F}^{T}
\boldsymbol{F}^{-T}
which may be alternatively formulated by expressions of the trace:
.. math::
\Delta\delta\psi(\boldsymbol{F}) = \mu \Delta\boldsymbol{F} :
\delta\boldsymbol{F} + (\lambda \ln(J) - \mu)
\text{tr}(\Delta\boldsymbol{F} \boldsymbol{F}^{-1}
\delta\boldsymbol{F} \boldsymbol{F}^{-1})
+ \lambda \text{tr}(\delta\boldsymbol{F} \boldsymbol{F}^{-1})
\text{tr}(\Delta\boldsymbol{F} \boldsymbol{F}^{-1})
"""
import numpy as np
from skfem import *
from skfem.helpers import grad, identity, ddot, det, transpose, inv, trace, mul
# note: rough mesh to make tests fast
mesh = MeshHex.init_tensor(
np.linspace(0, 1, 10),
np.linspace(-.1, .1, 3),
np.linspace(-.1, .1, 3),
).with_boundaries({
'left': lambda x: x[0] == 0.,
'right': lambda x: x[0] == 1.,
})
element = ElementVector(ElementHex1())
basis = Basis(mesh, element, intorder=1)
mu, lmbda = 1., 2.
def deformation_gradient(w):
dudX = grad(w["displacement"])
F = dudX + identity(dudX)
return F
import tensortrax as tr
import tensortrax.math as tm
def neo_hooke(F):
I1 = tm.ddot(F, F)
J = tm.linalg.det(F)
lnJ = tm.log(J)
return mu / 2 * (I1 - 3) - mu * lnJ + lmbda * lnJ ** 2 / 2
@LinearForm
def L(v, w):
F = deformation_gradient(w)
δF = grad(v)
return tr.gradient_vector_product(neo_hooke, ntrax=2)(F, δx=δF)
@BilinearForm
def a(u, v, w):
F = deformation_gradient(w)
ΔF = grad(u)
δF = grad(v)
return tr.hessian_vector_product(neo_hooke, ntrax=2)(F, δx=δF, Δx=ΔF)
u = basis.zeros()
du = basis.zeros()
uv = basis.interpolate(u)
dofs = basis.get_dofs('right')
dirichlet = basis.get_dofs({'right', 'left'})
f = L.assemble(basis, displacement=uv)
K = a.assemble(basis, displacement=uv)
right = -0.1
tol = 1e-10
nsteps = 7
for step in range(nsteps):
for iteration in range(10):
c = (step + 1.) / nsteps
du_D = u.copy()
x, y, z = basis.doflocs[:, dofs.nodal['u^1']]
du_D[dofs.nodal['u^1']] = c * right - du_D[dofs.nodal['u^1']]
du_D[dofs.nodal['u^2']] = (
y * np.cos(c * np.pi) - z * np.sin(c * np.pi) - y
- du_D[dofs.nodal['u^2']]
)
du_D[dofs.nodal['u^3']] = (
y * np.sin(c * np.pi) + z * np.cos(c * np.pi) - z
- du_D[dofs.nodal['u^3']]
)
du = solve(*condense(K, -f, x=du_D, D=dirichlet))
norm_du = np.linalg.norm(du)
u += du
uv = basis.interpolate(u)
f = L.assemble(basis, displacement=uv)
K = a.assemble(basis, displacement=uv)
print(1 + iteration, norm_du)
if norm_du < tol:
break
if __name__ == '__main__':
(mesh.translated(u[basis.nodal_dofs])
.draw('vedo', point_data={'uy': u[basis.nodal_dofs[1]]})
.show()) |
Beta Was this translation helpful? Give feedback.
-
|
Here is a picture I plan to use: The mesh is rough so that test suite runs fast. |
Beta Was this translation helpful? Give feedback.
-
|
See #841 |
Beta Was this translation helpful? Give feedback.
-
Hi,
as already mentioned in #703 I came across a solution how to drastically speed-up the assembly of bilinear forms for hyperelastic solids. The solution is to not explicitely calculate the deformation-dependent fourth-order elasticity tensor but instead formulate the problem in variational notation. Below follows an example - beside the speed-up this simplifies the analytic derivatives of the strain energy function. I really like that kind of coding 😎 👍🏻. In my opinion this is the way to go how to code hyperelastic forms in scikit-fem. Note that no fourth-order tensor is created inside the bilinear form. Please have a look at the following lines of theory + code and try it out for yourself:
Best regards,
Andreas
Edit: I've updated this post with the constitutive corrections from below (renamed
bulktolmbda, usemulhelper) so that if one only reads the original post at least the theory and the copy-paste code-block is up to date.Beta Was this translation helpful? Give feedback.
All reactions