/[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 153 - (hide annotations)
Tue Oct 25 01:51:20 2005 UTC (14 years, 5 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/GradTest.py
File MIME type: text/x-python
File size: 6911 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-10-25

1 jgs 102 # $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     """
9    
10     import sys
11     import os
12     import unittest
13 jgs 153 import time
14 jgs 82
15 jgs 149 from esys.escript import *
16     from esys.escript.linearPDEs import *
17     from esys import finley
18 jgs 82
19     from math import *
20     from numarray import array
21    
22 jgs 153 starttime = time.clock()
23    
24 jgs 82 numElements=10
25     max_error=0.
26     max_error_text=""
27 jgs 147 error_tol=pow(10,-10)
28 jgs 82
29     for dim in [2,3]:
30     for order in [1,2]:
31 jgs 102 for onFace in [True,False]:
32 jgs 82
33 jgs 102 print "\n"
34     print "-----------------------------------------------------------------------------------------------------"
35     print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace)
36     print "-----------------------------------------------------------------------------------------------------"
37 jgs 82
38     if onFace:
39     onFaceText=", on elements"
40     else:
41     onFaceText=""
42 jgs 102
43     case="dim=%d, order=%d%s" % (dim,order,onFaceText)
44 jgs 82
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 jgs 102 test="error gradient in interior (nodes)"
69    
70 jgs 82 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 jgs 102
75     text="%s: %55s = %e" % (case, test, error_norm)
76     print "%s" % text
77    
78 jgs 82 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 jgs 102 test="error gradient in interior (degrees of freedom)"
87    
88 jgs 82 x=n.getX()[0:2].interpolate(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 jgs 102 error_norm=Lsup(ref-g)/Lsup(ref)
92    
93     text="%s: %55s = %e" % (case, test, error_norm)
94     print "%s" % text
95    
96 jgs 82 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 jgs 102 test="error gradient in interior (reduced degrees of freedom)"
105    
106 jgs 82 x=n.getX()[0:2].interpolate(r)
107     g=grad(x+x[1]*[1,0])
108     ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
109 jgs 102 error_norm=Lsup(ref-g)/Lsup(ref)
110    
111     text="%s: %55s = %e" % (case, test, error_norm)
112     print "%s" % text
113    
114 jgs 82 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 jgs 102 test="error volume integration"
123    
124 jgs 82 x=e.getX()[0:2]
125 jgs 102 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 jgs 82 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 jgs 102 test="error gradient on boundary"
142    
143 jgs 82 x=n.getX()[0:2]
144 jgs 102 g=grad(x**order+x[1]*[1,0],where=f)
145 jgs 82 x=f.getX()[0:2]
146     ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
147 jgs 102 error_norm=Lsup(g-ref)
148    
149     text="%s: %55s = %e" % (case, test, error_norm)
150     print "%s" % text
151    
152 jgs 82 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 jgs 102 test="error gradient on boundary (degrees of freedom)"
161    
162 jgs 82 x=n.getX()[0:2].interpolate(d)
163 jgs 102 g=grad(x**order+x[1]*[1,0],where=f)
164 jgs 82 x=f.getX()[0:2]
165     ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
166 jgs 102 error_norm=Lsup(ref-g)/Lsup(ref)
167    
168     text="%s: %55s = %e" % (case, test, error_norm)
169     print "%s" % text
170    
171 jgs 82 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 jgs 102 test="error gradient on boundary (reduced degrees of freedom)"
180    
181 jgs 82 x=n.getX()[0:2].interpolate(r)
182 jgs 102 g=grad(x+x[1]*[1,0],where=f)
183 jgs 82 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
184 jgs 102 error_norm=Lsup(ref-g)/Lsup(ref)
185    
186     text="%s: %55s = %e" % (case, test, error_norm)
187     print "%s" % text
188    
189 jgs 82 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 jgs 102 test="error boundary integration"
198    
199 jgs 82 x=f.getX()[0:2]
200 jgs 102 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 jgs 82 if error_norm>max_error:
207     max_error_text=text
208     max_error=error_norm
209    
210     #
211     # normal test:
212     #
213    
214 jgs 102 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 jgs 82 if error_norm>max_error:
235     max_error_text=text
236     max_error=error_norm
237    
238 jgs 102 print "-----------------------------------------------------------------------------------------------------"
239 jgs 82
240 jgs 102 print "\n\n"
241     print "******************************************************************************************************************"
242     print "maximal error:", max_error_text
243     print "******************************************************************************************************************"
244 jgs 147
245 jgs 153 stoptime = time.clock()
246     elapsed = stoptime - starttime
247     print "\nElapsed time: ", elapsed, "\n"
248    
249 jgs 147 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