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

  ViewVC Help
Powered by ViewVC 1.1.26