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

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

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

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]
67         g=grad(x**order+x[1]*[1,0])         g=grad(x**order+x[1]*[1,0])
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)
85         g=grad(x**order+x[1]*[1,0])         g=grad(x**order+x[1]*[1,0])
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)
103         g=grad(x+x[1]*[1,0])         g=grad(x+x[1]*[1,0])
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]
140            #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])  
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)
159            #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])  
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)
178            #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])  
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

  ViewVC Help
Powered by ViewVC 1.1.26