revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC
# Line 1  Line 1
# \$Id\$
1  """  """
2  Test a grad, interpolate and integrate over the unit square.  Test a grad, interpolate and integrate over the unit square.
3
# Line 12  import sys Line 11  import sys
11  import os  import os
12  import unittest  import unittest
13
14  from esys.escript import *
15  from esys.linearPDEs import *  esys_root=os.getenv('ESYS_ROOT')
16  from esys import finley  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 *  from math import *
25  from numarray import array  from numarray import array
26
27
28  numElements=10  numElements=10
29
30  max_error=0.  max_error=0.
31  max_error_text=""  max_error_text=""
32
33  for dim in [2,3]:  for dim in [2,3]:
34    for order in [1,2]:    for order in [1,2]:
35      for onFace in [True,False]:      for onFace in [TRUE,FALSE]:
36
37         print "\n"         print "\ndim: %d order: %i onFace: %s" % (dim, order, onFace)
38         print "-----------------------------------------------------------------------------------------------------"         print "-------------------------------------------------------------------------------------"
print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace)
print "-----------------------------------------------------------------------------------------------------"
39
40         if onFace:         if onFace:
41            onFaceText=", on elements"            onFaceText=", on elements"
42         else:         else:
43            onFaceText=""            onFaceText=""
44           case="dim=%d, order=%d%s"%(dim,order,onFaceText)
case="dim=%d, order=%d%s" % (dim,order,onFaceText)
45
46         if dim==2:         if dim==2:
47           mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)           mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)
# Line 62  for dim in [2,3]: Line 66  for dim in [2,3]:
67         #         #
68

69         x=n.getX()[0:2]         x=n.getX()[0:2]
71         ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11         ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
72         error_norm=Lsup(ref-g)         error_norm=Lsup(ref-g)
73           text=" %s: error gradient in interior (nodes) = %e"%(case,error_norm)
74         text="%s: %55s = %e" % (case, test, error_norm)         print "\t%s"%text
print "%s" % text

75         if error_norm>max_error:         if error_norm>max_error:
76             max_error_text=text             max_error_text=text
77             max_error=error_norm             max_error=error_norm
# Line 80  for dim in [2,3]: Line 80  for dim in [2,3]:
80         # test gradient on degrees of freedom         # test gradient on degrees of freedom
81         #         #
82
test="error gradient in interior (degrees of freedom)"

83         x=n.getX()[0:2].interpolate(d)         x=n.getX()[0:2].interpolate(d)
85         ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11         ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
86         error_norm=Lsup(ref-g)/Lsup(ref)         error_norm=Lsup(ref-g)
87           text=" %s: error gradient in interior (degrees of freedom) = %e"%(case,error_norm)
88         text="%s: %55s = %e" % (case, test, error_norm)         print "\t%s"%text
print "%s" % text

89         if error_norm>max_error:         if error_norm>max_error:
90             max_error_text=text             max_error_text=text
91             max_error=error_norm             max_error=error_norm
# Line 98  for dim in [2,3]: Line 94  for dim in [2,3]:
94         # test gradient on reduced degrees of freedom         # test gradient on reduced degrees of freedom
95         #         #
96
test="error gradient in interior (reduced degrees of freedom)"

97         x=n.getX()[0:2].interpolate(r)         x=n.getX()[0:2].interpolate(r)
99         ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11         ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
100         error_norm=Lsup(ref-g)/Lsup(ref)         error_norm=Lsup(ref-g)
101           text=" %s: error gradient in interior (reduced degrees of freedom) = %e"%(case,error_norm)
102         text="%s: %55s = %e" % (case, test, error_norm)         print "\t%s"%text
print "%s" % text

103         if error_norm>max_error:         if error_norm>max_error:
104             max_error_text=text             max_error_text=text
105             max_error=error_norm             max_error=error_norm
# Line 115  for dim in [2,3]: Line 107  for dim in [2,3]:
107         #         #
108         # test integration over volume         # test integration over volume
109         #         #
110           # integrate causes: Fatal Python error: Call to numarray API function
111         test="error volume integration"         # without first calling import_libnumarray() in src/Data/Data.cpp
112
113         x=e.getX()[0:2]         x=e.getX()[0:2]
114         error=integrate(x**2+[0,2.]*x)-array([1./3.,1./3.+2*1./2.])         #error=integrate(x**2+[0,2.]*x)-array([1./3.,1./3.+2*1./2.])
115         error_norm=sqrt(numarray.innerproduct(error,error))         #error_norm=sqrt(numarray.innerproduct(error,error))
116           text=" %s: error volume integration= %e"%(case,error_norm)
117         text="%s: %55s = %e" % (case, test, error_norm)         print "\t%s"%text
print "%s" % text

