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

Diff of /trunk/doc/examples/cookbook/example06.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 35  import matplotlib Line 34  import matplotlib
34  matplotlib.use('agg') #It's just here for automated testing  matplotlib.use('agg') #It's just here for automated testing
35  from esys.pycad import *  from esys.pycad import *
36  from esys.pycad.gmsh import Design  from esys.pycad.gmsh import Design
 from esys.finley import MakeDomain  
37  from esys.escript import *  from esys.escript import *
38  import numpy as np  import numpy as np
39  import pylab as pl #Plotting package  import pylab as pl #Plotting package
# Line 44  from esys.escript.unitsSI import * Line 42  from esys.escript.unitsSI import *
42  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
43  import os, sys  import os, sys
44    
45    try:
46        # This imports the rectangle domain function
47        from esys.finley import MakeDomain#Converter for escript
48        HAVE_FINLEY = True
49    except ImportError:
50        print("Finley module not available")
51        HAVE_FINLEY = False
52  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
53  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
54          import sys          import sys
55          print("This example will not run in an MPI world.")          print("This example will not run in an MPI world.")
56          sys.exit(0)          sys.exit(0)
57  #################################################ESTABLISHING VARIABLES  
58  # where to put output files  if HAVE_FINLEY:
59  save_path= os.path.join("data","example06")      #################################################ESTABLISHING VARIABLES
60  mkDir(save_path)      # where to put output files
61        save_path= os.path.join("data","example06")
62  Ttop=20*Celsius # temperature at the top      mkDir(save_path)
63  qin=300.*Milli*W/(m*m) #our heat source temperature is now zero  
64        Ttop=20*Celsius # temperature at the top
65  ################################################ESTABLISHING PARAMETERS      qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
66  # Overall Domain  
67  p0=Point(0.0,        0.0, 0.0)      ################################################ESTABLISHING PARAMETERS
68  p1=Point(0.0,    -6000.0, 0.0)      # Overall Domain
69  p2=Point(5000.0, -6000.0, 0.0)      p0=Point(0.0,        0.0, 0.0)
70  p3=Point(5000.0,     0.0, 0.0)      p1=Point(0.0,    -6000.0, 0.0)
71        p2=Point(5000.0, -6000.0, 0.0)
72  ####################################################DOMAIN CONSTRUCTION      p3=Point(5000.0,     0.0, 0.0)
73  l01=Line(p0, p1)  
74  l12=Line(p1, p2)      ####################################################DOMAIN CONSTRUCTION
75  l23=Line(p2, p3)      l01=Line(p0, p1)
76  l30=Line(p3, p0)      l12=Line(p1, p2)
77        l23=Line(p2, p3)
78  # Generate Material Boundary      l30=Line(p3, p0)
79  p4=Point(0.0,    -2400.0, 0.0)  
80  p5=Point(2000.0, -2400.0, 0.0)      # Generate Material Boundary
81  p6=Point(3000.0, -6000.0, 0.0)      p4=Point(0.0,    -2400.0, 0.0)
82  p7=Point(5000.0, -2400.0, 0.0)      p5=Point(2000.0, -2400.0, 0.0)
83        p6=Point(3000.0, -6000.0, 0.0)
84  # Create TOP BLOCK      p7=Point(5000.0, -2400.0, 0.0)
85  tbl1=Line(p0,p4)  
86  tbl2=Line(p4,p5)      # Create TOP BLOCK
87  tbl3=Line(p5,p7)      tbl1=Line(p0,p4)
88  tbl4=Line(p7,p3)      tbl2=Line(p4,p5)
89  tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)      tbl3=Line(p5,p7)
90  tblock = PlaneSurface(tblockloop)      tbl4=Line(p7,p3)
91        tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)
92  # Create BOTTOM BLOCK LEFT      tblock = PlaneSurface(tblockloop)
93  bbll1=Line(p4,p1)  
94  bbll2=Line(p1,p6)      # Create BOTTOM BLOCK LEFT
95  bbll3=Line(p6,p5)      bbll1=Line(p4,p1)
96  bbll4=-tbl2      bbll2=Line(p1,p6)
97  bblockloopl = CurveLoop(bbll1,bbll2,bbll3,bbll4)      bbll3=Line(p6,p5)
98  bblockl = PlaneSurface(bblockloopl)      bbll4=-tbl2
99        bblockloopl = CurveLoop(bbll1,bbll2,bbll3,bbll4)
100  # Create BOTTOM BLOCK RIGHT      bblockl = PlaneSurface(bblockloopl)
101  bbrl1=Line(p6,p2)  
102  bbrl2=Line(p2,p7)      # Create BOTTOM BLOCK RIGHT
103  bbrl3=-tbl3      bbrl1=Line(p6,p2)
104  bbrl4=-bbll3      bbrl2=Line(p2,p7)
105  bblockloopr = CurveLoop(bbrl1,bbrl2,bbrl3,bbrl4)      bbrl3=-tbl3
106  bblockr = PlaneSurface(bblockloopr)      bbrl4=-bbll3
107        bblockloopr = CurveLoop(bbrl1,bbrl2,bbrl3,bbrl4)
108  #############################################EXPORTING MESH FOR ESCRIPT      bblockr = PlaneSurface(bblockloopr)
109  # Create a Design which can make the mesh  
110  d=Design(dim=2, element_size=200)      #############################################EXPORTING MESH FOR ESCRIPT
111  # Add the subdomains and flux boundaries.      # Create a Design which can make the mesh
112  d.addItems(PropertySet("top",tblock),\      d=Design(dim=2, element_size=200)
113                         PropertySet("bottomleft",bblockl),\      # Add the subdomains and flux boundaries.
114                         PropertySet("bottomright",bblockr),\      d.addItems(PropertySet("top",tblock),\
115                         PropertySet("linebottom",bbll2, bbrl1))                             PropertySet("bottomleft",bblockl),\
116                               PropertySet("bottomright",bblockr),\
117  # Create the geometry, mesh and Escript domain                             PropertySet("linebottom",bbll2, bbrl1))
118  d.setScriptFileName(os.path.join(save_path,"example06.geo"))  
119        # Create the geometry, mesh and Escript domain
120  d.setMeshFileName(os.path.join(save_path,"example06.msh"))      d.setScriptFileName(os.path.join(save_path,"example06.geo"))
121  domain=MakeDomain(d)  
122  print("Domain has been generated ...")      d.setMeshFileName(os.path.join(save_path,"example06.msh"))
123        domain=MakeDomain(d)
124  # set up kappa (thermal conductivity across domain) using tags      print("Domain has been generated ...")
125  kappa=Scalar(0,Function(domain))  
126  kappa.setTaggedValue("top",2.0)      # set up kappa (thermal conductivity across domain) using tags
127  kappa.setTaggedValue("bottomleft",10.0)      kappa=Scalar(0,Function(domain))
128  kappa.setTaggedValue("bottomright",6.0)      kappa.setTaggedValue("top",2.0)
129  ##############################################################SOLVE PDE      kappa.setTaggedValue("bottomleft",10.0)
130  mypde=LinearPDE(domain)      kappa.setTaggedValue("bottomright",6.0)
131  mypde.getSolverOptions().setVerbosityOn()      ##############################################################SOLVE PDE
132  mypde.setSymmetryOn()      mypde=LinearPDE(domain)
133  mypde.setValue(A=kappa*kronecker(domain))      mypde.getSolverOptions().setVerbosityOn()
134  x=Solution(domain).getX()      mypde.setSymmetryOn()
135  mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)      mypde.setValue(A=kappa*kronecker(domain))
136  qS=Scalar(0,FunctionOnBoundary(domain))      x=Solution(domain).getX()
137  qS.setTaggedValue("linebottom",qin)      mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
138  mypde.setValue(y=qS)      qS=Scalar(0,FunctionOnBoundary(domain))
139  print("PDE has been generated ...")      qS.setTaggedValue("linebottom",qin)
140  ###########################################################GET SOLUTION      mypde.setValue(y=qS)
141  T=mypde.getSolution()      print("PDE has been generated ...")
142  print("PDE has been solved ...")      ###########################################################GET SOLUTION
143  ###############################################################PLOTTING      T=mypde.getSolution()
144  # show temperature:      print("PDE has been solved ...")
145  xi, yi, zi = toRegGrid(T, nx=50, ny=50)      ###############################################################PLOTTING
146  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')      # show temperature:
147  pl.clabel(CS, inline=1, fontsize=8)      xi, yi, zi = toRegGrid(T, nx=50, ny=50)
148  # show sub domains:      CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
149  tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])      pl.clabel(CS, inline=1, fontsize=8)
150  pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)      # show sub domains:
151  bpgr=np.array([p.getCoordinates() for p in bblockloopr.getPolygon() ])      tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])
152  pl.fill(bpgr[:,0],bpgr[:,1],'orange',label='6 W/m/k',zorder=-1000)      pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)
153  bpgl=np.array([p.getCoordinates() for p in bblockloopl.getPolygon() ])      bpgr=np.array([p.getCoordinates() for p in bblockloopr.getPolygon() ])
154  pl.fill(bpgl[:,0],bpgl[:,1],'red',label='10 W/m/k',zorder=-1000)      pl.fill(bpgr[:,0],bpgr[:,1],'orange',label='6 W/m/k',zorder=-1000)
155  # show flux:      bpgl=np.array([p.getCoordinates() for p in bblockloopl.getPolygon() ])
156  xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)      pl.fill(bpgl[:,0],bpgl[:,1],'red',label='10 W/m/k',zorder=-1000)
157  pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")      # show flux:
158  # create plot      xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)
159  pl.title("Heat Refraction across an anisotropic structure\n with flux vectors.")      pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")
160  pl.xlabel("Horizontal Displacement (m)")      # create plot
161  pl.ylabel("Depth (m)")      pl.title("Heat Refraction across an anisotropic structure\n with flux vectors.")
162  pl.legend()      pl.xlabel("Horizontal Displacement (m)")
163  pl.savefig(os.path.join(save_path,"flux.png"))      pl.ylabel("Depth (m)")
164  print("Flux has been plotted  ...")      pl.legend()
165        pl.savefig(os.path.join(save_path,"flux.png"))
166        print("Flux has been plotted  ...")

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

  ViewVC Help
Powered by ViewVC 1.1.26