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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (hide annotations)
Mon Jan 20 03:37:18 2020 UTC (6 months, 2 weeks ago) by uqaeller
File MIME type: text/x-python
File size: 12458 byte(s)
Updated the copyright header.


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

  ViewVC Help
Powered by ViewVC 1.1.26