/[escript]/trunk/doc/examples/usersguide/slip.py
ViewVC logotype

Diff of /trunk/doc/examples/usersguide/slip.py

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

revision 5287 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 25  __url__="https://launchpad.net/escript-f Line 25  __url__="https://launchpad.net/escript-f
25  from esys.escript import *  from esys.escript import *
26  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
27  from esys.escript.models import FaultSystem  from esys.escript.models import FaultSystem
28  from esys.finley import Rectangle  try:
29        from esys.finley import Rectangle
30        HAVE_FINLEY = True
31    except ImportError:
32        HAVE_FINLEY = False
33  from esys.weipa import saveVTK  from esys.weipa import saveVTK
34  from esys.escript.unitsSI import DEG  from esys.escript.unitsSI import DEG
 #... set some parameters ...  
 lam=1.  
 mu=1  
 slip_max=1.  
35    
36  mydomain = Rectangle(l0=1.,l1=1.,n0=16, n1=16)  # n1 need to be multiple of 4!!!  if not HAVE_FINLEY:
37  # .. create the fault system      print("Finley module not available")
38  fs=FaultSystem(dim=2)  else:
39  fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)      #... set some parameters ...
40  # ... create a slip distribution on the fault:      lam=1.
41  p, m=fs.getParametrization(mydomain.getX(),tag=1)      mu=1
42  p0,p1= fs.getW0Range(tag=1)      slip_max=1.
43  s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]  
44  # ... calculate stress according to slip:      mydomain = Rectangle(l0=1.,l1=1.,n0=16, n1=16)  # n1 need to be multiple of 4!!!
45  D=symmetric(grad(s))      # .. create the fault system
46  chi, d=fs.getSideAndDistance(D.getFunctionSpace().getX(),tag=1)      fs=FaultSystem(dim=2)
47  sigma_s=(mu*D+lam*trace(D)*kronecker(mydomain))*chi      fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)
48  #... open symmetric PDE ...      # ... create a slip distribution on the fault:
49  mypde=LinearPDE(mydomain)      p, m=fs.getParametrization(mydomain.getX(),tag=1)
50  mypde.setSymmetryOn()      p0,p1= fs.getW0Range(tag=1)
51  #... set coefficients ...      s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]
52  C=Tensor4(0.,Function(mydomain))      # ... calculate stress according to slip:
53  for i in range(mydomain.getDim()):      D=symmetric(grad(s))
54    for j in range(mydomain.getDim()):      chi, d=fs.getSideAndDistance(D.getFunctionSpace().getX(),tag=1)
55       C[i,i,j,j]+=lam      sigma_s=(mu*D+lam*trace(D)*kronecker(mydomain))*chi
56       C[j,i,j,i]+=mu      #... open symmetric PDE ...
57       C[j,i,i,j]+=mu      mypde=LinearPDE(mydomain)
58  # ... fix displacement in normal direction      mypde.setSymmetryOn()
59  x=mydomain.getX()      #... set coefficients ...
60  msk=whereZero(x[0])*[1.,0.] + whereZero(x[0]-1.)*[1.,0.] \      C=Tensor4(0.,Function(mydomain))
61     +whereZero(x[1])*[0.,1.] + whereZero(x[1]-1.)*[0.,1.]      for i in range(mydomain.getDim()):
62  mypde.setValue(A=C,X=-0.5*sigma_s,q=msk)        for j in range(mydomain.getDim()):
63  #... solve pde ...           C[i,i,j,j]+=lam
64  mypde.getSolverOptions().setVerbosityOn()           C[j,i,j,i]+=mu
65  v=mypde.getSolution()           C[j,i,i,j]+=mu
66  # .. write the displacement to file:      # ... fix displacement in normal direction
67  D=symmetric(grad(v))      x=mydomain.getX()
68  sigma=(mu*D+lam*trace(D)*kronecker(mydomain))+0.5*sigma_s      msk=whereZero(x[0])*[1.,0.] + whereZero(x[0]-1.)*[1.,0.] \
69  saveVTK("slip.vtu",disp=v+0.5*chi*s, stress= sigma)         +whereZero(x[1])*[0.,1.] + whereZero(x[1]-1.)*[0.,1.]
70        mypde.setValue(A=C,X=-0.5*sigma_s,q=msk)
71        #... solve pde ...
72        mypde.getSolverOptions().setVerbosityOn()
73        v=mypde.getSolution()
74        # .. write the displacement to file:
75        D=symmetric(grad(v))
76        sigma=(mu*D+lam*trace(D)*kronecker(mydomain))+0.5*sigma_s
77        saveVTK("slip.vtu",disp=v+0.5*chi*s, stress= sigma)

Legend:
Removed from v.5287  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26