/[escript]/branches/3.4.1/modellib/py_src/mechanics.py
ViewVC logotype

Contents of /branches/3.4.1/modellib/py_src/mechanics.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4587 - (show annotations)
Wed Dec 11 06:17:09 2013 UTC (5 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 12003 byte(s)
Preparation begins

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

  ViewVC Help
Powered by ViewVC 1.1.26