test/python/GradTest.py

# $Id$
"""
Test a grad, interpolate and integrate over the unit square.

The tests are very basic.

by Lutz Gross, ACcESS, University of Queensland, Australia, 2003.
"""

import sys
import os
import unittest
import time

from esys.escript import *
from esys.escript.linearPDEs import *
from esys import finley

from math import *
from numarray import array

starttime = time.clock()

numElements=10
max_error=0.
max_error_text=""
error_tol=pow(10,-10)

for dim in [2,3]:
  for order in [1,2]:
    for onFace in [True,False]:

       print "\n"
       print "-----------------------------------------------------------------------------------------------------"
       print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace)
       print "-----------------------------------------------------------------------------------------------------"

       if onFace:
          onFaceText=", on elements"
       else:
          onFaceText=""

       case="dim=%d, order=%d%s" % (dim,order,onFaceText)
 
       if dim==2:
         mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)
         m00=[[1,0],[0,0]]
         m01=[[0,1],[0,0]]
         m11=[[0,0],[0,1]]
         h=5
       else:
         mydomain=finley.Brick(numElements,numElements,numElements,order=order,useElementsOnFace=onFace)
         m00=[[1,0,0],[0,0,0]]
         m01=[[0,1,0],[0,0,0]]
         m11=[[0,0,0],[0,1,0]]
         h=7

       n=ContinuousFunction(mydomain)
       e=Function(mydomain)
       f=FunctionOnBoundary(mydomain)
       d=Solution(mydomain)
       r=ReducedSolution(mydomain)

       #
       # test gradient 
       #

       test="error gradient in interior (nodes)"

       x=n.getX()[0:2]
       g=grad(x**order+x[1]*[1,0])
       ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
       error_norm=Lsup(ref-g)

       text="%s: %55s = %e" % (case, test, error_norm)
       print "%s" % text

       if error_norm>max_error: 
           max_error_text=text
           max_error=error_norm

       #
       # test gradient on degrees of freedom
       #

       test="error gradient in interior (degrees of freedom)"

       x=interpolate(n.getX()[0:2],d)
       g=grad(x**order+x[1]*[1,0])
       ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
       error_norm=Lsup(ref-g)/Lsup(ref)

       text="%s: %55s = %e" % (case, test, error_norm)
       print "%s" % text

       if error_norm>max_error: 
           max_error_text=text
           max_error=error_norm

       #
       # test gradient on reduced degrees of freedom
       #

       test="error gradient in interior (reduced degrees of freedom)"

       x=interpolate(n.getX()[0:2],r)
       g=grad(x+x[1]*[1,0])
       ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
       error_norm=Lsup(ref-g)/Lsup(ref)

       text="%s: %55s = %e" % (case, test, error_norm)
       print "%s" % text

       if error_norm>max_error: 
           max_error_text=text
           max_error=error_norm

       #
       # test integration over volume
       #

       test="error volume integration"

       x=e.getX()[0:2]
       error=integrate(x**2+[0,2.]*x)-array([1./3.,1./3.+2*1./2.])
       error_norm=sqrt(numarray.innerproduct(error,error))

       text="%s: %55s = %e" % (case, test, error_norm)
       print "%s" % text

       if error_norm>max_error: 
           max_error_text=text
           max_error=error_norm

       if onFace:

          #
          # gradient on the boundary:
          #

          test="error gradient on boundary"

          x=n.getX()[0:2]
          g=grad(x**order+x[1]*[1,0],where=f)
          x=f.getX()[0:2]
          ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
          error_norm=Lsup(g-ref)

          text="%s: %55s = %e" % (case, test, error_norm)
          print "%s" % text

          if error_norm>max_error: 
              max_error_text=text
              max_error=error_norm

          #
          # test gradient on degrees of freedom
          #

          test="error gradient on boundary (degrees of freedom)"

          x=interpolate(n.getX()[0:2],d)
          g=grad(x**order+x[1]*[1,0],where=f)
          x=f.getX()[0:2]
          ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11
          error_norm=Lsup(ref-g)/Lsup(ref)

          text="%s: %55s = %e" % (case, test, error_norm)
          print "%s" % text

          if error_norm>max_error: 
              max_error_text=text
              max_error=error_norm

          #
          # test gradient on reduced degrees of freedom
          #

          test="error gradient on boundary (reduced degrees of freedom)"

          x=interpolate(n.getX()[0:2],r)
          g=grad(x+x[1]*[1,0],where=f)
          ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11
          error_norm=Lsup(ref-g)/Lsup(ref)

          text="%s: %55s = %e" % (case, test, error_norm)
          print "%s" % text

          if error_norm>max_error: 
              max_error_text=text
              max_error=error_norm

       #
       # test integration over boundary
       #

       test="error boundary integration"

       x=f.getX()[0:2]
       error=integrate(x**2+[0,2.]*x)-array([h/3.,h/3.+2*(h-1)/2.])
       error_norm=sqrt(numarray.innerproduct(error,error))

       text="%s: %55s = %e" % (case, test, error_norm)
       print "%s" % text

       if error_norm>max_error: 
           max_error_text=text
           max_error=error_norm

       #
       # normal test:
       #

       test="error normals"

       refNormal=Vector(0,what=f)
       if dim==3:
           refNormal.setTaggedValue(2,[1,0,0])
           refNormal.setTaggedValue(1,[-1,0,0])
           refNormal.setTaggedValue(20,[0,1,0])
           refNormal.setTaggedValue(10,[0,-1,0])
           refNormal.setTaggedValue(200,[0,0,1])
           refNormal.setTaggedValue(100,[0,0,-1])
       else:
           refNormal.setTaggedValue(2,[1,0])
           refNormal.setTaggedValue(1,[-1,0])
           refNormal.setTaggedValue(20,[0,1])
           refNormal.setTaggedValue(10,[0,-1])
       error_norm=Lsup(f.getNormal()-refNormal)

       text="%s: %55s = %e" % (case, test, error_norm)
       print "%s" % text

       if error_norm>max_error: 
           max_error_text=text
           max_error=error_norm

       print "-----------------------------------------------------------------------------------------------------"

print "\n\n"
print "******************************************************************************************************************"
print "maximal error:", max_error_text
print "******************************************************************************************************************"

stoptime = time.clock()
elapsed = stoptime - starttime
print "\nElapsed time: ", elapsed, "\n"

if max_error > error_tol:
  print "max error exceeds tolerance"
  sys.exit(1)

sys.exit(0)
1	jgs	102	# $Id$
2	jgs	82	"""
3			Test a grad, interpolate and integrate over the unit square.
4
5			The tests are very basic.
6
7			by Lutz Gross, ACcESS, University of Queensland, Australia, 2003.
8			"""
9
10			import sys
11			import os
12			import unittest
13	jgs	153	import time
14	jgs	82
15	jgs	149	from esys.escript import *
16			from esys.escript.linearPDEs import *
17			from esys import finley
18	jgs	82
19			from math import *
20			from numarray import array
21
22	jgs	153	starttime = time.clock()
23
24	jgs	82	numElements=10
25			max_error=0.
26			max_error_text=""
27	jgs	147	error_tol=pow(10,-10)
28	jgs	82
29			for dim in [2,3]:
30			for order in [1,2]:
31	jgs	102	for onFace in [True,False]:
32	jgs	82
33	jgs	102	print "\n"
34			print "-----------------------------------------------------------------------------------------------------"
35			print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace)
36			print "-----------------------------------------------------------------------------------------------------"
37	jgs	82
38			if onFace:
39			onFaceText=", on elements"
40			else:
41			onFaceText=""
42	jgs	102
43			case="dim=%d, order=%d%s" % (dim,order,onFaceText)
44	jgs	82
45			if dim==2:
46			mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace)
47			m00=[[1,0],[0,0]]
48			m01=[[0,1],[0,0]]
49			m11=[[0,0],[0,1]]
50			h=5
51			else:
52			mydomain=finley.Brick(numElements,numElements,numElements,order=order,useElementsOnFace=onFace)
53			m00=[[1,0,0],[0,0,0]]
54			m01=[[0,1,0],[0,0,0]]
55			m11=[[0,0,0],[0,1,0]]
56			h=7
57
58			n=ContinuousFunction(mydomain)
59			e=Function(mydomain)
60			f=FunctionOnBoundary(mydomain)
61			d=Solution(mydomain)
62			r=ReducedSolution(mydomain)
63
64			#
65			# test gradient
66			#
67
68	jgs	102	test="error gradient in interior (nodes)"
69
70	jgs	82	x=n.getX()[0:2]
71			g=grad(x*order+x[1][1,0])
72			ref=orderx[0](order-1)m00+m01+orderx[1](order-1)m11
73			error_norm=Lsup(ref-g)
74	jgs	102
75			text="%s: %55s = %e" % (case, test, error_norm)
76			print "%s" % text
77
78	jgs	82	if error_norm>max_error:
79			max_error_text=text
80			max_error=error_norm
81
82			#
83			# test gradient on degrees of freedom
84			#
85
86	jgs	102	test="error gradient in interior (degrees of freedom)"
87
88	gross	428	x=interpolate(n.getX()[0:2],d)
89	jgs	82	g=grad(x*order+x[1][1,0])
90			ref=orderx[0](order-1)m00+m01+orderx[1](order-1)m11
91	jgs	102	error_norm=Lsup(ref-g)/Lsup(ref)
92
93			text="%s: %55s = %e" % (case, test, error_norm)
94			print "%s" % text
95
96	jgs	82	if error_norm>max_error:
97			max_error_text=text
98			max_error=error_norm
99
100			#
101			# test gradient on reduced degrees of freedom
102			#
103
104	jgs	102	test="error gradient in interior (reduced degrees of freedom)"
105
106	gross	428	x=interpolate(n.getX()[0:2],r)
107	jgs	82	g=grad(x+x[1]*[1,0])
108			ref=Scalar(1,what=r)m00+m01+Scalar(1,what=r)m11
109	jgs	102	error_norm=Lsup(ref-g)/Lsup(ref)
110
111			text="%s: %55s = %e" % (case, test, error_norm)
112			print "%s" % text
113
114	jgs	82	if error_norm>max_error:
115			max_error_text=text
116			max_error=error_norm
117
118			#
119			# test integration over volume
120			#
121
122	jgs	102	test="error volume integration"
123
124	jgs	82	x=e.getX()[0:2]
125	jgs	102	error=integrate(x*2+[0,2.]x)-array([1./3.,1./3.+2*1./2.])
126			error_norm=sqrt(numarray.innerproduct(error,error))
127
128			text="%s: %55s = %e" % (case, test, error_norm)
129			print "%s" % text
130
131	jgs	82	if error_norm>max_error:
132			max_error_text=text
133			max_error=error_norm
134
135			if onFace:
136
137			#
138			# gradient on the boundary:
139			#
140
141	jgs	102	test="error gradient on boundary"
142
143	jgs	82	x=n.getX()[0:2]
144	jgs	102	g=grad(x*order+x[1][1,0],where=f)
145	jgs	82	x=f.getX()[0:2]
146			ref=orderx[0](order-1)m00+m01+orderx[1](order-1)m11
147	jgs	102	error_norm=Lsup(g-ref)
148
149			text="%s: %55s = %e" % (case, test, error_norm)
150			print "%s" % text
151
152	jgs	82	if error_norm>max_error:
153			max_error_text=text
154			max_error=error_norm
155
156			#
157			# test gradient on degrees of freedom
158			#
159
160	jgs	102	test="error gradient on boundary (degrees of freedom)"
161
162	gross	428	x=interpolate(n.getX()[0:2],d)
163	jgs	102	g=grad(x*order+x[1][1,0],where=f)
164	jgs	82	x=f.getX()[0:2]
165			ref=orderx[0](order-1)m00+m01+orderx[1](order-1)m11
166	jgs	102	error_norm=Lsup(ref-g)/Lsup(ref)
167
168			text="%s: %55s = %e" % (case, test, error_norm)
169			print "%s" % text
170
171	jgs	82	if error_norm>max_error:
172			max_error_text=text
173			max_error=error_norm
174
175			#
176			# test gradient on reduced degrees of freedom
177			#
178
179	jgs	102	test="error gradient on boundary (reduced degrees of freedom)"
180
181	gross	428	x=interpolate(n.getX()[0:2],r)
182	jgs	102	g=grad(x+x[1]*[1,0],where=f)
183	jgs	82	ref=Scalar(1,what=r)m00+m01+Scalar(1,what=r)m11
184	jgs	102	error_norm=Lsup(ref-g)/Lsup(ref)
185
186			text="%s: %55s = %e" % (case, test, error_norm)
187			print "%s" % text
188
189	jgs	82	if error_norm>max_error:
190			max_error_text=text
191			max_error=error_norm
192
193			#
194			# test integration over boundary
195			#
196
197	jgs	102	test="error boundary integration"
198
199	jgs	82	x=f.getX()[0:2]
200	jgs	102	error=integrate(x*2+[0,2.]x)-array([h/3.,h/3.+2*(h-1)/2.])
201			error_norm=sqrt(numarray.innerproduct(error,error))
202
203			text="%s: %55s = %e" % (case, test, error_norm)
204			print "%s" % text
205
206	jgs	82	if error_norm>max_error:
207			max_error_text=text
208			max_error=error_norm
209
210			#
211			# normal test:
212			#
213
214	jgs	102	test="error normals"
215
216			refNormal=Vector(0,what=f)
217			if dim==3:
218			refNormal.setTaggedValue(2,[1,0,0])
219			refNormal.setTaggedValue(1,[-1,0,0])
220			refNormal.setTaggedValue(20,[0,1,0])
221			refNormal.setTaggedValue(10,[0,-1,0])
222			refNormal.setTaggedValue(200,[0,0,1])
223			refNormal.setTaggedValue(100,[0,0,-1])
224			else:
225			refNormal.setTaggedValue(2,[1,0])
226			refNormal.setTaggedValue(1,[-1,0])
227			refNormal.setTaggedValue(20,[0,1])
228			refNormal.setTaggedValue(10,[0,-1])
229			error_norm=Lsup(f.getNormal()-refNormal)
230
231			text="%s: %55s = %e" % (case, test, error_norm)
232			print "%s" % text
233
234	jgs	82	if error_norm>max_error:
235			max_error_text=text
236			max_error=error_norm
237
238	jgs	102	print "-----------------------------------------------------------------------------------------------------"
239	jgs	82
240	jgs	102	print "\n\n"
241			print "******************************************************************************************************************"
242			print "maximal error:", max_error_text
243			print "******************************************************************************************************************"
244	jgs	147
245	jgs	153	stoptime = time.clock()
246			elapsed = stoptime - starttime
247			print "\nElapsed time: ", elapsed, "\n"
248
249	jgs	147	if max_error > error_tol:
250			print "max error exceeds tolerance"
251			sys.exit(1)
252
253			sys.exit(0)
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision