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

Diff of /trunk/doc/examples/cookbook/example05c.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 37  import matplotlib Line 36  import matplotlib
36  matplotlib.use('agg') #It's just here for automated testing  matplotlib.use('agg') #It's just here for automated testing
37  from esys.pycad import * #domain constructor  from esys.pycad import * #domain constructor
38  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  
39  import os #file path tool  import os #file path tool
40  from math import * # math package  from math import * # math package
41  from esys.escript import *  from esys.escript import *
# Line 48  from cblib import toRegGrid, subsample Line 46  from cblib import toRegGrid, subsample
46  import pylab as pl #Plotting package  import pylab as pl #Plotting package
47  import numpy as np  import numpy as np
48    
49    try:
50        # This imports the rectangle domain function
51        from esys.finley import MakeDomain#Converter for escript
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  #set modal to 1 for a syncline or -1 for an anticline structural      #################################################ESTABLISHING VARIABLES
64  #configuration      #set modal to 1 for a syncline or -1 for an anticline structural
65  modal=-1      #configuration
66        modal=-1
67  # the folder to put our outputs in, leave blank "" for script path -  
68  # note this folder path must exist to work      # the folder to put our outputs in, leave blank "" for script path -
69  save_path= os.path.join("data","example05")      # note this folder path must exist to work
70  mkDir(save_path)      save_path= os.path.join("data","example05")
71        mkDir(save_path)
72  ################################################ESTABLISHING PARAMETERS  
73  #Model Parameters      ################################################ESTABLISHING PARAMETERS
74  width=5000.0*m   #width of model      #Model Parameters
75  depth=-6000.0*m  #depth of model      width=5000.0*m   #width of model
76  Ttop=20*K       # top temperature      depth=-6000.0*m  #depth of model
77  qin=70*Milli*W/(m*m) # bottom heat influx      Ttop=20*K       # top temperature
78        qin=70*Milli*W/(m*m) # bottom heat influx
79  sspl=51 #number of discrete points in spline  
80  dsp=width/(sspl-1) #dx of spline steps for width      sspl=51 #number of discrete points in spline
81  dep_sp=2500.0*m #avg depth of spline      dsp=width/(sspl-1) #dx of spline steps for width
82  h_sp=1500.0*m #heigh of spline      dep_sp=2500.0*m #avg depth of spline
83  orit=-1.0 #orientation of spline 1.0=>up -1.0=>down      h_sp=1500.0*m #heigh of spline
84        orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
85  ####################################################DOMAIN CONSTRUCTION  
86  # Domain Corners      ####################################################DOMAIN CONSTRUCTION
87  p0=Point(0.0,      0.0, 0.0)      # Domain Corners
88  p1=Point(0.0,    depth, 0.0)      p0=Point(0.0,      0.0, 0.0)
89  p2=Point(width, depth, 0.0)      p1=Point(0.0,    depth, 0.0)
90  p3=Point(width,   0.0, 0.0)      p2=Point(width, depth, 0.0)
91        p3=Point(width,   0.0, 0.0)
92  # Generate Material Boundary  
93  x=[ Point(i*dsp\      # Generate Material Boundary
94      ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\      x=[ Point(i*dsp\
95       for i in range(0,sspl)\          ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
96    ]           for i in range(0,sspl)\
97  mysp = Spline(*tuple(x))        ]
98  # Start and end of material boundary.      mysp = Spline(*tuple(x))
99  x1=mysp.getStartPoint()      # Start and end of material boundary.
100  x2=mysp.getEndPoint()      x1=mysp.getStartPoint()
101        x2=mysp.getEndPoint()
102  #  Create TOP BLOCK  
103  # lines      #  Create TOP BLOCK
104  tbl1=Line(p0,x1)      # lines
105  tbl2=mysp      tbl1=Line(p0,x1)
106  tbl3=Line(x2,p3)      tbl2=mysp
107  l30=Line(p3, p0)      tbl3=Line(x2,p3)
108  # curve      l30=Line(p3, p0)
109  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)      # curve
110  # surface      tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
111  tblock = PlaneSurface(tblockloop)      # surface
112  # Create BOTTOM BLOCK      tblock = PlaneSurface(tblockloop)
113  # lines      # Create BOTTOM BLOCK
114  bbl1=Line(x1,p1)      # lines
115  bbl3=Line(p2,x2)      bbl1=Line(x1,p1)
116  bbl4=-mysp      bbl3=Line(p2,x2)
117  l12=Line(p1, p2)      bbl4=-mysp
118  # curve      l12=Line(p1, p2)
119  bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)      # curve
120        bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)
121  # surface  
122  bblock = PlaneSurface(bblockloop)      # surface
123        bblock = PlaneSurface(bblockloop)
124  ################################################CREATE MESH FOR ESCRIPT  
125  # Create a Design which can make the mesh      ################################################CREATE MESH FOR ESCRIPT
126  d=Design(dim=2, element_size=200)      # Create a Design which can make the mesh
127  # Add the subdomains and flux boundaries.      d=Design(dim=2, element_size=200)
128  d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\      # Add the subdomains and flux boundaries.
129                                       PropertySet("linebottom",l12))      d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
130  # Create the geometry, mesh and Escript domain                                           PropertySet("linebottom",l12))
131  d.setScriptFileName(os.path.join(save_path,"example05.geo"))      # Create the geometry, mesh and Escript domain
132  d.setMeshFileName(os.path.join(save_path,"example05.msh"))      d.setScriptFileName(os.path.join(save_path,"example05.geo"))
133  domain=MakeDomain(d, optimizeLabeling=True)      d.setMeshFileName(os.path.join(save_path,"example05.msh"))
134  print("Domain has been generated ...")      domain=MakeDomain(d, optimizeLabeling=True)
135  ##############################################################SOLVE PDE      print("Domain has been generated ...")
136  mypde=LinearPDE(domain)      ##############################################################SOLVE PDE
137  mypde.getSolverOptions().setVerbosityOn()      mypde=LinearPDE(domain)
138  mypde.setSymmetryOn()      mypde.getSolverOptions().setVerbosityOn()
139  kappa=Scalar(0,Function(domain))      mypde.setSymmetryOn()
140  kappa.setTaggedValue("top",2.0*W/m/K)      kappa=Scalar(0,Function(domain))
141  kappa.setTaggedValue("bottom",4.0*W/m/K)      kappa.setTaggedValue("top",2.0*W/m/K)
142  mypde.setValue(A=kappa*kronecker(domain))      kappa.setTaggedValue("bottom",4.0*W/m/K)
143  x=Solution(domain).getX()      mypde.setValue(A=kappa*kronecker(domain))
144  mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)      x=Solution(domain).getX()
145  qS=Scalar(0,FunctionOnBoundary(domain))      mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
146  qS.setTaggedValue("linebottom",qin)      qS=Scalar(0,FunctionOnBoundary(domain))
147  mypde.setValue(y=qS)      qS.setTaggedValue("linebottom",qin)
148  print("PDE has been generated ...")      mypde.setValue(y=qS)
149  ###########################################################GET SOLUTION      print("PDE has been generated ...")
150  T=mypde.getSolution()      ###########################################################GET SOLUTION
151  print("PDE has been solved  ...")      T=mypde.getSolution()
152  ###############################################################PLOTTING      print("PDE has been solved  ...")
153  # show temperature:      ###############################################################PLOTTING
154  xi, yi, zi = toRegGrid(T, nx=50, ny=50)      # show temperature:
155  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')      xi, yi, zi = toRegGrid(T, nx=50, ny=50)
156  pl.clabel(CS, inline=1, fontsize=8)      CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
157  # show sub domains:      pl.clabel(CS, inline=1, fontsize=8)
158  tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])      # show sub domains:
159  pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)      tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])
160  bpg=np.array([p.getCoordinates() for p in bblockloop.getPolygon() ])      pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)
161  pl.fill(bpg[:,0],bpg[:,1],'red',label='4 W/m/k',zorder=-1000)      bpg=np.array([p.getCoordinates() for p in bblockloop.getPolygon() ])
162  # show flux:      pl.fill(bpg[:,0],bpg[:,1],'red',label='4 W/m/k',zorder=-1000)
163  xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)      # show flux:
164  pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")      xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)
165  # create plot      pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")
166  pl.title("Heat Refraction across a clinal structure\n with heat flux.")      # create plot
167  pl.xlabel("Horizontal Displacement (m)")      pl.title("Heat Refraction across a clinal structure\n with heat flux.")
168  pl.ylabel("Depth (m)")      pl.xlabel("Horizontal Displacement (m)")
169  pl.legend()      pl.ylabel("Depth (m)")
170  pl.savefig(os.path.join(save_path,"flux.png"))      pl.legend()
171  print("Flux has been plotted  ...")      pl.savefig(os.path.join(save_path,"flux.png"))
172        print("Flux has been plotted  ...")

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

  ViewVC Help
Powered by ViewVC 1.1.26