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)
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)')
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')
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')
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')
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)
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')
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)')
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')
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')
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')
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')
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')
In [ ]: