In [1]:
import sys
from sympy import *
from numpy import identity
import textwrap

In [2]:
WRITE_TO_FILE = False

In [3]:
Ident = Matrix(3,3, lambda i,j: 1 if i==j else 0)
ij_from_I = {0:(0,0),1:(1,1),2:(2,2),3:(0,1),4:(1,2),5:(0,2)}
I_from_ij = lambda m,n: dict([(v,k) for k,v in ij_from_I.items()])[m,n]
ij_from_I_9 = {0:(0,0),1:(0,1),2:(0,2),3:(1,0),4:(1,1),5:(1,2),6:(2,0),7:(2,1),8:(2,2)}
I_from_ij_9 = lambda m,n: dict([(v,k) for k,v in ij_from_I_9.items()])[m,n]
def _Matrix(M,ij):
    return Matrix(3,3,[Symbol('{}[{}]'.format(M,i)) for i in ij])
def SymMatrix(M, full=0):
    return _Matrix(M, (0,3,5,3,1,4,5,4,2))
def UnsymMatrix(M):
    return _Matrix(M, (0,1,2,3,4,5,6,7,8))

In [4]:
if WRITE_TO_FILE:
    STREAM = open('../materials/tensor.py', 'w')
else:
    STREAM = sys.stdout

In [5]:
def wrap(lines):
    kwds = {'width': 78,
            'break_long_words': False,
            'replace_whitespace': False}
    lines = lines.split('\n')
    for (i, line) in enumerate(lines):
        if len(line) <= kwds['width']:
            continue
        kwds['subsequent_indent'] = None
        if line.strip()[0].strip() and ' = ' in line:
            var, expr = line.split('=')
            kwds['subsequent_indent'] = ' ' * (len(var)+3)
            line = var + '= ({})'.format(expr.strip())
        lines[i] = textwrap.fill(line, **kwds)
    return '\n'.join(lines)

def write(string, wrapit=1):
    if wrapit:
        string = wrap(string)
    STREAM.write(string)
    STREAM.write('\n')
    STREAM.flush()

In [6]:
imports = """\
import numpy as np
from numpy import array, zeros"""
write(imports)
man = """tensor.py

THIS FILE IS AUTOMATICALLY GENERATED - DO NOT EDIT IT DIRECTLY!

Constants and functions having to deal with 3D symmetric and nonsymmetric 
second order tensors fourth-order tensors with minor and major symmetries.  
All of the functions are written with the following assumptions:

o Symmetric second-order tensors are stored as arrays of length 6 with the 
  following component ordering
      
  [XX, YY, ZZ, XY, YZ, ZZ]
      
o Nonsymmetric second-order tensors are stored as arrays of length 9 with the 
  following component ordering
    
  [XX, XY, XZ, YX, YY, YZ, ZX, ZY, ZZ]
      
o Fourth-order tensors are stored as 6x6 matrices using the same component 
  transformations as second-order symmetric tensors
      
"""
write('"""{}"""'.format(man), wrapit=0)


import numpy as np
from numpy import array, zeros
"""tensor.py

THIS FILE IS AUTOMATICALLY GENERATED - DO NOT EDIT IT DIRECTLY!

Constants and functions having to deal with 3D symmetric and nonsymmetric 
second order tensors fourth-order tensors with minor and major symmetries.  
All of the functions are written with the following assumptions:

o Symmetric second-order tensors are stored as arrays of length 6 with the 
  following component ordering
      
  [XX, YY, ZZ, XY, YZ, ZZ]
      
o Nonsymmetric second-order tensors are stored as arrays of length 9 with the 
  following component ordering
    
  [XX, XY, XZ, YX, YY, YZ, ZX, ZY, ZZ]
      
o Fourth-order tensors are stored as 6x6 matrices using the same component 
  transformations as second-order symmetric tensors
      
"""

In [7]:
d = identity(3, dtype=float)
II1 = []
II2 = []
II3 = []
II5 = []
for I in range(6):
    for J in range(6):
        i,j = ij_from_I[I]
        k,l = ij_from_I[J]
        II1.append('{0:.1f}'.format(d[i,j] * d[k,l]))
        II2.append('{0:.1f}'.format(d[i,k] * d[j,l]))
        II3.append('{0:.1f}'.format(d[i,l] * d[j,k]))
        II5.append('{0:.1f}'.format((d[i,k] * d[j,l] + d[i,l]*d[j,k])/2.))

write('# Fourth-order "identities"')
II1 = ', '.join(II1)
write('# II1: I[i,j] I[k,l]')
write('II1 = array([{}]).reshape(6,6)'.format(II1))

II2 = ', '.join(II2)
write('# II2: I[i,k] I[j,l]')
write('II2 = array([{}]).reshape(6,6)'.format(II2))

II3 = ', '.join(II3)
write('# II3: I[i,l] I[j,k]')
write('II3 = array([{}]).reshape(6,6)'.format(II3))

write('II4 = (II2 + II3) / 2')
II5 = ', '.join(II5)
write('# II5 = (I[i,k] I[j,l] + I[i,l] I[j,k]) / 2')
write('II5 = array([{}]).reshape(6,6)'.format(II5))

write('# Second-order identities')
write('I6 = array([1., 1., 1., 0., 0., 0.])')
write('I9 = array([1., 0., 0., 0., 1., 0., 0., 0., 1.])')
write('I3x3 = I9.reshape(3,3)')


# Fourth-order "identities"
# II1: I[i,j] I[k,l]
II1 = (array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0,
       1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]).reshape(6,6))
# II2: I[i,k] I[j,l]
II2 = (array([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]).reshape(6,6))
# II3: I[i,l] I[j,k]
II3 = (array([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]).reshape(6,6))
II4 = (II2 + II3) / 2
# II5 = (I[i,k] I[j,l] + I[i,l] I[j,k]) / 2
II5 = (array([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5]).reshape(6,6))
# Second-order identities
I6 = array([1., 1., 1., 0., 0., 0.])
I9 = array([1., 0., 0., 0., 1., 0., 0., 0., 1.])
I3x3 = I9.reshape(3,3)

In [8]:
F = UnsymMatrix('F')
FtF = F.T * F
write('def symsq(F):')
write('    """ Computes dot(F.T, F)"""')

write('    X = zeros(6)')
for I in range(6):
    i,j = ij_from_I[I]
    write('    X[{}] = {}'.format(I, FtF[i,j]))
