/[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 82 - (show annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years, 5 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/GradTest.py
File MIME type: text/x-python
File size: 6791 byte(s)
Initial revision

1 """
2 Test a grad, interpolate and integrate over the unit square.
3
4 The tests are very basic.
5
6 by Lutz Gross, ACcESS, University of Queensland, Australia, 2003.
7 Version $Id$
8 """
9
10 import sys
11 import os
12 import unittest
13
14
15 esys_root=os.getenv('ESYS_ROOT')
16 sys.path.append(esys_root+'/finley/lib')
17 sys.path.append(esys_root+'/escript/lib')
18 sys.path.append(esys_root+'/escript/py_src')
19
20 from escript import *
21 from util import *
22 import finley
23
24 from math import *
25 from numarray import array
26
27
28 numElements=10
29
30 max_error=0.
31 max_error_text=""
32
33 for dim in [2,3]:
34 for order in [1,2]:
35 for onFace in [TRUE,FALSE]:
36
37 print "\ndim: %d order: %i onFace: %s" % (dim, order, onFace)
38 print "-------------------------------------------------------------------------------------"
39
40 if onFace:
41 onFaceText=", on elements"
42 else:
43 onFaceText=""
44 case="dim=%d, order=%d%s"%(dim,order,onFaceText)
45
46 if dim==2:
47 mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)
48 m00=[[1,0],[0,0]]
49 m01=[[0,1],[0,0]]
50 m11=[[0,0],[0,1]]
51 h=5
52 else:
53 mydomain=finley.Brick(numElements,numElements,numElements,order=order,useElementsOnFace=onFace)
54 m00=[[1,0,0],[0,0,0]]
55 m01=[[0,1,0],[0,0,0]]
56 m11=[[0,0,0],[0,1,0]]
57 h=7
58
59 n=ContinuousFunction(mydomain)
60 e=Function(mydomain)
61 f=FunctionOnBoundary(mydomain)
62 d=Solution(mydomain)
63 r=ReducedSolution(mydomain)
64
65 #
66 # test gradient
67 #
68
69 x=n.getX()[0:2]
70 g=grad(x**order+x[1]*[1,0])
71 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
72 error_norm=Lsup(ref-g)
73 text=" %s: error gradient in interior (nodes) = %e"%(case,error_norm)
74 print "\t%s"%text
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 x=n.getX()[0:2].interpolate(d)
84 g=grad(x**order+x[1]*[1,0])
85 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
86 error_norm=Lsup(ref-g)
87 text=" %s: error gradient in interior (degrees of freedom) = %e"%(case,error_norm)
88 print "\t%s"%text
89 if error_norm>max_error:
90 max_error_text=text
91 max_error=error_norm
92
93 #
94 # test gradient on reduced degrees of freedom
95 #
96
97 x=n.getX()[0:2].interpolate(r)
98 g=grad(x+x[1]*[1,0])
99 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
100 error_norm=Lsup(ref-g)
101 text=" %s: error gradient in interior (reduced degrees of freedom) = %e"%(case,error_norm)
102 print "\t%s"%text
103 if error_norm>max_error:
104 max_error_text=text
105 max_error=error_norm
106
107 #
108 # test integration over volume
109 #
110 # integrate causes: Fatal Python error: Call to numarray API function
111 # without first calling import_libnumarray() in src/Data/Data.cpp
112
113 x=e.getX()[0:2]
114 #error=integrate(x**2+[0,2.]*x)-array([1./3.,1./3.+2*1./2.])
115 #error_norm=sqrt(numarray.innerproduct(error,error))
116 text=" %s: error volume integration= %e"%(case,error_norm)
117 print "\t%s"%text
118 if error_norm>max_error:
119 max_error_text=text
120 max_error=error_norm
121
122 if onFace:
123
124 #
125 # gradient on the boundary:
126 #
127 # Lsup fails - perhaps grad(what=f) is needed?
128
129 x=n.getX()[0:2]
130 #g=grad(x**order+x[1]*[1,0],what=f)
131 g=grad(x**order+x[1]*[1,0])
132 x=f.getX()[0:2]
133 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
134 #error_norm=Lsup(g-ref)
135 error_norm=0
136 text=" %s: error gradient on boundary = %e"%(case,error_norm)
137 print "\t%s"%text
138 if error_norm>max_error:
139 max_error_text=text
140 max_error=error_norm
141
142 #
143 # test gradient on degrees of freedom
144 #
145 # Lsup fails - perhaps grad(what=f) is needed?
146
147 x=n.getX()[0:2].interpolate(d)
148 #g=grad(x**order+x[1]*[1,0],what=f)
149 g=grad(x**order+x[1]*[1,0])
150 x=f.getX()[0:2]
151 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
152 #error_norm=Lsup(ref-g)
153 error_norm=0
154 text=" %s: error gradient on boundary (degrees of freedom) = %e"%(case,error_norm)
155 print "\t%s"%text
156 if error_norm>max_error:
157 max_error_text=text
158 max_error=error_norm
159
160 #
161 # test gradient on reduced degrees of freedom
162 #
163 # Lsup fails - perhaps grad(what=f) is needed?
164
165 x=n.getX()[0:2].interpolate(r)
166 #g=grad(x+x[1]*[1,0],what=f)
167 g=grad(x+x[1]*[1,0])
168 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
169 #error_norm=Lsup(ref-g)
170 error_norm=0
171 text=" %s: error gradient on boundary (reduced degrees of freedom) = %e"%(case,error_norm)
172 print "\t%s"%text
173 if error_norm>max_error:
174 max_error_text=text
175 max_error=error_norm
176
177 #
178 # test integration over boundary
179 #
180 # integrate causes: Fatal Python error: Call to numarray API function
181 # without first calling import_libnumarray() in src/Data/Data.cpp
182
183 x=f.getX()[0:2]
184 #error=integrate(x**2+[0,2.]*x)-array([h/3.,h/3.+2*(h-1)/2.])
185 #error_norm=sqrt(numarray.innerproduct(error,error))
186 text=" %s: error boundary integration= %e"%(case,error_norm)
187 print "\t%s"%text
188 if error_norm>max_error:
189 max_error_text=text
190 max_error=error_norm
191
192 #
193 # normal test:
194 #
195 # need to add wrapper for DataTagged::addTaggedValue
196 #
197
198 #refNormal=Vector(0,what=f)
199 #if dim==3:
200 # refNormal.addTaggedValue(2,[1,0,0])
201 # refNormal.addTaggedValue(1,[-1,0,0])
202 # refNormal.addTaggedValue(20,[0,1,0])
203 # refNormal.addTaggedValue(10,[0,-1,0])
204 # refNormal.addTaggedValue(200,[0,0,1])
205 # refNormal.addTaggedValue(100,[0,0,-1])
206 #else:
207 # refNormal.addTaggedValue(2,[1,0])
208 # refNormal.addTaggedValue(1,[-1,0])
209 # refNormal.addTaggedValue(20,[0,1])
210 # refNormal.addTaggedValue(10,[0,-1])
211 #error_norm=Lsup(f.getNormal()-refNormal)
212 text=" %s: error normals= %e"%(case,error_norm)
213 print "\t%s"%text
214 if error_norm>max_error:
215 max_error_text=text
216 max_error=error_norm
217
218 print "-------------------------------------------------------------------------------------"
219
220 print "\n\nmaximal error for",max_error_text

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26