/[escript]/trunk/finley/test/python/slip_stress.py
ViewVC logotype

Diff of /trunk/finley/test/python/slip_stress.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 892 by gross, Wed Nov 8 06:14:29 2006 UTC revision 893 by gross, Wed Nov 8 08:20:19 2006 UTC
# Line 33  rho=0. Line 33  rho=0.
33  lam_lmbd=1.7e11  lam_lmbd=1.7e11
34  lam_mu=1.7e11  lam_mu=1.7e11
35  g=9.81  g=9.81
36  fstart =  [50000.0, 44444.444444444445, 11666.666666666666]  fstart =  [50000.0, 40000.0, 10909.09090909091]
37  fend =  [50000.0, 50000.0, -13333.333333333336]  fend =  [50000.0, 60000.0, 19090.909090909092]
38  fstart =  [50000.0, 44444.444444444445, 11666.666666666666]  
 fend =  [50000.0, 55555.555555555555, -11666.666666666668]  
39    
40    
41    
# Line 66  class SlippingFault(SaddlePointProblem): Line 65  class SlippingFault(SaddlePointProblem):
65    
66        def solve_f(self,u,p,tol=1.e-8):        def solve_f(self,u,p,tol=1.e-8):
67           self.__pde_u.setTolerance(tol)           self.__pde_u.setTolerance(tol)
68           self.__pde_u.setValue(y_contact=p)           self.__pde_u.setValue(y_contact=-p)
69           print "p:",inf(p),sup(p)           # print "p:",inf(p),sup(p)
70           print "u:",inf(u),sup(u)           # print "u:",inf(u),sup(u)
71           self.__pde_u.setValue(y_contact=p)           self.__pde_u.setValue(y_contact=-p)
72           return  self.__pde_u.getSolution()           return  self.__pde_u.getSolution()
73    
74        def solve_g(self,u,tol=1.e-8):        def solve_g(self,u,tol=1.e-8):
75           dp=-(self.slip-jump(u))*lam_lmbd/FunctionOnContactZero(self.domain).getSize()           dp=Vector(0.,FunctionOnContactZero(self.domain))
76           print dp           h=FunctionOnContactZero(self.domain).getSize()
77             # print jump(u)-self.slip
78             dp[0]=(self.slip[0]-jump(u[0]))*lam_mu/h
79             dp[1]=(self.slip[1]-jump(u[1]))*lam_mu/h
80             dp[2]=(self.slip[2]-jump(u[2]))*lam_mu/h
81           return  dp           return  dp
82    
83    
# Line 84  d=dom.getDim() Line 87  d=dom.getDim()
87  x=dom.getX()[0]  x=dom.getX()[0]
88  # x=dom.getX()[d-1]  # x=dom.getX()[d-1]
89  mask=whereZero(x-inf(x))*numarray.ones((d,))  mask=whereZero(x-inf(x))*numarray.ones((d,))
 s=numarray.array([0.,1.,1.])  
90  x=FunctionOnContactZero(dom).getX()  x=FunctionOnContactZero(dom).getX()
91  s=numarray.array([0.,1.,1.])  s=numarray.array([-100000.,1.,1.])
92  for i in range(3):  for i in range(3):
93       d=fend[i]-fstart[i]       d=fend[i]-fstart[i]
94       if d>0:       if d>0:
# Line 98  for i in range(3): Line 100  for i in range(3):
100  u0=Vector(0.,Solution(dom))  u0=Vector(0.,Solution(dom))
101  p0=Vector(1.,FunctionOnContactZero(dom))  p0=Vector(1.,FunctionOnContactZero(dom))
102  prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu)  prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu)
103  u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.1,accepted_reduction=0.99)  u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.13,accepted_reduction=1.)
104  saveVTK("dis.xml",u=u)  saveVTK("dis.xml",u=u)
105  saveVTK("fault.xml",sigma=p,s=jump(u))  saveVTK("fault.xml",sigma=p,s=jump(u))
Legend:
Removed from v.892  
changed lines
  Added in v.893
  ViewVC Help
Powered by ViewVC 1.1.26