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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 141 - (show annotations)
Mon Jul 25 05:15:41 2005 UTC (13 years, 8 months ago) by jgs
File MIME type: text/x-python
File size: 6663 byte(s)
fixed to work with new python source/install directory structure

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
14 from escript.escript import *
15 from escript.linearPDEs import *
16 from finley import finley
17
18 from math import *
19 from numarray import array
20
21 numElements=10
22 max_error=0.
23 max_error_text=""
24
25 for dim in [2,3]:
26 for order in [1,2]:
27 for onFace in [True,False]:
28
29 print "\n"
30 print "-----------------------------------------------------------------------------------------------------"
31 print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace)
32 print "-----------------------------------------------------------------------------------------------------"
33
34 if onFace:
35 onFaceText=", on elements"
36 else:
37 onFaceText=""
38
39 case="dim=%d, order=%d%s" % (dim,order,onFaceText)
40
41 if dim==2:
42 mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)
43 m00=[[1,0],[0,0]]
44 m01=[[0,1],[0,0]]
45 m11=[[0,0],[0,1]]
46 h=5
47 else:
48 mydomain=finley.Brick(numElements,numElements,numElements,order=order,useElementsOnFace=onFace)
49 m00=[[1,0,0],[0,0,0]]
50 m01=[[0,1,0],[0,0,0]]
51 m11=[[0,0,0],[0,1,0]]
52 h=7
53
54 n=ContinuousFunction(mydomain)
55 e=Function(mydomain)
56 f=FunctionOnBoundary(mydomain)
57 d=Solution(mydomain)
58 r=ReducedSolution(mydomain)
59
60 #
61 # test gradient
62 #
63
64 test="error gradient in interior (nodes)"
65
66 x=n.getX()[0:2]
67 g=grad(x**order+x[1]*[1,0])
68 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
69 error_norm=Lsup(ref-g)
70
71 text="%s: %55s = %e" % (case, test, error_norm)
72 print "%s" % text
73
74 if error_norm>max_error:
75 max_error_text=text
76 max_error=error_norm
77
78 #
79 # test gradient on degrees of freedom
80 #
81
82 test="error gradient in interior (degrees of freedom)"
83
84 x=n.getX()[0:2].interpolate(d)
85 g=grad(x**order+x[1]*[1,0])
86 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
87 error_norm=Lsup(ref-g)/Lsup(ref)
88
89 text="%s: %55s = %e" % (case, test, error_norm)
90 print "%s" % text
91
92 if error_norm>max_error:
93 max_error_text=text
94 max_error=error_norm
95
96 #
97 # test gradient on reduced degrees of freedom
98 #
99
100 test="error gradient in interior (reduced degrees of freedom)"
101
102 x=n.getX()[0:2].interpolate(r)
103 g=grad(x+x[1]*[1,0])
104 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
105 error_norm=Lsup(ref-g)/Lsup(ref)
106
107 text="%s: %55s = %e" % (case, test, error_norm)
108 print "%s" % text
109
110 if error_norm>max_error:
111 max_error_text=text
112 max_error=error_norm
113
114 #
115 # test integration over volume
116 #
117
118 test="error volume integration"
119
120 x=e.getX()[0:2]
121 error=integrate(x**2+[0,2.]*x)-array([1./3.,1./3.+2*1./2.])
122 error_norm=sqrt(numarray.innerproduct(error,error))
123
124 text="%s: %55s = %e" % (case, test, error_norm)
125 print "%s" % text
126
127 if error_norm>max_error:
128 max_error_text=text
129 max_error=error_norm
130
131 if onFace:
132
133 #
134 # gradient on the boundary:
135 #
136
137 test="error gradient on boundary"
138
139 x=n.getX()[0:2]
140 g=grad(x**order+x[1]*[1,0],where=f)
141 x=f.getX()[0:2]
142 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
143 error_norm=Lsup(g-ref)
144
145 text="%s: %55s = %e" % (case, test, error_norm)
146 print "%s" % text
147
148 if error_norm>max_error:
149 max_error_text=text
150 max_error=error_norm
151
152 #
153 # test gradient on degrees of freedom
154 #
155
156 test="error gradient on boundary (degrees of freedom)"
157
158 x=n.getX()[0:2].interpolate(d)
159 g=grad(x**order+x[1]*[1,0],where=f)
160 x=f.getX()[0:2]
161 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
162 error_norm=Lsup(ref-g)/Lsup(ref)
163
164 text="%s: %55s = %e" % (case, test, error_norm)
165 print "%s" % text
166
167 if error_norm>max_error:
168 max_error_text=text
169 max_error=error_norm
170
171 #
172 # test gradient on reduced degrees of freedom
173 #
174
175 test="error gradient on boundary (reduced degrees of freedom)"
176
177 x=n.getX()[0:2].interpolate(r)
178 g=grad(x+x[1]*[1,0],where=f)
179 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
180 error_norm=Lsup(ref-g)/Lsup(ref)
181
182 text="%s: %55s = %e" % (case, test, error_norm)
183 print "%s" % text
184
185 if error_norm>max_error:
186 max_error_text=text
187 max_error=error_norm
188
189 #
190 # test integration over boundary
191 #
192
193 test="error boundary integration"
194
195 x=f.getX()[0:2]
196 error=integrate(x**2+[0,2.]*x)-array([h/3.,h/3.+2*(h-1)/2.])
197 error_norm=sqrt(numarray.innerproduct(error,error))
198
199 text="%s: %55s = %e" % (case, test, error_norm)
200 print "%s" % text
201
202 if error_norm>max_error:
203 max_error_text=text
204 max_error=error_norm
205
206 #
207 # normal test:
208 #
209
210 test="error normals"
211
212 refNormal=Vector(0,what=f)
213 if dim==3:
214 refNormal.setTaggedValue(2,[1,0,0])
215 refNormal.setTaggedValue(1,[-1,0,0])
216 refNormal.setTaggedValue(20,[0,1,0])
217 refNormal.setTaggedValue(10,[0,-1,0])
218 refNormal.setTaggedValue(200,[0,0,1])
219 refNormal.setTaggedValue(100,[0,0,-1])
220 else:
221 refNormal.setTaggedValue(2,[1,0])
222 refNormal.setTaggedValue(1,[-1,0])
223 refNormal.setTaggedValue(20,[0,1])
224 refNormal.setTaggedValue(10,[0,-1])
225 error_norm=Lsup(f.getNormal()-refNormal)
226
227 text="%s: %55s = %e" % (case, test, error_norm)
228 print "%s" % text
229
230 if error_norm>max_error:
231 max_error_text=text
232 max_error=error_norm
233
234 print "-----------------------------------------------------------------------------------------------------"
235
236 print "\n\n"
237 print "******************************************************************************************************************"
238 print "maximal error:", max_error_text
239 print "******************************************************************************************************************"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26