/[escript]/trunk/finley/test/python/GradTest.py
ViewVC logotype

Contents of /trunk/finley/test/python/GradTest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 428 - (show annotations)
Wed Jan 11 01:38:16 2006 UTC (14 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 6911 byte(s)
interplotared is changed from method to function call
1 # $Id$
2 """
3 Test a grad, interpolate and integrate over the unit square.
4
5 The tests are very basic.
6
7 by Lutz Gross, ACcESS, University of Queensland, Australia, 2003.
8 """
9
10 import sys
11 import os
12 import unittest
13 import time
14
15 from esys.escript import *
16 from esys.escript.linearPDEs import *
17 from esys import finley
18
19 from math import *
20 from numarray import array
21
22 starttime = time.clock()
23
24 numElements=10
25 max_error=0.
26 max_error_text=""
27 error_tol=pow(10,-10)
28
29 for dim in [2,3]:
30 for order in [1,2]:
31 for onFace in [True,False]:
32
33 print "\n"
34 print "-----------------------------------------------------------------------------------------------------"
35 print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace)
36 print "-----------------------------------------------------------------------------------------------------"
37
38 if onFace:
39 onFaceText=", on elements"
40 else:
41 onFaceText=""
42
43 case="dim=%d, order=%d%s" % (dim,order,onFaceText)
44
45 if dim==2:
46 mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)
47 m00=[[1,0],[0,0]]
48 m01=[[0,1],[0,0]]
49 m11=[[0,0],[0,1]]
50 h=5
51 else:
52 mydomain=finley.Brick(numElements,numElements,numElements,order=order,useElementsOnFace=onFace)
53 m00=[[1,0,0],[0,0,0]]
54 m01=[[0,1,0],[0,0,0]]
55 m11=[[0,0,0],[0,1,0]]
56 h=7
57
58 n=ContinuousFunction(mydomain)
59 e=Function(mydomain)
60 f=FunctionOnBoundary(mydomain)
61 d=Solution(mydomain)
62 r=ReducedSolution(mydomain)
63
64 #
65 # test gradient
66 #
67
68 test="error gradient in interior (nodes)"
69
70 x=n.getX()[0:2]
71 g=grad(x**order+x[1]*[1,0])
72 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
73 error_norm=Lsup(ref-g)
74
75 text="%s: %55s = %e" % (case, test, error_norm)
76 print "%s" % text
77
78 if error_norm>max_error:
79 max_error_text=text
80 max_error=error_norm
81
82 #
83 # test gradient on degrees of freedom
84 #
85
86 test="error gradient in interior (degrees of freedom)"
87
88 x=interpolate(n.getX()[0:2],d)
89 g=grad(x**order+x[1]*[1,0])
90 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
91 error_norm=Lsup(ref-g)/Lsup(ref)
92
93 text="%s: %55s = %e" % (case, test, error_norm)
94 print "%s" % text
95
96 if error_norm>max_error:
97 max_error_text=text
98 max_error=error_norm
99
100 #
101 # test gradient on reduced degrees of freedom
102 #
103
104 test="error gradient in interior (reduced degrees of freedom)"
105
106 x=interpolate(n.getX()[0:2],r)
107 g=grad(x+x[1]*[1,0])
108 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
109 error_norm=Lsup(ref-g)/Lsup(ref)
110
111 text="%s: %55s = %e" % (case, test, error_norm)
112 print "%s" % text
113
114 if error_norm>max_error:
115 max_error_text=text
116 max_error=error_norm
117
118 #
119 # test integration over volume
120 #
121
122 test="error volume integration"
123
124 x=e.getX()[0:2]
125 error=integrate(x**2+[0,2.]*x)-array([1./3.,1./3.+2*1./2.])
126 error_norm=sqrt(numarray.innerproduct(error,error))
127
128 text="%s: %55s = %e" % (case, test, error_norm)
129 print "%s" % text
130
131 if error_norm>max_error:
132 max_error_text=text
133 max_error=error_norm
134
135 if onFace:
136
137 #
138 # gradient on the boundary:
139 #
140
141 test="error gradient on boundary"
142
143 x=n.getX()[0:2]
144 g=grad(x**order+x[1]*[1,0],where=f)
145 x=f.getX()[0:2]
146 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
147 error_norm=Lsup(g-ref)
148
149 text="%s: %55s = %e" % (case, test, error_norm)
150 print "%s" % text
151
152 if error_norm>max_error:
153 max_error_text=text
154 max_error=error_norm
155
156 #
157 # test gradient on degrees of freedom
158 #
159
160 test="error gradient on boundary (degrees of freedom)"
161
162 x=interpolate(n.getX()[0:2],d)
163 g=grad(x**order+x[1]*[1,0],where=f)
164 x=f.getX()[0:2]
165 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
166 error_norm=Lsup(ref-g)/Lsup(ref)
167
168 text="%s: %55s = %e" % (case, test, error_norm)
169 print "%s" % text
170
171 if error_norm>max_error:
172 max_error_text=text
173 max_error=error_norm
174
175 #
176 # test gradient on reduced degrees of freedom
177 #
178
179 test="error gradient on boundary (reduced degrees of freedom)"
180
181 x=interpolate(n.getX()[0:2],r)
182 g=grad(x+x[1]*[1,0],where=f)
183 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
184 error_norm=Lsup(ref-g)/Lsup(ref)
185
186 text="%s: %55s = %e" % (case, test, error_norm)
187 print "%s" % text
188
189 if error_norm>max_error:
190 max_error_text=text
191 max_error=error_norm
192
193 #
194 # test integration over boundary
195 #
196
197 test="error boundary integration"
198
199 x=f.getX()[0:2]
200 error=integrate(x**2+[0,2.]*x)-array([h/3.,h/3.+2*(h-1)/2.])
201 error_norm=sqrt(numarray.innerproduct(error,error))
202
203 text="%s: %55s = %e" % (case, test, error_norm)
204 print "%s" % text
205
206 if error_norm>max_error:
207 max_error_text=text
208 max_error=error_norm
209
210 #
211 # normal test:
212 #
213
214 test="error normals"
215
216 refNormal=Vector(0,what=f)
217 if dim==3:
218 refNormal.setTaggedValue(2,[1,0,0])
219 refNormal.setTaggedValue(1,[-1,0,0])
220 refNormal.setTaggedValue(20,[0,1,0])
221 refNormal.setTaggedValue(10,[0,-1,0])
222 refNormal.setTaggedValue(200,[0,0,1])
223 refNormal.setTaggedValue(100,[0,0,-1])
224 else:
225 refNormal.setTaggedValue(2,[1,0])
226 refNormal.setTaggedValue(1,[-1,0])
227 refNormal.setTaggedValue(20,[0,1])
228 refNormal.setTaggedValue(10,[0,-1])
229 error_norm=Lsup(f.getNormal()-refNormal)
230
231 text="%s: %55s = %e" % (case, test, error_norm)
232 print "%s" % text
233
234 if error_norm>max_error:
235 max_error_text=text
236 max_error=error_norm
237
238 print "-----------------------------------------------------------------------------------------------------"
239
240 print "\n\n"
241 print "******************************************************************************************************************"
242 print "maximal error:", max_error_text
243 print "******************************************************************************************************************"
244
245 stoptime = time.clock()
246 elapsed = stoptime - starttime
247 print "\nElapsed time: ", elapsed, "\n"
248
249 if max_error > error_tol:
250 print "max error exceeds tolerance"
251 sys.exit(1)
252
253 sys.exit(0)

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26