first draft of a script to compute yield surfaces with the spectral solver

processing/misc/yieldSurface.py

@@ -0,0 +1,165 @@
+#!/usr/bin/python
+# -*- coding: UTF-8 no BOM -*-
+
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.linalg import svd
+import threading
+import time
+
+def asFullTensor(voigt):
+  return np.array([[voigt[0],voigt[3],voigt[5]],\
+                   [voigt[3],voigt[1],voigt[4]],\
+                   [voigt[5],voigt[4],voigt[2]]])
+
+def Hill48(x, F,G,H,L,M,N):
+  return F*(x[1]-x[2])**2 + G*(x[2]-x[0])**2 + H*(x[0]-x[1])** + \
+         2*L*x[4]**2 + 2*M*x[5]**2 + 2*N*x[3]**2 -1.
+
+def vonMises(x, S_y):
+  p = svd(asFullTensor(x)) 
+  return (s[2]-s[1])**2+(s[1]-s[0])**2+(s[0]-s[2])**2-2*S_y**2
+
+#---------------------------------------------------------------------------------------------------
+class Loadcase():
+#---------------------------------------------------------------------------------------------------
+  '''
+     Class for generating load cases for the spectral solver
+  '''
+
+# ------------------------------------------------------------------
+  def __init__(self):
+    print('using the random load case generator')
+
+  def getNext(self,N=0):
+    defgrad=['*']*9
+    stress =[0]*9
+    values=(np.random.random_sample(9)-.5)*scale*2
+
+    main=np.array([0,4,8])
+    np.random.shuffle(main)
+    for i in main[:2]:
+      defgrad[i]=1.+values[i]
+      stress[i]='*'
+    for off in [[1,3,0],[2,6,0],[5,7,0]]:
+      off=np.array(off)
+      np.random.shuffle(off)
+      print off
+      if off[0] != 0: 
+        defgrad[off[0]]=values[off[0]]
+        stress[off[0]]='*'
+    return 'f '+' '.join(str(c) for c in defgrad)+\
+          ' p '+' '.join(str(c) for c in stress)+\
+          ' incs %s'%incs+\
+          ' time %s'%duration
+
+#---------------------------------------------------------------------------------------------------
+class Criterion(object):
+#---------------------------------------------------------------------------------------------------
+  '''
+     Fitting to certain criterion
+  '''
+  def __init__(self,name):
+    self.name = name.lower()
+    if self.name not in ['hill48','vonmises']: print('Mist')
+    print('using the %s criterion'%self.name)
+    self.popt = 0.0
+
+  def fit(self,stress):
+    print len(stress)
+    if self.name == 'hill48':
+      try:
+        self.popt, pcov = curve_fit(Hill48, stress, np.zeros(np.shape(stress)[0]))
+        print self.popt
+      except:
+        pass
+    elif self.name == 'vonmises':
+      try:
+        self.popt, pcov = curve_fit(vonMises, stress.transpose(), np.shape(stress)[0])
+      except:
+        pass 
+
+
+#---------------------------------------------------------------------------------------------------
+
+#---------------------------------------------------------------------------------------------------
+  '''
+     Runner class
+  '''
+class myThread (threading.Thread):
+  def __init__(self, threadID):
+    threading.Thread.__init__(self)
+    self.threadID = threadID
+  def run(self):
+    s.acquire()
+    conv=converged()
+    s.release()
+    while not conv:
+      doSim(4.,self.name)
+      s.acquire()
+      conv=converged()
+      s.release()
+
+def doSim(delay,thread):
+  s.acquire()
+  me=getLoadcase()+1
+  print('starting sim %i from thread %s'%(me,thread))
+  f=open('%s.load'%me,'w')
+  f.write(myLoad.getNext(me))
+  f.close()
+  #dummy
+  voigt = np.random.random_sample(6)*90e6
+  global stressAll
+  stressAll=np.append(stressAll,asFullTensor(voigt).reshape(1,9))
+  stressAll=stressAll.reshape(len(stressAll)//9,9)
+  myFit.fit(stressAll)
+  s.release()
+  time.sleep(delay)
+  s.acquire()
+  print('doing postprocessing sim %i from thread %s'%(me,thread))
+  s.release()
+
+def getLoadcase():
+  global N_simulations
+  N_simulations+=1
+  return N_simulations
+
+def converged():
+  global N_simulations
+  global maxN_simulations
+  if N_simulations < maxN_simulations:
+    return False
+  else:
+    return True
+
+# main
+
+minN_simulations=20
+maxN_simulations=10
+N_simulations=0
+s=threading.Semaphore(1)
+scale = 0.02
+incs = 10
+duration = 10
+stressAll=np.zeros(0,'d')
+
+myLoad = Loadcase()
+myFit = Criterion('Hill48')
+
+N_threads=3
+t=[]
+
+for i in range(N_threads):
+  t.append(myThread(i))
+  t[i].start()
+
+for i in range(N_threads):
+  t[i].join()
+
+a, b = curve_fit(Hill48, stressAll.transpose(), np.zeros(10))
+print a
+print "Exiting Main Thread"
+
+
+
+