write('    return X')


def symsq(F):
    """ Computes dot(F.T, F)"""
    X = zeros(6)
    X[0] = F[0]**2 + F[3]**2 + F[6]**2
    X[1] = F[1]**2 + F[4]**2 + F[7]**2
    X[2] = F[2]**2 + F[5]**2 + F[8]**2
    X[3] = F[0]*F[1] + F[3]*F[4] + F[6]*F[7]
    X[4] = F[1]*F[2] + F[4]*F[5] + F[7]*F[8]
    X[5] = F[0]*F[2] + F[3]*F[5] + F[6]*F[8]
    return X

In [9]:
F = UnsymMatrix('F')
A = SymMatrix('A')
dA = A.det()
dF = F.det()
write('def det(A):')
write('    """ Computes the determininant of A"""')
write('    if A.size == 6:')
write('        X = {}'.format(dA))
write('    else:')
write('        X = {}'.format(dF).replace('F', 'A'))
write('    return X')


def det(A):
    """ Computes the determininant of A"""
    if A.size == 6:
        X = (A[0]*A[1]*A[2] - A[0]*A[4]**2 - A[1]*A[5]**2 - A[2]*A[3]**2 +
             2*A[3]*A[4]*A[5])
    else:
        X = (A[0]*A[4]*A[8] - A[0]*A[5]*A[7] - A[1]*A[3]*A[8] + A[1]*A[5]*A[6]
             + A[2]*A[3]*A[7] - A[2]*A[4]*A[6])
    return X

In [10]:
As = SymMatrix('A')
Bs = SymMatrix('B')
Au = UnsymMatrix('A')
Bu = UnsymMatrix('B')
AsBs = ' + '.join([str(As[i,j]*Bs[i,j]) for i in range(3) for j in range(3)])
AsBu = ' + '.join([str(As[i,j]*Bu[i,j]) for i in range(3) for j in range(3)])
AuBu = ' + '.join([str(Au[i,j]*Bu[i,j]) for i in range(3) for j in range(3)])
write('def ddot(A, B):')
write('    """ Computes A:B"""')
write('    if B.size == 6 and A.size == 9:')
write('        A, B = B, A')
write('    if A.size == 6 and B.size == 6:')
write('        X = {}'.format(AsBs))
write('    elif A.size == 6 and B.size == 9:')
write('        X = {}'.format(AsBu))
write('    elif A.size == 9 and B.size == 9:')
write('        X = {}'.format(AuBu))
write('    else:\n        raise NotImplementedError')
write('    return X')


def ddot(A, B):
    """ Computes A:B"""
    if B.size == 6 and A.size == 9:
        A, B = B, A
    if A.size == 6 and B.size == 6:
        X = (A[0]*B[0] + A[3]*B[3] + A[5]*B[5] + A[3]*B[3] + A[1]*B[1] +
             A[4]*B[4] + A[5]*B[5] + A[4]*B[4] + A[2]*B[2])
    elif A.size == 6 and B.size == 9:
        X = (A[0]*B[0] + A[3]*B[1] + A[5]*B[2] + A[3]*B[3] + A[1]*B[4] +
             A[4]*B[5] + A[5]*B[6] + A[4]*B[7] + A[2]*B[8])
    elif A.size == 9 and B.size == 9:
        X = (A[0]*B[0] + A[1]*B[1] + A[2]*B[2] + A[3]*B[3] + A[4]*B[4] +
             A[5]*B[5] + A[6]*B[6] + A[7]*B[7] + A[8]*B[8])
    else:
        raise NotImplementedError
    return X

In [11]:
write('def trace(A, metric=None):')
write('    """Computes the trace of a tensor"""', wrapit=0)
write('    if A.size not in (6, 9):')
write('        raise NotImplementedError')
write('    if A.size == 6 and metric is None:')
write('        metric, X = I6, I6')
write('    elif A.size == 9 and metric is None:')
write('        metric, X = I9, I9')
write('    else:')
write('        X = inv(metric)')
write('    return ddot(A, metric)')
write('def iso(A, metric=None):')
write('    """Computes the isotropic part of a tensor"""', wrapit=0)
write('    if A.size not in (6, 9):')
write('        raise NotImplementedError')
write('    if A.size == 6 and metric is None:')
write('        metric, X = I6, I6')
write('    elif A.size == 9 and metric is None:')
write('        metric, X = I9, I9')
write('    else:')
write('        X = inv(metric)')
write('    return trace(A, metric) / 3. * X')
write('def dev(A, metric=None):')
write('    """Computes the isotropic part of a tensor"""', wrapit=0)
write('    return A - iso(A, metric=metric)')
write('def mag(A):\n    """Computes the magnitude of a tensor"""\n    return np.sqrt(ddot(A, A))', wrapit=0)


def trace(A, metric=None):
    """Computes the trace of a tensor"""
    if A.size not in (6, 9):
        raise NotImplementedError
    if A.size == 6 and metric is None:
        metric, X = I6, I6
    elif A.size == 9 and metric is None:
        metric, X = I9, I9
    else:
        X = inv(metric)
    return ddot(A, metric)
def iso(A, metric=None):
    """Computes the isotropic part of a tensor"""
    if A.size not in (6, 9):
        raise NotImplementedError
    if A.size == 6 and metric is None:
        metric, X = I6, I6
    elif A.size == 9 and metric is None:
        metric, X = I9, I9
    else:
        X = inv(metric)
    return trace(A, metric) / 3. * X
def dev(A, metric=None):
    """Computes the isotropic part of a tensor"""
    return A - iso(A, metric=metric)
def mag(A):
    """Computes the magnitude of a tensor"""
    return np.sqrt(ddot(A, A))

In [12]:
A = SymMatrix('A')
write('def invariants(A, itype=0):', wrapit=0)
write('''\
    """ Computes the invariants of A

        Parameters
        ----------
        itype : int
            Type of invariants to compute
            
        Returns
        -------
        A, B, C : float
        
        Notes
        -----
        The values A, B, C depend on itype:
            =0 -> I1, I2, I3
            =1 -> I1, J2, J3
    """
''', wrapit=0)
I1 = A.trace()
I2 = (I1 ** 2 - (A * A).trace()) / 2
S = A - A.trace() / 3 * Ident
J2 = (S * S).trace() / 2
J3 = (S * S * S).trace() / 3
write('    if itype == 0:')
write('        I1 = {}'.format(I1))
write('        I2 = {}'.format(I2))
write('        I3 = det(A)')
write('        A, B, C = I1, I2, I3')
write('    elif itype == 1:')
write('        I1 = {}'.format(I1))
write('        J2 = {}'.format(J2))
write('        J3 = {}'.format(J3))
write('        A, B, C = I1, J2, J3')
write('    else:\n        raise NotImplementedError')
write('    return A, B, C')


