Assembling forms with internal variables

[bhaveshshrimali]

I am trying to implement a rubber-viscoelastic model in scikit-fem. This can be thought of as an extension of hyperelasticity (ex36) but now with internal variables at each quadrature point. Physically, this would correspond to dissipative phenomena with history dependence such as plasticity, viscoelasticity, viscoplasticity. With a DiscreteField object, I feel this should be straightforward, namely

@njit(cache=True, nogil=True, parallel=True)
def F1(dv, du, p, Cv):
    out = np.zeros_like(du[0, 0])
    F = np.zeros_like(du)
    for a in prange(du.shape[2]):
        for b in prange(du.shape[3]):
            F[0, 0, a, b] += 1.
            F[1, 1, a, b] += 1.
            F[2, 2, a, b] += 1.
    
    F += du
    Finv, J = vinv(F)
    Cvinv, _ = vinv(Cv)
    FCviv = np.zeros_like(Cvinv)
    for i in range(F.shape[0]):
        for j in range(F.shape[1]):
            for k in range(F.shape[0]):
                for a in prange(F.shape[-2]):
                    for b in prange(F.shape[-1]):
                        FCviv[i, j, a, b] += F[i, k, a, b] * Cvinv[k, j, a, b]

    for i in range(du.shape[0]):
        for j in range(du.shape[1]):
            for a in prange(J.shape[0]):
                for b in prange(J.shape[1]):
                    out[a,b] += (
                        mu * F[i, j, a, b] + p[a,b] * J[a,b] * Finv[j, i, a, b] + nu * FCviv[i, j, a, b]
                    ) * dv[i, j, a, b]
    return out
.....
.....
@LinearForm(nthreads=0)
def a1(v, w):
    return F1(*(grad(v), grad(w["disp"]), w["press"].value, w["Cv"].value))

....
....
asm(a1, basis["u"], disp=uv, press = pv, Cv=Cv)

Here Cv is a numpy array of shape like du (in this case 3, 3, num_elems, num_quad_points) which can be updated at each time step at quadrature point by some update rule.

Am I thinking this correctly? Or is there something that I am missing. I will test things tomorrow. The complete code is quite long. Let me see if I can shorten it up a bit to illustrate the main points.

I'll also put together a non-numba version as the code above is excessively convoluted.

[gdmcbain]

Am I thinking this correctly? Or is there something that I am missing.

I think this is correct. This is the kind of thing that I was learning about back in #24, which led to ex16 ‘Legendre's equation’.

A recent update is that the coefficients no longer need to be interpolated onto the abscissæ before passing to the assembler, in case that helps. #697

[bhaveshshrimali]

Perfect, thanks @gdmcbain !!

I'll be back on this tomorrow.

[kinnala]

Migrated to Discussions.

[bhaveshshrimali]

Here is a simple proof of concept. It is still surprising to me on how easy it was to incorporate this into existing code. The material model implemented is that of a Neo-Hookean solid with constant viscous dissipation, eqn-(8)-(9). And the boundary value problem is simple uniaxial loading and unloading.

The convergence isn't Newton-like, but that is likely due to the fact that my tangent modulus isn't correct (on what is commonly referred to as the "consistent-tangent"). The code isn't super clean but I wanted to just check if it worked at all, and indeed it does.

Code

import numpy as np, os, meshio
from scipy.sparse import bmat
from skfem.helpers import grad, transpose, det, inv, identity
from skfem import *
from numba import prange, njit
import numpy.linalg as nla

# material parameters
mu, lmbda = 1., 1.e3
nu, eta = 1., 1.

def doubleInner(A, B):
    return np.einsum("ij...,ij...", A, B)

def F1(w):
    u = w["disp"]
    p = w["press"]
    Cv = w["Cv"]
    Cvinv = inv(Cv.value)
    F = grad(u) + identity(u)
    Finv = inv(F)
    J = det(F)
    FCvinv = np.einsum("ik...,kj...->ij...", F, Cvinv)
    return mu * F + p * J * transpose(Finv) + nu * FCvinv

def sPiola(w):
    u = w["disp"]
    p = w["press"]
    Cv = w["Cv"]
    Cvinv = inv(Cv)
    F = grad(u) + identity(u)
    Finv = inv(F)
    J = det(F)
    FCvinv = np.einsum("ik...,kj...->ij...", F, Cvinv)
    return mu * F + p * J * transpose(Finv) + nu * FCvinv

