modellib/py_src/mechanics.py

# $Id:$

__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""

from esys.escript import *
from esys.escript.modelframe import Model,IterationDivergenceError
from esys.escript.linearPDEs import LinearPDE

class Mechanics(Model):
      """
      base class for mechanics models in updated lagrangean framework

      @ivar domain: domain (in)
      @ivar internal_force: =Data()
      @ivar external_force: =Data()
      @ivar prescribed_velocity: =Data()
      @ivar location_prescribed_velocity: =Data()
      @ivar temperature:  = None
      @ivar expansion_coefficient:  = 0.
      @ivar bulk_modulus: =1.
      @ivar shear_modulus: =1.
      @ivar rel_tol: =1.e-3
      @ivar abs_tol: =1.e-15
      @ivar max_iter: =10
      @ivar displacement: =None
      @ivar stress: =None
      """
      SAFTY_FACTOR_ITERATION=1./100.
      def __init__(self,**kwargs):
         """
         set up the model
         
         @param debug: debug flag
         @type debug: C{bool}
         """
         super(Mechanics, self).__init__(self,**kwargs)
         self.declareParameter(domain=None, \
                               displacement=None, \
                               stress=None, \
                               velocity=None, \
                               internal_force=None, \
                               external_force=None, \
                               prescribed_velocity=None, \
                               location_prescribed_velocity=None, \
                               temperature = None, \
                               expansion_coefficient = 0., \
                               bulk_modulus=2., \
                               shear_modulus=1., \
                               rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
         self.__iter=0

      def doInitialization(self):
           """
           initialize model
           """
           if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
           if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
           if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
           if not self.internal_force: self.internal_force = Data()
           if not self.external_force: self.external_force = Data()
           if not self.prescribed_velocity: self.prescribed_velocity = Data()
           if not self.location_prescribed_velocity: self.location_prescribed_velocity =Data()
           # save the old values:
           self.__stress_safe=self.stress
           self.__temperature_safe=self.temperature
           self.__displacement_safe=self.displacement
           self.__velocity_safe=self.velocity
           self.__velocity_old=None
           self.__old_dt=None
           self.__very_old_dt=None
           # get node cooridnates and apply initial displacement
           self.__x=self.domain.getX()
           self.domain.setX(self.__x+self.displacement)
           # open PDE:
           self.__pde=LinearPDE(self.domain)
           self.__pde.setSolverMethod(self.__pde.DIRECT)
           # self.__pde.setSymmetryOn()

      def doStepPreprocessing(self,dt):
            """
            step up pressure iteration

            if run within a time dependend problem extrapolation of pressure from previous time steps is used to
            get an initial guess (that needs some work!!!!!!!)
            """
            # reset iteration counters:
            self.__iter=0
            self.__diff=self.UNDEF_DT
            # set initial guesses for the iteration:
            self.displacement=self.__displacement_safe
            self.stress=self.__stress_safe
            self.velocity=self.__velocity_safe
            # update geometry
            self.domain.setX(self.__x+self.displacement)

      def doStep(self,dt):
          """
          """
          self.__iter+=1
          k3=kronecker(self.domain)
          # set new thermal stress increment
          if self.temperature == None:
             self.deps_th=0.
          else:
             self.deps_th=self.expansion_coefficient*(self.temperature-self.__temperature_safe)
          # set PDE coefficients:
          self.__pde.setValue(A=self.S)
          self.__pde.setValue(X=-self.stress-self.bulk_modulus*self.deps_th*k3)
          if self.internal_force: self.__pde.setValue(Y=self.internal_force)
          if self.external_force: self.__pde.setValue(y=self.external_force)
          self.__pde.setValue(q=self.location_prescribed_velocity, \
                              r=Data())
          if not self.prescribed_velocity.isEmpty() and self.__iter==1:
               self.__pde.setValue(r=dt*self.prescribed_velocity)
          # solve the PDE:
          self.__pde.setTolerance(self.rel_tol**2)
          self.du=self.__pde.getSolution(verbose=self.debug)
          # update geometry
          self.displacement=self.displacement+self.du
          self.domain.setX(self.__x+self.displacement)
          self.velocity=(self.displacement-self.__displacement_safe)/dt

          if self.debug:
             for i in range(self.domain.getDim()):
                self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i])))
             for i in range(self.domain.getDim()):
                self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))
             for i in range(self.domain.getDim()):
                self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
          self.__stress_last=self.stress

      def terminateIteration(self):
          """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
          if self.__iter>self.max_iter:
              raise IterationDivergenceError,"Maximum number of iterations steps reached"
          if self.__iter==0:
             self.__diff=self.UNDEF_DT
          else:
             self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff
             s_sup=Lsup(self.stress)
             self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))
             if self.__iter>2 and diff_safe<self.__diff:
                 raise IterationDivergenceError,"no improvement in stress iteration"
             return self.__diff<=self.rel_tol*self.SAFTY_FACTOR_ITERATION*s_sup+self.abs_tol

      def doStepPostprocessing(self,dt):
           """
           accept all the values:
           """
           self.__displacement_safe=self.displacement
           self.__temperature_safe=self.temperature
           self.__stress_safe=self.stress
           self.__velocity_safe=self.velocity

      def getSafeTimeStepSize(self,dt):
           """
           returns new step size
           """
           a=sup(length(self.velocity)/self.domain.getSize())
           if a>0:
              return 1./a
           else:
              return self.UNDEF_DT


class DruckerPrager(Mechanics):
      """

      """

      def __init__(self,**kwargs):
           """
           set up model
           """
           super(DruckerPrager, self).__init__(**kwargs)
           self.declareParameter(plastic_stress=0.,
                                 hardening=0.,
                                 friction_parameter=0.,
                                 dilatancy_parameter=0.,
                                 shear_length=1.e15)
      def doInitialization(self):
          """
          """
          super(DruckerPrager, self).doInitialization()
          self.__plastic_stress_safe=self.plastic_stress
          self.__shear_length_safe=self.shear_length
          self.__hardening_safe=self.hardening
          self.__chi_safe=0
          self.__tau_safe=0

      def doStepPreprocessing(self,dt):
          """
          """
          super(DruckerPrager, self).doStepPreprocessing(dt)
          # set initial guess for iteration:
          self.shear_length=self.__shear_length_safe
          self.plastic_stress=self.__plastic_stress_safe
          self.hardening=self.__hardening_safe
          self.__chi=self.__chi_safe
          self.__tau=self.__tau_safe

      def doStep(self,dt):
          # set new tangential operator:
          self.setTangentialTensor()
          # do the update step:
          super(DruckerPrager, self).doStep(dt)
          # update stresses:
          self.setStress()

      def doStepPostprocessing(self,dt):
          super(DruckerPrager, self).doStepPostprocessing(dt)
          self.__plastic_stress_safe=self.plastic_stress
          self.__shear_length_safe=self.shear_length
          self.__hardening_safe=self.hardening
          self.__chi_safe=self.__chi
          self.__tau_safe=self.__tau

      def setStress(self):
           d=self.domain.getDim()
           G=self.shear_modulus
           K=self.bulk_modulus
           alpha=self.friction_parameter
           beta=self.dilatancy_parameter
           h=self.hardening
           k3=kronecker(self.domain)
           # elastic trial stress:
           g=grad(self.du)
           D=symmetric(g)
           W=nonsymmetric(g)
           s_e=self.stress+K*self.deps_th*k3+ \
                      2*G*D+(K-2./3*G)*trace(D)*k3 \
                      +2*symmetric(matrix_mult(W,self.stress))
           p_e=-1./d*trace(s_e)
           s_e_dev=s_e+p_e*k3
           tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
           # yield conditon for elastic trial stress:
           F=tau_e-alpha*p_e-self.shear_length
           self.__chi=whereNonNegative(F+(self.rel_tol*(self.SAFTY_FACTOR_ITERATION)**2)*self.shear_length)
           # plastic stress increment:
           l=self.__chi*F/(h+G+beta*K)
           self.__tau=tau_e-G*l
           self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev-(p_e+l*beta*K)*k3
           self.plastic_stress=self.plastic_stress+l
           # update hardening
           self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l))

      def setTangentialTensor(self):
           d=self.domain.getDim()
           G=self.shear_modulus
           K=self.bulk_modulus
           alpha=self.friction_parameter
           beta=self.dilatancy_parameter
           tau_Y=self.shear_length
           chi=self.__chi
           tau=self.__tau
           h=self.hardening
           k3=kronecker(Function(self.domain))

           sXk3=outer(self.stress,k3)
           k3Xk3=outer(k3,k3)
           s_dev=self.stress-trace(self.stress)*(k3/d)
           tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))

           self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \
                 + (K-2./3*G)*k3Xk3 \
                 + (sXk3-swap_axes(swap_axes(sXk3,1,2),2,3))   \
                 + 1./2*(swap_axes(swap_axes(sXk3,0,2),2,3)   \
                        -swap_axes(swap_axes(sXk3,0,3),2,3)   \
                        -swap_axes(sXk3,1,2)                  \
                        +swap_axes(sXk3,1,3)                ) \
                 - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)
1	# $Id:$
2
3	__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
4	http://www.access.edu.au
5	Primary Business: Queensland, Australia"""
6	__license__="""Licensed under the Open Software License version 3.0
7	http://www.opensource.org/licenses/osl-3.0.php"""
8
9	from esys.escript import *
10	from esys.escript.modelframe import Model,IterationDivergenceError
11	from esys.escript.linearPDEs import LinearPDE
12
13	class Mechanics(Model):
14	"""
15	base class for mechanics models in updated lagrangean framework
16
17	@ivar domain: domain (in)
18	@ivar internal_force: =Data()
19	@ivar external_force: =Data()
20	@ivar prescribed_velocity: =Data()
21	@ivar location_prescribed_velocity: =Data()
22	@ivar temperature: = None
23	@ivar expansion_coefficient: = 0.
24	@ivar bulk_modulus: =1.
25	@ivar shear_modulus: =1.
26	@ivar rel_tol: =1.e-3
27	@ivar abs_tol: =1.e-15
28	@ivar max_iter: =10
29	@ivar displacement: =None
30	@ivar stress: =None
31	"""
32	SAFTY_FACTOR_ITERATION=1./100.
33	def __init__(self,**kwargs):
34	"""
35	set up the model
36
37	@param debug: debug flag
38	@type debug: C{bool}
39	"""
40	super(Mechanics, self).__init__(self,**kwargs)
41	self.declareParameter(domain=None, \
42	displacement=None, \
43	stress=None, \
44	velocity=None, \
45	internal_force=None, \
46	external_force=None, \
47	prescribed_velocity=None, \
48	location_prescribed_velocity=None, \
49	temperature = None, \
50	expansion_coefficient = 0., \
51	bulk_modulus=2., \
52	shear_modulus=1., \
53	rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
54	self.__iter=0
55
56	def doInitialization(self):
57	"""
58	initialize model
59	"""
60	if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
61	if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
62	if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
63	if not self.internal_force: self.internal_force = Data()
64	if not self.external_force: self.external_force = Data()
65	if not self.prescribed_velocity: self.prescribed_velocity = Data()
66	if not self.location_prescribed_velocity: self.location_prescribed_velocity =Data()
67	# save the old values:
68	self.__stress_safe=self.stress
69	self.__temperature_safe=self.temperature
70	self.__displacement_safe=self.displacement
71	self.__velocity_safe=self.velocity
72	self.__velocity_old=None
73	self.__old_dt=None
74	self.__very_old_dt=None
75	# get node cooridnates and apply initial displacement
76	self.__x=self.domain.getX()
77	self.domain.setX(self.__x+self.displacement)
78	# open PDE:
79	self.__pde=LinearPDE(self.domain)
80	self.__pde.setSolverMethod(self.__pde.DIRECT)
81	# self.__pde.setSymmetryOn()
82
83	def doStepPreprocessing(self,dt):
84	"""
85	step up pressure iteration
86
87	if run within a time dependend problem extrapolation of pressure from previous time steps is used to
88	get an initial guess (that needs some work!!!!!!!)
89	"""
90	# reset iteration counters:
91	self.__iter=0
92	self.__diff=self.UNDEF_DT
93	# set initial guesses for the iteration:
94	self.displacement=self.__displacement_safe
95	self.stress=self.__stress_safe
96	self.velocity=self.__velocity_safe
97	# update geometry
98	self.domain.setX(self.__x+self.displacement)
99
100	def doStep(self,dt):
101	"""
102	"""
103	self.__iter+=1
104	k3=kronecker(self.domain)
105	# set new thermal stress increment
106	if self.temperature == None:
107	self.deps_th=0.
108	else:
109	self.deps_th=self.expansion_coefficient*(self.temperature-self.__temperature_safe)
110	# set PDE coefficients:
111	self.__pde.setValue(A=self.S)
112	self.__pde.setValue(X=-self.stress-self.bulk_modulusself.deps_thk3)
113	if self.internal_force: self.__pde.setValue(Y=self.internal_force)
114	if self.external_force: self.__pde.setValue(y=self.external_force)
115	self.__pde.setValue(q=self.location_prescribed_velocity, \
116	r=Data())
117	if not self.prescribed_velocity.isEmpty() and self.__iter==1:
118	self.__pde.setValue(r=dt*self.prescribed_velocity)
119	# solve the PDE:
120	self.__pde.setTolerance(self.rel_tol**2)
121	self.du=self.__pde.getSolution(verbose=self.debug)
122	# update geometry
123	self.displacement=self.displacement+self.du
124	self.domain.setX(self.__x+self.displacement)
125	self.velocity=(self.displacement-self.__displacement_safe)/dt
126
127	if self.debug:
128	for i in range(self.domain.getDim()):
129	self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i])))
130	for i in range(self.domain.getDim()):
131	self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))
132	for i in range(self.domain.getDim()):
133	self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
134	self.__stress_last=self.stress
135
136	def terminateIteration(self):
137	"""iteration is terminateIterationd if relative pressure change is less then rel_tol"""
138	if self.__iter>self.max_iter:
139	raise IterationDivergenceError,"Maximum number of iterations steps reached"
140	if self.__iter==0:
141	self.__diff=self.UNDEF_DT
142	else:
143	self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff
144	s_sup=Lsup(self.stress)
145	self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))
146	if self.__iter>2 and diff_safe<self.__diff:
147	raise IterationDivergenceError,"no improvement in stress iteration"
148	return self.__diff<=self.rel_tolself.SAFTY_FACTOR_ITERATIONs_sup+self.abs_tol
149
150	def doStepPostprocessing(self,dt):
151	"""
152	accept all the values:
153	"""
154	self.__displacement_safe=self.displacement
155	self.__temperature_safe=self.temperature
156	self.__stress_safe=self.stress
157	self.__velocity_safe=self.velocity
158
159	def getSafeTimeStepSize(self,dt):
160	"""
161	returns new step size
162	"""
163	a=sup(length(self.velocity)/self.domain.getSize())
164	if a>0:
165	return 1./a
166	else:
167	return self.UNDEF_DT
168
169
170
171	class DruckerPrager(Mechanics):
172	"""
173
174	"""
175
176	def __init__(self,**kwargs):
177	"""
178	set up model
179	"""
180	super(DruckerPrager, self).__init__(**kwargs)
181	self.declareParameter(plastic_stress=0.,
182	hardening=0.,
183	friction_parameter=0.,
184	dilatancy_parameter=0.,
185	shear_length=1.e15)
186	def doInitialization(self):
187	"""
188	"""
189	super(DruckerPrager, self).doInitialization()
190	self.__plastic_stress_safe=self.plastic_stress
191	self.__shear_length_safe=self.shear_length
192	self.__hardening_safe=self.hardening
193	self.__chi_safe=0
194	self.__tau_safe=0
195
196	def doStepPreprocessing(self,dt):
197	"""
198	"""
199	super(DruckerPrager, self).doStepPreprocessing(dt)
200	# set initial guess for iteration:
201	self.shear_length=self.__shear_length_safe
202	self.plastic_stress=self.__plastic_stress_safe
203	self.hardening=self.__hardening_safe
204	self.__chi=self.__chi_safe
205	self.__tau=self.__tau_safe
206
207	def doStep(self,dt):
208	# set new tangential operator:
209	self.setTangentialTensor()
210	# do the update step:
211	super(DruckerPrager, self).doStep(dt)
212	# update stresses:
213	self.setStress()
214
215	def doStepPostprocessing(self,dt):
216	super(DruckerPrager, self).doStepPostprocessing(dt)
217	self.__plastic_stress_safe=self.plastic_stress
218	self.__shear_length_safe=self.shear_length
219	self.__hardening_safe=self.hardening
220	self.__chi_safe=self.__chi
221	self.__tau_safe=self.__tau
222
223	def setStress(self):
224	d=self.domain.getDim()
225	G=self.shear_modulus
226	K=self.bulk_modulus
227	alpha=self.friction_parameter
228	beta=self.dilatancy_parameter
229	h=self.hardening
230	k3=kronecker(self.domain)
231	# elastic trial stress:
232	g=grad(self.du)
233	D=symmetric(g)
234	W=nonsymmetric(g)
235	s_e=self.stress+Kself.deps_thk3+ \
236	2GD+(K-2./3G)trace(D)*k3 \
237	+2*symmetric(matrix_mult(W,self.stress))
238	p_e=-1./d*trace(s_e)
239	s_e_dev=s_e+p_e*k3
240	tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
241	# yield conditon for elastic trial stress:
242	F=tau_e-alpha*p_e-self.shear_length
243	self.__chi=whereNonNegative(F+(self.rel_tol(self.SAFTY_FACTOR_ITERATION)2)self.shear_length)
244	# plastic stress increment:
245	l=self.__chiF/(h+G+betaK)
246	self.__tau=tau_e-G*l
247	self.stress=self.__tau/(tau_e+self.abs_tolwhereZero(tau_e,self.abs_tol))s_e_dev-(p_e+lbetaK)*k3
248	self.plastic_stress=self.plastic_stress+l
249	# update hardening
250	self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l))
251
252	def setTangentialTensor(self):
253	d=self.domain.getDim()
254	G=self.shear_modulus
255	K=self.bulk_modulus
256	alpha=self.friction_parameter
257	beta=self.dilatancy_parameter
258	tau_Y=self.shear_length
259	chi=self.__chi
260	tau=self.__tau
261	h=self.hardening
262	k3=kronecker(Function(self.domain))
263
264	sXk3=outer(self.stress,k3)
265	k3Xk3=outer(k3,k3)
266	s_dev=self.stress-trace(self.stress)*(k3/d)
267	tmp=Gs_dev/(tau+self.abs_tolwhereZero(tau,self.abs_tol))
268
269	self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \
270	+ (K-2./3G)k3Xk3 \
271	+ (sXk3-swap_axes(swap_axes(sXk3,1,2),2,3)) \
272	+ 1./2*(swap_axes(swap_axes(sXk3,0,2),2,3) \
273	-swap_axes(swap_axes(sXk3,0,3),2,3) \
274	-swap_axes(sXk3,1,2) \
275	+swap_axes(sXk3,1,3) ) \
276	- outer(chi/(h+G+alphabetaK)(tmp+betaKk3),tmp+alphaK*k3)