def invariants(A, itype=0):
    """ Computes the invariants of A

        Parameters
        ----------
        itype : int
            Type of invariants to compute
            
        Returns
        -------
        A, B, C : float
        
        Notes
        -----
        The values A, B, C depend on itype:
            =0 -> I1, I2, I3
            =1 -> I1, J2, J3
    """

    if itype == 0:
        I1 = A[0] + A[1] + A[2]
        I2 = (-A[0]**2/2 - A[1]**2/2 - A[2]**2/2 - A[3]**2 - A[4]**2 - A[5]**2
              + (A[0] + A[1] + A[2])**2/2)
        I3 = det(A)
        A, B, C = I1, I2, I3
    elif itype == 1:
        I1 = A[0] + A[1] + A[2]
        J2 = (A[3]**2 + A[4]**2 + A[5]**2 + (-A[0]/3 - A[1]/3 + 2*A[2]/3)**2/2
              + (-A[0]/3 + 2*A[1]/3 - A[2]/3)**2/2 + (2*A[0]/3 - A[1]/3 -
              A[2]/3)**2/2)
        J3 = (2*A[3]*(A[3]*(-A[0]/3 + 2*A[1]/3 - A[2]/3) + A[3]*(2*A[0]/3 -
              A[1]/3 - A[2]/3) + A[4]*A[5])/3 + 2*A[4]*(A[3]*A[5] +
              A[4]*(-A[0]/3 - A[1]/3 + 2*A[2]/3) + A[4]*(-A[0]/3 + 2*A[1]/3 -
              A[2]/3))/3 + 2*A[5]*(A[3]*A[4] + A[5]*(-A[0]/3 - A[1]/3 +
              2*A[2]/3) + A[5]*(2*A[0]/3 - A[1]/3 - A[2]/3))/3 + (-A[0]/3 -
              A[1]/3 + 2*A[2]/3)*(A[4]**2 + A[5]**2 + (-A[0]/3 - A[1]/3 +
              2*A[2]/3)**2)/3 + (-A[0]/3 + 2*A[1]/3 - A[2]/3)*(A[3]**2 +
              A[4]**2 + (-A[0]/3 + 2*A[1]/3 - A[2]/3)**2)/3 + (2*A[0]/3 -
              A[1]/3 - A[2]/3)*(A[3]**2 + A[5]**2 + (2*A[0]/3 - A[1]/3 -
              A[2]/3)**2)/3)
        A, B, C = I1, J2, J3
    else:
        raise NotImplementedError
    return A, B, C

In [13]:
F = UnsymMatrix('F')
A = SymMatrix('A')
X = F * A * F.T
write('def push(F, A):')
write('    """Computes the push operation F A F.T / J"""', wrapit=0)
write('    if A.size == 6:\n        return push6(F, A)')
write('    elif A.size == 36:\n        return push66(F, A)')
write('    raise NotImplementedError')
write('def push6(F, A):')
write('    X = zeros(6)')
for I in range(6):
    i,j = ij_from_I[I]
    write('    X[{}] = {}'.format(I, X[i,j]))
write('    return X / det(F)')
write('def push66(F, A):')
write('    Q = symleaf(F)')
write('    X = np.dot(np.dot(Q, A), Q.T)')
write('    return X / det(F)')


def push(F, A):
    """Computes the push operation F A F.T / J"""
    if A.size == 6:
        return push6(F, A)
    elif A.size == 36:
        return push66(F, A)
    raise NotImplementedError
def push6(F, A):
    X = zeros(6)
    X[0] = (F[0]*(A[0]*F[0] + A[3]*F[1] + A[5]*F[2]) + F[1]*(A[1]*F[1] +
            A[3]*F[0] + A[4]*F[2]) + F[2]*(A[2]*F[2] + A[4]*F[1] + A[5]*F[0]))
    X[1] = (F[3]*(A[0]*F[3] + A[3]*F[4] + A[5]*F[5]) + F[4]*(A[1]*F[4] +
            A[3]*F[3] + A[4]*F[5]) + F[5]*(A[2]*F[5] + A[4]*F[4] + A[5]*F[3]))
    X[2] = (F[6]*(A[0]*F[6] + A[3]*F[7] + A[5]*F[8]) + F[7]*(A[1]*F[7] +
            A[3]*F[6] + A[4]*F[8]) + F[8]*(A[2]*F[8] + A[4]*F[7] + A[5]*F[6]))
    X[3] = (F[3]*(A[0]*F[0] + A[3]*F[1] + A[5]*F[2]) + F[4]*(A[1]*F[1] +
            A[3]*F[0] + A[4]*F[2]) + F[5]*(A[2]*F[2] + A[4]*F[1] + A[5]*F[0]))
    X[4] = (F[6]*(A[0]*F[3] + A[3]*F[4] + A[5]*F[5]) + F[7]*(A[1]*F[4] +
            A[3]*F[3] + A[4]*F[5]) + F[8]*(A[2]*F[5] + A[4]*F[4] + A[5]*F[3]))
    X[5] = (F[6]*(A[0]*F[0] + A[3]*F[1] + A[5]*F[2]) + F[7]*(A[1]*F[1] +
            A[3]*F[0] + A[4]*F[2]) + F[8]*(A[2]*F[2] + A[4]*F[1] + A[5]*F[0]))
    return X / det(F)
def push66(F, A):
    Q = symleaf(F)
    X = np.dot(np.dot(Q, A), Q.T)
    return X / det(F)

In [14]:
F = UnsymMatrix('F')
A = SymMatrix('A')
Ai = A.inv()
Fi = F.inv()
write('def inv(A):')
write('    """Computes the inverse of A"""', wrapit=0)
write('    if A.size == 6:')
write('        X = zeros(6)')
for I in range(6):
    i,j = ij_from_I[I]
    write('        X[{}] = {}'.format(I, Ai[i,j]))
