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

Revision 2061 - (show annotations)
Wed Nov 19 03:40:21 2008 UTC (11 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 11748 byte(s)
```Fixing some warnings from epydoc.
Disabled c++ signatures in python docstrings.
Removed references to Bruce in epydoc and the users guide.
```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2003-2008 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2003-2008 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 20 __url__= 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: C{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.__pde.setSymmetryOn() 95 96 def doStepPreprocessing(self,dt): 97 """ 98 step up pressure iteration 99 100 if run within a time dependend problem extrapolation of pressure from previous time steps is used to 101 get an initial guess (that needs some work!!!!!!!) 102 """ 103 # reset iteration counters: 104 self.__iter=0 105 self.__diff=self.UNDEF_DT 106 # set initial guesses for the iteration: 107 self.displacement=self.__displacement_safe 108 self.stress=self.__stress_safe 109 self.velocity=self.__velocity_safe 110 # update geometry 111 self.domain.setX(self.__x+self.displacement) 112 113 def doStep(self,dt): 114 """ 115 """ 116 self.__iter+=1 117 k3=kronecker(self.domain) 118 # set new thermal stress increment 119 if self.temperature == None: 120 self.deps_th=0. 121 else: 122 self.deps_th=self.expansion_coefficient*(self.temperature-self.__temperature_safe) 123 # set PDE coefficients: 124 self.__pde.setValue(A=self.S) 125 self.__pde.setValue(X=-self.stress-self.bulk_modulus*self.deps_th*k3) 126 if self.internal_force: self.__pde.setValue(Y=self.internal_force) 127 if self.external_force: self.__pde.setValue(y=self.external_force) 128 self.__pde.setValue(q=self.location_prescribed_velocity, \ 129 r=Data()) 130 if not self.prescribed_velocity.isEmpty() and self.__iter==1: 131 self.__pde.setValue(r=dt*self.prescribed_velocity) 132 # solve the PDE: 133 self.__pde.setTolerance(self.rel_tol**2) 134 self.du=self.__pde.getSolution(verbose=self.debug) 135 # update geometry 136 self.displacement=self.displacement+self.du 137 self.domain.setX(self.__x+self.displacement) 138 self.velocity=(self.displacement-self.__displacement_safe)/dt 139 140 if self.debug: 141 for i in range(self.domain.getDim()): 142 self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i]))) 143 for i in range(self.domain.getDim()): 144 self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i]))) 145 for i in range(self.domain.getDim()): 146 self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i]))) 147 self.__stress_last=self.stress 148 149 def terminateIteration(self): 150 """iteration is terminateIterationd if relative pressure change is less then rel_tol""" 151 if self.__iter>self.max_iter: 152 raise IterationDivergenceError,"Maximum number of iterations steps reached" 153 if self.__iter==0: 154 self.__diff=self.UNDEF_DT 155 else: 156 self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff 157 s_sup=Lsup(self.stress) 158 self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff)) 159 if self.__iter>2 and diff_safe0: 178 return 1./a 179 else: 180 return self.UNDEF_DT 181 182 183 184 class DruckerPrager(Mechanics): 185 """ 186 187 """ 188 189 def __init__(self,**kwargs): 190 """ 191 set up model 192 """ 193 super(DruckerPrager, self).__init__(**kwargs) 194 self.declareParameter(plastic_stress=0., 195 hardening=0., 196 friction_parameter=0., 197 dilatancy_parameter=0., 198 shear_length=1.e15) 199 def doInitialization(self): 200 """ 201 """ 202 super(DruckerPrager, self).doInitialization() 203 self.__plastic_stress_safe=self.plastic_stress 204 self.__shear_length_safe=self.shear_length 205 self.__hardening_safe=self.hardening 206 self.__chi_safe=0 207 self.__tau_safe=0 208 209 def doStepPreprocessing(self,dt): 210 """ 211 """ 212 super(DruckerPrager, self).doStepPreprocessing(dt) 213 # set initial guess for iteration: 214 self.shear_length=self.__shear_length_safe 215 self.plastic_stress=self.__plastic_stress_safe 216 self.hardening=self.__hardening_safe 217 self.__chi=self.__chi_safe 218 self.__tau=self.__tau_safe 219 220 def doStep(self,dt): 221 # set new tangential operator: 222 self.setTangentialTensor() 223 # do the update step: 224 super(DruckerPrager, self).doStep(dt) 225 # update stresses: 226 self.setStress() 227 228 def doStepPostprocessing(self,dt): 229 super(DruckerPrager, self).doStepPostprocessing(dt) 230 self.__plastic_stress_safe=self.plastic_stress 231 self.__shear_length_safe=self.shear_length 232 self.__hardening_safe=self.hardening 233 self.__chi_safe=self.__chi 234 self.__tau_safe=self.__tau 235 236 def setStress(self): 237 d=self.domain.getDim() 238 G=self.shear_modulus 239 K=self.bulk_modulus 240 alpha=self.friction_parameter 241 beta=self.dilatancy_parameter 242 h=self.hardening 243 k3=kronecker(self.domain) 244 # elastic trial stress: 245 g=grad(self.du) 246 D=symmetric(g) 247 W=nonsymmetric(g) 248 s_e=self.stress+K*self.deps_th*k3+ \ 249 2*G*D+(K-2./3*G)*trace(D)*k3 \ 250 +2*symmetric(matrix_mult(W,self.stress)) 251 p_e=-1./d*trace(s_e) 252 s_e_dev=s_e+p_e*k3 253 tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev)) 254 # yield conditon for elastic trial stress: 255 F=tau_e-alpha*p_e-self.shear_length 256 self.__chi=whereNonNegative(F+(self.rel_tol*(self.SAFTY_FACTOR_ITERATION)**2)*self.shear_length) 257 # plastic stress increment: 258 l=self.__chi*F/(h+G+beta*K) 259 self.__tau=tau_e-G*l 260 self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev-(p_e+l*beta*K)*k3 261 self.plastic_stress=self.plastic_stress+l 262 # update hardening 263 self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l)) 264 265 def setTangentialTensor(self): 266 d=self.domain.getDim() 267 G=self.shear_modulus 268 K=self.bulk_modulus 269 alpha=self.friction_parameter 270 beta=self.dilatancy_parameter 271 tau_Y=self.shear_length 272 chi=self.__chi 273 tau=self.__tau 274 h=self.hardening 275 k3=kronecker(Function(self.domain)) 276 277 sXk3=outer(self.stress,k3) 278 k3Xk3=outer(k3,k3) 279 s_dev=self.stress-trace(self.stress)*(k3/d) 280 tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol)) 281 282 self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \ 283 + (K-2./3*G)*k3Xk3 \ 284 + (sXk3-swap_axes(swap_axes(sXk3,1,2),2,3)) \ 285 + 1./2*(swap_axes(swap_axes(sXk3,0,2),2,3) \ 286 -swap_axes(swap_axes(sXk3,0,3),2,3) \ 287 -swap_axes(sXk3,1,2) \ 288 +swap_axes(sXk3,1,3) ) \ 289 - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)

 ViewVC Help Powered by ViewVC 1.1.26