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

  ViewVC Help
Powered by ViewVC 1.1.26