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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years, 7 months ago) by jgs
File MIME type: text/x-python
File size: 6791 byte(s)
Initial revision

1 jgs 82 """
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