/[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 814 by gross, Wed Aug 23 09:31:10 2006 UTC revision 821 by gross, Mon Aug 28 08:36:30 2006 UTC
# Line 15  class Mechanics(Model): Line 15  class Mechanics(Model):
15        base class for mechanics models in updated lagrangean framework        base class for mechanics models in updated lagrangean framework
16    
17        @ivar domain: domain (in)        @ivar domain: domain (in)
18        @ivar displacement: current displacements        @ivar internal_force: =Data()
19          @ivar external_force: =Data()
20          @ivar prescribed_velocity: =Data()
21          @ivar location_prescribed_velocity: =Data()
22          @ivar temperature:  = None
23          @ivar expansion_coefficient:  = 0.
24          @ivar bulk_modulus: =1.
25          @ivar shear_modulus: =1.
26          @ivar rel_tol: =1.e-3
27          @ivar abs_tol: =1.e-15
28          @ivar max_iter: =10
29          @ivar displacement: =None
30          @ivar stress: =None
31          @ivar velocity: =None
32        """        """
33        def __init__(self,debug=False):        def __init__(self,debug=False):
34           """           """
# Line 36  class Mechanics(Model): Line 48  class Mechanics(Model):
48                                 location_prescribed_velocity=Data(), \                                 location_prescribed_velocity=Data(), \
49                                 temperature = None, \                                 temperature = None, \
50                                 expansion_coefficient = 0., \                                 expansion_coefficient = 0., \
51                                 bulk_modulus=1., \                                 bulk_modulus=2., \
52                                 shear_modulus=1., \                                 shear_modulus=1., \
53                                 rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)                                 rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
54           self.__iter=0           self.__iter=0
# Line 49  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             self.__pde=LinearPDE(self.domain)             self.__pde=LinearPDE(self.domain)
64               self.__pde.setSymmetryOn()
65             self.__displacement_old=self.displacement             self.__displacement_old=self.displacement
66             self.stress_old=self.stress             self.stress_old=self.stress
67             self.__velocity_old=self.velocity             self.__velocity_old=self.velocity
# Line 95  class Mechanics(Model): Line 108  class Mechanics(Model):
108                                q=self.location_prescribed_velocity)                                q=self.location_prescribed_velocity)
109            # solve the PDE:            # solve the PDE:
110            self.__pde.setTolerance(self.rel_tol/100.)            self.__pde.setTolerance(self.rel_tol/100.)
111            self.velocity=self.__pde.getSolution()            self.velocity=self.__pde.getSolution(verbose=True)
112            # calculate convergence indicators:            # calculate convergence indicators:
113            self.__diff,diff_old=Lsup(self.velocity-self.__velocity_last),self.__diff            self.__diff,diff_old=Lsup(self.velocity-self.__velocity_last),self.__diff
114            self.__velocity_last=self.velocity            self.__velocity_last=self.velocity
115            self.displacement=self.__displacement_old+dt*self.velocity            self.displacement=self.__displacement_old+dt*self.velocity
116            self.__iter+=1            self.__iter+=1
117            self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))            if self.debug:
118                 for i in range(self.domain.getDim()):
119                    self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
120                 self.trace("velocity increment %s"%self.__diff)
121            if self.__iter>2 and diff_old<self.__diff:            if self.__iter>2 and diff_old<self.__diff:
122                raise IterationDivergenceError,"no improvement in stress iteration"                raise IterationDivergenceError,"no improvement in stress iteration"
123            if self.__iter>self.max_iter:            if self.__iter>self.max_iter:
# Line 145  class DruckerPrager(Mechanics): Line 161  class DruckerPrager(Mechanics):
161             self.declareParameter(plastic_stress=0.,             self.declareParameter(plastic_stress=0.,
162                                   friction_parameter=0.,                                   friction_parameter=0.,
163                                   dilatancy_parameter=0.,                                   dilatancy_parameter=0.,
164                                   shear_length=0.)                                   shear_length=1.e15)
165    
166        def doInitialization(self):        def doInitialization(self):
167             """             """
# Line 153  class DruckerPrager(Mechanics): Line 169  class DruckerPrager(Mechanics):
169             """             """
170             super(DruckerPrager, self).doInitialization()             super(DruckerPrager, self).doInitialization()
171             self.__plastic_stress_old=self.plastic_stress             self.__plastic_stress_old=self.plastic_stress
172             self.__tau_y_old=tau_y=self.shear_length             self.__tau_y_old=self.shear_length
173    
174        def doStepPreprocessing(self,dt):        def doStepPreprocessing(self,dt):
175              """              """
# Line 164  class DruckerPrager(Mechanics): Line 180  class DruckerPrager(Mechanics):
180              """              """
181              super(DruckerPrager, self).doStepPreprocessing(dt)              super(DruckerPrager, self).doStepPreprocessing(dt)
182              self.plastic_stress=self.__plastic_stress_old              self.plastic_stress=self.__plastic_stress_old
             self.tau_y=self.__tau_y_old  
183    
184        def doStep(self,dt):        def doStep(self,dt):
185             G=self.shear_modulus             G=self.shear_modulus
# Line 174  class DruckerPrager(Mechanics): Line 189  class DruckerPrager(Mechanics):
189             tau_Y=self.shear_length             tau_Y=self.shear_length
190             if self.__plastic_stress_old:             if self.__plastic_stress_old:
191                dps=self.plastic_stress-self.__plastic_stress_old                dps=self.plastic_stress-self.__plastic_stress_old
192                h=(tau_Y-self.tau_Y_last)/(dps+self.abs_tol*whereZero(dps))                h=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))
193             else:             else:
194                h=0                h=0
195             # set new tangential operator:             # set new tangential operator:
# Line 191  class DruckerPrager(Mechanics): Line 206  class DruckerPrager(Mechanics):
206        def doStepPostprocessing(self,dt):        def doStepPostprocessing(self,dt):
207            super(DruckerPrager, self).doStepPostprocessing(dt)            super(DruckerPrager, self).doStepPostprocessing(dt)
208            self.plastic_stress=self.__plastic_stress_old=self.plastic_stress            self.plastic_stress=self.__plastic_stress_old=self.plastic_stress
           self.shear_length=self.__tau_y_old  
209    
210        def getNewStress(self,s,gamma_p,du,ds_therm,tau_Y,G,K,alpha,beta,h):        def getNewStress(self,s,gamma_p,du,ds_therm,tau_Y,G,K,alpha,beta,h):
211              k3=kronecker(self.domain)              k3=kronecker(self.domain)
# Line 225  class DruckerPrager(Mechanics): Line 239  class DruckerPrager(Mechanics):
239             S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \             S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \
240                 + (K-2./3*G)*k3Xk3 \                 + (K-2./3*G)*k3Xk3 \
241                 + sXk3-swap_axes(sXk3,1,3) \                 + sXk3-swap_axes(sXk3,1,3) \
242                 + 1./2*(swap_axes(sXk3,0,3)+swap_axes(sXk3,1,2) \                 + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \
243                        -swap_axes(sXk3,1,3)-swap_axes(sXk3,0,2))                        -swap_axes(sXk3,1,3)+swap_axes(sXk3,0,2))
244                 # - chi/(h+G+alpha*beta*K)*outer(tmp+beta*K*k3,tmp+alpha*K*k3)\                 # - chi/(h+G+alpha*beta*K)*outer(tmp+beta*K*k3,tmp+alpha*K*k3)\
245             return S             return S
Legend:
Removed from v.814  
changed lines
  Added in v.821
  ViewVC Help
Powered by ViewVC 1.1.26