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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 820 - (hide annotations)
Mon Aug 28 06:55:36 2006 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 9566 byte(s)
make the modelframe test run
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_modulus*self.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*(2*G*D+(K-2./3*G)*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=chi*F/(h+G+beta*K)
222     s=(1.-l*G/tau_e)*s_e_dev+(p_e+l*beta*K)*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./3*G)*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+alpha*beta*K)*outer(tmp+beta*K*k3,tmp+alpha*K*k3)\
243     return S

  ViewVC Help
Powered by ViewVC 1.1.26