# Diff of /trunk/esys2/finley/test/python/GradTest.py

revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC revision 102 by jgs, Wed Dec 15 07:08:39 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
5  The tests are very basic.  The tests are very basic.
6
7  by Lutz Gross, ACcESS, University of Queensland, Australia, 2003.  by Lutz Gross, ACcESS, University of Queensland, Australia, 2003.
Version \$Id\$
8  """  """
9
10  import sys  import sys
11  import os  import os
12  import unittest  import unittest
13
14    from esys.escript import *
15  esys_root=os.getenv('ESYS_ROOT')  from esys.linearPDEs import *
16  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
17
18  from math import *  from math import *
19  from numarray import array  from numarray import array
20

21  numElements=10  numElements=10

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

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