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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 97 - (hide annotations)
Tue Dec 14 05:39:33 2004 UTC (15 years, 3 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/GradTest.py
File MIME type: text/x-python
File size: 6668 byte(s)
*** empty log message ***

1 jgs 97 # $Id$
2 jgs 82 """
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     Version $Id$
9     """
10    
11     import sys
12     import os
13     import unittest
14    
15 jgs 97 from esys.escript import *
16     from esys.linearPDEs import *
17     from esys import finley
18 jgs 82
19     from math import *
20     from numarray import array
21    
22     numElements=10
23     max_error=0.
24     max_error_text=""
25    
26     for dim in [2,3]:
27     for order in [1,2]:
28 jgs 97 for onFace in [True,False]:
29 jgs 82
30 jgs 97 print "\n"
31     print "-----------------------------------------------------------------------------------------------------"
32     print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace)
33     print "-----------------------------------------------------------------------------------------------------"
34 jgs 82
35     if onFace:
36     onFaceText=", on elements"
37     else:
38     onFaceText=""
39 jgs 97
40     case="dim=%d, order=%d%s" % (dim,order,onFaceText)
41 jgs 82
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 jgs 97 test="error gradient in interior (nodes)"
66    
67 jgs 82 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 jgs 97
72     text="%s: %55s = %e" % (case, test, error_norm)
73     print "%s" % text
74    
75 jgs 82 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 jgs 97 test="error gradient in interior (degrees of freedom)"
84    
85 jgs 82 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 jgs 97 error_norm=Lsup(ref-g)/Lsup(ref)
89    
90     text="%s: %55s = %e" % (case, test, error_norm)
91     print "%s" % text
92    
93 jgs 82 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 jgs 97 test="error gradient in interior (reduced degrees of freedom)"
102    
103 jgs 82 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 jgs 97 error_norm=Lsup(ref-g)/Lsup(ref)
107    
108     text="%s: %55s = %e" % (case, test, error_norm)
109     print "%s" % text
110    
111 jgs 82 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 jgs 97 test="error volume integration"
120    
121 jgs 82 x=e.getX()[0:2]
122 jgs 97 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 jgs 82 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 jgs 97 test="error gradient on boundary"
139    
140 jgs 82 x=n.getX()[0:2]
141 jgs 97 g=grad(x**order+x[1]*[1,0],where=f)
142 jgs 82 x=f.getX()[0:2]
143     ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
144 jgs 97 error_norm=Lsup(g-ref)
145    
146     text="%s: %55s = %e" % (case, test, error_norm)
147     print "%s" % text
148    
149 jgs 82 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 jgs 97 test="error gradient on boundary (degrees of freedom)"
158    
159 jgs 82 x=n.getX()[0:2].interpolate(d)
160 jgs 97 g=grad(x**order+x[1]*[1,0],where=f)
161 jgs 82 x=f.getX()[0:2]
162     ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
163 jgs 97 error_norm=Lsup(ref-g)/Lsup(ref)
164    
165     text="%s: %55s = %e" % (case, test, error_norm)
166     print "%s" % text
167    
168 jgs 82 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 jgs 97 test="error gradient on boundary (reduced degrees of freedom)"
177    
178 jgs 82 x=n.getX()[0:2].interpolate(r)
179 jgs 97 g=grad(x+x[1]*[1,0],where=f)
180 jgs 82 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
181 jgs 97 error_norm=Lsup(ref-g)/Lsup(ref)
182    
183     text="%s: %55s = %e" % (case, test, error_norm)
184     print "%s" % text
185    
186 jgs 82 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 jgs 97 test="error boundary integration"
195    
196 jgs 82 x=f.getX()[0:2]
197 jgs 97 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 jgs 82 if error_norm>max_error:
204     max_error_text=text
205     max_error=error_norm
206    
207     #
208     # normal test:
209     #
210    
211 jgs 97 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 jgs 82 if error_norm>max_error:
232     max_error_text=text
233     max_error=error_norm
234    
235 jgs 97 print "-----------------------------------------------------------------------------------------------------"
236 jgs 82
237 jgs 97 print "\n\n"
238     print "******************************************************************************************************************"
239     print "maximal error:", max_error_text
240     print "******************************************************************************************************************"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26