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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (hide annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 3696 byte(s)
Merging dudley and scons updates from branches

1 gross 2948
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 gross 2951 # example04b.py
28 gross 2948 # 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 gross 2951 save_path= os.path.join("data","example04")
52 gross 2948 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 jfenwick 3259 #############################################MAKE THE DOMAIN
84 gross 2948 domain=MakeDomain(d, optimizeLabeling=True)
85     print "Domain has been generated ..."
86 ahallam 2977 ##############################################################SOLVE PDE
87 gross 2948 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 gross 2951 pl.savefig(os.path.join(save_path,"example04.png"))
107 gross 2948 print "Solution has been plotted ..."

  ViewVC Help
Powered by ViewVC 1.1.26