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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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]:
62         # test gradient         # test gradient
63         #         #
64    
65           test="error gradient in interior (nodes)"
66    
67         x=n.getX()[0:2]         x=n.getX()[0:2]
68         g=grad(x**order+x[1]*[1,0])         g=grad(x**order+x[1]*[1,0])
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)
86         g=grad(x**order+x[1]*[1,0])         g=grad(x**order+x[1]*[1,0])
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)
104         g=grad(x+x[1]*[1,0])         g=grad(x+x[1]*[1,0])
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            #            #
135            # gradient on the boundary:            # gradient on the boundary:
136            #            #
137            # Lsup fails - perhaps grad(what=f) is needed?  
138              test="error gradient on boundary"
139    
140            x=n.getX()[0:2]            x=n.getX()[0:2]
141            #g=grad(x**order+x[1]*[1,0],what=f)            g=grad(x**order+x[1]*[1,0],where=f)
           g=grad(x**order+x[1]*[1,0])  
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)
160            #g=grad(x**order+x[1]*[1,0],what=f)            g=grad(x**order+x[1]*[1,0],where=f)
           g=grad(x**order+x[1]*[1,0])  
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)
179            #g=grad(x+x[1]*[1,0],what=f)            g=grad(x+x[1]*[1,0],where=f)
           g=grad(x+x[1]*[1,0])  
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         #         #
        # need to add wrapper for DataTagged::addTaggedValue  
        #  
210    
211         #refNormal=Vector(0,what=f)         test="error normals"
212         #if dim==3:  
213         #    refNormal.addTaggedValue(2,[1,0,0])         refNormal=Vector(0,what=f)
214         #    refNormal.addTaggedValue(1,[-1,0,0])         if dim==3:
215         #    refNormal.addTaggedValue(20,[0,1,0])             refNormal.setTaggedValue(2,[1,0,0])
216         #    refNormal.addTaggedValue(10,[0,-1,0])             refNormal.setTaggedValue(1,[-1,0,0])
217         #    refNormal.addTaggedValue(200,[0,0,1])             refNormal.setTaggedValue(20,[0,1,0])
218         #    refNormal.addTaggedValue(100,[0,0,-1])             refNormal.setTaggedValue(10,[0,-1,0])
219         #else:             refNormal.setTaggedValue(200,[0,0,1])
220         #    refNormal.addTaggedValue(2,[1,0])             refNormal.setTaggedValue(100,[0,0,-1])
221         #    refNormal.addTaggedValue(1,[-1,0])         else:
222         #    refNormal.addTaggedValue(20,[0,1])             refNormal.setTaggedValue(2,[1,0])
223         #    refNormal.addTaggedValue(10,[0,-1])             refNormal.setTaggedValue(1,[-1,0])
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

  ViewVC Help
Powered by ViewVC 1.1.26