modellib/py_src/mechanics.py

# $Id:$

__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""

from esys.escript import *
from esys.escript.modelframe import Model,IterationDivergenceError
from esys.escript.linearPDEs import LinearPDE

class Mechanics(Model):
      """
      base class for mechanics models in updated lagrangean framework

      @ivar domain: domain (in)
      @ivar internal_force=Data()
      @ivar external_force=Data()
      @ivar prescribed_velocity=Data()
      @ivar location_prescribed_velocity=Data()
      @ivar temperature = None
      @ivar expansion_coefficient = 0.
      @ivar bulk_modulus=1.
      @ivar shear_modulus=1.
      @ivar rel_tol=1.e-3
      @ivar abs_tol=1.e-15
      @ivar max_iter=10
      @ivar displacement=None
      @ivar stress=None
      @ivar velocity=None
      """
      def __init__(self,debug=False):
         """
         set up the model
         
         @param debug: debug flag
         @type debug: C{bool}
         """
         super(Mechanics, self).__init__(self,debug=debug)
         self.declareParameter(domain=None, \
                               displacement=None, \
                               stress=None, \
                               velocity=None, \
                               internal_force=Data(), \
                               external_force=Data(), \
                               prescribed_velocity=Data(), \
                               location_prescribed_velocity=Data(), \
                               temperature = None, \
                               expansion_coefficient = 0., \
                               bulk_modulus=1., \
                               shear_modulus=1., \
                               rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
         self.__iter=0

      def doInitialization(self):
           """
           initialize model
           """
           if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
           if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
           if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
           self.__pde=LinearPDE(self.domain)
           self.__displacement_old=self.displacement
           self.stress_old=self.stress
           self.__velocity_old=self.velocity
           self.__temperature_old=self.temperature
           
      def doStepPreprocessing(self,dt):
            """
            step up pressure iteration

            if run within a time dependend problem extrapolation of pressure from previous time steps is used to
            get an initial guess (that needs some work!!!!!!!)
            """
            self.__iter=0
            self.__diff=self.UNDEF_DT
            # set new values:
            self.displacement=self.__displacement_old
            self.stress=self.stress_old
            self.velocity=self.__velocity_old
            self.temperature=self.__temperature_old
            self.__velocity_last=self.velocity

      def doStep(self,dt):
          """

          performs an iteration step of the penalty method.
          IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.
     
          requires self.S to be set 
          updates the thermal stress increment

          """
          k3=kronecker(self.domain)
          # set new thermal stress increment
          if self.temperature:
             self.dthermal_stress=self.bulk_modulus*self.self.expansion_coefficient*(self.temperature-self.__temperature_old)
          else:
             self.dthermal_stress=0.
          # set PDE coefficients:
          self.__pde.setValue(A=self.S)
          self.__pde.setValue(X=self.stress_old-self.dthermal_stress*k3)
          if self.internal_force: self.__pde.setValue(Y=self.internal_force)
          if self.external_force: self.__pde.setValue(y=self.external_force)
          print self.prescribed_velocity
          print self.location_prescribed_velocity
          self.__pde.setValue(r=self.prescribed_velocity, \
                              q=self.location_prescribed_velocity)
          # solve the PDE:
          self.__pde.setTolerance(self.rel_tol/100.)
          self.velocity=self.__pde.getSolution()
          # calculate convergence indicators:
          self.__diff,diff_old=Lsup(self.velocity-self.__velocity_last),self.__diff
          self.__velocity_last=self.velocity
          self.displacement=self.__displacement_old+dt*self.velocity
          self.__iter+=1
          self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
          if self.__iter>2 and diff_old<self.__diff:
              raise IterationDivergenceError,"no improvement in stress iteration"
          if self.__iter>self.max_iter:
              raise IterationDivergenceError,"Maximum number of iterations steps reached"

      def terminateIteration(self):
          """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
          return self.__diff<=self.rel_tol*Lsup(self.velocity)+self.abs_tol

      def doStepPostprocessing(self,dt):
           """
           accept all the values:
           """
           self.displacement=self.__displacement_old
           self.stress=self.stress_old
           self.velocity=self.__velocity_old
           self.temperature=self.__temperature_old

      def getSafeTimeStepSize(self,dt):
           """
           returns new step size
           """
           d=Lsup(self.velocity-self.__velocity_old)/dt
           if d>0:
               return Lsup(self.displacement)/d
           else:
               return self.UNDEF_DT


class DruckerPrager(Mechanics):
      """

      """

      def __init__(self,debug=False):
           """
           set up model
           """
           super(DruckerPrager, self).__init__(debug=debug)
           self.declareParameter(plastic_stress=0.,
                                 friction_parameter=0.,
                                 dilatancy_parameter=0.,
                                 shear_length=0.)

      def doInitialization(self):
           """
           initialize model
           """
           super(DruckerPrager, self).doInitialization()
           self.__plastic_stress_old=self.plastic_stress
           self.__tau_y_old=self.shear_length

      def doStepPreprocessing(self,dt):
            """
            step up pressure iteration

            if run within a time dependend problem extrapolation of pressure from previous time steps is used to
            get an initial guess (that needs some work!!!!!!!)
            """
            super(DruckerPrager, self).doStepPreprocessing(dt)
            self.plastic_stress=self.__plastic_stress_old

      def doStep(self,dt):
           G=self.shear_modulus
           K=self.bulk_modulus
           alpha=self.friction_parameter
           beta=self.dilatancy_parameter
           tau_Y=self.shear_length
           if self.__plastic_stress_old:
              dps=self.plastic_stress-self.__plastic_stress_old
              h=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))
           else:
              h=0
           # set new tangential operator:
           self.S=self.getTangentialTensor(self.stress,
                                           tau_Y,G,K,alpha,beta,h)
           # do the update step:
           super(DruckerPrager, self).doStep(dt)

           # update stresses:
           self.stress,self.plastic_stress=self.getNewStress(self.stress_old,self.__plastic_stress_old,
                                                        self.velocity*dt,
                                                        self.dthermal_stress,tau_Y,G,K,alpha,beta,h)

      def doStepPostprocessing(self,dt):
          super(DruckerPrager, self).doStepPostprocessing(dt)
          self.plastic_stress=self.__plastic_stress_old=self.plastic_stress

      def getNewStress(self,s,gamma_p,du,ds_therm,tau_Y,G,K,alpha,beta,h):
            k3=kronecker(self.domain)
            dt=1.
            g=grad(du)
            D=symmetric(g)
            W=nonsymmetric(g)
            s_e=s+ds_therm+dt*(2*G*D+(K-2./3*G)*trace(D)*k3 \
                               +2*nonsymmetric(matrix_mult(W,s)))
            p_e=-1./3*trace(s_e)
            s_e_dev=s_e+p_e*k3
            tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
            F=tau_e-alpha*p_e-tau_Y
            chi=whereNonNegative(F)
            l=chi*F/(h+G+beta*K)
            s=(1.-l*G/tau_e)*s_e_dev+(p_e+l*beta*K)*k3
            gamma_p=gamma_p+l
            return s, gamma_p
            

      def getTangentialTensor(self,s,tau_Y,G,K,alpha,beta,h):
           d=self.domain.getDim()
           k3=kronecker(Function(self.domain))
           p=-1./d*trace(s)
           s_dev=s+p*k3 
           tau=sqrt(1./2*inner(s_dev,s_dev))
           chi=whereNonNegative(tau-alpha*p-tau_Y)
           sXk3=outer(s,k3)
           k3Xk3=outer(k3,k3)
           tmp=G*s_dev/tau
           S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \
               + (K-2./3*G)*k3Xk3 \
               + sXk3-swap_axes(sXk3,1,3) \
               + 1./2*(swap_axes(sXk3,0,3)+swap_axes(sXk3,1,2) \
                      -swap_axes(sXk3,1,3)-swap_axes(sXk3,0,2))
               # - chi/(h+G+alpha*beta*K)*outer(tmp+beta*K*k3,tmp+alpha*K*k3)\
           return S
1	gross	814	# $Id:$
2
3			__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
4			http://www.access.edu.au
5			Primary Business: Queensland, Australia"""
6			__license__="""Licensed under the Open Software License version 3.0
7			http://www.opensource.org/licenses/osl-3.0.php"""
8
9			from esys.escript import *
10			from esys.escript.modelframe import Model,IterationDivergenceError
11			from esys.escript.linearPDEs import LinearPDE
12
13			class Mechanics(Model):
14			"""
15			base class for mechanics models in updated lagrangean framework
16
17			@ivar domain: domain (in)
18	gross	820	@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	gross	814	"""
33			def __init__(self,debug=False):
34			"""
35			set up the model
36
37			@param debug: debug flag
38			@type debug: C{bool}
39			"""
40			super(Mechanics, self).__init__(self,debug=debug)
41			self.declareParameter(domain=None, \
42			displacement=None, \
43			stress=None, \
44			velocity=None, \
45			internal_force=Data(), \
46			external_force=Data(), \
47			prescribed_velocity=Data(), \
48			location_prescribed_velocity=Data(), \
49			temperature = None, \
50			expansion_coefficient = 0., \
51			bulk_modulus=1., \
52			shear_modulus=1., \
53			rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
54			self.__iter=0
55
56			def doInitialization(self):
57			"""
58			initialize model
59			"""
60			if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
61			if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
62			if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
63			self.__pde=LinearPDE(self.domain)
64			self.__displacement_old=self.displacement
65			self.stress_old=self.stress
66			self.__velocity_old=self.velocity
67			self.__temperature_old=self.temperature
68
69			def doStepPreprocessing(self,dt):
70			"""
71			step up pressure iteration
72
73			if run within a time dependend problem extrapolation of pressure from previous time steps is used to
74			get an initial guess (that needs some work!!!!!!!)
75			"""
76			self.__iter=0
77			self.__diff=self.UNDEF_DT
78			# set new values:
79			self.displacement=self.__displacement_old
80			self.stress=self.stress_old
81			self.velocity=self.__velocity_old
82			self.temperature=self.__temperature_old
83			self.__velocity_last=self.velocity
84
85			def doStep(self,dt):
86			"""
87
88			performs an iteration step of the penalty method.
89			IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.
90
91			requires self.S to be set
92			updates the thermal stress increment
93
94			"""
95			k3=kronecker(self.domain)
96			# set new thermal stress increment
97			if self.temperature:
98			self.dthermal_stress=self.bulk_modulusself.self.expansion_coefficient(self.temperature-self.__temperature_old)
99			else:
100			self.dthermal_stress=0.
101			# set PDE coefficients:
102			self.__pde.setValue(A=self.S)
103			self.__pde.setValue(X=self.stress_old-self.dthermal_stress*k3)
104			if self.internal_force: self.__pde.setValue(Y=self.internal_force)
105			if self.external_force: self.__pde.setValue(y=self.external_force)
106	gross	819	print self.prescribed_velocity
107			print self.location_prescribed_velocity
108	gross	814	self.__pde.setValue(r=self.prescribed_velocity, \
109			q=self.location_prescribed_velocity)
110			# solve the PDE:
111			self.__pde.setTolerance(self.rel_tol/100.)
112			self.velocity=self.__pde.getSolution()
113			# calculate convergence indicators:
114			self.__diff,diff_old=Lsup(self.velocity-self.__velocity_last),self.__diff
115			self.__velocity_last=self.velocity
116			self.displacement=self.__displacement_old+dt*self.velocity
117			self.__iter+=1
118			self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
119			if self.__iter>2 and diff_old<self.__diff:
120			raise IterationDivergenceError,"no improvement in stress iteration"
121			if self.__iter>self.max_iter:
122			raise IterationDivergenceError,"Maximum number of iterations steps reached"
123
124			def terminateIteration(self):
125			"""iteration is terminateIterationd if relative pressure change is less then rel_tol"""
126			return self.__diff<=self.rel_tol*Lsup(self.velocity)+self.abs_tol
127
128			def doStepPostprocessing(self,dt):
129			"""
130			accept all the values:
131			"""
132			self.displacement=self.__displacement_old
133			self.stress=self.stress_old
134			self.velocity=self.__velocity_old
135			self.temperature=self.__temperature_old
136
137			def getSafeTimeStepSize(self,dt):
138			"""
139			returns new step size
140			"""
141			d=Lsup(self.velocity-self.__velocity_old)/dt
142			if d>0:
143			return Lsup(self.displacement)/d
144			else:
145			return self.UNDEF_DT
146
147
148
149			class DruckerPrager(Mechanics):
150			"""
151
152			"""
153
154			def __init__(self,debug=False):
155			"""
156			set up model
157			"""
158			super(DruckerPrager, self).__init__(debug=debug)
159			self.declareParameter(plastic_stress=0.,
160			friction_parameter=0.,
161			dilatancy_parameter=0.,
162			shear_length=0.)
163
164			def doInitialization(self):
165			"""
166			initialize model
167			"""
168			super(DruckerPrager, self).doInitialization()
169			self.__plastic_stress_old=self.plastic_stress
170	gross	816	self.__tau_y_old=self.shear_length
171	gross	814
172			def doStepPreprocessing(self,dt):
173			"""
174			step up pressure iteration
175
176			if run within a time dependend problem extrapolation of pressure from previous time steps is used to
177			get an initial guess (that needs some work!!!!!!!)
178			"""
179			super(DruckerPrager, self).doStepPreprocessing(dt)
180			self.plastic_stress=self.__plastic_stress_old
181
182			def doStep(self,dt):
183			G=self.shear_modulus
184			K=self.bulk_modulus
185			alpha=self.friction_parameter
186			beta=self.dilatancy_parameter
187			tau_Y=self.shear_length
188			if self.__plastic_stress_old:
189			dps=self.plastic_stress-self.__plastic_stress_old
190	gross	816	h=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))
191	gross	814	else:
192			h=0
193			# set new tangential operator:
194			self.S=self.getTangentialTensor(self.stress,
195			tau_Y,G,K,alpha,beta,h)
196			# do the update step:
197			super(DruckerPrager, self).doStep(dt)
198
199			# update stresses:
200			self.stress,self.plastic_stress=self.getNewStress(self.stress_old,self.__plastic_stress_old,
201			self.velocity*dt,
202			self.dthermal_stress,tau_Y,G,K,alpha,beta,h)
203
204			def doStepPostprocessing(self,dt):
205			super(DruckerPrager, self).doStepPostprocessing(dt)
206			self.plastic_stress=self.__plastic_stress_old=self.plastic_stress
207
208			def getNewStress(self,s,gamma_p,du,ds_therm,tau_Y,G,K,alpha,beta,h):
209			k3=kronecker(self.domain)
210			dt=1.
211			g=grad(du)
212			D=symmetric(g)
213			W=nonsymmetric(g)
214			s_e=s+ds_therm+dt(2GD+(K-2./3G)trace(D)k3 \
215			+2*nonsymmetric(matrix_mult(W,s)))
216			p_e=-1./3*trace(s_e)
217			s_e_dev=s_e+p_e*k3
218			tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
219			F=tau_e-alpha*p_e-tau_Y
220			chi=whereNonNegative(F)
221			l=chiF/(h+G+betaK)
222			s=(1.-lG/tau_e)s_e_dev+(p_e+lbetaK)*k3
223			gamma_p=gamma_p+l
224			return s, gamma_p
225
226
227			def getTangentialTensor(self,s,tau_Y,G,K,alpha,beta,h):
228			d=self.domain.getDim()
229			k3=kronecker(Function(self.domain))
230			p=-1./d*trace(s)
231			s_dev=s+p*k3
232			tau=sqrt(1./2*inner(s_dev,s_dev))
233			chi=whereNonNegative(tau-alpha*p-tau_Y)
234			sXk3=outer(s,k3)
235			k3Xk3=outer(k3,k3)
236			tmp=G*s_dev/tau
237			S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \
238			+ (K-2./3G)k3Xk3 \
239			+ sXk3-swap_axes(sXk3,1,3) \
240			+ 1./2*(swap_axes(sXk3,0,3)+swap_axes(sXk3,1,2) \
241			-swap_axes(sXk3,1,3)-swap_axes(sXk3,0,2))
242			# - chi/(h+G+alphabetaK)outer(tmp+betaKk3,tmp+alphaK*k3)\
243			return S