/[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 831 by gross, Thu Aug 31 08:47:40 2006 UTC revision 832 by gross, Fri Sep 1 02:04:46 2006 UTC
# Line 61  class Mechanics(Model): Line 61  class Mechanics(Model):
61             if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))             if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
62             if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))             if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
63             # save the old values:             # save the old values:
64             self.__stress_old=self.stress             self.__stress_safe=self.stress
65             self.__temperature_old=self.temperature             self.__temperature_safe=self.temperature
66             self.__displacement_old=self.displacement             self.__displacement_safe=self.displacement
67               self.__velocity_safe=self.velocity
68               self.__velocity_last=None
69             # get node cooridnates and apply initial displacement             # get node cooridnates and apply initial displacement
70             self.__x=self.domain.getX()             self.__x=self.domain.getX()
71             self.domain.setX(self.__x+self.displacement)             self.domain.setX(self.__x+self.displacement)
# Line 83  class Mechanics(Model): Line 85  class Mechanics(Model):
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 values for the iteration:
88              self.temperature=self.__temperature_old              self.displacement=self.__displacement_safe
89              self.displacement=self.__displacement_old              self.stress=self.__stress_safe
90              self.stress=self.__stress_old              self.velocity=self.__velocity_safe
91                self.__velocity_last=self.__velocity_safe
92                print "set self.__velocity_last",dt
93              # update geometry              # update geometry
94              self.domain.setX(self.__x+self.displacement)              self.domain.setX(self.__x+self.displacement)
95    
# Line 96  class Mechanics(Model): Line 100  class Mechanics(Model):
100            k3=kronecker(self.domain)            k3=kronecker(self.domain)
101            # set new thermal stress increment            # set new thermal stress increment
102            if self.temperature:            if self.temperature:
103               self.deps_th=self.self.expansion_coefficient*(self.temperature-self.__temperature_old)               self.deps_th=self.self.expansion_coefficient*(self.temperature-self.__temperature_safe)
104            else:            else:
105               self.deps_th=0.               self.deps_th=0.
106            # set PDE coefficients:            # set PDE coefficients:
# Line 113  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
120              self.velocity=(self.displacement-self.__displacement_safe)/dt
121            self.domain.setX(self.__x+self.displacement)            self.domain.setX(self.__x+self.displacement)
122    
123            if self.debug:            if self.debug:
# Line 120  class Mechanics(Model): Line 125  class Mechanics(Model):
125                  self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i])))                  self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i])))
126               for i in range(self.domain.getDim()):               for i in range(self.domain.getDim()):
127                  self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))                  self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))
128                 for i in range(self.domain.getDim()):
129                    self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
130            self.__stress_last=self.stress            self.__stress_last=self.stress
131    
132        def terminateIteration(self):        def terminateIteration(self):
133            """iteration is terminateIterationd if relative pressure change is less then rel_tol"""            """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
134            if self.__iter>self.max_iter:            if self.__iter>self.max_iter:
135                raise IterationDivergenceError,"Maximum number of iterations steps reached"                raise IterationDivergenceError,"Maximum number of iterations steps reached"
           print self.__iter  
136            if self.__iter==0:            if self.__iter==0:
137               self.__diff=self.UNDEF_DT               self.__diff=self.UNDEF_DT
138            else:            else:
139               self.__diff,diff_old=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_old<self.__diff:               if 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 141  class Mechanics(Model): Line 147  class Mechanics(Model):
147             """             """
148             accept all the values:             accept all the values:
149             """             """
150             self.__displacement_old=self.displacement             self.__displacement_safe=self.displacement
151             self.__temperature_old=self.temperature             self.__temperature_safe=self.temperature
152             self.__stress_old=self.stress             self.__stress_safe=self.stress
153               self.__velocity_safe=self.velocity
154    
155        def getSafeTimeStepSize(self,dt):        def getSafeTimeStepSize(self,dt):
156             """             """
157             returns new step size             returns new step size
            d=Lsup(self.velocity-self.__velocity_old)*Lsup(self.displacement)  
            if d>0:  
               return dt/d*self.rel_tol  
            else:  
158             """             """
159             return self.UNDEF_DT             if self.__velocity_last:
160                  a=Lsup(self.velocity-self.__velocity_last)/dt
161                  if a>0:
162                     print "new dt:= ",Lsup(self.displacement)/a*self.rel_tol,dt,a
163                     # return Lsup(self.displacement)/a*self.rel_tol
164                     return self.UNDEF_DT
165                  else:
166                     return self.UNDEF_DT
167               else:
168                  return self.UNDEF_DT
169    
170    
171    
# Line 176  class DruckerPrager(Mechanics): Line 188  class DruckerPrager(Mechanics):
188            """            """
189            """            """
190            super(DruckerPrager, self).doInitialization()            super(DruckerPrager, self).doInitialization()
191            self.__plastic_stress_old=self.plastic_stress            self.__plastic_stress_safe=self.plastic_stress
192            self.__shear_length_old=self.shear_length            self.__shear_length_safe=self.shear_length
193            self.__hardening_old=self.hardening            self.__hardening_safe=self.hardening
194            self.__chi_old=0            self.__chi_safe=0
195            self.__tau_old=0            self.__tau_safe=0
196    
197        def doStepPreprocessing(self,dt):        def doStepPreprocessing(self,dt):
198            """            """
199            """            """
200            super(DruckerPrager, self).doStepPreprocessing(dt)            super(DruckerPrager, self).doStepPreprocessing(dt)
201            # set initial guess for iteration:            # set initial guess for iteration:
202            self.shear_length=self.__shear_length_old            self.shear_length=self.__shear_length_safe
203            self.plastic_stress=self.__plastic_stress_old            self.plastic_stress=self.__plastic_stress_safe
204            self.hardening=self.__hardening_old            self.hardening=self.__hardening_safe
205            self.__chi=self.__chi_old            self.__chi=self.__chi_safe
206            self.__tau=self.__tau_old            self.__tau=self.__tau_safe
207    
208        def doStep(self,dt):        def doStep(self,dt):
209            # set new tangential operator:            # set new tangential operator:
# Line 203  class DruckerPrager(Mechanics): Line 215  class DruckerPrager(Mechanics):
215    
216        def doStepPostprocessing(self,dt):        def doStepPostprocessing(self,dt):
217            super(DruckerPrager, self).doStepPostprocessing(dt)            super(DruckerPrager, self).doStepPostprocessing(dt)
218            self.__plastic_stress_old=self.plastic_stress            self.__plastic_stress_safe=self.plastic_stress
219            self.__shear_length_old=self.shear_length            self.__shear_length_safe=self.shear_length
220            self.__hardening_old=self.hardening            self.__hardening_safe=self.hardening
221            self.__chi_old=self.__chi            self.__chi_safe=self.__chi
222            self.__tau_old=self.__tau            self.__tau_safe=self.__tau
223    
224        def setStress(self):        def setStress(self):
225             d=self.domain.getDim()             d=self.domain.getDim()
# Line 233  class DruckerPrager(Mechanics): Line 245  class DruckerPrager(Mechanics):
245             # plastic stress increment:             # plastic stress increment:
246             l=self.__chi*F/(h+G+beta*K)             l=self.__chi*F/(h+G+beta*K)
247             self.__tau=tau_e-G*l             self.__tau=tau_e-G*l
248             self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev+(p_e+l*beta*K)*k3             self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev-(p_e+l*beta*K)*k3
249             self.plastic_stress=self.plastic_stress+l             self.plastic_stress=self.plastic_stress+l
250             # update hardening             # update hardening
251             self.hardening=(self.shear_length-self.__shear_length_old)/(l+self.abs_tol*whereZero(l))             self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l))
252    
253        def setTangentialTensor(self):        def setTangentialTensor(self):
254             d=self.domain.getDim()             d=self.domain.getDim()
# Line 248  class DruckerPrager(Mechanics): Line 260  class DruckerPrager(Mechanics):
260             chi=self.__chi             chi=self.__chi
261             tau=self.__tau             tau=self.__tau
262             h=self.hardening             h=self.hardening
   
263             k3=kronecker(Function(self.domain))             k3=kronecker(Function(self.domain))
264    
265               ###
266               p_e=-1./d*trace(self.stress)
267               s_e_dev=self.stress+p_e*k3
268               tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
269               F=tau_e-alpha*p_e-self.shear_length
270               # print F
271               ###
272    
273    
274             sXk3=outer(self.stress,k3)             sXk3=outer(self.stress,k3)
275             k3Xk3=outer(k3,k3)             k3Xk3=outer(k3,k3)
276             s_dev=self.stress-trace(self.stress)*(k3/d)             s_dev=self.stress-trace(self.stress)*(k3/d)
Legend:
Removed from v.831  
changed lines
  Added in v.832
  ViewVC Help
Powered by ViewVC 1.1.26