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

Diff of /trunk/doc/examples/cookbook/example05a.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 1  Line 1 
1  from __future__ import division  from __future__ import division, print_function
 from __future__ import print_function  
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2009-2014 by University of Queensland  # Copyright (c) 2009-2014 by University of Queensland
# Line 36  import matplotlib Line 35  import matplotlib
35  matplotlib.use('agg') #It's just here for automated testing  matplotlib.use('agg') #It's just here for automated testing
36  from esys.pycad import * #domain constructor  from esys.pycad import * #domain constructor
37  from esys.pycad.gmsh import Design #Finite Element meshing package  from esys.pycad.gmsh import Design #Finite Element meshing package
 from esys.finley import MakeDomain #Converter for escript  
38  import os #file path tool  import os #file path tool
39  from math import * # math package  from math import * # math package
40  from esys.escript import *  from esys.escript import *
# Line 45  from esys.escript.linearPDEs import Line Line 43  from esys.escript.linearPDEs import Line
43  from cblib import toRegGrid  from cblib import toRegGrid
44  import pylab as pl #Plotting package  import pylab as pl #Plotting package
45    
46    try:
47        # This imports the rectangle domain function
48        from esys.finley import MakeDomain#Converter for escript
49        HAVE_FINLEY = True
50    except ImportError:
51        print("Finley module not available")
52        HAVE_FINLEY = False
53  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
54  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
55          import sys          import sys
56          print("This example will not run in an MPI world.")          print("This example will not run in an MPI world.")
57          sys.exit(0)          sys.exit(0)
58    
59  #################################################ESTABLISHING VARIABLES  if HAVE_FINLEY:
60  #set modal to 1 for a syncline or -1 for an anticline structural      #################################################ESTABLISHING VARIABLES
61  #configuration      #set modal to 1 for a syncline or -1 for an anticline structural
62  modal=-1      #configuration
63        modal=-1
64  # the folder to put our outputs in, leave blank "" for script path -  
65  # note this folder path must exist to work      # the folder to put our outputs in, leave blank "" for script path -
66  save_path= os.path.join("data","example05")      # note this folder path must exist to work
67  mkDir(save_path)      save_path= os.path.join("data","example05")
68        mkDir(save_path)
69  ################################################ESTABLISHING PARAMETERS  
70  #Model Parameters      ################################################ESTABLISHING PARAMETERS
71  width=5000.0*m   #width of model      #Model Parameters
72  depth=-6000.0*m  #depth of model      width=5000.0*m   #width of model
73  Ttop=20*K       # top temperature      depth=-6000.0*m  #depth of model
74  qin=70*Milli*W/(m*m) # bottom heat influx      Ttop=20*K       # top temperature
75        qin=70*Milli*W/(m*m) # bottom heat influx
76  sspl=51 #number of discrete points in spline  
77  dsp=width/(sspl-1) #dx of spline steps for width      sspl=51 #number of discrete points in spline
78  dep_sp=2500.0*m #avg depth of spline      dsp=width/(sspl-1) #dx of spline steps for width
79  h_sp=1500.0*m #heigh of spline      dep_sp=2500.0*m #avg depth of spline
80  orit=-1.0 #orientation of spline 1.0=>up -1.0=>down      h_sp=1500.0*m #heigh of spline
81        orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
82  ####################################################DOMAIN CONSTRUCTION  
83  # Domain Corners      ####################################################DOMAIN CONSTRUCTION
84  p0=Point(0.0,      0.0, 0.0)      # Domain Corners
85  p1=Point(0.0,    depth, 0.0)      p0=Point(0.0,      0.0, 0.0)
86  p2=Point(width, depth, 0.0)      p1=Point(0.0,    depth, 0.0)
87  p3=Point(width,   0.0, 0.0)      p2=Point(width, depth, 0.0)
88        p3=Point(width,   0.0, 0.0)
89  # Generate Material Boundary  
90  x=[ Point(i*dsp\      # Generate Material Boundary
91      ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\      x=[ Point(i*dsp\
92       for i in range(0,sspl)\          ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
93    ]           for i in range(0,sspl)\
94  mysp = Spline(*tuple(x))        ]
95  # Start and end of material boundary.      mysp = Spline(*tuple(x))
96  x1=mysp.getStartPoint()      # Start and end of material boundary.
97  x2=mysp.getEndPoint()      x1=mysp.getStartPoint()
98        x2=mysp.getEndPoint()
99  #  Create TOP BLOCK  
100  # lines      #  Create TOP BLOCK
101  tbl1=Line(p0,x1)      # lines
102  tbl2=mysp      tbl1=Line(p0,x1)
103  tbl3=Line(x2,p3)      tbl2=mysp
104  l30=Line(p3, p0)      tbl3=Line(x2,p3)
105  # curve      l30=Line(p3, p0)
106  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)      # curve
107  # surface      tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
108  tblock = PlaneSurface(tblockloop)      # surface
109        tblock = PlaneSurface(tblockloop)
110    
111  # Create BOTTOM BLOCK  
112  # lines      # Create BOTTOM BLOCK
113  bbl1=Line(x1,p1)      # lines
114  bbl3=Line(p2,x2)      bbl1=Line(x1,p1)
115  bbl4=-mysp      bbl3=Line(p2,x2)
116  l12=Line(p1, p2)      bbl4=-mysp
117  # curve      l12=Line(p1, p2)
118  bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)      # curve
119  # surface      bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)
120  bblock = PlaneSurface(bblockloop)      # surface
121        bblock = PlaneSurface(bblockloop)
122  ################################################CREATE MESH FOR ESCRIPT  
123  # Create a Design which can make the mesh      ################################################CREATE MESH FOR ESCRIPT
124  d=Design(dim=2, element_size=200)      # Create a Design which can make the mesh
125  # Add the subdomains and flux boundaries.      d=Design(dim=2, element_size=200)
126  d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\      # Add the subdomains and flux boundaries.
127                                       PropertySet("linebottom",l12))      d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
128  # Create the geometry, mesh and Escript domain                                           PropertySet("linebottom",l12))
129  d.setScriptFileName(os.path.join(save_path,"example05.geo"))      # Create the geometry, mesh and Escript domain
130  d.setMeshFileName(os.path.join(save_path,"example05.msh"))      d.setScriptFileName(os.path.join(save_path,"example05.geo"))
131  domain=MakeDomain(d, optimizeLabeling=True)      d.setMeshFileName(os.path.join(save_path,"example05.msh"))
132  print("Domain has been generated ...")      domain=MakeDomain(d, optimizeLabeling=True)
133  ##############################################################SOLVE PDE      print("Domain has been generated ...")
134  mypde=LinearPDE(domain)      ##############################################################SOLVE PDE
135  mypde.getSolverOptions().setVerbosityOn()      mypde=LinearPDE(domain)
136  mypde.setSymmetryOn()      mypde.getSolverOptions().setVerbosityOn()
137  kappa=Scalar(0,Function(domain))      mypde.setSymmetryOn()
138  kappa.setTaggedValue("top",2.0*W/m/K)      kappa=Scalar(0,Function(domain))
139  kappa.setTaggedValue("bottom",4.0*W/m/K)      kappa.setTaggedValue("top",2.0*W/m/K)
140  mypde.setValue(A=kappa*kronecker(domain))      kappa.setTaggedValue("bottom",4.0*W/m/K)
141  x=Solution(domain).getX()      mypde.setValue(A=kappa*kronecker(domain))
142  mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)      x=Solution(domain).getX()
143  qS=Scalar(0,FunctionOnBoundary(domain))      mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
144  qS.setTaggedValue("linebottom",qin)      qS=Scalar(0,FunctionOnBoundary(domain))
145  mypde.setValue(y=qS)      qS.setTaggedValue("linebottom",qin)
146  print("PDE has been generated ...")      mypde.setValue(y=qS)
147  ###########################################################GET SOLUTION      print("PDE has been generated ...")
148  T=mypde.getSolution()      ###########################################################GET SOLUTION
149  print("PDE has been solved  ...")      T=mypde.getSolution()
150        print("PDE has been solved  ...")
151  #######################################################################  
152  xi, yi, zi = toRegGrid(T, nx=50, ny=50)      #######################################################################
153  pl.matplotlib.pyplot.autumn()      xi, yi, zi = toRegGrid(T, nx=50, ny=50)
154  pl.contourf(xi,yi,zi,10)      pl.matplotlib.pyplot.autumn()
155  pl.xlabel("Horizontal Displacement (m)")      pl.contourf(xi,yi,zi,10)
156  pl.ylabel("Depth (m)")      pl.xlabel("Horizontal Displacement (m)")
157  pl.savefig(os.path.join(save_path,"Tcontour.png"))      pl.ylabel("Depth (m)")
158  print("Solution has been plotted  ...")      pl.savefig(os.path.join(save_path,"Tcontour.png"))
159        print("Solution has been plotted  ...")

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

  ViewVC Help
Powered by ViewVC 1.1.26