/[escript]/trunk/modellib/py_src/mechanics.py
ViewVC logotype

Contents of /trunk/modellib/py_src/mechanics.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26