write('    elif A.size == 9:')
write('        X = zeros(9)')
for I in range(9):
    i,j = ij_from_I_9[I]
    write('        X[{}] = {}'.format(I, Fi[i,j]))
#write('        X = np.linalg.inv(A.reshape(3,3))')
#write('        X = X.reshape(9)')
write('    else:\n        raise NotImplementedError')
write('    return X')


def inv(A):
    """Computes the inverse of A"""
    if A.size == 6:
        X = zeros(6)
        X[0] = (-(A[0]*A[1] - A[3]**2)*(-A[3]*(A[4] -
                A[3]*A[5]/A[0])/(A[0]*(A[1] - A[3]**2/A[0])) +
                A[5]/A[0])*(A[3]*(A[4] - A[3]*A[5]/A[0])/(A[0]*(A[1] -
                A[3]**2/A[0])) - A[5]/A[0])/(A[0]*A[1]*A[2] - A[0]*A[4]**2 -
                A[1]*A[5]**2 - A[2]*A[3]**2 + 2*A[3]*A[4]*A[5]) + 1/A[0] +
                A[3]**2/(A[0]**2*(A[1] - A[3]**2/A[0])))
        X[1] = (1/(A[1] - A[3]**2/A[0]) + (A[4] -
                A[3]*A[5]/A[0])**2*(A[0]*A[1] - A[3]**2)/((A[1] -
                A[3]**2/A[0])**2*(A[0]*A[1]*A[2] - A[0]*A[4]**2 - A[1]*A[5]**2
                - A[2]*A[3]**2 + 2*A[3]*A[4]*A[5])))
        X[2] = ((A[0]*A[1] - A[3]**2)/(A[0]*A[1]*A[2] - A[0]*A[4]**2 -
                A[1]*A[5]**2 - A[2]*A[3]**2 + 2*A[3]*A[4]*A[5]))
        X[3] = ((A[4] - A[3]*A[5]/A[0])*(A[0]*A[1] - A[3]**2)*(-A[3]*(A[4] -
                A[3]*A[5]/A[0])/(A[0]*(A[1] - A[3]**2/A[0])) +
                A[5]/A[0])/((A[1] - A[3]**2/A[0])*(A[0]*A[1]*A[2] -
                A[0]*A[4]**2 - A[1]*A[5]**2 - A[2]*A[3]**2 +
                2*A[3]*A[4]*A[5])) - A[3]/(A[0]*(A[1] - A[3]**2/A[0])))
        X[4] = (-(A[4] - A[3]*A[5]/A[0])*(A[0]*A[1] - A[3]**2)/((A[1] -
                A[3]**2/A[0])*(A[0]*A[1]*A[2] - A[0]*A[4]**2 - A[1]*A[5]**2 -
                A[2]*A[3]**2 + 2*A[3]*A[4]*A[5])))
        X[5] = (-(A[0]*A[1] - A[3]**2)*(-A[3]*(A[4] -
                A[3]*A[5]/A[0])/(A[0]*(A[1] - A[3]**2/A[0])) +
                A[5]/A[0])/(A[0]*A[1]*A[2] - A[0]*A[4]**2 - A[1]*A[5]**2 -
                A[2]*A[3]**2 + 2*A[3]*A[4]*A[5]))
    elif A.size == 9:
        X = zeros(9)
        X[0] = (-(F[0]*F[4] - F[1]*F[3])*(-F[1]*(F[5] -
                F[2]*F[3]/F[0])/(F[0]*(F[4] - F[1]*F[3]/F[0])) +
                F[2]/F[0])*(F[3]*(F[7] - F[1]*F[6]/F[0])/(F[0]*(F[4] -
                F[1]*F[3]/F[0])) - F[6]/F[0])/(F[0]*F[4]*F[8] - F[0]*F[5]*F[7]
                - F[1]*F[3]*F[8] + F[1]*F[5]*F[6] + F[2]*F[3]*F[7] -
                F[2]*F[4]*F[6]) + 1/F[0] + F[1]*F[3]/(F[0]**2*(F[4] -
                F[1]*F[3]/F[0])))
        X[1] = ((F[7] - F[1]*F[6]/F[0])*(F[0]*F[4] - F[1]*F[3])*(-F[1]*(F[5] -
                F[2]*F[3]/F[0])/(F[0]*(F[4] - F[1]*F[3]/F[0])) +
                F[2]/F[0])/((F[4] - F[1]*F[3]/F[0])*(F[0]*F[4]*F[8] -
                F[0]*F[5]*F[7] - F[1]*F[3]*F[8] + F[1]*F[5]*F[6] +
                F[2]*F[3]*F[7] - F[2]*F[4]*F[6])) - F[1]/(F[0]*(F[4] -
                F[1]*F[3]/F[0])))
        X[2] = (-(F[0]*F[4] - F[1]*F[3])*(-F[1]*(F[5] -
                F[2]*F[3]/F[0])/(F[0]*(F[4] - F[1]*F[3]/F[0])) +
                F[2]/F[0])/(F[0]*F[4]*F[8] - F[0]*F[5]*F[7] - F[1]*F[3]*F[8] +
                F[1]*F[5]*F[6] + F[2]*F[3]*F[7] - F[2]*F[4]*F[6]))
        X[3] = (-(F[5] - F[2]*F[3]/F[0])*(F[0]*F[4] - F[1]*F[3])*(F[3]*(F[7] -
                F[1]*F[6]/F[0])/(F[0]*(F[4] - F[1]*F[3]/F[0])) -
                F[6]/F[0])/((F[4] - F[1]*F[3]/F[0])*(F[0]*F[4]*F[8] -
                F[0]*F[5]*F[7] - F[1]*F[3]*F[8] + F[1]*F[5]*F[6] +
                F[2]*F[3]*F[7] - F[2]*F[4]*F[6])) - F[3]/(F[0]*(F[4] -
                F[1]*F[3]/F[0])))
        X[4] = (1/(F[4] - F[1]*F[3]/F[0]) + (F[5] - F[2]*F[3]/F[0])*(F[7] -
                F[1]*F[6]/F[0])*(F[0]*F[4] - F[1]*F[3])/((F[4] -
                F[1]*F[3]/F[0])**2*(F[0]*F[4]*F[8] - F[0]*F[5]*F[7] -
                F[1]*F[3]*F[8] + F[1]*F[5]*F[6] + F[2]*F[3]*F[7] -
                F[2]*F[4]*F[6])))
        X[5] = (-(F[5] - F[2]*F[3]/F[0])*(F[0]*F[4] - F[1]*F[3])/((F[4] -
                F[1]*F[3]/F[0])*(F[0]*F[4]*F[8] - F[0]*F[5]*F[7] -
                F[1]*F[3]*F[8] + F[1]*F[5]*F[6] + F[2]*F[3]*F[7] -
                F[2]*F[4]*F[6])))
        X[6] = ((F[0]*F[4] - F[1]*F[3])*(F[3]*(F[7] -
                F[1]*F[6]/F[0])/(F[0]*(F[4] - F[1]*F[3]/F[0])) -
                F[6]/F[0])/(F[0]*F[4]*F[8] - F[0]*F[5]*F[7] - F[1]*F[3]*F[8] +
                F[1]*F[5]*F[6] + F[2]*F[3]*F[7] - F[2]*F[4]*F[6]))
        X[7] = (-(F[7] - F[1]*F[6]/F[0])*(F[0]*F[4] - F[1]*F[3])/((F[4] -
                F[1]*F[3]/F[0])*(F[0]*F[4]*F[8] - F[0]*F[5]*F[7] -
                F[1]*F[3]*F[8] + F[1]*F[5]*F[6] + F[2]*F[3]*F[7] -
                F[2]*F[4]*F[6])))
        X[8] = ((F[0]*F[4] - F[1]*F[3])/(F[0]*F[4]*F[8] - F[0]*F[5]*F[7] -
                F[1]*F[3]*F[8] + F[1]*F[5]*F[6] + F[2]*F[3]*F[7] -
                F[2]*F[4]*F[6]))
    else:
        raise NotImplementedError
    return X

