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

Contents of /trunk/doc/examples/cookbook/example04b.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: 3703 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
26 ############################################################FILE HEADER
27 # example04b.py
28 # Create either a 2D mesh for a rectangle model using pycad meshing
29 # tools.
30 #
31 #######################################################EXTERNAL MODULES
32 from esys.pycad import * #domain constructor
33 from esys.pycad.gmsh import Design #Finite Element meshing package
34 from esys.finley import MakeDomain #Converter for escript
35 from esys.escript import *
36 from esys.escript.unitsSI import *
37 from esys.escript.linearPDEs import LinearPDE
38 import matplotlib
39 matplotlib.use('agg') #It's just here for automated testing
40 import pylab as pl #Plotting package
41 from cblib import toRegGrid
42 import os
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
50 # make sure path exists
51 save_path= os.path.join("data","example04")
52 mkDir(save_path)
53
54 ################################################ESTABLISHING PARAMETERS
55 #Model Parameters
56 width=5000.0*m #width of model
57 depth=-6000.0*m #depth of model
58 kappa=2.0*W/m/K # watts/m.Kthermal conductivity
59 Ttop=20*K # top temperature
60 qin=70*Milli*W/(m*m) # bottom heat influx
61
62 ####################################################DOMAIN CONSTRUCTION
63 # Domain Corners
64 p0=Point(0.0, 0.0, 0.0)
65 p1=Point(0.0, depth, 0.0)
66 p2=Point(width, depth, 0.0)
67 p3=Point(width, 0.0, 0.0)
68 # Join corners in anti-clockwise manner.
69 l01=Line(p0, p1)
70 l12=Line(p1, p2)
71 l23=Line(p2, p3)
72 l30=Line(p3, p0)
73 # Join line segments to create domain boundary.
74 c=CurveLoop(l01, l12, l23, l30)
75 # surface
76 rec = PlaneSurface(c)
77
78 #############################################EXPORTING MESH FOR ESCRIPT
79 # Create a Design which can make the mesh
80 d=Design(dim=2, element_size=200*m)
81 # Add the subdomains and flux boundaries.
82 d.addItems(rec, PropertySet("linebottom",l12))
83 #############################################MAKE THE FINLEY DOMAIN
84 domain=MakeDomain(d, optimizeLabeling=True)
85 print "Domain has been generated ..."
86 ##############################################################SOLVE PDE
87 mypde=LinearPDE(domain)
88 mypde.getSolverOptions().setVerbosityOn()
89 mypde.setSymmetryOn()
90 mypde.setValue(A=kappa*kronecker(domain))
91 x=Solution(domain).getX()
92 mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
93 qS=Scalar(0,FunctionOnBoundary(domain))
94 qS.setTaggedValue("linebottom",qin)
95 mypde.setValue(y=-qS)
96 print "PDE has been generated ..."
97 ###########################################################GET SOLUTION
98 T=mypde.getSolution()
99 print "PDE has been solved ..."
100 ###########################################################
101 xi, yi, zi = toRegGrid(T, nx=50, ny=50)
102 pl.matplotlib.pyplot.autumn()
103 pl.contourf(xi,yi,zi,10)
104 pl.xlabel("Horizontal Displacement (m)")
105 pl.ylabel("Depth (m)")
106 pl.savefig(os.path.join(save_path,"example04.png"))
107 print "Solution has been plotted ..."

  ViewVC Help
Powered by ViewVC 1.1.26