/[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 6798 - (hide annotations)
Wed Mar 6 06:25:41 2019 UTC (7 weeks, 1 day ago) by aellery
File MIME type: text/x-python
File size: 4335 byte(s)
Removed the natgrid references from the cookbook and changed calls to matplotlib.griddata to scipy.interpolate.griddata (as the matplotlib griddata function is deprecated). Also, updated the cookbook to reflect these changes.

This closes bug 448 (problems with calls to matplotlib griddata.)



1 jfenwick 3981 ##############################################################################
2 gross 2948 #
3 jfenwick 6651 # Copyright (c) 2009-2018 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 gross 2948 #
6     # Primary Business: Queensland, Australia
7 jfenwick 6112 # Licensed under the Apache License, version 2.0
8     # http://www.apache.org/licenses/LICENSE-2.0
9 gross 2948 #
10 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 gross 2948
16 sshaw 5775 from __future__ import division, print_function
17 jfenwick 6651 __copyright__="""Copyright (c) 2009-2018 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 gross 2948 Primary Business: Queensland, Australia"""
20 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
21     http://www.apache.org/licenses/LICENSE-2.0"""
22 gross 2948 __url__="https://launchpad.net/escript-finley"
23    
24     """
25     Author: Antony Hallam antony.hallam@uqconnect.edu.au
26     """
27    
28     ############################################################FILE HEADER
29 gross 2951 # example04b.py
30 gross 2948 # Create either a 2D mesh for a rectangle model using pycad meshing
31     # tools.
32     #
33     #######################################################EXTERNAL MODULES
34 caltinay 4087 import matplotlib
35     matplotlib.use('agg') #It's just here for automated testing
36 gross 2948 from esys.pycad import * #domain constructor
37     from esys.pycad.gmsh import Design #Finite Element meshing package
38     from esys.escript import *
39     from esys.escript.unitsSI import *
40     from esys.escript.linearPDEs import LinearPDE
41     import pylab as pl #Plotting package
42 aellery 6798 from cblib import toRegGrid
43 gross 2948 import os
44 sshaw 5288 try:
45     # This imports the rectangle domain function
46     from esys.finley import MakeDomain#Converter for escript
47     HAVE_FINLEY = True
48     except ImportError:
49     print("Finley module not available")
50     HAVE_FINLEY = False
51 gross 2948
52     ########################################################MPI WORLD CHECK
53     if getMPISizeWorld() > 1:
54 sshaw 5705 import sys
55     print("This example will not run in an MPI world.")
56     sys.exit(0)
57 gross 2948
58 aellery 6798 if HAVE_FINLEY:
59 sshaw 5288 # make sure path exists
60     save_path= os.path.join("data","example04")
61     mkDir(save_path)
62 gross 2948
63 sshaw 5288 ################################################ESTABLISHING PARAMETERS
64     #Model Parameters
65     width=5000.0*m #width of model
66     depth=-6000.0*m #depth of model
67     kappa=2.0*W/m/K # watts/m.Kthermal conductivity
68     Ttop=20*K # top temperature
69     qin=70*Milli*W/(m*m) # bottom heat influx
70 gross 2948
71 sshaw 5288 ####################################################DOMAIN CONSTRUCTION
72     # Domain Corners
73     p0=Point(0.0, 0.0, 0.0)
74     p1=Point(0.0, depth, 0.0)
75     p2=Point(width, depth, 0.0)
76     p3=Point(width, 0.0, 0.0)
77     # Join corners in anti-clockwise manner.
78     l01=Line(p0, p1)
79     l12=Line(p1, p2)
80     l23=Line(p2, p3)
81     l30=Line(p3, p0)
82     # Join line segments to create domain boundary.
83     c=CurveLoop(l01, l12, l23, l30)
84     # surface
85     rec = PlaneSurface(c)
86 gross 2948
87 sshaw 5288 #############################################EXPORTING MESH FOR ESCRIPT
88     # Create a Design which can make the mesh
89     d=Design(dim=2, element_size=200*m)
90     # Add the subdomains and flux boundaries.
91     d.addItems(rec, PropertySet("linebottom",l12))
92     d.addItems(l01, l23, l30) # just in case we need them
93     #############################################MAKE THE DOMAIN
94     domain=MakeDomain(d, optimizeLabeling=True)
95     print("Domain has been generated ...")
96     ##############################################################SOLVE PDE
97     mypde=LinearPDE(domain)
98     mypde.getSolverOptions().setVerbosityOn()
99     mypde.setSymmetryOn()
100     mypde.setValue(A=kappa*kronecker(domain))
101     x=Solution(domain).getX()
102     mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
103     qS=Scalar(0,FunctionOnBoundary(domain))
104     qS.setTaggedValue("linebottom",qin)
105     mypde.setValue(y=-qS)
106     print("PDE has been generated ...")
107     ###########################################################GET SOLUTION
108     T=mypde.getSolution()
109     print("PDE has been solved ...")
110     ###########################################################
111     xi, yi, zi = toRegGrid(T, nx=50, ny=50)
112     pl.matplotlib.pyplot.autumn()
113     pl.contourf(xi,yi,zi,10)
114     pl.xlabel("Horizontal Displacement (m)")
115     pl.ylabel("Depth (m)")
116     pl.savefig(os.path.join(save_path,"example04.png"))
117     print("Solution has been plotted ...")

  ViewVC Help
Powered by ViewVC 1.1.26