In [15]:
write('def dyad(A, B):')
write('    """Computes the outer product of A and B"""', wrapit=0)
write('    if A.size == 6 and B.size == 6:')
write('        X = zeros((6,6))')
for I in range(6):
    for J in range(6):
        ij = I_from_ij(*sorted(ij_from_I[I]))
        kl = I_from_ij(*sorted(ij_from_I[J]))
        write('        X[{},{}] = A[{}] * B[{}]'.format(I,J,ij,kl))
write('    elif A.size == 3 and B.size == 3:')
write('        X = zeros(6)')

for I in range(6):
    i,j = ij_from_I[I]
    write('        X[{}] = A[{}] * B[{}]'.format(I,i,j))
write('    else:\n        raise NotImplementedError')
write('    return X')


def dyad(A, B):
    """Computes the outer product of A and B"""
    if A.size == 6 and B.size == 6:
        X = zeros((6,6))
        X[0,0] = A[0] * B[0]
        X[0,1] = A[0] * B[1]
        X[0,2] = A[0] * B[2]
        X[0,3] = A[0] * B[3]
        X[0,4] = A[0] * B[4]
        X[0,5] = A[0] * B[5]
        X[1,0] = A[1] * B[0]
        X[1,1] = A[1] * B[1]
        X[1,2] = A[1] * B[2]
        X[1,3] = A[1] * B[3]
        X[1,4] = A[1] * B[4]
        X[1,5] = A[1] * B[5]
        X[2,0] = A[2] * B[0]
        X[2,1] = A[2] * B[1]
        X[2,2] = A[2] * B[2]
        X[2,3] = A[2] * B[3]
        X[2,4] = A[2] * B[4]
        X[2,5] = A[2] * B[5]
        X[3,0] = A[3] * B[0]
        X[3,1] = A[3] * B[1]
        X[3,2] = A[3] * B[2]
        X[3,3] = A[3] * B[3]
        X[3,4] = A[3] * B[4]
        X[3,5] = A[3] * B[5]
        X[4,0] = A[4] * B[0]
        X[4,1] = A[4] * B[1]
        X[4,2] = A[4] * B[2]
        X[4,3] = A[4] * B[3]
        X[4,4] = A[4] * B[4]
        X[4,5] = A[4] * B[5]
        X[5,0] = A[5] * B[0]
        X[5,1] = A[5] * B[1]
        X[5,2] = A[5] * B[2]
        X[5,3] = A[5] * B[3]
        X[5,4] = A[5] * B[4]
        X[5,5] = A[5] * B[5]
    elif A.size == 3 and B.size == 3:
        X = zeros(6)
        X[0] = A[0] * B[0]
        X[1] = A[1] * B[1]
        X[2] = A[2] * B[2]
        X[3] = A[0] * B[1]
        X[4] = A[1] * B[2]
        X[5] = A[0] * B[2]
    else:
        raise NotImplementedError
    return X

In [16]:
write('def symshuffle(A, B):')
write('    """ Computes the product Xijkl = .5 (Aik Bjl + Ail Bjk)"""', wrapit=0)
write('    X = zeros((6,6))')
write('    if A.size == 6 and B.size == 6:')
for I in range(6):
    for J in range(6):
        i,j = ij_from_I[I]
        k,l = ij_from_I[J]
        ik = I_from_ij(*sorted((i,k)))
        jl = I_from_ij(*sorted((j,l)))
        il = I_from_ij(*sorted((i,l)))
        jk = I_from_ij(*sorted((j,k)))
        write('        X[{},{}] = (A[{}] * B[{}] + A[{}] * B[{}]) / 2.'.format(I,J,ik,jl,il,jk))
write('    else:')
write('        raise NotImplementedError')
write('    return X')
write('odot = symshuffle')