def sPiola33(w):
    return F1(w)[2, 2]

def F2(w):
    u = w["disp"]
    p = w["press"].value
    F = grad(u) + identity(u)
    J = det(F)
    Js = .5 * (lmbda + p + 2. * np.sqrt(lmbda * mu + .25 * (lmbda + p) ** 2)) / lmbda
    dJsdp = ((.25 * lmbda + .25 * p + .5 * np.sqrt(lmbda * mu + .25 * (lmbda + p) ** 2))
             / (lmbda * np.sqrt(lmbda * mu + .25 * (lmbda + p) ** 2)))
    return J - (Js + (p + mu / Js - lmbda * (Js - 1)) * dJsdp)

def A11(w):
    u = w["disp"]
    p = w["press"]
    Cv = w["Cv"]
    eye = identity(u)
    F = grad(u) + eye
    J = det(F)
    Cvinv = inv(Cv.value)

    Finv = inv(F)
    L = (p * J * np.einsum("lk...,ji...->ijkl...", Finv, Finv)
         - p * J * np.einsum("jk...,li...->ijkl...", Finv, Finv)
         + mu * np.einsum("ik...,jl...->ijkl...", eye, eye) + 
         + nu * np.einsum("ik...,lj...->ijkl...", eye, Cvinv))
    return L

def A12(w):
    u = w["disp"]
    F = grad(u) + identity(u)
    J = det(F)
    Finv = inv(F)
    return J * transpose(Finv)

def A22(w):
    u = w["disp"]
    p = w["press"].value
    Js = .5 * (lmbda + p + 2. * np.sqrt(lmbda * mu + .25 * (lmbda + p) ** 2)) / lmbda
    dJsdp = ((.25 * lmbda + .25 * p + .5 * np.sqrt(lmbda * mu + .25 * (lmbda + p) ** 2))
             / (lmbda * np.sqrt(lmbda * mu + .25 * (lmbda + p) ** 2)))
    d2Jdp2 = .25 * mu / (lmbda * mu + .25 * (lmbda + p) ** 2) ** (3/2)
    L = (-2. * dJsdp - p * d2Jdp2 + mu / Js ** 2 * dJsdp ** 2 - mu / Js * d2Jdp2
         + lmbda * (Js - 1.) * d2Jdp2 + lmbda * dJsdp ** 2)
    return L

def volume(w):
    dw = w["disp"].grad
    F = dw + identity(dw)
    J = det(F)
    return J

def Cvdot(Cv, C, nu_m, eta_m):
    Cvinv = inv(Cv)
    CCvinv = doubleInner(C, Cvinv)
    Cvdot = nu_m/eta_m * (C - 1./3 * CCvinv * Cv)
    return Cvdot

def calculateCStrain(gradu):
    F = gradu + identity(gradu)
    C = np.einsum("ki...,kj...->ij...", F, F)
    return C

def evolEq(dt, Cvn, C, Cn, nu_m, eta_m):
    C_half = Cn + 0.5 * (C - Cn)
    C_quart = Cn + 0.25 * (C - Cn)
    C_three_quart = Cn + 0.75 * (C - Cn)
    G1 = Cvdot(Cvn, Cn, nu_m, eta_m)
    G2 = Cvdot(Cvn + G1 * dt/2., C_half, nu_m, eta_m)
    G3 = Cvdot(Cvn + dt/16. * (3. * G1 + G2), C_quart, nu_m, eta_m)
    G4 = Cvdot(Cvn + G3 * dt/2., C_half, nu_m, eta_m)
    G5 = Cvdot(Cvn + 3. * dt/16. * (-G2 + 2. * G3 + 3. * G4), C_three_quart, nu_m, eta_m)
    G6 = Cvdot(Cvn + dt/7. * (G1 + 4. * G2 + 6. * G3 - 12. * G4 + 8. * G5), C, nu_m, eta_m)
    return Cvn + dt / 90.0 * (7.0 * G1 + 32.0 * G3 + 12.0 * G4 + 32.0 * G5 + 7.0 * G6)


