/[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 149 - (show annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 10 months ago) by jgs
File MIME type: text/x-python
File size: 6777 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26