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 |
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()): |
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 |
|
|
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 |
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) |