/[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 886 by gross, Mon Oct 30 06:37:30 2006 UTC revision 887 by gross, Thu Nov 2 07:17:07 2006 UTC
# Line 30  from esys.finley import ReadMesh Line 30  from esys.finley import ReadMesh
30    
31    
32  rho=0.  rho=0.
33  lam_lmbd=1.  lam_lmbd=1.7e11
34  lam_mu=1.  lam_mu=1.7e11
35  g=9.81  g=9.81
36    fstart =  [50000.0, 44444.444444444445, 11666.666666666666]
37    fend =  [50000.0, 50000.0, -13333.333333333336]
38    fstart =  [50000.0, 44444.444444444445, 11666.666666666666]
39    fend =  [50000.0, 55555.555555555555, -11666.666666666668]
40    
41    
42    
43  class SlippingFault(SaddlePointProblem):  class SlippingFault(SaddlePointProblem):
44        """        """
# Line 56  class SlippingFault(SaddlePointProblem): Line 62  class SlippingFault(SaddlePointProblem):
62           self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction)           self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction)
63    
64        def inner(self,p0,p1):        def inner(self,p0,p1):
65           return integrate(p0*p1,FunctionOnContactZero(self.domain))           return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))
66    
67        def solve_f(self,u,p,tol=1.e-8):        def solve_f(self,u,p,tol=1.e-8):
68           self.__pde_u.setTolerance(tol)           self.__pde_u.setTolerance(tol)
69           self.__pde_u.setValue(y_contact=-p)           self.__pde_u.setValue(y_contact=p)
70           return  self.__pde_u.getSolution()           return  self.__pde_u.getSolution()
71    
72        def solve_g(self,u,tol=1.e-8):        def solve_g(self,u,tol=1.e-8):
73           dp=-(self.slip-jump(u))           dp=(self.slip-jump(u))*lam_lmbd/FunctionOnContactZero(self.domain).getX()
74           return  dp           return  dp
75    
76    
 s=numarray.array([0.,1.,1.])  
77  dom=ReadMesh("meshfault3D.fly")  dom=ReadMesh("meshfault3D.fly")
78  prop=SlippingFault(dom)  prop=SlippingFault(dom)
79  d=dom.getDim()  d=dom.getDim()
80  x=dom.getX()[d-1]  x=dom.getX()[d-1]
81  mask=whereZero(x-inf(x))*numarray.ones((d,))  mask=whereZero(x-inf(x))*numarray.ones((d,))
82    s=numarray.array([0.,1.,1.])
83    x=FunctionOnContactZero(dom).getX()
84    s=numarray.array([0.,1.,1.])
85    for i in range(3):
86         d=fend[i]-fstart[i]
87         if d>0:
88             q=(x[i]-fstart[i])/d
89             s=q*(1-q)*4*s
90         elif d<0:
91             q=(x[i]-fend[i])/d
92             s=q*(1-q)*4*s
93  u0=Vector(0.,Solution(dom))  u0=Vector(0.,Solution(dom))
94  p0=Vector(1.,FunctionOnContactZero(dom))  p0=Vector(1.,FunctionOnContactZero(dom))
95  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)
96  u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.01)  u,p=prop.solve(u0,p0,iter_max=100,tolerance=0.5,accepted_reduction=1.1 )
97  saveVTK("dis.xml",u=u)  saveVTK("dis.xml",u=u)
98  saveVTK("fault.xml",sigma=p,s=jump(u))  saveVTK("fault.xml",sigma=p,s=jump(u))
Legend:
Removed from v.886  
changed lines
  Added in v.887
  ViewVC Help
Powered by ViewVC 1.1.26