223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

add the calculation of Jacobian of isotropic and anisotropic Barlat1991, shows better performance.

This commit is contained in:
Haiming Zhang 2015-02-12 11:04:00 +00:00
parent e4874550dc
commit 69af205721
1 changed files with 60 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

86
processing/misc/yieldSurface.py
View File

@ -172,7 +172,7 @@ class Barlat1991iso(object):
def fun(self, (sigma0, m), ydata, sigmas): def fun(self, (sigma0, m), ydata, sigmas):
return Barlat1991Basis(sigma0, 1.0,1.0,1.0,1.0,1.0,1.0, m, sigmas) return Barlat1991Basis(sigma0, 1.0,1.0,1.0,1.0,1.0,1.0, m, sigmas)
def jac(self, (sigma0, m), ydata, sigmas): def jac(self, (sigma0, m), ydata, sigmas):
pass return Barlat1991Basis(sigma0, 1.0,1.0,1.0,1.0,1.0,1.0, m, sigmas, Jac=True, nParas=2)
class Barlat1991aniso(object): class Barlat1991aniso(object):
''' '''
@ -181,7 +181,7 @@ class Barlat1991aniso(object):
def fun(self, (sigma0, a,b,c,f,g,h, m), ydata, sigmas): def fun(self, (sigma0, a,b,c,f,g,h, m), ydata, sigmas):
return Barlat1991Basis(sigma0, a,b,c,f,g,h, m, sigmas) return Barlat1991Basis(sigma0, a,b,c,f,g,h, m, sigmas)
def jac(self, (sigma0, a,b,c,f,g,h, m), ydata, sigmas): def jac(self, (sigma0, a,b,c,f,g,h, m), ydata, sigmas):
pass return Barlat1991Basis(sigma0, a,b,c,f,g,h, m, sigmas, Jac=True, nParas=8)
def Cazacu_Barlat3D(sigma0,a1,a2,a3,a4,a5,a6, b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11, c, def Cazacu_Barlat3D(sigma0,a1,a2,a3,a4,a5,a6, b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11, c,
ydata, sigmas): ydata, sigmas):
@ -303,37 +303,71 @@ def HosfordBasis(sigma0, F,G,H, a, sigmas, Jac=False, nParas=1):
for a,b,c,d,e in zip(j1, j2,j3,j4,j5): jaco.append([a,b,c,d,e]) for a,b,c,d,e in zip(j1, j2,j3,j4,j5): jaco.append([a,b,c,d,e])
return np.array(jaco) return np.array(jaco)
def Barlat1991Basis(sigma0, a, b, c, f, g, h, order, sigmas): def Barlat1991Basis(sigma0, a, b, c, f, g, h, m, sigmas, Jac=False, nParas=2):
''' '''
residuum of Barlat 1997 yield criterion residuum of Barlat 1997 yield criterion
''' '''
cos = np.cos; pi = np.pi; abs = np.abs cos = np.cos; sin = np.sin; pi = np.pi; abs = np.abs
A = a*(sigmas[1] - sigmas[2]) dAda = sigmas[1] - sigmas[2]; A = a*dAda
B = b*(sigmas[2] - sigmas[0]) dBdb = sigmas[2] - sigmas[0]; B = b*dBdb
C = c*(sigmas[0] - sigmas[1]) dCdc = sigmas[0] - sigmas[1]; C = c*dCdc
F = f* sigmas[4] dFdf = sigmas[4]; F = f*dFdf
G = g* sigmas[5] dGdg = sigmas[5]; G = g*dGdg
H = h* sigmas[3] dHdh = sigmas[3]; H = h*dHdh
I2 = (F*F + G*G + H*H)/3.0 + ((A-C)**2+(C-B)**2+(B-A)**2)/54.0 I2 = (F*F + G*G + H*H)/3.0 + ((A-C)**2+(C-B)**2+(B-A)**2)/54.0
I3 = (C-B)*(A-C) * (B-A)/54.0 + F*G*H - \ I3 = (C-B)*(A-C)*(B-A)/54.0 + F*G*H - \
( (C-B)*F*F + (A-C)*G*G + (B-A)*H*H )/6.0 ( (C-B)*F*F + (A-C)*G*G + (B-A)*H*H )/6.0
theta = np.arccos(I3/I2**1.5) theta = np.arccos(I3/I2**1.5)
Phi = np.sqrt(3.0*I2)* ( phi1 = (2.0*theta + pi)/6.0; cos1 = 2.0*cos(phi1); absc1 = abs(cos1)
(abs(2.0*cos((2.0*theta + pi)/6.0)))**order + phi2 = (2.0*theta + pi*3.0)/6.0; cos2 = 2.0*cos(phi2); absc2 = abs(cos2)
(abs(2.0*cos((2.0*theta + pi*3.0)/6.0)))**order + phi3 = (2.0*theta + pi*5.0)/6.0; cos3 = 2.0*cos(phi3); absc3 = abs(cos3)
(abs(2.0*cos(( 2.0*theta + pi*5.0)/6.0)))**order ratio= np.sqrt(3.0*I2)/sigma0; expo = 1.0/m
)**(1.0/order) left = ( absc1**m + absc2**m + absc3**m )/2.0
# r = Phi/2.0**(1.0/order) - sigma0 leftNorm = left**expo
r = Phi/2.0**(1.0/order)/sigma0 - 1.0 r = ratio*leftNorm - 1.0
# Phi = (3.0*I2)**(order/2.0) * (
# (abs(2.0*cos((2.0*theta + pi)/6.0))) **order + if not Jac:
# (abs(2.0*cos((2.0*theta + pi*3.0)/6.0)))**order + return r.ravel()
# (abs(2.0*cos((2.0*theta + pi*5.0)/6.0)))**order else:
# ) ln = lambda x : np.log(x + 1.0e-32)
# r = (Phi - 2.0*sigma0**order)**(1.0/order) jaco = []
dfdl = expo*leftNorm/left
js = -(r + 1.0)/sigma0
jm = (r+1.0)*ln(left)*(-expo*expo) + ratio*dfdl*0.5*(
absc1**m*ln(absc1) + absc2**m*ln(absc2) + absc3**m*ln(absc3) )
if nParas == 2:
for j1,j2 in zip(js, jm): jaco.append([j1,j2])
return np.array(jaco)
else:
dI2da = (2.0*A-B-C)*dAda/27.0; dI2df = 2.0*F*dFdf/3.0
dI2db = (2.0*B-C-A)*dBdb/27.0; dI2dg = 2.0*G*dGdg/3.0
dI2dc = (2.0*C-A-B)*dCdc/27.0; dI2dh = 2.0*H*dHdh/3.0
dI3da = dI2da*(B-C)/2.0 + (H**2 - G**2)*dAda/6.0
dI3db = dI2db*(C-A)/2.0 + (F**2 - H**2)*dBdb/6.0
dI3dc = dI2dc*(A-B)/2.0 + (G**2 - F**2)*dCdc/6.0
dI3df = ( (H*G + (B-C)) * F/3.0 )*dFdf
dI3dg = ( (F*H + (C-A)) * G/3.0 )*dGdg
dI3dh = ( (G*F + (A-B)) * H/3.0 )*dHdh
darccos = -(1.0 - I3**2/I2**3)**(-0.5)
dthedI2 = darccos*I3*(-1.5)*I2**(-2.5)
dthedI3 = darccos*I2**(-1.5)
dc1dthe = -sin(phi1)/3.0; dc2dthe = -sin(phi2)/3.0; dc3dthe = -sin(phi3)/3.0
dfdc = ratio * dfdl * 0.5 * m
dfdc1 = dfdc* absc1**(expo-1.0)*np.sign(cos1)
dfdc2 = dfdc* absc2**(expo-1.0)*np.sign(cos2)
dfdc3 = dfdc* absc3**(expo-1.0)*np.sign(cos3)
dfdthe= (dfdc1*dc1dthe + dfdc2*dc2dthe + dfdc2*dc2dthe)*2.0
dfdI2 = dfdthe*dthedI2; dfdI3 = dfdthe*dthedI3
ja = dfdI2*dI2da + dfdI3*dI3da; jf = dfdI2*dI2df + dfdI3*dI3df
jb = dfdI2*dI2db + dfdI3*dI3db; jg = dfdI2*dI2dg + dfdI3*dI3dg
jc = dfdI2*dI2dc + dfdI3*dI3dc; jh = dfdI2*dI2dh + dfdI3*dI3dh
for j1,j2,j3,j4,j5,j6,j7,j8 in zip(js,ja,jb,jc,jf,jg,jh,jm):
jaco.append([j1,j2,j3,j4,j5,j6,j7,j8])
return np.array(jaco)
return r.ravel()
fittingCriteria = { fittingCriteria = {
'tresca' :{'func' : Tresca, 'tresca' :{'func' : Tresca,