/[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 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]:
66         # test gradient         # test gradient
67         #         #
68    
        test="error gradient in interior (nodes)"  
   
69         x=n.getX()[0:2]         x=n.getX()[0:2]
70         g=grad(x**order+x[1]*[1,0])         g=grad(x**order+x[1]*[1,0])
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)
84         g=grad(x**order+x[1]*[1,0])         g=grad(x**order+x[1]*[1,0])
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)
98         g=grad(x+x[1]*[1,0])         g=grad(x+x[1]*[1,0])
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            #            #
125            # gradient on the boundary:            # gradient on the boundary:
126            #            #
127              # Lsup fails - perhaps grad(what=f) is needed?
           test="error gradient on boundary"  
128    
129            x=n.getX()[0:2]            x=n.getX()[0:2]
130            g=grad(x**order+x[1]*[1,0],where=f)            #g=grad(x**order+x[1]*[1,0],what=f)
131              g=grad(x**order+x[1]*[1,0])
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)
148            g=grad(x**order+x[1]*[1,0],where=f)            #g=grad(x**order+x[1]*[1,0],what=f)
149              g=grad(x**order+x[1]*[1,0])
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)
166            g=grad(x+x[1]*[1,0],where=f)            #g=grad(x+x[1]*[1,0],what=f)
167              g=grad(x+x[1]*[1,0])
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         #         #
195           # need to add wrapper for DataTagged::addTaggedValue
196           #
197    
198         test="error normals"         #refNormal=Vector(0,what=f)
199           #if dim==3:
200         refNormal=Vector(0,what=f)         #    refNormal.addTaggedValue(2,[1,0,0])
201         if dim==3:         #    refNormal.addTaggedValue(1,[-1,0,0])
202             refNormal.setTaggedValue(2,[1,0,0])         #    refNormal.addTaggedValue(20,[0,1,0])
203             refNormal.setTaggedValue(1,[-1,0,0])         #    refNormal.addTaggedValue(10,[0,-1,0])
204             refNormal.setTaggedValue(20,[0,1,0])         #    refNormal.addTaggedValue(200,[0,0,1])
205             refNormal.setTaggedValue(10,[0,-1,0])         #    refNormal.addTaggedValue(100,[0,0,-1])
206             refNormal.setTaggedValue(200,[0,0,1])         #else:
207             refNormal.setTaggedValue(100,[0,0,-1])         #    refNormal.addTaggedValue(2,[1,0])
208         else:         #    refNormal.addTaggedValue(1,[-1,0])
209             refNormal.setTaggedValue(2,[1,0])         #    refNormal.addTaggedValue(20,[0,1])
210             refNormal.setTaggedValue(1,[-1,0])         #    refNormal.addTaggedValue(10,[0,-1])
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

  ViewVC Help
Powered by ViewVC 1.1.26