def solveSystem(du, dp, Cv, I, basis):
    uv = basis["u"].interpolate(du)
    pv = basis["p"].interpolate(dp)
    K11 = asm(b11, basis["u"], basis["u"], disp=uv, press = pv, Cv=Cv)
    K12 = asm(b12, basis["p"], basis["u"], disp=uv, press = pv)
    K22 = asm(b22, basis["p"], basis["p"], disp=uv, press = pv)
    f = np.concatenate((
        asm(a1, basis["u"], disp=uv, press = pv, Cv=Cv),
        asm(a2, basis["p"], disp=uv, press = pv)
    ))
    K = bmat(
        [[K11, K12],
         [K12.T, K22]], "csr"
    )
    uvp = solve(*condense(K, -f, I=I))
    delu, delp = np.split(uvp, [du.shape[0]])
    du += delu
    dp += delp
    normu = nla.norm(delu)
    normp = nla.norm(delp)
    norm_res = nla.norm(f[I])
    # print(f"{itr+1}, norm_du: {normu}, norm_dp: {normp}, norm_res = {norm_res}")
    sbar = stress.assemble(basis["u"], disp=uv, press=pv, Cv=Cv)
    Cstrain = calculateCStrain(uv.grad)
    return (norm_res, np.copy(du), np.copy(dp), sbar, Cstrain)
    

# mesh = from_meshio(meshio.xdmf.read("fenicsmesh.xdmf"))
mesh = MeshTet()
uelem = ElementVectorH1(ElementTetP2())
pelem = ElementTetP1()
elems = {
    "u": uelem,
    "p": pelem
}
basis = {
    field: Basis(mesh, e, intorder=2)
    for field, e in elems.items()
}

du = np.zeros(basis["u"].N)
dp = -1. * (mu + nu) * np.ones(basis["p"].N)

stretch_ = 0.0

dofsold = [
    basis["u"].find_dofs({"left":mesh.facets_satisfying(lambda x: x[0] < 1.e-6)}, skip=["u^2", "u^3"]),
    basis["u"].find_dofs({"bottom":mesh.facets_satisfying(lambda x: x[1] < 1.e-6)}, skip=["u^1", "u^3"]),
    basis["u"].find_dofs({"back":mesh.facets_satisfying(lambda x: x[2] < 1.e-6)}, skip=["u^1", "u^2"]),
    basis["u"].find_dofs({"front":mesh.facets_satisfying(lambda x: np.abs(x[2]-1.) < 1.e-6 )}, skip=["u^1", "u^2"])
]

dofs = {}
for dof in dofsold:
    dofs.update(dof)

du[dofs["left"].all()] = 0.
du[dofs["bottom"].all()] = 0.
du[dofs["back"].all()] = 0.
du[dofs["front"].all()] = stretch_

I = np.hstack((
    basis["u"].complement_dofs(dofs),
    basis["u"].N + np.arange(basis["p"].N)
))

@LinearForm
def a1(v, w):
    return np.einsum("ij...,ij...", F1(w), grad(v))

@LinearForm
def a2(v, w):
    return F2(w) * v

@BilinearForm
def b11(u, v, w):
    return np.einsum("ijkl...,ij...,kl...", A11(w), grad(u), grad(v))

@BilinearForm
def b12(u, v, w):
    return np.einsum("ij...,ij...", A12(w), grad(v)) * u

@BilinearForm
def b22(u, v, w):
    return A22(w) * u * v

@Functional
def vol(w):
    return volume(w)

@Functional
def stress(w):
    return sPiola33(w)

ldot = 1.0
dt = 0.05
final_stretch = 2.
Tf = 2. * final_stretch/ldot
num_steps = Tf/dt
timevals = np.linspace(0., Tf, int(num_steps + 1))
t1 = timevals[timevals <= Tf/2.]
t2 = timevals[timevals > Tf/2.]
stretchVals = np.hstack((
    np.linspace(1., final_stretch, t1.shape[0]),
    np.linspace(final_stretch - ldot * dt, 1., t2.shape[0])
))

meshioPoints = mesh.p.T
meshioCells = [("tetra", mesh.t.T)]
results_dir = os.path.join(os.getcwd(), "resultsViscoTimeSteps")
os.makedirs(results_dir, exist_ok=True)
filename = os.path.join(results_dir, f"disps_{ldot}.xdmf")
avgStress = np.zeros_like(stretchVals)
func_dict = {
    "disp": duv,
    "press": basis["p"].interpolate(dp),
    "Cv": Cv
}

