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

Contents of /trunk/doc/examples/cookbook/example06.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2977 - (show annotations)
Tue Mar 9 00:33:28 2010 UTC (9 years, 9 months ago) by ahallam
File MIME type: text/x-python
File size: 5369 byte(s)
cookbook review final 3.1
1
2 ########################################################
3 #
4 # Copyright (c) 2009-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2009-2010 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 """
23 Author: Antony Hallam antony.hallam@uqconnect.edu.au
24 """
25 ############################################################FILE HEADER
26 # example06.py
27 # Create a 3 block fault and overburden style model in 2d using pycad
28 # meshing tools.
29
30 #######################################################EXTERNAL MODULES
31 import matplotlib
32 matplotlib.use('agg') #It's just here for automated testing
33 from esys.pycad import *
34 from esys.pycad.gmsh import Design
35 from esys.finley import MakeDomain
36 from esys.escript import *
37 import numpy as np
38 import pylab as pl #Plotting package
39 from cblib import toRegGrid, subsample
40 from esys.escript.unitsSI import *
41 from esys.escript.linearPDEs import LinearPDE
42 import os, sys
43
44 ########################################################MPI WORLD CHECK
45 if getMPISizeWorld() > 1:
46 import sys
47 print "This example will not run in an MPI world."
48 sys.exit(0)
49 #################################################ESTABLISHING VARIABLES
50 # where to put output files
51 save_path= os.path.join("data","example06")
52 mkDir(save_path)
53
54 Ttop=20*Celsius # temperature at the top
55 qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
56
57 ################################################ESTABLISHING PARAMETERS
58 # Overall Domain
59 p0=Point(0.0, 0.0, 0.0)
60 p1=Point(0.0, -6000.0, 0.0)
61 p2=Point(5000.0, -6000.0, 0.0)
62 p3=Point(5000.0, 0.0, 0.0)
63
64 ####################################################DOMAIN CONSTRUCTION
65 l01=Line(p0, p1)
66 l12=Line(p1, p2)
67 l23=Line(p2, p3)
68 l30=Line(p3, p0)
69
70 # Generate Material Boundary
71 p4=Point(0.0, -2400.0, 0.0)
72 p5=Point(2000.0, -2400.0, 0.0)
73 p6=Point(3000.0, -6000.0, 0.0)
74 p7=Point(5000.0, -2400.0, 0.0)
75
76 # Create TOP BLOCK
77 tbl1=Line(p0,p4)
78 tbl2=Line(p4,p5)
79 tbl3=Line(p5,p7)
80 tbl4=Line(p7,p3)
81 tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)
82 tblock = PlaneSurface(tblockloop)
83
84 # Create BOTTOM BLOCK LEFT
85 bbll1=Line(p4,p1)
86 bbll2=Line(p1,p6)
87 bbll3=Line(p6,p5)
88 bbll4=-tbl2
89 bblockloopl = CurveLoop(bbll1,bbll2,bbll3,bbll4)
90 bblockl = PlaneSurface(bblockloopl)
91
92 # Create BOTTOM BLOCK RIGHT
93 bbrl1=Line(p6,p2)
94 bbrl2=Line(p2,p7)
95 bbrl3=-tbl3
96 bbrl4=-bbll3
97 bblockloopr = CurveLoop(bbrl1,bbrl2,bbrl3,bbrl4)
98 bblockr = PlaneSurface(bblockloopr)
99
100 #############################################EXPORTING MESH FOR ESCRIPT
101 # Create a Design which can make the mesh
102 d=Design(dim=2, element_size=200)
103 # Add the subdomains and flux boundaries.
104 d.addItems(PropertySet("top",tblock),\
105 PropertySet("bottomleft",bblockl),\
106 PropertySet("bottomright",bblockr),\
107 PropertySet("linebottom",bbll2, bbrl1))
108
109 # Create the geometry, mesh and Escript domain
110 d.setScriptFileName(os.path.join(save_path,"example06.geo"))
111
112 d.setMeshFileName(os.path.join(save_path,"example06.msh"))
113 domain=MakeDomain(d)
114 print "Domain has been generated ..."
115
116 # set up kappa (thermal conductivity across domain) using tags
117 kappa=Scalar(0,Function(domain))
118 kappa.setTaggedValue("top",2.0)
119 kappa.setTaggedValue("bottomleft",10.0)
120 kappa.setTaggedValue("bottomright",6.0)
121 ##############################################################SOLVE PDE
122 mypde=LinearPDE(domain)
123 mypde.getSolverOptions().setVerbosityOn()
124 mypde.setSymmetryOn()
125 mypde.setValue(A=kappa*kronecker(domain))
126 x=Solution(domain).getX()
127 mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
128 qS=Scalar(0,FunctionOnBoundary(domain))
129 qS.setTaggedValue("linebottom",qin)
130 mypde.setValue(y=qS)
131 print "PDE has been generated ..."
132 ###########################################################GET SOLUTION
133 T=mypde.getSolution()
134 print "PDE has been solved ..."
135 ###############################################################PLOTTING
136 # show temperature:
137 xi, yi, zi = toRegGrid(T, nx=50, ny=50)
138 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
139 pl.clabel(CS, inline=1, fontsize=8)
140 # show sub domains:
141 tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])
142 pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)
143 bpgr=np.array([p.getCoordinates() for p in bblockloopr.getPolygon() ])
144 pl.fill(bpgr[:,0],bpgr[:,1],'orange',label='6 W/m/k',zorder=-1000)
145 bpgl=np.array([p.getCoordinates() for p in bblockloopl.getPolygon() ])
146 pl.fill(bpgl[:,0],bpgl[:,1],'red',label='10 W/m/k',zorder=-1000)
147 # show flux:
148 xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)
149 pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")
150 # create plot
151 pl.title("Heat Refraction across an anisotropic structure\n with flux vectors.")
152 pl.xlabel("Horizontal Displacement (m)")
153 pl.ylabel("Depth (m)")
154 pl.legend()
155 pl.savefig(os.path.join(save_path,"flux.png"))
156 print "Flux has been plotted ..."

  ViewVC Help
Powered by ViewVC 1.1.26