118         if error_norm>max_error:         if error_norm>max_error:
119             max_error_text=text             max_error_text=text
120             max_error=error_norm             max_error=error_norm
# Line 134  for dim in [2,3]: Line 124  for dim in [2,3]:
124            #            #
126            #            #
127              # Lsup fails - perhaps grad(what=f) is needed?
128
129            x=n.getX()[0:2]            x=n.getX()[0:2]
132            x=f.getX()[0:2]            x=f.getX()[0:2]
133            ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11            ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
134            error_norm=Lsup(g-ref)            #error_norm=Lsup(g-ref)
135              error_norm=0
136            text="%s: %55s = %e" % (case, test, error_norm)            text=" %s: error gradient on boundary = %e"%(case,error_norm)
137            print "%s" % text            print "\t%s"%text

138            if error_norm>max_error:            if error_norm>max_error:
139                max_error_text=text                max_error_text=text
140                max_error=error_norm                max_error=error_norm
# Line 153  for dim in [2,3]: Line 142  for dim in [2,3]:
142            #            #
143            # test gradient on degrees of freedom            # test gradient on degrees of freedom
144            #            #
145              # Lsup fails - perhaps grad(what=f) is needed?
test="error gradient on boundary (degrees of freedom)"
146
147            x=n.getX()[0:2].interpolate(d)            x=n.getX()[0:2].interpolate(d)
150            x=f.getX()[0:2]            x=f.getX()[0:2]
151            ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11            ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
152            error_norm=Lsup(ref-g)/Lsup(ref)            #error_norm=Lsup(ref-g)
153              error_norm=0
154            text="%s: %55s = %e" % (case, test, error_norm)            text=" %s: error gradient on boundary (degrees of freedom) = %e"%(case,error_norm)
155            print "%s" % text            print "\t%s"%text

156            if error_norm>max_error:            if error_norm>max_error:
157                max_error_text=text                max_error_text=text
158                max_error=error_norm                max_error=error_norm
# Line 172  for dim in [2,3]: Line 160  for dim in [2,3]:
160            #            #
161            # test gradient on reduced degrees of freedom            # test gradient on reduced degrees of freedom
162            #            #
163              # Lsup fails - perhaps grad(what=f) is needed?
test="error gradient on boundary (reduced degrees of freedom)"
164
165            x=n.getX()[0:2].interpolate(r)            x=n.getX()[0:2].interpolate(r)
168            ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11            ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
169            error_norm=Lsup(ref-g)/Lsup(ref)            #error_norm=Lsup(ref-g)
170              error_norm=0
171            text="%s: %55s = %e" % (case, test, error_norm)            text=" %s: error gradient on boundary (reduced degrees of freedom) = %e"%(case,error_norm)
172            print "%s" % text            print "\t%s"%text

173            if error_norm>max_error:            if error_norm>max_error:
174                max_error_text=text                max_error_text=text
175                max_error=error_norm                max_error=error_norm
# Line 190  for dim in [2,3]: Line 177  for dim in [2,3]:
177         #         #
178         # test integration over boundary         # test integration over boundary
179         #         #
180           # integrate causes: Fatal Python error: Call to numarray API function
181         test="error boundary integration"         # without first calling import_libnumarray() in src/Data/Data.cpp
182
183         x=f.getX()[0:2]         x=f.getX()[0:2]
184         error=integrate(x**2+[0,2.]*x)-array([h/3.,h/3.+2*(h-1)/2.])         #error=integrate(x**2+[0,2.]*x)-array([h/3.,h/3.+2*(h-1)/2.])
185         error_norm=sqrt(numarray.innerproduct(error,error))         #error_norm=sqrt(numarray.innerproduct(error,error))
186           text=" %s: error boundary integration= %e"%(case,error_norm)
187         text="%s: %55s = %e" % (case, test, error_norm)         print "\t%s"%text
print "%s" % text

188         if error_norm>max_error:         if error_norm>max_error:
189             max_error_text=text             max_error_text=text
190             max_error=error_norm             max_error=error_norm
# Line 207  for dim in [2,3]: Line 192  for dim in [2,3]:
192         #         #
193         # normal test:         # normal test:
194         #         #
196           #
197
198         test="error normals"         #refNormal=Vector(0,what=f)
199           #if dim==3:
206             refNormal.setTaggedValue(200,[0,0,1])         #else:
211             refNormal.setTaggedValue(20,[0,1])         #error_norm=Lsup(f.getNormal()-refNormal)
212             refNormal.setTaggedValue(10,[0,-1])         text=" %s: error normals= %e"%(case,error_norm)
213         error_norm=Lsup(f.getNormal()-refNormal)         print "\t%s"%text

text="%s: %55s = %e" % (case, test, error_norm)
print "%s" % text

214         if error_norm>max_error:         if error_norm>max_error:
215             max_error_text=text             max_error_text=text
216             max_error=error_norm             max_error=error_norm
217
218         print "-----------------------------------------------------------------------------------------------------"         print "-------------------------------------------------------------------------------------"
219
220  print "\n\n"  print "\n\nmaximal error for",max_error_text
print "******************************************************************************************************************"
print "maximal error:", max_error_text
print "******************************************************************************************************************"

Legend:
 Removed from v.97 changed lines Added in v.100