with meshio.xdmf.TimeSeriesWriter(filename) as writer:
    writer.write_points_cells(meshioPoints, meshioCells)
    for idx, tv in enumerate(timevals):
        print(f"Time: {tv} of {Tf}")
        dispVal = stretchVals[idx] - 1.
        du[dofs["front"].all()] = dispVal
        maxiters = 10
        for itr in range(maxiters):
            norm_res, du, dp, sbar, Cstrain = solveSystem(du, dp, Cv, I, basis)
            print(f"res: {norm_res}")
            Cv = evolEq(dt, Cvn, Cstrain, Cn, nu, eta)
            if norm_res < 1.e-6:
                converged = True
                break
        Cvn = Cv.copy()
        Cn = Cstrain.copy()
        vol_def = vol.assemble(basis["u"], disp=basis["u"].interpolate(du))
        writer.write_data(
            tv, point_data={
                "u": du[basis["u"].nodal_dofs].T,
                "p": dp[basis["p"].nodal_dofs].T,
            }
        )
        avgStress[idx] = sbar

np.savetxt(os.path.join(results_dir, f"stresses_{ldot}.txt"), avgStress)
np.savetxt(os.path.join(results_dir, f"stretches_{ldot}.txt"), stretchVals)
np.savetxt(os.path.join(results_dir, f"timevals_{ldot}.txt"), timevals)

Results

I reproduced a simple uniaxial calculation (solving an ODE instead of the entire full blown 3D elasticity) and it is correct.

[bhaveshshrimali]

As a side note, I think computing the deformed volume in ex36 is also a more meaningful check. Given that it is "nearly incompressible hyperelasticity".

[bhaveshshrimali]

Let me know what you think. I can change/append the tests for example 36 accordingly

[kinnala]

Sounds good. Any improvements to the tests are welcome.

[kinnala]

No need to remove any tests though, just append to it if you feel like it.

[bhaveshshrimali]

OK, will do it later today. Thanks

[bhaveshshrimali]

Hmmm, switching to master now doesn't allow one to use data directly at the quadrature points in the form of np.ndarray. I think this is what @gdmcbain alluded to above

Changed: Form.assemble allows now only DiscreteField or 1d/2d ndarray objects, and will automatically pass 1d ndarray objects to Basis.interpolate for convenience

The above examples raises

Traceback (most recent call last):
  File "D:\viscoelasticityMixed.py", line 296, in <module>
    norm_res, du, dp, sbar, Cstrain = solveSystem(du, dp, Cv, I, basis)
  File "D:\viscoelasticityMixed.py", line 128, in solveSystem
    K11 = asm(b11, basis["u"], basis["u"], disp=uv, press = pv, Cv=Cv)
  File "C:\Users\bshri\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\assembly\__init__.py", line 70, in asm
    return form.assemble(*args, **kwargs)
  File "C:\Users\bshri\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\assembly\form\bilinear_form.py", line 142, in assemble
    return COOData._assemble_scipy_csr(*self._assemble(*args, **kwargs))
  File "C:\Users\bshri\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\assembly\form\bilinear_form.py", line 73, in _assemble
    **self.dictify(kwargs, ubasis),
  File "C:\Users\bshri\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\assembly\form\form.py", line 60, in 
dictify
    raise ValueError("The given type '{}' for the list of extra "
ValueError: The given type '<class 'numpy.ndarray'>' for the list of extra form parameters w cannot be converted to DiscreteField.

I am not sure what would be the best way around here. One way could be to create something like an ElementTensor and interpolate at the quadrature points to get a tensorial DiscreteField, or pass the components of the entire tensor as an array of DiscreteFields. Will try to see if I can get the latter to work somehow

To add this wouldn't be a problem with 1D arrays like a variable material (scalar) property as those can be interpolated at quadrature points using existing machinery. Happy to be proven wrong though.

[bhaveshshrimali]

Oh wait, I can simply do DiscreteField(my_array), where my_array is of shape (3,3,num_cells,num_qps) in my case, to get it at the quadrature points. It seems to work with 3.2.0. Let me test it on master. Got this from #564

Edit: It works on master :)