def symshuffle(A, B):
    """ Computes the product Xijkl = .5 (Aik Bjl + Ail Bjk)"""
    X = zeros((6,6))
    if A.size == 6 and B.size == 6:
        X[0,0] = (A[0] * B[0] + A[0] * B[0]) / 2.
        X[0,1] = (A[3] * B[3] + A[3] * B[3]) / 2.
        X[0,2] = (A[5] * B[5] + A[5] * B[5]) / 2.
        X[0,3] = (A[0] * B[3] + A[3] * B[0]) / 2.
        X[0,4] = (A[3] * B[5] + A[5] * B[3]) / 2.
        X[0,5] = (A[0] * B[5] + A[5] * B[0]) / 2.
        X[1,0] = (A[3] * B[3] + A[3] * B[3]) / 2.
        X[1,1] = (A[1] * B[1] + A[1] * B[1]) / 2.
        X[1,2] = (A[4] * B[4] + A[4] * B[4]) / 2.
        X[1,3] = (A[3] * B[1] + A[1] * B[3]) / 2.
        X[1,4] = (A[1] * B[4] + A[4] * B[1]) / 2.
        X[1,5] = (A[3] * B[4] + A[4] * B[3]) / 2.
        X[2,0] = (A[5] * B[5] + A[5] * B[5]) / 2.
        X[2,1] = (A[4] * B[4] + A[4] * B[4]) / 2.
        X[2,2] = (A[2] * B[2] + A[2] * B[2]) / 2.
        X[2,3] = (A[5] * B[4] + A[4] * B[5]) / 2.
        X[2,4] = (A[4] * B[2] + A[2] * B[4]) / 2.
        X[2,5] = (A[5] * B[2] + A[2] * B[5]) / 2.
        X[3,0] = (A[0] * B[3] + A[0] * B[3]) / 2.
        X[3,1] = (A[3] * B[1] + A[3] * B[1]) / 2.
        X[3,2] = (A[5] * B[4] + A[5] * B[4]) / 2.
        X[3,3] = (A[0] * B[1] + A[3] * B[3]) / 2.
        X[3,4] = (A[3] * B[4] + A[5] * B[1]) / 2.
        X[3,5] = (A[0] * B[4] + A[5] * B[3]) / 2.
        X[4,0] = (A[3] * B[5] + A[3] * B[5]) / 2.
        X[4,1] = (A[1] * B[4] + A[1] * B[4]) / 2.
        X[4,2] = (A[4] * B[2] + A[4] * B[2]) / 2.
        X[4,3] = (A[3] * B[4] + A[1] * B[5]) / 2.
        X[4,4] = (A[1] * B[2] + A[4] * B[4]) / 2.
        X[4,5] = (A[3] * B[2] + A[4] * B[5]) / 2.
        X[5,0] = (A[0] * B[5] + A[0] * B[5]) / 2.
        X[5,1] = (A[3] * B[4] + A[3] * B[4]) / 2.
        X[5,2] = (A[5] * B[2] + A[5] * B[2]) / 2.
        X[5,3] = (A[0] * B[4] + A[3] * B[5]) / 2.
        X[5,4] = (A[3] * B[2] + A[5] * B[4]) / 2.
        X[5,5] = (A[0] * B[2] + A[5] * B[5]) / 2.
    else:
        raise NotImplementedError
    return X
odot = symshuffle

Push transformation

$$A_{ij} = F_{im}C_{mn}F_{jn} = L_{imjn}C_{mn}$$

where

$$ L_{imjn} = F_{im}F_{jn}$$

In [17]:
L = {}
for i in range(3):
    for j in range(3):
        p, q = sorted((i,j))
        I = I_from_ij(p,q)
        row = []
        for m in range(3):
            for n in range(3):
                p, q = sorted((m,n))
                J = I_from_ij(p,q)
                im = I_from_ij_9(i,m)
                jn = I_from_ij_9(j,n)
                L.setdefault((I,J), []).append('F[{}] * F[{}]'.format(im,jn))

for (k, v) in L.items():
    L[k] = ' + '.join(v[:2])

write('def symleaf(F):')
write('''\
    """ COMPUTE A 6X6 MANDEL MATRIX THAT IS THE SYM-LEAF TRANSFORMATION OF THE
    INPUT 3X3 MATRIX F.

    Parameters
    ----------
    F : ANY 3X3 MATRIX (IN CONVENTIONAL 3X3 STORAGE)

    Returns
    -------
    X : 6X6 MANDEL MATRIX FOR THE SYM-LEAF TRANSFORMATION MATRIX

    Notes
    -----
    IF A IS ANY SYMMETRIC TENSOR, AND IF {A} IS ITS 6X1 MANDEL ARRAY, THEN THE
    6X1 MANDEL ARRAY FOR THE TENSOR B=F.A.TRANSPOSE[F] MAY BE COMPUTED BY

                          {B}=[FF]{A}

    IF F IS A DEFORMATION F, THEN B IS THE "PUSH" (SPATIAL) TRANSFORMATION OF
    THE REFERENCE TENSOR A IF F IS Inverse[F], THEN B IS THE "PULL"
    (REFERENCE) TRANSFORMATION OF THE SPATIAL TENSOR A, AND THEREFORE B WOULD
    BE Inverse[FF]{A}.

    IF F IS A ROTATION, THEN B IS THE ROTATION OF A, AND
    FF WOULD BE BE A 6X6 ORTHOGONAL MATRIX, JUST AS IS F

    """''', wrapit=0)
write('    X = zeros((6,6))')
for i in range(6):
    for j in range(6):
        write('    X[{},{}] = {}'.format(i,j,L[i,j]))
write('    return X')


