revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC
# Line 1  Line 1
1    # \$Id\$
2  """  """
3  Test a grad, interpolate and integrate over the unit square.  Test a grad, interpolate and integrate over the unit square.
4
# Line 11  import sys Line 12  import sys
12  import os  import os
13  import unittest  import unittest
14
15    from esys.escript import *
16  esys_root=os.getenv('ESYS_ROOT')  from esys.linearPDEs import *
17  sys.path.append(esys_root+'/finley/lib')  from esys import finley
sys.path.append(esys_root+'/escript/lib')
sys.path.append(esys_root+'/escript/py_src')

from escript import *
from util import *
import finley
18
19  from math import *  from math import *
20  from numarray import array  from numarray import array
21

22  numElements=10  numElements=10

23  max_error=0.  max_error=0.
24  max_error_text=""  max_error_text=""
25
26  for dim in [2,3]:  for dim in [2,3]:
27    for order in [1,2]:    for order in [1,2]:
28      for onFace in [TRUE,FALSE]:      for onFace in [True,False]:
29
30         print "\ndim: %d order: %i onFace: %s" % (dim, order, onFace)         print "\n"
31         print "-------------------------------------------------------------------------------------"         print "-----------------------------------------------------------------------------------------------------"
32           print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace)
33           print "-----------------------------------------------------------------------------------------------------"
34
35         if onFace:         if onFace:
36            onFaceText=", on elements"            onFaceText=", on elements"
37         else:         else:
38            onFaceText=""            onFaceText=""
39         case="dim=%d, order=%d%s"%(dim,order,onFaceText)
40           case="dim=%d, order=%d%s" % (dim,order,onFaceText)
41
42         if dim==2:         if dim==2:
43           mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)           mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)
# Line 66  for dim in [2,3]: Line 62  for dim in [2,3]:
63         #         #
64
65           test="error gradient in interior (nodes)"
66
67         x=n.getX()[0:2]         x=n.getX()[0:2]
69         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
70         error_norm=Lsup(ref-g)         error_norm=Lsup(ref-g)
71         text=" %s: error gradient in interior (nodes) = %e"%(case,error_norm)
72         print "\t%s"%text         text="%s: %55s = %e" % (case, test, error_norm)
73           print "%s" % text
74
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
83           test="error gradient in interior (degrees of freedom)"
84
85         x=n.getX()[0:2].interpolate(d)         x=n.getX()[0:2].interpolate(d)
87         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
88         error_norm=Lsup(ref-g)         error_norm=Lsup(ref-g)/Lsup(ref)
89         text=" %s: error gradient in interior (degrees of freedom) = %e"%(case,error_norm)
90         print "\t%s"%text         text="%s: %55s = %e" % (case, test, error_norm)
91           print "%s" % text
92
93         if error_norm>max_error:         if error_norm>max_error:
94             max_error_text=text             max_error_text=text
95             max_error=error_norm             max_error=error_norm
# Line 94  for dim in [2,3]: Line 98  for dim in [2,3]:
98         # test gradient on reduced degrees of freedom         # test gradient on reduced degrees of freedom
99         #         #
100
101           test="error gradient in interior (reduced degrees of freedom)"
102
103         x=n.getX()[0:2].interpolate(r)         x=n.getX()[0:2].interpolate(r)
105         ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11         ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
106         error_norm=Lsup(ref-g)         error_norm=Lsup(ref-g)/Lsup(ref)
107         text=" %s: error gradient in interior (reduced degrees of freedom) = %e"%(case,error_norm)
108         print "\t%s"%text         text="%s: %55s = %e" % (case, test, error_norm)
109           print "%s" % text
110
111         if error_norm>max_error:         if error_norm>max_error:
112             max_error_text=text             max_error_text=text
113             max_error=error_norm             max_error=error_norm
# Line 107  for dim in [2,3]: Line 115  for dim in [2,3]:
115         #         #
116         # test integration over volume         # test integration over volume
117         #         #
118         # integrate causes: Fatal Python error: Call to numarray API function
119         # without first calling import_libnumarray() in src/Data/Data.cpp         test="error volume integration"
120
121         x=e.getX()[0:2]         x=e.getX()[0:2]
122         #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.])
123         #error_norm=sqrt(numarray.innerproduct(error,error))         error_norm=sqrt(numarray.innerproduct(error,error))
124         text=" %s: error volume integration= %e"%(case,error_norm)
125         print "\t%s"%text         text="%s: %55s = %e" % (case, test, error_norm)
126           print "%s" % text
127
128         if error_norm>max_error:         if error_norm>max_error:
129             max_error_text=text             max_error_text=text
130             max_error=error_norm             max_error=error_norm
# Line 124  for dim in [2,3]: Line 134  for dim in [2,3]:
134            #            #
136            #            #
137            # Lsup fails - perhaps grad(what=f) is needed?
139
140            x=n.getX()[0:2]            x=n.getX()[0:2]
142            x=f.getX()[0:2]            x=f.getX()[0:2]
143            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
144            #error_norm=Lsup(g-ref)            error_norm=Lsup(g-ref)
145            error_norm=0
146            text=" %s: error gradient on boundary = %e"%(case,error_norm)            text="%s: %55s = %e" % (case, test, error_norm)
147            print "\t%s"%text            print "%s" % text
148
149            if error_norm>max_error:            if error_norm>max_error:
150                max_error_text=text                max_error_text=text
151                max_error=error_norm                max_error=error_norm
# Line 142  for dim in [2,3]: Line 153  for dim in [2,3]:
153            #            #
154            # test gradient on degrees of freedom            # test gradient on degrees of freedom
155            #            #
156            # Lsup fails - perhaps grad(what=f) is needed?
157              test="error gradient on boundary (degrees of freedom)"
158
159            x=n.getX()[0:2].interpolate(d)            x=n.getX()[0:2].interpolate(d)
161            x=f.getX()[0:2]            x=f.getX()[0:2]
162            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
163            #error_norm=Lsup(ref-g)            error_norm=Lsup(ref-g)/Lsup(ref)
164            error_norm=0
165            text=" %s: error gradient on boundary (degrees of freedom) = %e"%(case,error_norm)            text="%s: %55s = %e" % (case, test, error_norm)
166            print "\t%s"%text            print "%s" % text
167
168            if error_norm>max_error:            if error_norm>max_error:
169                max_error_text=text                max_error_text=text
170                max_error=error_norm                max_error=error_norm
# Line 160  for dim in [2,3]: Line 172  for dim in [2,3]:
172            #            #
173            # test gradient on reduced degrees of freedom            # test gradient on reduced degrees of freedom
174            #            #
175            # Lsup fails - perhaps grad(what=f) is needed?
176              test="error gradient on boundary (reduced degrees of freedom)"
177
178            x=n.getX()[0:2].interpolate(r)            x=n.getX()[0:2].interpolate(r)
180            ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11            ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
181            #error_norm=Lsup(ref-g)            error_norm=Lsup(ref-g)/Lsup(ref)
182            error_norm=0
183            text=" %s: error gradient on boundary (reduced degrees of freedom) = %e"%(case,error_norm)            text="%s: %55s = %e" % (case, test, error_norm)
184            print "\t%s"%text            print "%s" % text
185
186            if error_norm>max_error:            if error_norm>max_error:
187                max_error_text=text                max_error_text=text
188                max_error=error_norm                max_error=error_norm
# Line 177  for dim in [2,3]: Line 190  for dim in [2,3]:
190         #         #
191         # test integration over boundary         # test integration over boundary
192         #         #
193         # integrate causes: Fatal Python error: Call to numarray API function
194         # without first calling import_libnumarray() in src/Data/Data.cpp         test="error boundary integration"
195
196         x=f.getX()[0:2]         x=f.getX()[0:2]
197         #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.])
198         #error_norm=sqrt(numarray.innerproduct(error,error))         error_norm=sqrt(numarray.innerproduct(error,error))
199         text=" %s: error boundary integration= %e"%(case,error_norm)
200         print "\t%s"%text         text="%s: %55s = %e" % (case, test, error_norm)
201           print "%s" % text
202
203         if error_norm>max_error:         if error_norm>max_error:
204             max_error_text=text             max_error_text=text
205             max_error=error_norm             max_error=error_norm
# Line 192  for dim in [2,3]: Line 207  for dim in [2,3]:
207         #         #
208         # normal test:         # normal test:
209         #         #
#
210
211         #refNormal=Vector(0,what=f)         test="error normals"
212         #if dim==3:
219         #else:             refNormal.setTaggedValue(200,[0,0,1])
224         #error_norm=Lsup(f.getNormal()-refNormal)             refNormal.setTaggedValue(20,[0,1])
225         text=" %s: error normals= %e"%(case,error_norm)             refNormal.setTaggedValue(10,[0,-1])
226         print "\t%s"%text         error_norm=Lsup(f.getNormal()-refNormal)
227
228           text="%s: %55s = %e" % (case, test, error_norm)
229           print "%s" % text
230
231         if error_norm>max_error:         if error_norm>max_error:
232             max_error_text=text             max_error_text=text
233             max_error=error_norm             max_error=error_norm
234
235         print "-------------------------------------------------------------------------------------"         print "-----------------------------------------------------------------------------------------------------"
236
237  print "\n\nmaximal error for",max_error_text  print "\n\n"
238    print "******************************************************************************************************************"
239    print "maximal error:", max_error_text
240    print "******************************************************************************************************************"

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