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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 833 by gross, Fri Sep 1 02:04:46 2006 UTC revision 834 by gross, Fri Sep 1 09:05:16 2006 UTC
# Line 84  class Mechanics(Model): Line 84  class Mechanics(Model):
84              # reset iteration counters:              # reset iteration counters:
85              self.__iter=0              self.__iter=0
86              self.__diff=self.UNDEF_DT              self.__diff=self.UNDEF_DT
87              # set initial values for the iteration:              # set initial guesses for the iteration:
88              self.displacement=self.__displacement_safe              self.displacement=self.__displacement_safe
89              self.stress=self.__stress_safe              self.stress=self.__stress_safe
90              self.velocity=self.__velocity_safe              self.velocity=self.__velocity_safe
# Line 117  class Mechanics(Model): Line 117  class Mechanics(Model):
117            self.du=self.__pde.getSolution(verbose=True)            self.du=self.__pde.getSolution(verbose=True)
118            # update geometry            # update geometry
119            self.displacement=self.displacement+self.du            self.displacement=self.displacement+self.du
           self.velocity=(self.displacement-self.__displacement_safe)/dt  
120            self.domain.setX(self.__x+self.displacement)            self.domain.setX(self.__x+self.displacement)
121              self.velocity=(self.displacement-self.__displacement_safe)/dt
122    
123            if self.debug:            if self.debug:
124               for i in range(self.domain.getDim()):               for i in range(self.domain.getDim()):
# Line 139  class Mechanics(Model): Line 139  class Mechanics(Model):
139               self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff               self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff
140               s_sup=Lsup(self.stress)               s_sup=Lsup(self.stress)
141               self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))               self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))
142               if diff_safe<self.__diff:               if self.__iter>2 and diff_safe<self.__diff:
143                   raise IterationDivergenceError,"no improvement in stress iteration"                   raise IterationDivergenceError,"no improvement in stress iteration"
144               return self.__diff<=self.rel_tol*self.SAFTY_FACTOR_ITERATION*s_sup+self.abs_tol               return self.__diff<=self.rel_tol*self.SAFTY_FACTOR_ITERATION*s_sup+self.abs_tol
145    
# Line 232  class DruckerPrager(Mechanics): Line 232  class DruckerPrager(Mechanics):
232             # elastic trial stress:             # elastic trial stress:
233             g=grad(self.du)             g=grad(self.du)
234             D=symmetric(g)             D=symmetric(g)
235               print trace(D)
236             W=nonsymmetric(g)             W=nonsymmetric(g)
237             s_e=self.stress+K*self.deps_th*k3+ \             s_e=self.stress+K*self.deps_th*k3+ \
238                        2*G*D+(K-2./3*G)*trace(D)*k3 \                        2*G*D+(K-2./3*G)*trace(D)*k3 \
239                        +2*symmetric(matrix_mult(W,self.stress))                        +2*symmetric(matrix_mult(W,self.stress))
240               print s_e
241             p_e=-1./d*trace(s_e)             p_e=-1./d*trace(s_e)
242             s_e_dev=s_e+p_e*k3             s_e_dev=s_e+p_e*k3
243             tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))             tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
244             # yield conditon for elastic trial stress:             # yield conditon for elastic trial stress:
245             F=tau_e-alpha*p_e-self.shear_length             F=tau_e-alpha*p_e-self.shear_length
246             self.__chi=whereNonNegative(F)             self.__chi=whereNonNegative(F)
247               print "chi =",Lsup(self.__chi)
248             # plastic stress increment:             # plastic stress increment:
249             l=self.__chi*F/(h+G+beta*K)             l=self.__chi*F/(h+G+beta*K)
250             self.__tau=tau_e-G*l             self.__tau=tau_e-G*l
# Line 275  class DruckerPrager(Mechanics): Line 278  class DruckerPrager(Mechanics):
278             k3Xk3=outer(k3,k3)             k3Xk3=outer(k3,k3)
279             s_dev=self.stress-trace(self.stress)*(k3/d)             s_dev=self.stress-trace(self.stress)*(k3/d)
280             tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))             tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))
281             self.S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \             self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \
282                   + (K-2./3*G)*k3Xk3 \                   + (K-2./3*G)*k3Xk3 \
283                   + (sXk3-swap_axes(sXk3,1,3)) \                   + (sXk3-swap_axes(sXk3,1,3)) \
284                   + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \                   + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \
285                        -swap_axes(sXk3,1,3)+swap_axes(sXk3,0,2)) \                        +swap_axes(sXk3,1,3)-swap_axes(sXk3,0,2)) \
286                   - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)                   - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)
Legend:
Removed from v.833  
changed lines
  Added in v.834
  ViewVC Help
Powered by ViewVC 1.1.26