/[escript]/release/5.2/modellib/py_src/mechanics.py
ViewVC logotype

Contents of /release/5.2/modellib/py_src/mechanics.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6692 - (show annotations)
Mon Jun 25 02:31:06 2018 UTC (3 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 12388 byte(s)
Fix

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

  ViewVC Help
Powered by ViewVC 1.1.26