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

trying to simplify

This commit is contained in:
Martin Diehl 2015-10-13 21:00:12 +00:00
parent 3a998ba114
commit 180d4625c1
1 changed files with 213 additions and 225 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

438
processing/misc/yieldSurface.py
View File

@ -1,9 +1,8 @@
#!/usr/bin/python #!/usr/bin/python
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import threading,time,os,subprocess,shlex,string import threading,time,os,subprocess,string,sys
import numpy as np import numpy as np
from scipy.linalg import svd
from optparse import OptionParser from optparse import OptionParser
import damask import damask
from damask.util import leastsqBound from damask.util import leastsqBound
@ -12,21 +11,23 @@ scriptID = string.replace('$Id$','\n','\\n')
scriptName = os.path.splitext(scriptID.split()[1])[0] scriptName = os.path.splitext(scriptID.split()[1])[0]
def runFit(exponent, eqStress, dimension, criterion): def runFit(exponent, eqStress, dimension, criterion):
global s, threads, myFit global s, threads, myFit, myLoad
global fitResults, fitErrors, fitResidual, stressAll, strainAll global fitResults, fitErrors, fitResidual, stressAll, strainAll
global N_simulations, Guess, dDim, numParas global N_simulations, Guess, dDim
fitResults = []; fitErrors = []; fitResidual = []; Guess = []; threads=[] fitResults = []; fitErrors = []; fitResidual = []; Guess = []; threads=[]
dDim = dimension - 3 dDim = dimension - 3
numParas = len(fitCriteria[criterion]['bound'][dDim]) Nparas = len(fitCriteria[criterion]['bound'][dDim])
nExpo = fitCriteria[criterion]['nExpo'] nExpo = fitCriteria[criterion]['nExpo']
if exponent > 0.0: if exponent > 0.0: # User defined exponents
numParas = numParas-nExpo # User defines the exponents Nparas = Nparas-nExpo
fitCriteria[criterion]['bound'][dDim] = fitCriteria[criterion]['bound'][dDim][:numParas] fitCriteria[criterion]['bound'][dDim] = fitCriteria[criterion]['bound'][dDim][:Nparas]
for i in xrange(numParas):
for i in xrange(Nparas):
temp = fitCriteria[criterion]['bound'][dDim][i] temp = fitCriteria[criterion]['bound'][dDim][i]
if fitCriteria[criterion]['bound'][dDim][i] == (None,None): Guess.append(1.0) if fitCriteria[criterion]['bound'][dDim][i] == (None,None):
Guess.append(1.0)
else: else:
g = (temp[0]+temp[1])/2.0 g = (temp[0]+temp[1])/2.0
if g == 0: g = temp[1]*0.5 if g == 0: g = temp[1]*0.5
@ -36,17 +37,19 @@ def runFit(exponent, eqStress, dimension, criterion):
s=threading.Semaphore(1) s=threading.Semaphore(1)
myLoad = Loadcase(options.load[0],options.load[1],options.load[2], myLoad = Loadcase(options.load[0],options.load[1],options.load[2],
nSet = 10, dimension = dimension, vegter = options.vegter) nSet = 10, dimension = dimension, vegter = options.vegter)
stressAll= [np.zeros(0,'d').reshape(0,0) for i in xrange(int(options.yieldValue[2]))] stressAll= [np.zeros(0,'d').reshape(0,0) for i in xrange(int(options.yieldValue[2]))]
strainAll= [np.zeros(0,'d').reshape(0,0) for i in xrange(int(options.yieldValue[2]))] strainAll= [np.zeros(0,'d').reshape(0,0) for i in xrange(int(options.yieldValue[2]))]
myFit = Criterion(exponent,eqStress, dimension, criterion) myFit = Criterion(exponent,eqStress, dimension, criterion)
for i in range(options.threads): for t in range(options.threads):
threads.append(myThread(i)) threads.append(myThread(t))
threads[i].start() threads[t].start()
for i in range(options.threads): for t in range(options.threads):
threads[i].join() threads[t].join()
print fitResidual damask.util.croak('Residuals')
damask.util.croak(fitResidual)
def principalStresses(sigmas): def principalStresses(sigmas):
''' '''
@ -61,32 +64,33 @@ def principalStresses(sigmas):
return lambdas return lambdas
def principalStress(p): def principalStress(p):
sin = np.sin; cos = np.cos I = invariant(p)
I1,I2,I3 = invariant(p)
I1s3I2= (I1**2 - 3.0*I2)**0.5 I1s3I2= (I[0]**2 - 3.0*I[1])**0.5
numer = 2.0*I1**3 - 9.0*I1*I2 + 27.0*I3 numer = 2.0*I[0]**3 - 9.0*I[0]*I[1] + 27.0*I[2]
denom = 2.0*I1s3I2**3 denom = 2.0*I1s3I2**3
cs = numer/denom cs = numer/denom
phi = np.arccos(cs)/3.0 phi = np.arccos(cs)/3.0
t1 = I1/3.0; t2 = 2.0/3.0*I1s3I2 t1 = I[0]/3.0; t2 = 2.0/3.0*I1s3I2
return np.array( [t1 + t2*cos(phi), t1+t2*cos(phi+np.pi*2.0/3.0), t1+t2*cos(phi+np.pi*4.0/3.0)]) return np.array( [t1 + t2*np.cos(phi),
t1 + t2*np.cos(phi+np.pi*2.0/3.0),
t1 + t2*np.cos(phi+np.pi*4.0/3.0)])
def principalStrs_Der(p, (s1, s2, s3, s4, s5, s6), dim, Karafillis=False): def principalStrs_Der(p, (s1, s2, s3, s4, s5, s6), dim, Karafillis=False):
''' '''
The derivative of principal stress with respect to stress Derivative of principal stress with respect to stress
''' '''
sin = np.sin; cos = np.cos third = 1.0/3.0
I1,I2,I3 = invariant(p) third2 = 2.0*third
third = 1.0/3.0 I = invariant(p)
I1s3I2= np.sqrt(I1**2 - 3.0*I2) I1s3I2= np.sqrt(I[0]**2 - 3.0*I[1])
numer, denom = 2.0*I1**3 - 9.0*I1*I2 + 27.0*I3, 2.0*I1s3I2**3 numer = 2.0*I1**3 - 9.0*I[0]*I[1] + 27.0*I[2]
denom = 2.0*I1s3I2**3
cs = numer/denom cs = numer/denom
phi = np.arccos(cs)/3.0 phi = np.arccos(cs)/3.0
dphidcs = -third/np.sqrt(1.0 - cs**2) dphidcs = -third/np.sqrt(1.0 - cs**2)
dcsddenom = 0.5*numer*(-1.5)*I1s3I2**(-5.0) dcsddenom = 0.5*numer*(-1.5)*I1s3I2**(-5.0)
dcsdI1 = (6.0*I1**2 - 9.0*I2)*denom + dcsddenom*(2.0*I1) dcsdI1 = (6.0*I1**2 - 9.0*I2)*denom + dcsddenom*(2.0*I1)
@ -94,12 +98,13 @@ def principalStrs_Der(p, (s1, s2, s3, s4, s5, s6), dim, Karafillis=False):
dcsdI3 = 27.0*denom dcsdI3 = 27.0*denom
dphidI1, dphidI2, dphidI3 = dphidcs*dcsdI1, dphidcs*dcsdI2, dphidcs*dcsdI3 dphidI1, dphidI2, dphidI3 = dphidcs*dcsdI1, dphidcs*dcsdI2, dphidcs*dcsdI3
dI1s3I2dI1= I1/I1s3I2; dI1s3I2dI2 = -1.5/I1s3I2 dI1s3I2dI1 = I1/I1s3I2
third2 = 2.0*third; tcoeff = third2*I1s3I2 dI1s3I2dI2 = -1.5/I1s3I2
tcoeff = third2*I1s3I2
dSidIj = lambda theta : ( tcoeff*(-sin(theta))*dphidI1 + third2*dI1s3I2dI1*cos(theta) + third, dSidIj = lambda theta : ( tcoeff*(-np.sin(theta))*dphidI1 + third2*dI1s3I2dI1*np.cos(theta) + third,
tcoeff*(-sin(theta))*dphidI2 + third2*dI1s3I2dI2*cos(theta), tcoeff*(-np.sin(theta))*dphidI2 + third2*dI1s3I2dI2*np.cos(theta),
tcoeff*(-sin(theta))*dphidI3) tcoeff*(-np.sin(theta))*dphidI3)
dSdI = np.array([dSidIj(phi),dSidIj(phi+np.pi*2.0/3.0),dSidIj(phi+np.pi*4.0/3.0)]) # i=1,2,3; j=1,2,3 dSdI = np.array([dSidIj(phi),dSidIj(phi+np.pi*2.0/3.0),dSidIj(phi+np.pi*4.0/3.0)]) # i=1,2,3; j=1,2,3
# calculate the derivation of principal stress with regards to the anisotropic coefficients # calculate the derivation of principal stress with regards to the anisotropic coefficients
@ -111,8 +116,11 @@ def principalStrs_Der(p, (s1, s2, s3, s4, s5, s6), dim, Karafillis=False):
if Karafillis: if Karafillis:
dpdc = np.array([[zero,s1-s3,s1-s2], [s2-s3,zero,s2-s1], [s3-s2,s3-s1,zero]])/3.0 dpdc = np.array([[zero,s1-s3,s1-s2], [s2-s3,zero,s2-s1], [s3-s2,s3-s1,zero]])/3.0
dSdp = np.array([np.dot(dSdI[:,:,i],dIdp[:,:,i]).T for i in xrange(num)]).T dSdp = np.array([np.dot(dSdI[:,:,i],dIdp[:,:,i]).T for i in xrange(num)]).T
if dim == 2: temp = np.vstack([dSdp[:,3]*s4]).T.reshape(num,1,3).T if dim == 2:
else: temp = np.vstack([dSdp[:,3]*s4,dSdp[:,4]*s5,dSdp[:,5]*s6]).T.reshape(num,3,3).T temp = np.vstack([dSdp[:,3]*s4]).T.reshape(num,1,3).T
else:
temp = np.vstack([dSdp[:,3]*s4,dSdp[:,4]*s5,dSdp[:,5]*s6]).T.reshape(num,3,3).T
return np.concatenate((np.array([np.dot(dSdp[:,0:3,i], dpdc[:,:,i]).T for i in xrange(num)]).T, return np.concatenate((np.array([np.dot(dSdp[:,0:3,i], dpdc[:,:,i]).T for i in xrange(num)]).T,
temp), axis=1) temp), axis=1)
else: else:
@ -127,14 +135,12 @@ def principalStrs_Der(p, (s1, s2, s3, s4, s5, s6), dim, Karafillis=False):
return np.array([np.dot(dSdI[:,:,i],dIdc[:,:,i]).T for i in xrange(num)]).T return np.array([np.dot(dSdI[:,:,i],dIdc[:,:,i]).T for i in xrange(num)]).T
def invariant(sigmas): def invariant(sigmas):
I=np.zeros(3)
s11,s22,s33,s12,s23,s31 = sigmas s11,s22,s33,s12,s23,s31 = sigmas
I1 = s11 + s22 + s33 I[0] = s11 + s22 + s33
I2 = s11*s22 + s22*s33 + s33*s11 - s12**2 - s23**2 - s31**2 I[1] = s11*s22 + s22*s33 + s33*s11 - s12**2 - s23**2 - s31**2
I3 = s11*s22*s33 + 2.0*s12*s23*s31 - s12**2*s33 - s23**2*s11 - s31**2*s22 I[2] = s11*s22*s33 + 2.0*s12*s23*s31 - s12**2*s33 - s23**2*s11 - s31**2*s22
return (I1,I2,I3) return I
def formatOutput(n, type='%-14.6f'):
return ''.join([type for i in xrange(n)])
def math_ln(x): def math_ln(x):
return np.log(x + 1.0e-32) return np.log(x + 1.0e-32)
@ -155,9 +161,9 @@ class Criteria(object):
''' '''
def __init__(self, criterion, uniaxialStress,exponent, dimension): def __init__(self, criterion, uniaxialStress,exponent, dimension):
self.stress0 = uniaxialStress self.stress0 = uniaxialStress
if exponent < 0.0: # The exponent m is undetermined if exponent < 0.0: # Fitting exponent m
self.mFix = [False, exponent] self.mFix = [False, exponent]
else: # The exponent m is fixed else: # fixed exponent m
self.mFix = [True, exponent] self.mFix = [True, exponent]
self.func = fitCriteria[criterion]['func'] self.func = fitCriteria[criterion]['func']
self.criteria = criterion self.criteria = criterion
@ -234,7 +240,7 @@ def VetgerCriterion(stress,lankford, rhoBi0, theta=0.0):
nps = len(stress) nps = len(stress)
if nps%4 != 0: if nps%4 != 0:
print ('Warning: the number of stress points is uncorrect, stress points of %s are missing in set %i'%( damask.util.croak('Warning: the number of stress points is uncorrect, stress points of %s are missing in set %i'%(
['eq-biaxial, plane strain & uniaxial', 'eq-biaxial & plane strain','eq-biaxial'][nps%4-1],nps/4+1)) ['eq-biaxial, plane strain & uniaxial', 'eq-biaxial & plane strain','eq-biaxial'][nps%4-1],nps/4+1))
else: else:
nset = nps/4 nset = nps/4
@ -338,15 +344,12 @@ def Drucker(eqStress, paras, sigmas, mFix, criteria, dim, Jac = False):
sigma0, C_D = paras[0:2] sigma0, C_D = paras[0:2]
if mFix[0]: p = mFix[1] if mFix[0]: p = mFix[1]
else: p = paras[-1] else: p = paras[-1]
I1,I2,I3 = invariant(sigmas) I = invariant(sigmas)
#I = invariant(sigmas) J = np.zeros([3])
#J = np.zeros([3]) J[1] = I[0]**2/3.0 - I[1]
J2 = I1**2/3.0 - I2 J[2] = I[0]**3/13.5 - I[0]*I[1]/3.0 + I[2]
#J[1] = I[0]**2/3.0 - I[1] J2_3p = J[1]**(3.0*p)
J3 = I1**3/13.5 - I1*I2/3.0 + I3 J3_2p = J[2]**(2.0*p)
#J[2] = I[0]**3/13.5 - I[0]*I[1]/3.0 + I[2] etc.
J2_3p = J2**(3.0*p)
J3_2p = J3**(2.0*p)
left = J2_3p - C_D*J3_2p left = J2_3p - C_D*J3_2p
r = left**(1.0/(6.0*p))*3.0**0.5/sigma0 r = left**(1.0/(6.0*p))*3.0**0.5/sigma0
@ -501,7 +504,6 @@ def Barlat1991(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False):
if mFix[0]: m = mFix[1] if mFix[0]: m = mFix[1]
else: m = paras[-1] else: m = paras[-1]
cos = np.cos; sin = np.sin; pi = np.pi; abs = np.abs
s11,s22,s33,s12,s23,s31 = sigmas s11,s22,s33,s12,s23,s31 = sigmas
if dim == 2: if dim == 2:
dXdx = np.array([s22,-s11,s11-s22,s12]) dXdx = np.array([s22,-s11,s11-s22,s12])
@ -512,9 +514,9 @@ def Barlat1991(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False):
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 - ((C-B)*F*F + (A-C)*G*G + (B-A)*H*H)/6.0 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
phi1 = np.arccos(I3/I2**1.5)/3.0 + pi/6.0; absc1 = 2.0*abs(cos(phi1)) phi1 = np.arccos(I3/I2**1.5)/3.0 + np.pi/6.0; absc1 = 2.0*np.abs(np.cos(phi1))
phi2 = phi1 + pi/3.0; absc2 = 2.0*abs(cos(phi2)) phi2 = phi1 + np.pi/3.0; absc2 = 2.0*np.abs(np.cos(phi2))
phi3 = phi2 + pi/3.0; absc3 = 2.0*abs(cos(phi3)) phi3 = phi2 + np.pi/3.0; absc3 = 2.0*np.abs(np.cos(phi3))
left = ( absc1**m + absc2**m + absc3**m ) left = ( absc1**m + absc2**m + absc3**m )
r = (0.5*left)**(1.0/m)*np.sqrt(3.0*I2)/eqStress r = (0.5*left)**(1.0/m)*np.sqrt(3.0*I2)/eqStress
@ -542,7 +544,7 @@ def Barlat1991(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False):
(G*F*3.0 + (A-B))*H ])/3.0*dXdx (G*F*3.0 + (A-B))*H ])/3.0*dXdx
darccos = -1.0/np.sqrt(1.0 - I3**2/I2**3) darccos = -1.0/np.sqrt(1.0 - I3**2/I2**3)
dfdcos = lambda phi : dfdl*m*(2.0*abs(cos(phi)))**(m-1.0)*np.sign(cos(phi))*(-sin(phi)/1.5) dfdcos = lambda phi : dfdl*m*(2.0*abs(np.cos(phi)))**(m-1.0)*np.sign(np.cos(phi))*(-np.sin(phi)/1.5)
dfdthe= (dfdcos(phi1) + dfdcos(phi2) + dfdcos(phi3)) dfdthe= (dfdcos(phi1) + dfdcos(phi2) + dfdcos(phi3))
dfdI2, dfdI3 = dfdthe*darccos*I3*(-1.5)*I2**(-2.5)+r/2.0/I2, dfdthe*darccos*I2**(-1.5) dfdI2, dfdI3 = dfdthe*darccos*I3*(-1.5)*I2**(-2.5)+r/2.0/I2, dfdthe*darccos*I2**(-1.5)
@ -842,150 +844,132 @@ def KarafillisBoyce(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False):
fitCriteria = { fitCriteria = {
'tresca' :{'func' : Tresca, 'tresca' :{'name': 'Tresca',
'nExpo': 0,'err':np.inf, 'func': Tresca,
'dimen': 3, 'nExpo': 0,'err':np.inf,
'bound': [ [(None,None)] ], 'dimen': 3,
'paras': [ 'sigma0' ], 'bound': [[(None,None)]],
'text' : '\nCoefficient of Tresca criterion: ', 'labels': [['sigma0']],
'error': 'The standard deviation error is: '
}, },
'vonmises' :{'func' : Hosford, 'vonmises' :{'name': 'Huber-Mises-Hencky',
'nExpo': 0,'err':np.inf, 'func' : Hosford,
'dimen': 3, 'nExpo': 0,'err':np.inf,
'bound': [ [(None,None)] ], 'dimen': 3,
'paras': [ 'sigma0' ], 'bound': [[(None,None)]],
'text' : '\nCoefficient of Huber-Mises-Hencky criterion: ', 'labels': [['sigma0']],
'error': 'The standard deviation error is: '
}, },
'hershey' :{'func' : Hosford, 'hershey' :{'name': 'Hershey',
'nExpo': 1,'err':np.inf, 'func': Hosford,
'dimen': 3, 'nExpo': 1,'err':np.inf,
'bound': [ [(None,None)]+[(1.0,8.0)] ], 'dimen': 3,
'paras': [ 'sigma0, a' ], 'bound': [[(None,None)]+[(1.0,8.0)]],
'text' : '\nCoefficients of Hershey criterion: ', 'labels': [['sigma0','a']],
'error': 'The standard deviation errors are: '
}, },
'ghosford' :{'func' : Hosford, 'hosford' :{'name': 'General Hosford',
'nExpo': 1,'err':np.inf, 'func': Hosford,
'dimen': 3, 'nExpo': 1,'err':np.inf,
'bound': [ [(0.0,2.0)]*3+[(1.0,8.0)] ], 'dimen': 3,
'paras': [ 'F, G, H, a' ], 'bound': [[(0.0,2.0)]*3+[(1.0,8.0)] ],
'text' : '\nCoefficients of Hosford criterion: ', 'labels': [['F','G','H','a']],
'error': 'The standard deviation errors are: '
}, },
'hill1948' :{'func' : Hill1948, 'hill1948' :{'name': 'Hill 1948',
'nExpo': 0,'err':np.inf, 'func': Hill1948,
'dimen': 3, 'nExpo': 0,'err':np.inf,
'bound': [ [(None,None)]*6, [(None,None)]*4 ], 'dimen': 3,
'paras': [ 'F, G, H, L, M, N', 'F, G, H, N'], 'bound': [[(None,None)]*6, [(None,None)]*4 ],
'text' : '\nCoefficients of Hill1948 criterion: ', 'labels': [['F','G','H','L','M','N'],['F','G','H','N']],
'error': 'The standard deviation errors are: '
}, },
'hill1979' :{'func' : Hill1979, 'hill1979' :{'name': 'Hill 1979',
'nExpo': 1,'err':np.inf, 'func': Hill1979,
'dimen': 3, 'nExpo': 1,'err':np.inf,
'bound': [ [(-2.0,2.0)]*6+[(1.0,8.0)] ], 'dimen': 3,
'paras': [ 'f,g,h,a,b,c,m' ], 'bound': [[(-2.0,2.0)]*6+[(1.0,8.0)] ],
'text' : '\nCoefficients of Hill1979 criterion: ' , 'labels': [['f','g','h','a','b','c','m']],
'error': 'The standard deviation errors are: '
}, },
'drucker' :{'func' : Drucker, 'drucker' :{'name': 'Drucker',
'nExpo': 0,'err':np.inf, 'func': Drucker,
'dimen': 3, 'nExpo': 0,'err':np.inf,
'bound': [ [(None,None)]+[(-3.375, 2.25)] ], 'dimen': 3,
'paras': [ 'sigma0, C_D' ], 'bound': [[(None,None)]+[(-3.375, 2.25)]],
'text' : '\nCoefficients of Drucker criterion: ', 'labels': [['sigma0','C_D']],
'error': 'The standard deviation errors are: '
}, },
'gdrucker' :{'func' : Drucker, 'gdrucker' :{'name': 'General Drucker',
'nExpo': 1,'err':np.inf, 'func': Drucker,
'dimen': 3, 'nExpo': 1,'err':np.inf,
'bound': [ [(None,None)]+[(-3.375, 2.25)]+[(1.0,8.0)] ], 'dimen': 3,
'paras': [ 'sigma0, C_D, p' ], 'bound': [[(None,None)]+[(-3.375, 2.25)]+[(1.0,8.0)] ],
'text' : '\nCoefficients of general Drucker criterion: ', 'labels': [['sigma0','C_D', 'p']],
'error': 'The standard deviation errors are: '
}, },
'barlat1989' :{'func' : Barlat1989, 'barlat1989' :{'name': 'Barlat 1989',
'nExpo': 1,'err':np.inf, 'func': Barlat1989,
'dimen': 2, 'nExpo': 1,'err':np.inf,
'bound': [ [(-3.0,3.0)]*4+[(1.0,8.0)] ], 'dimen': 2,
'paras': [ 'a,c,h,f, m' ], 'bound': [[(-3.0,3.0)]*4+[(1.0,8.0)] ],
'text' : '\nCoefficients of isotropic Barlat 1989 criterion: ', 'labels': [['a','c','h','f', 'm']],
'error': 'The standard deviation errors are: '
}, },
'barlat1991' :{'func' : Barlat1991, 'barlat1991' :{'name': 'Barlat 1991',
'nExpo': 1,'err':np.inf, 'func': Barlat1991,
'dimen': 3, 'nExpo': 1,'err':np.inf,
'bound': [ [(-2,2)]*6+[(1.0,8.0)], [(-2,2)]*4+[(1.0,8.0)] ], 'dimen': 3,
'paras': ['a, b, c, f, g, h, m', 'a, b, c, f, m'], 'bound': [[(-2,2)]*6+[(1.0,8.0)], [(-2,2)]*4+[(1.0,8.0)]],
'text' : '\nCoefficients of anisotropic Barlat 1991 criterion: ', 'labels': [['a','b','c','f','g','h','m'],['a','b','c','f','m']],
'error': 'The standard deviation errors are: '
}, },
'bbc2000' :{'func' : BBC2000, 'bbc2000' :{'name': 'Banabic-Balan-Comsa 2000',
'nExpo': 1,'err':np.inf, 'func': BBC2000,
'dimen': 2, 'nExpo': 1,'err':np.inf,
'bound': [ [(None,None)]*7+[(1.0,8.0)] ], #[(None,None)]*6+[(0.0,1.0)]+[(1.0,9.0)], 'dimen': 2,
'paras': [ 'd,e,f,g, b,c,a, k' ], 'bound': [[(None,None)]*7+[(1.0,8.0)]], #[(None,None)]*6+[(0.0,1.0)]+[(1.0,9.0)],
'text' : '\nCoefficients of Banabic-Balan-Comsa 2000 criterion: ', 'labels': [['d','e','f','g','b','c','a','k']],
'error': 'The standard deviation errors are: '
}, },
'bbc2003' :{'func' : BBC2003, 'bbc2003' :{'name': 'Banabic-Balan-Comsa 2003',
'nExpo': 1,'err':np.inf, 'func' : BBC2003,
'dimen': 2, 'nExpo': 1,'err':np.inf,
'bound': [ [(None,None)]*8+[(1.0,8.0)] ], #[(None,None)]*7+[(0.0,1.0)]+[(1.0,9.0)], 'dimen': 2,
'paras': [ 'M, N, P, Q, R, S, T, a, k' ], 'bound': [[(None,None)]*8+[(1.0,8.0)]], #[(None,None)]*7+[(0.0,1.0)]+[(1.0,9.0)],
'text' : '\nCoefficients of Banabic-Balan-Comsa 2003 criterion: ', 'labels': [['M','N','P','Q','R','S','T','a','k']],
'error': 'The standard deviation errors are: '
}, },
'bbc2005' :{'func' : BBC2005, 'bbc2005' :{'name': 'Banabic-Balan-Comsa 2005',
'nExpo': 1,'err':np.inf, 'func' : BBC2005,
'dimen': 2, 'nExpo': 1,'err':np.inf,
'bound': [ [(None,None)]*8+[(1.0,8.0)] ], #[(None,None)]*6+[(0.0,1.0)]*2+[(1.0,9.0)], 'dimen': 2,
'paras': [ 'L ,M, N, P, Q, R, a, b, k' ], 'bound': [[(None,None)]*8+[(1.0,8.0)] ], #[(None,None)]*6+[(0.0,1.0)]*2+[(1.0,9.0)],
'text' : '\nCoefficients of Banabic-Balan-Comsa 2005 criterion: ', 'labels': [['L','M','N','P','Q','R','a','b','k']],
'error': 'The standard deviation errors are: '
}, },
'cazacu' :{'func' : Cazacu_Barlat, 'cazacu' :{'name': 'Cazacu Barlat',
'nExpo': 0,'err':np.inf, 'func': Cazacu_Barlat,
'dimen': 3, 'nExpo': 0,'err':np.inf,
'bound': [ [(None,None)]*16+[(-2.5,2.5)]+[(None,None)] ], 'dimen': 3,
'paras': [ 'a1,a2,a3,a4,a5,a6; b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11; c','a1,a2,a3,a6; b1,b2,b3,b4,b5,b10; c'], 'bound': [[(None,None)]*16+[(-2.5,2.5)]+[(None,None)]],
'text' : '\nCoefficients of Cazacu Barlat yield criterion: ', 'labels': [['a1','a2','a3','a4','a5','a6', 'b1','b2','b3','b4','b5','b6','b7','b8','b9','b10','b11', 'c'],
'error': 'The standard deviation errors are: ' ['a1','a2','a3','a6', 'b1','b2','b3','b4','b5','b10', 'c']],
}, },
'yld2000' :{'func' : Yld2000, 'yld2000' :{'name': 'Yld2000-2D',
'nExpo': 1,'err':np.inf, 'func': Yld2000,
'dimen': 2, 'nExpo': 1,'err':np.inf,
'bound': [ [(None,None)]*8+[(1.0,8.0)] ], 'dimen': 2,
'paras': [ 'a1,a2,a7,a3,a4,a5,a6,a8,m' ], 'bound': [[(None,None)]*8+[(1.0,8.0)]],
'text' : '\nCoefficients of Yld2000-2D yield criterion: ', 'labels': [['a1','a2','a7','a3','a4','a5','a6','a8','m']],
'error': 'The standard deviation errors are: '
}, },
'yld200418p' :{'func' : Yld200418p, 'yld200418p' :{'name': 'Yld2004-18p',
'nExpo': 1,'err':np.inf, 'func' : Yld200418p,
'dimen': 3, 'nExpo': 1,'err':np.inf,
'bound': [ [(None,None)]*18+[(1.0,8.0)], [(None,None)]*14+[(1.0,8.0)] ], 'dimen': 3,
'paras': [ 'c12,c21,c23,c32,c31,c13,c44,c55,c66,d12,d21,d23,d32,d31,d13,d44,d55,d66,m', \ 'bound': [[(None,None)]*18+[(1.0,8.0)], [(None,None)]*14+[(1.0,8.0)]],
'c12,c21,c23,c32,c31,c13,c44,d12,d21,d23,d32,d31,d13,d44,m' ], 'labels': [['c12','c21','c23','c32','c31','c13','c44','c55','c66','d12','d21','d23','d32','d31','d13','d44','d55','d66','m'],
'text' : '\nCoefficients of Yld2004-18p yield criterion: ', ['c12','c21','c23','c32','c31','c13','c44','d12','d21','d23','d32','d31','d13','d44','m']],
'error': 'The standard deviation errors are: '
}, },
'karafillis' :{'func' : KarafillisBoyce, 'karafillis' :{'name': 'Karafillis-Boyce',
'nExpo': 1,'err':np.inf, 'func' : KarafillisBoyce,
'dimen': 3, 'nExpo': 1,'err':np.inf,
'bound': [ [(None,None)]*6+[(0.0,1.0)]+[(1.0,8.0)], [(None,None)]*4+[(0.0,1.0)]+[(1.0,8.0)]], 'dimen': 3,
'paras': [ 'c11,c12,c13,c14,c15,c16,c,m', \ 'bound': [[(None,None)]*6+[(0.0,1.0)]+[(1.0,8.0)], [(None,None)]*4+[(0.0,1.0)]+[(1.0,8.0)]],
'c11,c12,c13,c14,c15,c16,c,m' ], 'labels': [['c11','c12','c13','c14','c15','c16','c','m'],
'text' : '\nCoefficients of Karafillis-Boyce yield criterion: ', ['c11','c12','c13','c14','c15','c16','c','m']],
'error': 'The standard deviation errors are: '
}, },
'all' : 'fit all the criteria'
} }
thresholdParameter = ['totalshear','equivalentStrain'] thresholdParameter = ['totalshear','equivalentStrain']
#--------------------------------------------------------------------------------------------------- #---------------------------------------------------------------------------------------------------
class Loadcase(): class Loadcase():
#--------------------------------------------------------------------------------------------------- #---------------------------------------------------------------------------------------------------
@ -995,7 +979,7 @@ class Loadcase():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def __init__(self,finalStrain,incs,time,ND=3,RD=1,nSet=1,dimension=3,vegter=False): def __init__(self,finalStrain,incs,time,ND=3,RD=1,nSet=1,dimension=3,vegter=False):
print('using the random load case generator') damask.util.croak('Using the random load case generator')
self.finalStrain = finalStrain self.finalStrain = finalStrain
self.incs = incs self.incs = incs
self.time = time self.time = time
@ -1010,14 +994,14 @@ class Loadcase():
def getLoadcase(self,number): def getLoadcase(self,number):
if self.dimension == 3: if self.dimension == 3:
print 'generate random 3D load case' damask.util.croak('Generate random 3D load case')
return self._getLoadcase3D() return self._getLoadcase3D()
else: else:
if self.vegter is True: if self.vegter is True:
print 'generate load case for Vegter' damask.util.croak('Generate load case for Vegter')
return self._getLoadcase2dVegter(number) return self._getLoadcase2dVegter(number)
else: else:
print 'generate random 2D load case' damask.util.croak('Generate random 2D load case')
return self._getLoadcase2dRandom() return self._getLoadcase2dRandom()
def _getLoadcase3D(self): def _getLoadcase3D(self):
@ -1127,17 +1111,26 @@ class Criterion(object):
''' '''
Fitting to certain criterion Fitting to certain criterion
''' '''
def __init__(self, exponent, uniaxial, dimension, name='vonmises'): def __init__(self, exponent, uniaxial, dimension, label='vonmises'):
self.name = name self.name = label
self.expo = exponent self.expo = exponent
self.uniaxial= uniaxial self.uniaxial= uniaxial
self.dimen = dimension self.dimen = dimension
self.results = fitCriteria self.results = fitCriteria
if self.name.lower() not in map(str.lower, self.results.keys()): if self.name.lower() not in map(str.lower, self.results.keys()):
raise Exception('no suitable fitting criterion selected') raise Exception('No suitable fitting criterion selected')
else: else:
print('fitting to the %s criterion'%name) damask.util.croak('Fitting to the %s criterion'%fitCriteria[self.name]['name'])
def report_labels(self):
if len(fitCriteria[self.name]['labels']) > 1 and self.dimen == 2:
return fitCriteria[self.name]['labels'][1]
else:
return fitCriteria[self.name]['labels'][0]
def report_name(self):
return fitCriteria[self.name]['name']
def fit(self,stress): def fit(self,stress):
global fitResults; fitErrors; fitResidual global fitResults; fitErrors; fitResidual
@ -1145,23 +1138,19 @@ class Criterion(object):
else: nExponent = 0 else: nExponent = 0
nameCriterion = self.name.lower() nameCriterion = self.name.lower()
criteria = Criteria(nameCriterion,self.uniaxial,self.expo, self.dimen) criteria = Criteria(nameCriterion,self.uniaxial,self.expo, self.dimen)
textParas = fitCriteria[nameCriterion]['text']+fitCriteria[nameCriterion]['paras'][dDim]+':\n' + \ bounds = fitCriteria[nameCriterion]['bound'][dDim] # Default bounds, no bound
formatOutput(numParas+nExponent) guess0 = Guess # Default initial guess, depends on bounds
textError = fitCriteria[nameCriterion]['error']+ formatOutput(numParas+nExponent,'%-14.8f')
bounds = fitCriteria[nameCriterion]['bound'][dDim] # Default bounds, no bound
guess0 = Guess # Default initial guess, depends on bounds
if fitResults == [] : initialguess = guess0 if fitResults == [] : initialguess = guess0
else : initialguess = np.array(fitResults[-1]) else : initialguess = np.array(fitResults[-1])
weight = get_weight(np.shape(stress)[1])
ydata = np.zeros(np.shape(stress)[1]) ydata = np.zeros(np.shape(stress)[1])
try: try:
popt, pcov, infodict, errmsg, ierr = \ popt, pcov, infodict, errmsg, ierr = \
leastsqBound (criteria.fun, initialguess, args=(ydata,stress), leastsqBound (criteria.fun, initialguess, args=(ydata,stress),
bounds=bounds, Dfun=criteria.jac, full_output=True) bounds=bounds, Dfun=criteria.jac, full_output=True)
if ierr not in [1, 2, 3, 4]: if ierr not in [1, 2, 3, 4]:
raise RuntimeError("Optimal parameters not found: " + errmsg) raise RuntimeError("Optimal parameters not found: "+errmsg)
else: else:
residual = criteria.fun(popt, ydata, stress) residual = criteria.fun(popt, ydata, stress)
fitResidual.append(np.linalg.norm(residual)/np.sqrt(len(residual))) fitResidual.append(np.linalg.norm(residual)/np.sqrt(len(residual)))
@ -1175,12 +1164,11 @@ class Criterion(object):
popt = np.concatenate((np.array(popt), np.repeat(options.exponent,nExponent))) popt = np.concatenate((np.array(popt), np.repeat(options.exponent,nExponent)))
perr = np.concatenate((np.array(perr), np.repeat(0.0,nExponent))) perr = np.concatenate((np.array(perr), np.repeat(0.0,nExponent)))
print (textParas%array2tuple(popt)) damask.util.croak('Needed {} function calls for fitting'.format(infodict['nfev']))
print (textError%array2tuple(perr))
print('Number of function calls =', infodict['nfev'])
except Exception as detail: except Exception as detail:
print detail damask.util.croak(detail)
pass pass
return popt
#--------------------------------------------------------------------------------------------------- #---------------------------------------------------------------------------------------------------
@ -1205,9 +1193,10 @@ class myThread (threading.Thread):
def doSim(delay,thread): def doSim(delay,thread):
s.acquire() s.acquire()
global myLoad
loadNo=loadcaseNo() loadNo=loadcaseNo()
if not os.path.isfile('%s.load'%loadNo): if not os.path.isfile('%s.load'%loadNo):
print('generating loadcase for sim %s from %s'%(loadNo,thread)) damask.util.croak('Generating load case for simulation %s (%s)'%(loadNo,thread))
f=open('%s.load'%loadNo,'w') f=open('%s.load'%loadNo,'w')
f.write(myLoad.getLoadcase(loadNo)) f.write(myLoad.getLoadcase(loadNo))
f.close() f.close()
@ -1216,27 +1205,26 @@ def doSim(delay,thread):
s.acquire() s.acquire()
if not os.path.isfile('%s_%i.spectralOut'%(options.geometry,loadNo)): if not os.path.isfile('%s_%i.spectralOut'%(options.geometry,loadNo)):
print('starting simulation %s from %s'%(loadNo,thread)) damask.util.croak('Starting simulation %i (%s)'%(loadNo,thread))
s.release() s.release()
damask.util.execute('DAMASK_spectral -g %s -l %i'%(options.geometry,loadNo)) damask.util.execute('DAMASK_spectral -g %s -l %i'%(options.geometry,loadNo))
else: s.release() else: s.release()
s.acquire() s.acquire()
if not os.path.isfile('./postProc/%s_%i.txt'%(options.geometry,loadNo)): if not os.path.isfile('./postProc/%s_%i.txt'%(options.geometry,loadNo)):
print('starting post processing for sim %i from %s'%(loadNo,thread)) damask.util.croak('Starting post processing for simulation %i (%s)'%(loadNo,thread))
s.release() s.release()
try: try:
damask.util.execute('postResults --cr f,p --co totalshear %s_%i.spectralOut'%(options.geometry,loadNo)) damask.util.execute('postResults --cr f,p --co totalshear %s_%i.spectralOut'%(options.geometry,loadNo))
except: except:
damask.util.execute('postResults --cr f,p %s_%i.spectralOut'%(options.geometry,loadNo)) damask.util.execute('postResults --cr f,p %s_%i.spectralOut'%(options.geometry,loadNo))
damask.util.execute('addCauchy ./postProc/%s_%i.txt'%(options.geometry,loadNo)) damask.util.execute('addCauchy ./postProc/%s_%i.txt'%(options.geometry,loadNo))
damask.util.execute('addStrainTensors -l -v ./postProc/%s_%i.txt'%(options.geometry,loadNo)) damask.util.execute('addStrainTensors -0 -v ./postProc/%s_%i.txt'%(options.geometry,loadNo))
damask.util.execute('addMises -s Cauchy -e ln(V) ./postProc/%s_%i.txt'%(options.geometry,loadNo)) damask.util.execute('addMises -s Cauchy -e ln(V) ./postProc/%s_%i.txt'%(options.geometry,loadNo))
else: s.release() else: s.release()
s.acquire() s.acquire()
print('-'*10) damask.util.croak('Reading values from simulation %i (%s)'%(loadNo,thread))
print('reading values for sim %i from %s'%(loadNo,thread))
s.release() s.release()
refFile = './postProc/%s_%i.txt'%(options.geometry,loadNo) refFile = './postProc/%s_%i.txt'%(options.geometry,loadNo)
@ -1248,7 +1236,7 @@ def doSim(delay,thread):
thresholdKey = 'totalshear' thresholdKey = 'totalshear'
s.acquire() s.acquire()
for l in [thresholdKey,'1_Cauchy']: for l in [thresholdKey,'1_Cauchy']:
if l not in table.labels: print '%s not found'%l if l not in table.labels: damask.util.croak('%s not found'%l)
s.release() s.release()
table.data_readArray(['%i_Cauchy'%(i+1) for i in xrange(9)]+[thresholdKey]+['%i_ln(V)'%(i+1) for i in xrange(9)]) table.data_readArray(['%i_Cauchy'%(i+1) for i in xrange(9)]+[thresholdKey]+['%i_ln(V)'%(i+1) for i in xrange(9)])
@ -1282,23 +1270,25 @@ def doSim(delay,thread):
else: else:
line+=1 line+=1
if not validity[i]: if not validity[i]:
print ('The data of sim %i at the threshold %f is invalid, the fitting at this point is skipped'%(loadNo,threshold)) damask.util.croak('The data of result %i at the threshold %f is invalid,'\
+'the fitting at this point is skipped'%(loadNo,threshold))
s.acquire() s.acquire()
global stressAll, strainAll global stressAll, strainAll
print('number of yield points of sim %i: %i'%(loadNo,len(yieldStress))) f=open(options.geometry+'_'+options.criterion+'.txt','w')
print('starting fitting for sim %i from %s'%(loadNo,thread)) f.write(' '.join([options.fitting]+myFit.report_labels())+'\n')
try: try:
for i in xrange(int(options.yieldValue[2])): for i,threshold in enumerate(np.linspace(options.yieldValue[0],options.yieldValue[1],options.yieldValue[2])):
if validity[i]: if validity[i]:
a = (yieldStress[0][2]/stressUnit)**2 + (yieldStress[0][4]/stressUnit)**2 + (yieldStress[0][5]/stressUnit)**2 a = (yieldStress[0][2]/stressUnit)**2 + (yieldStress[0][4]/stressUnit)**2 + (yieldStress[0][5]/stressUnit)**2
stressAll[i]=np.append(stressAll[i], yieldStress[i]/stressUnit) stressAll[i]=np.append(stressAll[i], yieldStress[i]/stressUnit)
strainAll[i]=np.append(strainAll[i], deformationRate[i]) strainAll[i]=np.append(strainAll[i], deformationRate[i])
myFit.fit(stressAll[i].reshape(len(stressAll[i])//6,6).transpose()) f.write( str(threshold)+' '+
' '.join(map(str,myFit.fit(stressAll[i].reshape(len(stressAll[i])//6,6).transpose())))+'\n')
except Exception as detail: except Exception as detail:
print('could not fit for sim %i from %s'%(loadNo,thread)) damask.util.croak('Could not fit results of simulation (%s)'%thread)
print detail
s.release() s.release()
return return
damask.util.croak('\n')
s.release() s.release()
def loadcaseNo(): def loadcaseNo():
@ -1312,7 +1302,8 @@ def converged():
if N_simulations < options.max: if N_simulations < options.max:
if len(fitResidual) > 5: if len(fitResidual) > 5:
residualList = np.array(fitResidual[len(fitResidual)-5:]) residualList = np.array(fitResidual[len(fitResidual)-5:])
if np.std(residualList)/np.max(residualList) < 0.05: return True if np.std(residualList)/np.max(residualList) < 0.05:
return True
return False return False
else: else:
return True return True
@ -1324,8 +1315,8 @@ def converged():
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Performs calculations with various loads on given geometry file and fits yield surface. Performs calculations with various loads on given geometry file and fits yield surface.
""", version=string.replace(scriptID,'\n','\\n') """, version = scriptID)
)
# maybe make an option to specifiy if 2D/3D fitting should be done? # maybe make an option to specifiy if 2D/3D fitting should be done?
parser.add_option('-l','--load' , dest='load', type='float', nargs=3, parser.add_option('-l','--load' , dest='load', type='float', nargs=3,
@ -1387,14 +1378,11 @@ if options.yieldValue[2] != int(options.yieldValue[2]):
parser.error('count must be an integer') parser.error('count must be an integer')
if options.dimension not in [2,3]: if options.dimension not in [2,3]:
parser.error('Dimension is wrong, should be 2(plane stress state) or 3(general stress state)') parser.error('Dimension is wrong, should be 2(plane stress state) or 3(general stress state)')
#if options.criterion not in ['tresca', 'vonmises', 'hershey','drucker', 'gdrucker', 'hill1948']:
# if options.eqStress == None:
# parser.error("The equivalent stress is indispensable for the yield criterion '"+ options.criterion+"'")
if not os.path.isfile('numerics.config'): if not os.path.isfile('numerics.config'):
print('numerics.config file not found') damask.util.croak('numerics.config file not found')
if not os.path.isfile('material.config'): if not os.path.isfile('material.config'):
print('material.config file not found') damask.util.croak('material.config file not found')
if options.vegter is True: if options.vegter is True:
@ -1426,4 +1414,4 @@ else:
run = runFit(options.exponent, eqStress, fitCriteria[criter]['dimen'], criter) run = runFit(options.exponent, eqStress, fitCriteria[criter]['dimen'], criter)
fit_allResults[criter] = {'results': fitResults, 'errors': fitErrors, 'residual': fitResidual} fit_allResults[criter] = {'results': fitResults, 'errors': fitErrors, 'residual': fitResidual}
print 'Finished fitting to yield criteria' damask.util.croak('Finished fitting to yield criteria')