/[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 102 - (hide annotations)
Wed Dec 15 07:08:39 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: 6655 byte(s)
*** empty log message ***

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26