/[escript]/trunk/escript/py_src/rheologies.py
ViewVC logotype

Annotation of /trunk/escript/py_src/rheologies.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2344 - (hide annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 22557 byte(s)
Change __url__ to launchpad site

1 ksteube 1809 ########################################################
2 gross 1639 #
3 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
4     # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6 gross 1639 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10 gross 1639 #
11 ksteube 1809 ########################################################
12 gross 1639
13 ksteube 1809 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
14     Earth Systems Science Computational Center (ESSCC)
15     http://www.uq.edu.au/esscc
16     Primary Business: Queensland, Australia"""
17     __license__="""Licensed under the Open Software License version 3.0
18     http://www.opensource.org/licenses/osl-3.0.php"""
19 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
20 ksteube 1809
21 gross 1639 """
22     Some models for flow
23    
24     @var __author__: name of author
25     @var __copyright__: copyrights
26     @var __license__: licence agreement
27     @var __url__: url entry point on documentation
28     @var __version__: version
29     @var __date__: date of the version
30     """
31    
32     __author__="Lutz Gross, l.gross@uq.edu.au"
33    
34     from escript import *
35     import util
36     from linearPDEs import LinearPDE
37 gross 2100 from pdetools import Defect, NewtonGMRES, ArithmeticTuple
38 gross 1639
39 gross 2300 class PowerLaw(object):
40     """
41     this implements the power law for a composition of a set of materials where the viscosity eta of each material is given by a
42     power law relationship of the form
43    
44     M{eta=eta_N*(tau/tau_t)**(1.-1./power)}
45    
46     where tau is equivalent stress and eta_N, tau_t and power are given constant. Moreover an elastic component can be considered.
47     Moreover tau meets the Drucker-Prager type yield condition
48    
49     M{tau <= tau_Y + friction * pressure}
50    
51     where gamma_dot is the equivalent.
52     """
53     def __init__(self, numMaterials=1,verbose=False):
54     """
55     initializes a power law
56    
57     @param numMaterials: number of materials
58     @type numMaterials: C{int}
59     @param verbose: if C{True} some informations are printed.
60     @type verbose: C{bool}
61     """
62     if numMaterials<1:
63     raise ValueError,"at least one material must be defined."
64     self.__numMaterials=numMaterials
65     self.__eta_N=[None for i in xrange(self.__numMaterials)]
66     self.__tau_t=[1. for i in xrange(self.__numMaterials)]
67     self.__power=[1. for i in xrange(self.__numMaterials)]
68     self.__tau_Y=None
69 gross 2301 self.__friction=None
70 gross 2300 self.__mu=None
71     self.__verbose=verbose
72     self.setEtaTolerance()
73     #===========================================================================
74     def getNumMaterials(self):
75     """
76     returns the numebr of materials
77    
78     @return: number of materials
79     @rtype: C{int}
80     """
81     return self.__numMaterials
82     def validMaterialId(self,id=0):
83     """
84     checks if a given material id is valid
85    
86     @param id: a material id
87     @type id: C{int}
88     @return: C{True} is the id is valid
89     @rtype: C{bool}
90     """
91     return 0<=id and id<self.getNumMaterials()
92     def setEtaTolerance(self,rtol=util.sqrt(util.EPSILON)):
93     """
94     sets the relative tolerance for the effectice viscosity.
95    
96     @param rtol: relative tolerance
97     @type rtol: positive C{float}
98     """
99     if rtol<=0:
100     raise ValueError,"rtol needs to positive."
101     self.__rtol=rtol
102     def getEtaTolerance(self):
103     """
104     returns the relative tolerance for the effectice viscosity.
105    
106     @return: relative tolerance
107     @rtype rtol: positive C{float}
108     """
109     return self.__rtol
110     #===========================================================================
111 gross 2301 def setDruckerPragerLaw(self,tau_Y=None,friction=None):
112 gross 2300 """
113     Sets the parameters for the Drucker-Prager model.
114    
115     @param tau_Y: yield stress
116     @param friction: friction coefficient
117     """
118     self.__tau_Y=tau_Y
119     self.__friction=friction
120     def getFriction(self):
121     """
122     returns the friction coefficient
123    
124     @return: friction coefficient
125     """
126     return self.__friction
127     def getTauY(self):
128     """
129     returns the yield stress
130    
131     @return: the yield stress
132     """
133     return self.__tau_Y
134     #===========================================================================
135     def getElasticShearModulus(self):
136     """
137     returns the elastic shear modulus.
138    
139     @return: elastic shear modulus
140     """
141     return self.__mu
142     def setElasticShearModulus(self,mu=None):
143     """
144     Sets the elastic shear modulus.
145    
146     @param mu: elastic shear modulus
147     """
148     self.__mu=mu
149     #===========================================================================
150     def getPower(self, id=None):
151     """
152     returns the power in the power law
153    
154     @param id: if present, the power for material C{id} is returned.
155     @type id: C{int}
156     @return: the list of the powers for all matrials is returned. If C{id} is present only the power for material C{id} is returned.
157     """
158     if id == None:
159     return self.__power
160     else:
161     if self.validMaterialId(id):
162     return self.__power[id]
163     else:
164     raise ValueError,"Illegal material id %s."%id
165     def getEtaN(self, id=None):
166     """
167     returns the viscosity
168    
169     @param id: if present, the viscosity for material C{id} is returned.
170     @type id: C{int}
171     @return: the list of the viscosities for all matrials is returned. If C{id} is present only the viscosity for material C{id} is returned.
172     """
173     if id == None:
174     return self.__eta_N
175     else:
176     if self.validMaterialId(id):
177     return self.__eta_N[id]
178     else:
179     raise ValueError,"Illegal material id %s."%id
180     def getTauT(self, id=None):
181     """
182     returns the transition stress
183    
184     @param id: if present, the transition stress for material C{id} is returned.
185     @type id: C{int}
186     @return: the list of the transition stresses for all matrials is returned. If C{id} is present only the transition stress for material C{id} is returned.
187     """
188     if id == None:
189     return self.__tau_t
190     else:
191     if self.validMaterialId(id):
192     return self.__tau_t[id]
193     else:
194     raise ValueError,"Illegal material id %s."%id
195    
196     def setPowerLaw(self,eta_N, id=0, tau_t=1, power=1):
197     """
198     Sets the power-law parameters for material id
199    
200     @param id: material id
201     @type id: C{int}
202     @param eta_N: viscosity for tau=tau_t
203     @param tau_t: transition stress
204     @param power: power law coefficient
205     """
206     if self.validMaterialId(id):
207     self.__eta_N[id]=eta_N
208     self.__power[id]=power
209     self.__tau_t[id]=tau_t
210     else:
211     raise ValueError,"Illegal material id %s."%id
212    
213     def setPowerLaws(self,eta_N, tau_t, power):
214     """
215     Sets the parameters of the power-law for all materials.
216    
217     @param eta_N: list of viscosities for tau=tau_t
218     @param tau_t: list of transition stresses
219     @param power: list of power law coefficient
220     """
221     if len(eta_N)!=self.__numMaterials or len(tau_t)!=self.__numMaterials or len(power)!=self.__numMaterials:
222     raise ValueError,"%s materials are expected."%self.__numMaterials
223     for i in xrange(self.__numMaterials):
224     self.setPowerLaw(id=i, eta_N=eta_N[i],tau_t=tau_t[i],power=power[i])
225    
226     #===========================================================================
227 gross 2301 def getEtaEff(self,gamma_dot, eta0=None, pressure=None,dt=None, iter_max=10):
228 gross 2300 """
229     returns the effective viscosity eta_eff such that
230    
231     M{tau=eta_eff * gamma_dot}
232    
233     by solving a non-linear problem for tau.
234    
235     @param gamma_dot: equivalent strain gamma_dot
236     @param eta0: initial guess for the effective viscosity (e.g from a previous time step). If not present, an initial guess is calculated.
237     @param pressure: pressure used to calculate yield condition
238     @param dt: time step size. only needed if elastic component is considered.
239     @type dt: positive C{float} if present
240     @param iter_max: maximum number of iteration steps.
241     @type iter_max: C{int}
242     @return: effective viscosity.
243     """
244 gross 2301 SMALL=1./(util.DBLE_MAX/100.)
245 gross 2300 numMaterial=self.getNumMaterials()
246     s=[1.-1./p for p in self.getPower() ]
247     eta_N=self.getEtaN()
248     tau_t=self.getTauT()
249     mu=self.getElasticShearModulus()
250     fric=self.getFriction()
251     tau_Y=self.getTauY()
252     if eta0==None:
253     theta=0.
254     for i in xrange(numMaterial):
255     inv_eta_i=0**s[i]/eta_N[i]
256     theta=theta+inv_eta_i
257     if util.inf(theta)<=0:
258     raise ValueError,"unable to set positive initial guess for eta_eff. Most likely no power law with power 1 set."
259     eta_eff=1./theta
260     else:
261     if util.inf(eta0)<=0:
262     raise ValueError,"initial guess for eta_eff is not positive."
263     eta_eff=eta0
264    
265     if mu !=None and dt == None:
266     raise ValueError,"Time stepsize dt must be given."
267     if dt !=None:
268 gross 2301 if dt<=0: raise ValueError,"time step size must be positive."
269     if tau_Y==None and fric==None:
270 gross 2300 eta_max=None
271     else:
272 gross 2301 if fric == None:
273     eta_max=tau_Y/(gamma_dot+SMALL*util.whereZero(gamma_dot))
274     else:
275     if tau_Y==None: tau_Y==0
276     if util.inf(fric)<=0:
277     raise ValueError,"if friction present it needs to be positive."
278     eta_max=fric*util.clip(tau_Y/fric+pressure,minval=0)/(gamma_dot+SMALL*util.whereZero(gamma_dot))
279 gross 2300 rtol=self.getEtaTolerance()
280     iter =0
281     converged=False
282     tau=eta_eff*gamma_dot
283     if self.__verbose: print "Start calculation of eta_eff (tolerance = %s)\ninitial max eta_eff = %s, tau = %s."%(rtol,util.Lsup(eta_eff),util.Lsup(tau))
284     while not converged:
285     if iter>max(iter_max,1):
286 gross 2301 raise RuntimeError,"tolerance not reached after %s steps."%max(iter_max,1)
287 gross 2300 #===========================================
288     theta=0. # =1/eta
289     omega=0. # = tau*theta'= eta'*tau/eta**2
290     if mu !=None: theta=1./(dt*mu)
291     for i in xrange(numMaterial):
292     inv_eta_i=(tau/tau_t[i])**s[i]/eta_N[i]
293     theta=theta+inv_eta_i
294     omega=omega+s[i]*inv_eta_i
295     #===========================================
296     eta_eff, eta_eff_old=util.clip(eta_eff*(theta+omega)/(eta_eff*theta**2+omega),maxval=eta_max), eta_eff
297     tau=eta_eff*gamma_dot
298     d=util.Lsup(eta_eff-eta_eff_old)
299     l=util.Lsup(eta_eff)
300     iter+=1
301     if self.__verbose: print "step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))
302     converged= d<= rtol* l
303     return eta_eff
304    
305     #====================================================================================================================================
306    
307 gross 2100 class IncompressibleIsotropicKelvinFlow(Defect):
308 gross 1639 """
309 caltinay 2169 This class implements the rheology of an isotropic Kelvin material.
310 gross 2100
311 caltinay 2169 Typical usage::
312 gross 1639
313 caltinay 2169 sp = IncompressibleIsotropicKelvinFlow(domain, stress=0, v=0)
314     sp.setTolerance()
315     sp.initialize(...)
316     v,p = sp.solve(v0, p0)
317 gross 1639
318 caltinay 2169 @note: This model has been used in the self-consistent plate-mantle model
319     proposed in U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>}
320     and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
321     I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
322     see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract?CRETRY=1&SRETRY=0>}.
323    
324 gross 1639 """
325 gross 2100 def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, useJaumannStress=True, **kwargs):
326 gross 1639 """
327 caltinay 2169 Initializes the model.
328 gross 1639
329     @param domain: problem domain
330     @type domain: L{domain}
331 gross 2100 @param stress: initial deviatoric stress
332 gross 1639 @param v: initial velocity field
333 gross 2100 @param p: initial pressure
334 gross 1639 @param t: initial time
335 caltinay 2169 @param useJaumannStress: C{True} if Jaumann stress is used
336     (not supported yet)
337 gross 1639 """
338 gross 2100 super(IncompressibleIsotropicKelvinFlow, self).__init__(**kwargs)
339 gross 1639 self.__domain=domain
340     self.__t=t
341     self.__vol=util.integrate(1.,Function(self.__domain))
342     self.__useJaumannStress=useJaumannStress
343 gross 2100 self.__v_boundary=Vector(0,Solution(self.__domain))
344 gross 1639 #=======================
345     # state variables:
346     #
347     if isinstance(stress,Data):
348     self.__stress=util.interpolate(stress,Function(domain))
349     else:
350     self.__stress=Data(stress,(domain.getDim(),domain.getDim()),Function(domain))
351     self.__stress-=util.trace(self.__stress)*(util.kronecker(domain)/domain.getDim())
352     if isinstance(v,Data):
353     self.__v=util.interpolate(v,Solution(domain))
354     else:
355 caltinay 2169 self.__v=Data(v,(domain.getDim(),),Solution(domain))
356 gross 1639 if isinstance(p,Data):
357     self.__p=util.interpolate(p,ReducedSolution(domain))
358     else:
359     self.__p=Data(p,(),ReducedSolution(domain))
360     #=======================
361 caltinay 2169 # PDE related stuff
362 gross 2100 self.__pde_v=LinearPDE(domain,numEquations=self.getDomain().getDim(),numSolutions=self.getDomain().getDim())
363     self.__pde_v.setSymmetryOn()
364     self.__pde_v.setSolverMethod(preconditioner=LinearPDE.RILU)
365 caltinay 2169
366 gross 2100 self.__pde_p=LinearPDE(domain)
367     self.__pde_p.setReducedOrderOn()
368     self.__pde_p.setSymmetryOn()
369 gross 1639
370 gross 2100 self.setTolerance()
371     self.setSmall()
372 caltinay 2169
373 gross 1639 def useJaumannStress(self):
374     """
375 caltinay 2169 Returns C{True} if Jaumann stress is included.
376     """
377 gross 1639 return self.__useJaumannStress
378 gross 2100
379 gross 1639 def setSmall(self,small=util.sqrt(util.EPSILON)):
380     """
381 caltinay 2169 Sets a small value to be used.
382    
383 gross 1639 @param small: positive small value
384     """
385     self.__small=small
386 gross 2100
387 gross 1639 def getSmall(self):
388     """
389 caltinay 2169 Returns small value.
390 gross 1639 @rtype: positive float
391     """
392     return self.__small
393 gross 2100
394     def setTolerance(self,tol=1.e-4):
395 gross 1639 """
396 caltinay 2169 Sets the tolerance.
397 gross 1639 """
398 gross 2100 self.__pde_v.setTolerance(tol**2)
399     self.__pde_p.setTolerance(tol**2)
400     self.__tol=tol
401    
402     def getTolerance(self):
403 gross 1639 """
404 caltinay 2169 Returns the set tolerance.
405 gross 2100 @rtype: positive float
406 gross 1639 """
407 gross 2100 return self.__tol
408    
409     def getDomain(self):
410 gross 1639 """
411 caltinay 2169 Returns the domain.
412 gross 1639 """
413 gross 2100 return self.__domain
414 caltinay 2169
415 gross 1639 def getTime(self):
416     """
417 caltinay 2169 Returns current time.
418 gross 1639 """
419     return self.__t
420 gross 2100
421     def setExternals(self, F=None, f=None, q=None, v_boundary=None):
422 caltinay 2169 """
423     Sets externals.
424 gross 1639
425 caltinay 2169 @param F: external force
426 gross 2100 @param f: surface force
427 gross 1639 @param q: location of constraints
428     """
429 gross 2100 if F != None: self.__pde_v.setValue(Y=F)
430     if f != None: self.__pde_v.setValue(y=f)
431     if q != None: self.__pde_v.setValue(q=q)
432     if v_boundary != None: self.__v_boundary=v_boundary
433 gross 1639
434 gross 2100 def bilinearform(self,arg0,arg1):
435     s0=util.deviatoric(util.symmetric(util.grad(arg0[0])))
436     s1=util.deviatoric(util.symmetric(util.grad(arg1[0])))
437     # s0=util.interpolate(arg0[0],Function(self.getDomain()))
438     # s1=util.interpolate(arg1[0],Function(self.getDomain()))
439     p0=util.interpolate(arg0[1],Function(self.getDomain()))
440     p1=util.interpolate(arg1[1],Function(self.getDomain()))
441     a=util.integrate(self.__p_weight**2*util.inner(s0,s1))+util.integrate(p0*p1)
442     return a
443    
444     def getEtaEff(self,strain, pressure):
445 caltinay 2169 if self.__mu==None:
446 gross 2100 eps=util.length(strain)*util.sqrt(2)
447 gross 1639 else:
448 gross 2100 eps=util.length(strain+self.__stress/((2*self.__dt)*self.__mu))*util.sqrt(2)
449     p=util.interpolate(pressure,eps.getFunctionSpace())
450     if self.__tau_Y!= None:
451     tmp=self.__tau_Y+self.__friction*p
452     m=util.wherePositive(eps)*util.wherePositive(tmp)
453     eta_max=m*tmp/(eps+(1-m)*util.EPSILON)+(1-m)*util.DBLE_MAX
454 gross 1639 else:
455 gross 2100 eta_max=util.DBLE_MAX
456     # initial guess:
457     tau=util.length(self.__stress)/util.sqrt(2)
458     # start the iteration:
459     cc=0
460     TOL=1e-7
461     dtau=util.DBLE_MAX
462     print "tau = ", tau, "eps =",eps
463     while cc<10 and dtau>TOL*util.Lsup(tau):
464     eta_eff2,eta_eff_dash=self.evalEtaEff(tau,return_dash=True)
465     eta_eff=util.clip(eta_eff2-eta_eff_dash*tau/(1-eta_eff_dash*eps),maxval=eta_max)
466     tau, tau2=eta_eff*eps, tau
467     dtau=util.Lsup(tau2-tau)
468     print "step ",cc,dtau, util.Lsup(tau)
469     cc+=1
470     return eta_eff
471 gross 1639
472 gross 2100 def getEtaCharacteristic(self):
473     a=0
474     for i in xrange(self.__numMaterials):
475     a=a+1./self.__eta_N[i]
476     return 1/a
477 caltinay 2169
478 gross 2100 def evalEtaEff(self, tau, return_dash=False):
479     a=Scalar(0,tau.getFunctionSpace()) # =1/eta
480     if return_dash: a_dash=Scalar(0,tau.getFunctionSpace()) # =(1/eta)'
481     s=util.Lsup(tau)
482     if s>0:
483     m=util.wherePositive(tau)
484     tau2=s*util.EPSILON*(1.-m)+m*tau
485     for i in xrange(self.__numMaterials):
486     eta_N=self.__eta_N[i]
487     tau_t=self.__tau_t[i]
488     if tau_t==None:
489     a+=1./eta_N
490     else:
491     power=1.-1./self.__power[i]
492     c=1./(tau_t**power*eta_N)
493     a+=c*tau2**power
494     if return_dash: a_dash+=power*c*tau2**(power-1.)
495     else:
496     for i in xrange(self.__numMaterials):
497     eta_N=self.__eta_N[i]
498     power=1.-1./self.__power[i]
499     a+=util.whereZero(power)/eta_N
500     if self.__mu!=None: a+=1./(self.__dt*self.__mu)
501     out=1/a
502     if return_dash:
503     return out,-out**2*a_dash
504     else:
505     return out
506 caltinay 2169
507 gross 2100 def eval(self,arg):
508     v=arg[0]
509     p=arg[1]
510     D=self.getDeviatoricStrain(v)
511     eta_eff=self.getEtaEff(D,p)
512     print "eta_eff=",eta_eff
513     # solve for dv
514     self.__pde_v.setValue(A=Data()) # save memory!
515     k3=util.kronecker(Function(self.getDomain()))
516     k3Xk3=util.outer(k3,k3)
517     self.__pde_v.setValue(A=eta_eff*(util.swap_axes(k3Xk3,0,3)+util.swap_axes(k3Xk3,1,3)),X=-eta_eff*D+p*util.kronecker(self.getDomain()))
518     dv=self.__pde_v.getSolution(verbose=self.__verbose)
519     print "resistep dv =",dv
520     # solve for dp
521     v2=v+dv
522     self.__pde_p.setValue(D=1/eta_eff,Y=util.div(v2))
523     dp=self.__pde_p.getSolution(verbose=self.__verbose)
524     print "resistep dp =",dp
525     return ArithmeticTuple(dv,dp)
526    
527 caltinay 2169 def update(self, dt, iter_max=100, inner_iter_max=20, verbose=False):
528 gross 1639 """
529 caltinay 2169 Updates stress, velocity and pressure for time increment dt.
530 gross 1639
531     @param dt: time increment
532 caltinay 2169 @param inner_iter_max: maximum number of iteration steps in the
533     incompressible solver
534     @param verbose: prints some infos in the incompressible solver
535 gross 1639 """
536 gross 2100 self.__verbose=verbose
537     self.__dt=dt
538     tol=self.getTolerance()
539     # set the initial velocity:
540     m=util.wherePositive(self.__pde_v.getCoefficient("q"))
541     v_new=self.__v*(1-m)+self.__v_boundary*m
542     # and off we go:
543     x=ArithmeticTuple(v_new, self.__p)
544     # self.__p_weight=util.interpolate(1./self.getEtaCharacteristic(),Function(self.__domain))**2
545     self.__p_weight=self.getEtaCharacteristic()
546     # self.__p_weight=util.interpolate(1./self.getEtaCharacteristic()**2,self.__p.getFunctionSpace())
547     atol=self.norm(x)*self.__tol
548     x_new=NewtonGMRES(self, x, iter_max=iter_max,sub_iter_max=inner_iter_max, atol=atol,rtol=0., verbose=verbose)
549     self.__v=x_new[0]
550     self.__p=x_new[1]
551     1/0
552     # self.__stress=self.getUpdatedStress(...)
553     self.__t+=dt
554     return self.__v, self.__p
555 gross 1639
556 caltinay 2169 #========================================================================
557 gross 1639
558 gross 2100 def getNewDeviatoricStress(self,D,eta_eff=None):
559     if eta_eff==None: eta_eff=self.evalEtaEff(self.__stress,D,self.__p)
560     s=(2*eta_eff)*D
561     if self.__mu!=None: s+=eta_eff/(self.__dt*self.__mu)*self.__last_stress
562     return s
563 gross 1639
564 gross 2100 def getDeviatoricStress(self):
565     """
566 caltinay 2169 Returns current stress.
567 gross 2100 """
568     return self.__stress
569 caltinay 2169
570 gross 2100 def getDeviatoricStrain(self,velocity=None):
571     """
572 caltinay 2169 Returns strain.
573 gross 2100 """
574     if velocity==None: velocity=self.getVelocity()
575     return util.deviatoric(util.symmetric(util.grad(velocity)))
576 gross 1659
577 gross 2100 def getPressure(self):
578     """
579 caltinay 2169 Returns current pressure.
580 gross 2100 """
581     return self.__p
582 caltinay 2169
583 gross 2100 def getVelocity(self):
584     """
585 caltinay 2169 Returns current velocity.
586 gross 2100 """
587     return self.__v
588 gross 1639
589 gross 2100 def getTau(self,stress=None):
590 gross 1639 """
591 caltinay 2169 Returns current second stress deviatoric invariant.
592 gross 1639 """
593 gross 2100 if stress==None: stress=self.getDeviatoricStress()
594     return util.sqrt(0.5)*util.length(stress)
595 gross 1639
596 gross 2100 def getGammaDot(self,strain=None):
597     """
598 caltinay 2169 Returns current second stress deviatoric invariant.
599 gross 2100 """
600     if strain==None: strain=self.getDeviatoricStrain()
601     return util.sqrt(2)*util.length(strain)
602 gross 1639

  ViewVC Help
Powered by ViewVC 1.1.26