added call to DAMASK spectral

processing/misc/yieldSurface.py

@@ -4,8 +4,19 @@
 import numpy as np
 from scipy.optimize import curve_fit
 from scipy.linalg import svd
-import threading
+import threading,time,os,subprocess,shlex
-import time
+import damask
+
+def execute(cmd,dir='./'):
+      
+  initialPath=os.getcwd()
+  out = ''
+  line = True
+  process = subprocess.Popen(shlex.split(cmd),stdout=subprocess.PIPE,stderr = subprocess.STDOUT)
+  while line:
+    line = process.stdout.readline()
+    out += line
+  os.chdir(initialPath)
 
 def asFullTensor(voigt):
   return np.array([[voigt[0],voigt[3],voigt[5]],\
@@ -18,8 +29,13 @@ def Hill48(x, F,G,H,L,M,N):
   return a.ravel()
 
 def vonMises(x, S_y):
-  p = svd(asFullTensor(x)) 
+  sv=np.zeros(0,'d')
-  return (s[2]-s[1])**2+(s[1]-s[0])**2+(s[0]-s[2])**2-2*S_y**2
+  for i in xrange(np.shape(x)[1]):
+    U, l, Vh = svd(np.array(x[:,i]).reshape(3,3))
+    sv = np.append(sv,l)
+  sv = sv.reshape(3,np.shape(x)[1])
+  ooo= (sv[2,:]-sv[1,:])**2+(sv[1,:]-sv[0,:])**2+(sv[0,:]-sv[2,:])**2-2*S_y**2
+  return ooo
 
 #---------------------------------------------------------------------------------------------------
 class Loadcase():
@@ -39,15 +55,17 @@ class Loadcase():
 
     main=np.array([0,4,8])
     np.random.shuffle(main)
-    for i in main[:2]:
+    for i in main[:2]: # fill 2 out of 3 main entries
       defgrad[i]=1.+values[i]
       stress[i]='*'
-    for off in [[1,3,0],[2,6,0],[5,7,0]]:
+    for off in [[1,3,0],[2,6,0],[5,7,0]]: # fill 3 off-diagonal pairs of defgrad (1 or 2 entries)
       off=np.array(off)
       np.random.shuffle(off)
-      if off[0] != 0: 
+      for i in off[0:2]:
-        defgrad[off[0]]=values[off[0]]
+        if i != 0: 
-        stress[off[0]]='*'
+          defgrad[i]=values[i]
+          stress[i]='*'
+
     return 'f '+' '.join(str(c) for c in defgrad)+\
           ' p '+' '.join(str(c) for c in stress)+\
           ' incs %s'%incs+\
@@ -66,17 +84,14 @@ class Criterion(object):
     self.popt = 0.0
 
   def fit(self,stress):
-    if self.name == 'hill48':
     try:
-        self.popt, pcov = curve_fit(Hill48, stress, np.zeros(np.shape(stress)[1]))
+      popt, pcov = curve_fit(Hill48, stress, np.zeros(np.shape(stress)[1]))
-        print self.popt
+      print 'Hill 48', popt
-      except:
+    except Exception as detail:
-        pass
+      print detail
-    elif self.name == 'vonmises':
-      try:
-        self.popt, pcov = curve_fit(vonMises, stress.transpose(), np.shape(stress)[1])
-      except:
       pass
+    popt, pcov = curve_fit(vonMises, stress, np.zeros(np.shape(stress)[1]))
+    print 'von Mises', popt
 
 
 #---------------------------------------------------------------------------------------------------
@@ -100,23 +115,43 @@ class myThread (threading.Thread):
       s.release()
 
 def doSim(delay,thread):
+
   s.acquire()
   me=getLoadcase()
-  print('starting sim %i from %s'%(me,thread))
+  print('generating loadcase for sim %i from %s'%(me,thread))
   f=open('%s.load'%me,'w')
   f.write(myLoad.getNext(me))
   f.close()
-  #dummy
+
-  print('doing postprocessing sim %i from %s'%(me,thread))
+  print('starting simulation %i from %s'%(me,thread))
-  voigt = np.random.random_sample(6)*90e6
+  s.release()
+
+  execute('DAMASK_spectral -l %i -g 20grains16x16x16'%me)
+
+  s.acquire()
+  print('startin post processing for sim %i from %s'%(me,thread))
+  s.release()
+  execute('postResults --cr f,p 20grains16x16x16_%i.spectralOut'%me)
+  execute('addCauchy ./postProc/20grains16x16x16_%i.txt'%me)
+  execute('addStrainTensors -l -v ./postProc/20grains16x16x16_%i.txt'%me)
+  execute('addMises -s Cauchy -e ln(V) ./postProc/20grains16x16x16_%i.txt'%me)
+
+  refFile = open('./postProc/20grains16x16x16_%i.txt'%me)
+  table = damask.ASCIItable(refFile)
+  table.head_read()
+  for l in ['Mises(ln(V))','1_Cauchy']:
+    if l not in table.labels: print '%s not found'%l
+  while table.data_read():
+    if float(table.data[table.labels.index('Mises(ln(V))')]) > 0.002:
+      yieldStress = np.array(table.data[table.labels.index('1_Cauchy'):table.labels.index('9_Cauchy')+1],'d').reshape(3,3)
+
+  s.acquire()
+  print('startin fitting for sim %i from %s'%(me,thread))
   global stressAll
-  stressAll=np.append(stressAll,asFullTensor(voigt).reshape(9))
+  stressAll=np.append(stressAll,yieldStress.reshape(9))
   stressAll=stressAll.reshape(9,len(stressAll)//9)
   myFit.fit(stressAll)
   s.release()
-  time.sleep(delay)
-  s.acquire()
-  s.release()
 
 def getLoadcase():
   global N_simulations
@@ -134,7 +169,7 @@ def converged():
 # main
 
 minN_simulations=20
-maxN_simulations=10
+maxN_simulations=20
 N_simulations=0
 s=threading.Semaphore(1)
 scale = 0.02