/[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 918 - (show annotations)
Wed Jan 3 06:30:00 2007 UTC (12 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 11279 byte(s)
fixes on ESySXML to get runmodel going.

    * object ids are now local for read and write of XML
    * ParameterSets are now written with class name
    * ParameterSets linked by other ParameterSets (previously only Models) are written now 
    * list are now lists of str (rather than bools), lists with bool, int or float are mapped to numarray
    * numarray are writen with generic types Bool, Int, Float (portability!)


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_modulus*self.deps_th*k3)
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_tol*self.SAFTY_FACTOR_ITERATION*s_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+K*self.deps_th*k3+ \
236 2*G*D+(K-2./3*G)*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.__chi*F/(h+G+beta*K)
246 self.__tau=tau_e-G*l
247 self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev-(p_e+l*beta*K)*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=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))
268
269 self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \
270 + (K-2./3*G)*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+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)

  ViewVC Help
Powered by ViewVC 1.1.26