def symleaf(F):
    """ COMPUTE A 6X6 MANDEL MATRIX THAT IS THE SYM-LEAF TRANSFORMATION OF THE
    INPUT 3X3 MATRIX F.

    Parameters
    ----------
    F : ANY 3X3 MATRIX (IN CONVENTIONAL 3X3 STORAGE)

    Returns
    -------
    X : 6X6 MANDEL MATRIX FOR THE SYM-LEAF TRANSFORMATION MATRIX

    Notes
    -----
    IF A IS ANY SYMMETRIC TENSOR, AND IF {A} IS ITS 6X1 MANDEL ARRAY, THEN THE
    6X1 MANDEL ARRAY FOR THE TENSOR B=F.A.TRANSPOSE[F] MAY BE COMPUTED BY

                          {B}=[FF]{A}

    IF F IS A DEFORMATION F, THEN B IS THE "PUSH" (SPATIAL) TRANSFORMATION OF
    THE REFERENCE TENSOR A IF F IS Inverse[F], THEN B IS THE "PULL"
    (REFERENCE) TRANSFORMATION OF THE SPATIAL TENSOR A, AND THEREFORE B WOULD
    BE Inverse[FF]{A}.

    IF F IS A ROTATION, THEN B IS THE ROTATION OF A, AND
    FF WOULD BE BE A 6X6 ORTHOGONAL MATRIX, JUST AS IS F

    """
    X = zeros((6,6))
    X[0,0] = F[0] * F[0]
    X[0,1] = F[1] * F[1]
    X[0,2] = F[2] * F[2]
    X[0,3] = F[0] * F[1] + F[1] * F[0]
    X[0,4] = F[1] * F[2] + F[2] * F[1]
    X[0,5] = F[0] * F[2] + F[2] * F[0]
    X[1,0] = F[3] * F[3]
    X[1,1] = F[4] * F[4]
    X[1,2] = F[5] * F[5]
    X[1,3] = F[3] * F[4] + F[4] * F[3]
    X[1,4] = F[4] * F[5] + F[5] * F[4]
    X[1,5] = F[3] * F[5] + F[5] * F[3]
    X[2,0] = F[6] * F[6]
    X[2,1] = F[7] * F[7]
    X[2,2] = F[8] * F[8]
    X[2,3] = F[6] * F[7] + F[7] * F[6]
    X[2,4] = F[7] * F[8] + F[8] * F[7]
    X[2,5] = F[6] * F[8] + F[8] * F[6]
    X[3,0] = F[0] * F[3] + F[3] * F[0]
    X[3,1] = F[1] * F[4] + F[4] * F[1]
    X[3,2] = F[2] * F[5] + F[5] * F[2]
    X[3,3] = F[0] * F[4] + F[1] * F[3]
    X[3,4] = F[1] * F[5] + F[2] * F[4]
    X[3,5] = F[0] * F[5] + F[2] * F[3]
    X[4,0] = F[3] * F[6] + F[6] * F[3]
    X[4,1] = F[4] * F[7] + F[7] * F[4]
    X[4,2] = F[5] * F[8] + F[8] * F[5]
    X[4,3] = F[3] * F[7] + F[4] * F[6]
    X[4,4] = F[4] * F[8] + F[5] * F[7]
    X[4,5] = F[3] * F[8] + F[5] * F[6]
    X[5,0] = F[0] * F[6] + F[6] * F[0]
    X[5,1] = F[1] * F[7] + F[7] * F[1]
    X[5,2] = F[2] * F[8] + F[8] * F[2]
    X[5,3] = F[0] * F[7] + F[1] * F[6]
    X[5,4] = F[1] * F[8] + F[2] * F[7]
    X[5,5] = F[0] * F[8] + F[2] * F[6]
    return X

In [18]:
print('    @classmethod')
print('    def param_names(cls, n):')
N = 3
M = 3
ij = []
for n in range(N+2):
    ij.extend([(i,j) for i in range(n)[::-1] for j in range(n) if i+j==n-1])
ij = ij[1:] 
cij = ["'C{}{}'".format(i,j) for i,j in ij]
dij = ['D{0}'.format(i) for i in range(M)]
params = '        return [{}]'.format(','.join(cij + dij))
print(params)

NCI = [x[0] for x in ij]
NCJ = [x[1] for x in ij]
nc = len(NCI)
nd = len(dij)

u = ['u = zeros(2)']
for k in range(nc):
    i = NCI[k]
    j = NCJ[k]
    w = 'c[{0}] * ((IB1 - 3.) ** {1}) * ((IB2 - 3.) ** {2})'.format(k, i, j)
    u.append('u[1] += {0}'.format(w))
for k in range(nd):
    K = k + 1
    w = '1. / c[{0}] * (Jac - 1.) ** {1}'.format(nc+k, 2*K)
    u.append('if c[{0}] > 0.: u[0] += {1}'.format(nc+k, w))
u.append('u[0] = u[0] + u[1]')

du = []
for k in range(nc):
    i = NCI[k]
    j = NCJ[k]
    if i > 0:
        fac = i
        w1 = '(IB1 - 3.) ** {0}'.format(i-1)
        w2 = '(IB2 - 3.) ** {0}'.format(j)
        w = '{0:.1f} * c[{1}] * ({2}) * ({3})'.format(fac, k, w1, w2)
        du.append('du[0] += {0}'.format(w))

    if j > 0:
        fac = j
        w1 = '(IB1 - 3.) ** {0}'.format(i)
        w2 = '(IB2 - 3.) ** {0}'.format(j-1)
        w = '{0:.1f} * c[{1}] * ({2}) * ({3})'.format(fac, k, w1, w2)
        du.append('du[1] += {0}'.format(w))
du = ['du = zeros(3)'] + sorted(du, key=lambda x: int(x[3]))

for k in range(nd):
    K = k + 1
    fac = 2 * K
    if fac <= 0:
        continue
    w = '{0:.1f} / c[{1}] * (Jac - 1.) ** {2}'.format(fac, nc+k, 2*K-1)
    du.append('if c[{0}] > 0.: du[2] += {1}'.format(nc+k, w))

d2u = []
for k in range(nc):
    i = NCI[k]
    j = NCJ[k]

    fac = i * (i - 1)
    if fac > 0:
        w1 = '(IB1 - 3.) ** {0}'.format(i-2)
        w2 = '(IB2 - 3.) ** {0}'.format(j)
        w = '{0:.1f} * c[{1}] * ({2}) * ({3})'.format(fac, k, w1, w2)
        d2u.append('ddu[0] += {0}'.format(w))

    fac = j * (j - 1)
    if fac > 0:
        w1 = '(IB1 - 3.) ** {0}'.format(i)
        w2 = '(IB2 - 3.) ** {0}'.format(j-2)
        w = '{0:.1f} * c[{1}] * ({2}) * ({3})'.format(fac, k, w1, w2)
        d2u.append('ddu[1] += {0}'.format(w))

    fac = i * j
    if fac > 0:
        w1 = '(IB1 - 3.) ** {0}'.format(i-1)
        w2 = '(IB2 - 3.) ** {0}'.format(j-1)
        w = '{0:.1f} * c[{1}] * ({2}) * ({3})'.format(fac, k, w1, w2)
        d2u.append('ddu[3] += {0}'.format(w))

d2u = ['ddu = zeros(6)'] + sorted(d2u, key=lambda x: int(x[4]))
for k in range(nd):
    K = k + 1
    fac = 2 * K * (2 * K - 1)
    if fac <= 0:
        continue
    w = '{0:.1f} / c[{1}] * (Jac - 1.) ** {2}'.format(fac, nc+k, 2*K-2)
    d2u.append('if c[{0}] > 0.: ddu[2] += {1}'.format(nc+k, w))

