/[escript]/trunk/doc/examples/cookbook/example11b.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/example11b.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 43  from esys.escript import * # This import Line 43  from esys.escript import * # This import
43  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
44  from esys.escript.pdetools import Projector  from esys.escript.pdetools import Projector
45  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
 from esys.finley import Rectangle # This imports the rectangle domain function from finley  
46  from esys.weipa import saveVTK  from esys.weipa import saveVTK
47    
48  from cblib import toRegGrid  from cblib import toRegGrid
49    
50    try:
51        from esys.finley import Rectangle
52        HAVE_FINLEY = True
53    except ImportError:
54        print("Finley module not available")
55        HAVE_FINLEY = False
56  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
57  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
58      import sys      import sys
59      print("This example will not run in an MPI world.")      print("This example will not run in an MPI world.")
60      sys.exit(0)      sys.exit(0)
61    
62  #################################################ESTABLISHING VARIABLES  if HAVE_FINLEY:
63  #Domain related.      #################################################ESTABLISHING VARIABLES
64  mx = 1000*m #meters - model length      #Domain related.
65  my = -250*m #meters - model depth      mx = 1000*m #meters - model length
66  ndx = 200 # mesh steps in x direction      my = -250*m #meters - model depth
67  ndy = 200 # mesh steps in y direction - one dimension means one element      ndx = 200 # mesh steps in x direction
68  #PDE related      ndy = 200 # mesh steps in y direction - one dimension means one element
69  res1=500.0      #PDE related
70  res2=10.0      res1=500.0
71  res3=10000.0      res2=10.0
72  #con=1/res      res3=10000.0
73  cur=1000.      #con=1/res
74        cur=1000.
75  ################################################ESTABLISHING PARAMETERS  
76  #the folder to put our outputs in, leave blank "" for script path      ################################################ESTABLISHING PARAMETERS
77  save_path= os.path.join("data","example11")      #the folder to put our outputs in, leave blank "" for script path
78  #ensure the dir exists      save_path= os.path.join("data","example11")
79  mkDir(save_path)      #ensure the dir exists
80        mkDir(save_path)
81  ####################################################DOMAIN CONSTRUCTION  
82  domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)      ####################################################DOMAIN CONSTRUCTION
83  x=Solution(domain).getX()      domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
84        x=Solution(domain).getX()
85  kro=kronecker(domain)  
86  source1=[3.*mx/8.,0]; source2=[5.*mx/8.,0]      kro=kronecker(domain)
87        source1=[3.*mx/8.,0]; source2=[5.*mx/8.,0]
88  c1=length(exp(-length(x-source1)/(10.))); c1=c1/integrate(c1)  
89  c2=-length(exp(-length(x-source2)/(10.))); c2=c2/integrate(c2)      c1=length(exp(-length(x-source1)/(10.))); c1=c1/integrate(c1)
90  sourceg=cur*(c1-c2)      c2=-length(exp(-length(x-source2)/(10.))); c2=c2/integrate(c2)
91        sourceg=cur*(c1-c2)
92  res=res1*wherePositive(x[1]-my/3)+res2*whereNegative(x[1]-my/3)*wherePositive(x[1]-my*2/3)+res3*whereNegative(x[1]-my*2/3)  
93  con=1/res      res=res1*wherePositive(x[1]-my/3)+res2*whereNegative(x[1]-my/3)*wherePositive(x[1]-my*2/3)+res3*whereNegative(x[1]-my*2/3)
94  q=whereZero(x[1]-my)+whereZero(x[0])+whereZero(x[0]-mx)      con=1/res
95  ###############################################ESCRIPT PDE CONSTRUCTION      q=whereZero(x[1]-my)+whereZero(x[0])+whereZero(x[0]-mx)
96        ###############################################ESCRIPT PDE CONSTRUCTION
97  mypde=LinearPDE(domain)  
98  mypde.setValue(A=con*kro,Y=sourceg,q=q,r=0)      mypde=LinearPDE(domain)
99  #mypde.setSymmetryOn()      mypde.setValue(A=con*kro,Y=sourceg,q=q,r=0)
100  sol=mypde.getSolution()      #mypde.setSymmetryOn()
101        sol=mypde.getSolution()
102  # Save the output to file.  
103  saveVTK(os.path.join(save_path,"ex11b.vtu"),\      # Save the output to file.
104          source=sourceg,\      saveVTK(os.path.join(save_path,"ex11b.vtu"),\
105          res_pot=sol,\              source=sourceg,\
106          res=res,\              res_pot=sol,\
107          curden=-con*grad(sol),\              res=res,\
108          efield=-grad(sol))              curden=-con*grad(sol),\
109                efield=-grad(sol))

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

  ViewVC Help
Powered by ViewVC 1.1.26