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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3892 - (hide annotations)
Tue Apr 10 08:57:23 2012 UTC (7 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 5374 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



1 gross 2952
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 jfenwick 3892 print("This example will not run in an MPI world.")
48 gross 2952 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 jfenwick 3892 print("Domain has been generated ...")
115 gross 2952
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 ahallam 2977 ##############################################################SOLVE PDE
122 gross 2952 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 jfenwick 3892 print("PDE has been generated ...")
132 gross 2952 ###########################################################GET SOLUTION
133     T=mypde.getSolution()
134 jfenwick 3892 print("PDE has been solved ...")
135 ahallam 2977 ###############################################################PLOTTING
136 gross 2952 # 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 jfenwick 3892 print("Flux has been plotted ...")

  ViewVC Help
Powered by ViewVC 1.1.26