d3u = ['dddu = zeros(6)']
for k in range(nd):
    K = k + 1
    fac = 4 * K * (2 * K - 1) * (K - 1)
    if fac <= 0:
        continue
    w = '{0:.1f} / c[{1}] * (Jac - 1.) ** {2}'.format(fac, nc+k, 2*K-3)
    d3u.append('if c[{0}] > 0.: dddu[5] += {1}'.format(nc+k, w))

print()
print('def uhyper(c, IB1, IB2, Jac):')
print('    """Computes the energies and derivatives of the isotropic')
print('       hyperelastic energy"""\n')
print('       THIS FUNCTION WAS GENERATED AUTOMATICALLY, DO NOT EDIT IT DIRECTLY')
print('    """')

s = 'Energies'
print('    # {0}'.format(s))
print('    {0}'.format('\n    '.join(u)))

s = 'First Derivatives'
print('\n    # {0}'.format(s))
print('    {0}'.format('\n    '.join(du)))

s = 'Second Derivatives'
print('\n    # {0}'.format(s))
print('    {0}'.format('\n    '.join(d2u)))

s = 'Third Derivatives'
print('\n    # {0}'.format(s))
print('    {0}'.format('\n    '.join(d3u)))

print('    return u, du, ddu, dddu')


    @classmethod
    def param_names(cls, n):
        return ['C10','C01','C20','C11','C02','C30','C21','C12','C03',D0,D1,D2]

def uhyper(c, IB1, IB2, Jac):
    """Computes the energies and derivatives of the isotropic
       hyperelastic energy"""

       THIS FUNCTION WAS GENERATED AUTOMATICALLY, DO NOT EDIT IT DIRECTLY
    """
    # Energies
    u = zeros(2)
    u[1] += c[0] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 0)
    u[1] += c[1] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 1)
    u[1] += c[2] * ((IB1 - 3.) ** 2) * ((IB2 - 3.) ** 0)
    u[1] += c[3] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 1)
    u[1] += c[4] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 2)
    u[1] += c[5] * ((IB1 - 3.) ** 3) * ((IB2 - 3.) ** 0)
    u[1] += c[6] * ((IB1 - 3.) ** 2) * ((IB2 - 3.) ** 1)
    u[1] += c[7] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 2)
    u[1] += c[8] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 3)
    if c[9] > 0.: u[0] += 1. / c[9] * (Jac - 1.) ** 2
    if c[10] > 0.: u[0] += 1. / c[10] * (Jac - 1.) ** 4
    if c[11] > 0.: u[0] += 1. / c[11] * (Jac - 1.) ** 6
    u[0] = u[0] + u[1]

    # First Derivatives
    du = zeros(3)
    du[0] += 1.0 * c[0] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 0)
    du[0] += 2.0 * c[2] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 0)
    du[0] += 1.0 * c[3] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 1)
    du[0] += 3.0 * c[5] * ((IB1 - 3.) ** 2) * ((IB2 - 3.) ** 0)
    du[0] += 2.0 * c[6] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 1)
    du[0] += 1.0 * c[7] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 2)
    du[1] += 1.0 * c[1] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 0)
    du[1] += 1.0 * c[3] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 0)
    du[1] += 2.0 * c[4] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 1)
    du[1] += 1.0 * c[6] * ((IB1 - 3.) ** 2) * ((IB2 - 3.) ** 0)
    du[1] += 2.0 * c[7] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 1)
    du[1] += 3.0 * c[8] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 2)
    if c[9] > 0.: du[2] += 2.0 / c[9] * (Jac - 1.) ** 1
    if c[10] > 0.: du[2] += 4.0 / c[10] * (Jac - 1.) ** 3
    if c[11] > 0.: du[2] += 6.0 / c[11] * (Jac - 1.) ** 5

    # Second Derivatives
    ddu = zeros(6)
    ddu[0] += 2.0 * c[2] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 0)
    ddu[0] += 6.0 * c[5] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 0)
    ddu[0] += 2.0 * c[6] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 1)
    ddu[1] += 2.0 * c[4] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 0)
    ddu[1] += 2.0 * c[7] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 0)
    ddu[1] += 6.0 * c[8] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 1)
    ddu[3] += 1.0 * c[3] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 0)
    ddu[3] += 2.0 * c[6] * ((IB1 - 3.) ** 1) * ((IB2 - 3.) ** 0)
    ddu[3] += 2.0 * c[7] * ((IB1 - 3.) ** 0) * ((IB2 - 3.) ** 1)
    if c[9] > 0.: ddu[2] += 2.0 / c[9] * (Jac - 1.) ** 0
    if c[10] > 0.: ddu[2] += 12.0 / c[10] * (Jac - 1.) ** 2
    if c[11] > 0.: ddu[2] += 30.0 / c[11] * (Jac - 1.) ** 4

    # Third Derivatives
    dddu = zeros(6)
    if c[10] > 0.: dddu[5] += 24.0 / c[10] * (Jac - 1.) ** 1
    if c[11] > 0.: dddu[5] += 120.0 / c[11] * (Jac - 1.) ** 3
    return u, du, ddu, dddu
$$\begin{align} \dot{Y} &= Y_1 \dot{\epsilon}^p_{eq} \\ &= Y_1\sqrt{\frac{2}{3}}\lVert{\rm dev}\left(\dot{\pmb{\epsilon}}^p\right)\rVert \\ &= Y_1\sqrt{\frac{2}{3}}\sqrt{\pmb{m}{:}\pmb{m}}\dot{\lambda} \\ &= Y_1\sqrt{\frac{2}{3}}\sqrt{\frac{1}{\sqrt{2J_2}}\frac{1}{\sqrt{2J_2}}\pmb{s}{:}\pmb{s}}\dot{\lambda}\\ &= Y_1\sqrt{\frac{2}{3}}\frac{1}{\sqrt{2J_2}}\sqrt{\pmb{s}{:}\pmb{s}}\dot{\lambda} \\ &= Y_1\frac{1}{\sqrt{3J_2}}\sqrt{\pmb{s}{:}\pmb{s}}\dot{\lambda} \\ &= Y_1\frac{1}{\sqrt{3/2}\sqrt{2J_2}}\sqrt{\pmb{s}{:}\pmb{s}}\dot{\lambda}\\ &= Y_1\sqrt{\frac{2}{3}}\dot{\lambda}\\ &= H_Y\dot{\lambda} \end{align}$$

In [ ]: