/[escript]/trunk/ripley/test/cook_1.py
ViewVC logotype

Contents of /trunk/ripley/test/cook_1.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3479 - (show annotations)
Wed Mar 23 04:03:22 2011 UTC (8 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 3337 byte(s)


1 #!/usr/bin/env python
2 """
3 Illustrate simple contour plotting, contours on an image with
4 a colorbar for the contours, and labelled contours.
5
6 See also contour_image.py.
7 """
8 import matplotlib
9 matplotlib.use('Agg')
10 import numpy as np
11 import matplotlib.cm as cm
12 import matplotlib.mlab as mlab
13 import matplotlib.pyplot as plt
14
15 from esys.escript import *
16 from esys.finley import Rectangle
17 from esys.escript.linearPDEs import LinearSinglePDE
18
19
20 matplotlib.rcParams['xtick.direction'] = 'out'
21 matplotlib.rcParams['ytick.direction'] = 'out'
22
23 NE=100
24 H=1.
25 L=1.
26
27
28 class PoissonTest(object):
29 def __init__(self, domain, f0, r, xc):
30 self.domain=domain
31 DIM=self.domain.getDim()
32 self.pde=LinearSinglePDE(self.domain)
33 x=self.domain.getX()
34 self.pde.setValue(A=kronecker(domain), q=whereZero(x[DIM-1]-inf(x[DIM-1])))
35 self.u_ref=self.getU(f0, r, xc)
36 self.f0=f0
37 self.r=r
38 self.xc=xc
39 self.omega=1.
40
41 def getU(self, f0, r, xc):
42 x=Function(self.domain).getX()
43 self.pde.setValue(Y=exp(f0-length(x-xc)**2/r**2))
44 u=self.pde.getSolution()
45 print "evaluated for f0, r, xc, u=",f0, r, xc,inf(u),sup(u)
46 return u
47
48 def getCostFunction(self, f0, r, xc):
49 u=interpolate(self.getU(f0, r, xc), FunctionOnBoundary(self.domain))
50 u_ref=interpolate(self.u_ref, FunctionOnBoundary(self.domain))
51 C=0.5*integrate(self.omega*length(u-u_ref)**2)
52 return C
53
54 def plot(X, Y, Z, filename="test.png", title="None"):
55
56 plt.figure()
57 im = plt.imshow(Z, interpolation='bilinear', origin='lower',
58 cmap=cm.gray,
59 extent=(0,1.,0,1.))
60
61 levels = np.linspace(np.amin(Z), np.amax(Z), 20)
62 CS = plt.contour(X, Y, Z,levels,
63 origin='lower',
64 linewidths=2)
65
66 #Thicken the zero contour.
67 zc = CS.collections[6]
68 plt.setp(zc, linewidth=4)
69
70 plt.clabel(CS, levels[1::2], # label every second level
71 inline=1,
72 fmt='%2.2e',
73 fontsize=14)
74
75 # make a colorbar for the contour lines
76 #CB = plt.colorbar(CS, shrink=0.8, extend='both')
77
78 plt.title(title)
79 plt.xlabel('r')
80 plt.ylabel('f_0')
81 #plt.hot() # Now change the colormap for the contour lines and colorbar
82 plt.flag()
83
84 # We can still add a colorbar for the image, too.
85 # CBI = plt.colorbar(im, orientation='horizontal', shrink=0.8)
86
87 # This makes the original colorbar look a bit out of place,
88 # so let's improve its position.
89
90 # l,b,w,h = plt.gca().get_position().bounds
91 # ll,bb,ww,hh = CB.ax.get_position().bounds
92 # CB.ax.set_position([ll, b+0.1*h, ww, h*0.8])
93
94
95 plt.savefig(filename, dpi=100)
96
97 def runTestOverGrid(p):
98 delta=0.1
99 # r = np.arange(max(delta,2./NE), 1.0+delta, delta)
100 r = np.arange(0., 1.0+delta, delta)
101 f0 = np.arange(0, 1.0+delta, delta)
102 R, F0= np.meshgrid(r, f0)
103 C=R.copy()
104 for i in range(C.shape[0]):
105 for j in range(C.shape[1]):
106 # C[i,j]=p.getCostFunction(F0[i,j], R[i,j], p.xc)
107 C[i,j]=p.getCostFunction(p.f0, p.r, (R[i,j],F0[i,j]))
108 plot(R,F0,C, "test.png", "test")
109
110 p1=PoissonTest(Rectangle(NE,int(NE*H/L),l0=L,l1=H), f0=0.3, r=0.1, xc=(0.5,H*0.5))
111 runTestOverGrid(p1)

  ViewVC Help
Powered by ViewVC 1.1.26