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

Contents of /trunk/finley/test/python/run_generators.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1341 - (show annotations)
Thu Nov 8 23:55:20 2007 UTC (11 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 18414 byte(s)
tests for numerical integration added.
1 #
2 # $Id$
3 #
4 #######################################################
5 #
6 # Copyright 2003-2007 by ACceSS MNRF
7 # Copyright 2007 by University of Queensland
8 #
9 # http://esscc.uq.edu.au
10 # Primary Business: Queensland, Australia
11 # Licensed under the Open Software License version 3.0
12 # http://www.opensource.org/licenses/osl-3.0.php
13 #
14 #######################################################
15 #
16
17 """
18 checks the mesh generators against the reference meshes in test_meshes and test the finley integration schemes.
19 """
20
21 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
22 http://www.access.edu.au
23 Primary Business: Queensland, Australia"""
24 __license__="""Licensed under the Open Software License version 3.0
25 http://www.opensource.org/licenses/osl-3.0.php"""
26
27 import sys
28 import os
29 import unittest
30 from esys.escript import *
31 from esys.finley import Rectangle,Brick,JoinFaces, ReadGmsh, ReadMesh
32
33 try:
34 FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
35 except KeyError:
36 FINLEY_TEST_DATA='.'
37
38 try:
39 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
40 except KeyError:
41 FINLEY_WORKDIR='.'
42
43 FINLEY_TEST_MESH_PATH=FINLEY_TEST_DATA+"/data_meshes/"
44 if os.name == "nt":
45 FINLEY_TEST_MESH_PATH = FINLEY_TEST_MESH_PATH+"win32/"
46 FINLEY_WORKDIR_PATH=FINLEY_WORKDIR+"/"
47
48 TEST_FILE_EXT=".test"
49 class Test_Generators(unittest.TestCase):
50
51 def checker(self,dom,reference):
52 dom_file=FINLEY_WORKDIR_PATH+TEST_FILE_EXT
53 dom.write(dom_file)
54 # Uncomment this section to dump the files for regression testing
55 # if True:
56 # dom.write(FINLEY_TEST_MESH_PATH+reference)
57 dom_string=open(dom_file).read().splitlines()
58 ref_string=open(FINLEY_TEST_MESH_PATH+reference).read().splitlines()
59 self.failUnlessEqual(len(dom_string),len(ref_string),"number of lines in mesh files does not match reference")
60 for l in range(1,len(ref_string)):
61 self.failUnlessEqual(dom_string[l].strip(),ref_string[l].strip(),"line %d (%s) in mesh files does not match reference (%s)"%(l,ref_string[l].strip(),dom_string[l].strip()))
62
63 def test_hex_2D_order1(self):
64 file="hex_2D_order1.msh"
65 my_dom=Rectangle(1,1,1)
66 self.checker(my_dom,file)
67
68 def test_hex_2D_order1_onFace(self):
69 file="hex_2D_order1_onFace.msh"
70 my_dom=Rectangle(1,1,1,useElementsOnFace=1)
71 self.checker(my_dom,file)
72
73 def test_hex_2D_order2(self):
74 file="hex_2D_order2.msh"
75 my_dom=Rectangle(1,1,2)
76 self.checker(my_dom,file)
77
78 def test_hex_2D_order2_onFace(self):
79 file="hex_2D_order2_onFace.msh"
80 my_dom=Rectangle(1,1,2,useElementsOnFace=1)
81 self.checker(my_dom,file)
82
83 def test_hex_3D_order1(self):
84 file="hex_3D_order1.msh"
85 my_dom=Brick(1,1,1,1)
86 self.checker(my_dom,file)
87
88 def test_hex_3D_order1_onFace(self):
89 file="hex_3D_order1_onFace.msh"
90 my_dom=Brick(1,1,1,1,useElementsOnFace=1)
91 self.checker(my_dom,file)
92
93 def test_hex_3D_order2(self):
94 file="hex_3D_order2.msh"
95 my_dom=Brick(1,1,1,2)
96 self.checker(my_dom,file)
97
98 def test_hex_3D_order2_onFace(self):
99 file="hex_3D_order2_onFace.msh"
100 my_dom=Brick(1,1,1,2,useElementsOnFace=1)
101 self.checker(my_dom,file)
102
103 def test_hex_contact_2D_order1(self):
104 file="hex_contact_2D_order1.msh"
105 ms1=Rectangle(1,1,1,l1=0.5,useElementsOnFace=False)
106 ms2=Rectangle(1,1,1,l1=0.5,useElementsOnFace=False)
107 ms2.setX(ms2.getX()+[0,0.5])
108 my_dom=JoinFaces([ms1,ms2],optimize=False)
109 self.checker(my_dom,file)
110
111 def test_hex_contact_2D_order1_onFace(self):
112 file="hex_contact_2D_order1_onFace.msh"
113 ms1=Rectangle(1,1,1,l1=0.5,useElementsOnFace=True)
114 ms2=Rectangle(1,1,1,l1=0.5,useElementsOnFace=True)
115 ms2.setX(ms2.getX()+[0,0.5])
116 my_dom=JoinFaces([ms1,ms2],optimize=False)
117 self.checker(my_dom,file)
118
119 def test_hex_contact_2D_order2(self):
120 file="hex_contact_2D_order2.msh"
121 ms1=Rectangle(1,1,2,l1=0.5,useElementsOnFace=False)
122 ms2=Rectangle(1,1,2,l1=0.5,useElementsOnFace=False)
123 ms2.setX(ms2.getX()+[0,0.5])
124 my_dom=JoinFaces([ms1,ms2],optimize=False)
125 self.checker(my_dom,file)
126
127 def test_hex_contact_2D_order2_onFace(self):
128 file="hex_contact_2D_order2_onFace.msh"
129 ms1=Rectangle(1,1,2,l1=0.5,useElementsOnFace=True)
130 ms2=Rectangle(1,1,2,l1=0.5,useElementsOnFace=True)
131 ms2.setX(ms2.getX()+[0,0.5])
132 my_dom=JoinFaces([ms1,ms2],optimize=False)
133 self.checker(my_dom,file)
134
135 def test_hex_contact_3D_order1(self):
136 file="hex_contact_3D_order1.msh"
137 ms1=Brick(1,1,1,1,l2=0.5,useElementsOnFace=False)
138 ms2=Brick(1,1,1,1,l2=0.5,useElementsOnFace=False)
139 ms2.setX(ms2.getX()+[0,0,0.5])
140 my_dom=JoinFaces([ms1,ms2],optimize=False)
141 self.checker(my_dom,file)
142
143 def test_hex_contact_3D_order1_onFace(self):
144 file="hex_contact_3D_order1_onFace.msh"
145 ms1=Brick(1,1,1,1,l2=0.5,useElementsOnFace=True)
146 ms2=Brick(1,1,1,1,l2=0.5,useElementsOnFace=True)
147 ms2.setX(ms2.getX()+[0,0,0.5])
148 my_dom=JoinFaces([ms1,ms2],optimize=False)
149 self.checker(my_dom,file)
150
151 def test_hex_contact_3D_order2(self):
152 file="hex_contact_3D_order2.msh"
153 ms1=Brick(1,1,1,2,l2=0.5,useElementsOnFace=False)
154 ms2=Brick(1,1,1,2,l2=0.5,useElementsOnFace=False)
155 ms2.setX(ms2.getX()+[0,0,0.5])
156 my_dom=JoinFaces([ms1,ms2],optimize=False)
157 self.checker(my_dom,file)
158
159 def test_hex_contact_3D_order2_onFace(self):
160 file="hex_contact_3D_order2_onFace.msh"
161 ms1=Brick(1,1,1,2,l2=0.5,useElementsOnFace=True)
162 ms2=Brick(1,1,1,2,l2=0.5,useElementsOnFace=True)
163 ms2.setX(ms2.getX()+[0,0,0.5])
164 my_dom=JoinFaces([ms1,ms2],optimize=False)
165 self.checker(my_dom,file)
166
167 class Test_GMSHReader(unittest.TestCase):
168 def compare(self, test_file, reference_file):
169 dom_string=open(test_file).read().splitlines()
170 ref_string=open(reference_file).read().splitlines()
171 self.failUnlessEqual(len(dom_string),len(ref_string),"number of lines in mesh files does not match reference")
172 for l in range(1,len(ref_string)):
173 self.failUnlessEqual(dom_string[l].strip(),ref_string[l].strip(),"line %d (%s) in mesh files does not match reference (%s)"%(l,ref_string[l].strip(),dom_string[l].strip()))
174
175 def test_Tri3(self):
176 file="tri3_gmsh.msh"
177 ref ="tri3.fly"
178 test = FINLEY_WORKDIR+os.sep+"tri3_test.fly"
179 dom = ReadGmsh(FINLEY_TEST_MESH_PATH+os.sep+file,2,optimize=False)
180 dom.write(test)
181 self.compare(test, FINLEY_TEST_MESH_PATH+os.sep+ref)
182
183 def test_Tri6(self):
184 file="tri6_gmsh.msh"
185 ref="tri6.fly"
186 test = FINLEY_WORKDIR+os.sep+"tri8_test.fly"
187 dom = ReadGmsh(FINLEY_TEST_MESH_PATH+os.sep+file,2,optimize=False)
188 dom.write(test)
189 self.compare(test, FINLEY_TEST_MESH_PATH+os.sep+ref)
190
191 def test_Tet4(self):
192 file="tet4_gmsh.msh"
193 ref="tet4.fly"
194 test = FINLEY_WORKDIR+os.sep+"tet4_test.fly"
195 dom = ReadGmsh(FINLEY_TEST_MESH_PATH+os.sep+file,3,optimize=False)
196 dom.write(test)
197 self.compare(test, FINLEY_TEST_MESH_PATH+os.sep+ref)
198
199 def test_Tet10(self):
200 file="tet10_gmsh.msh"
201 ref="tet10.fly"
202 test = FINLEY_WORKDIR+os.sep+"tet10_test.fly"
203 dom = ReadGmsh(FINLEY_TEST_MESH_PATH+os.sep+file,3,optimize=False)
204 dom.write(test)
205 self.compare(test, FINLEY_TEST_MESH_PATH+os.sep+ref)
206
207 class Test_Reader(unittest.TestCase):
208 def test_ReadWriteTagNames(self):
209 file="hex_2D_order2.msh"
210 test = FINLEY_WORKDIR+os.sep+"test.fly"
211 dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+file,3,optimize=False)
212 insertTagNames(dom,A=1,B=2)
213 dom.write(test)
214 dom2 = ReadMesh(test,3,optimize=False)
215 t=getTagNames(dom)
216 self.failUnless(len(t)==6)
217 self.failUnless("A" in t)
218 self.failUnless("B" in t)
219 self.failUnless(dom2.getTag("A") == 1)
220 self.failUnless(dom2.getTag("B") == 2)
221 self.failUnless(dom2.isValidTagName("A"))
222 self.failUnless(dom2.isValidTagName("B"))
223
224 class Test_Integration(unittest.TestCase):
225 TOL=EPSILON*100.
226 def __test_2DQ(self,dom,order):
227 x=Function(dom).getX()
228 x_bound=FunctionOnBoundary(dom).getX()
229 for i in xrange(order+1):
230 for j in xrange(order+1):
231 res=integrate(x[0]**i*x[1]**j)
232 ref=1./((i+1)*(j+1))
233 error=abs(res-ref)/abs(ref)
234 self.failUnless(error<=self.TOL,"integration for order (%s,%s) failed. True value = %s, calculated = %s"%(i,j,ref,res))
235
236 res=integrate(x_bound[0]**i*x_bound[1]**j)
237 ref=0
238 if i == 0:
239 ref += 2./(j+1)
240 else:
241 ref += 1./(j+1)
242 if j == 0:
243 ref += 2./(i+1)
244 else:
245 ref += 1./(i+1)
246 error=abs(res-ref)/abs(ref)
247 self.failUnless(error<=self.TOL,"surface integration for order (%s,%s) failed. True value = %s, calculated = %s"%(i,j,ref,res))
248
249 def __test_2DT(self,dom,order,raise_tol=1.):
250 x=Function(dom).getX()
251 x_bound=FunctionOnBoundary(dom).getX()
252 for i in xrange(order+1):
253 for j in xrange(order+1):
254 if i+j<=order:
255 res=integrate(x[0]**i*x[1]**j)
256 ref=1./((i+1)*(j+1))
257 error=abs(res-ref)/abs(ref)
258 # print order,error
259 self.failUnless(error<=self.TOL*raise_tol,"integration for order (%s,%s) failed. True value = %s, calculated = %s"%(i,j,ref,res))
260
261 res=integrate(x_bound[0]**i*x_bound[1]**j)
262 ref=0
263 if i == 0:
264 ref += 2./(j+1)
265 else:
266 ref += 1./(j+1)
267 if j == 0:
268 ref += 2./(i+1)
269 else:
270 ref += 1./(i+1)
271 error=abs(res-ref)/abs(ref)
272 self.failUnless(error<=self.TOL*raise_tol,"surface integration for order (%s,%s) failed. True value = %s, calculated = %s"%(i,j,ref,res))
273
274
275 def __test_3DQ(self,dom,order):
276 x=Function(dom).getX()
277 x_bound=FunctionOnBoundary(dom).getX()
278 for i in xrange(order+1):
279 for j in xrange(order+1):
280 for k in xrange(order+1):
281 res=integrate(x[0]**i*x[1]**j*x[2]**k)
282 ref=1./((i+1)*(j+1)*(k+1))
283 error=abs(res-ref)/abs(ref)
284 self.failUnless(error<=self.TOL,"integration for order (%s,%s,%s) failed. True value = %s, calculated = %s"%(i,j,k,ref,res))
285
286 res=integrate(x_bound[0]**i*x_bound[1]**j*x_bound[2]**k)
287 ref=0
288 if i == 0:
289 ref += 2./((j+1)*(k+1))
290 else:
291 ref += 1./((j+1)*(k+1))
292 if j == 0:
293 ref += 2./((i+1)*(k+1))
294 else:
295 ref += 1./((i+1)*(k+1))
296 if k == 0:
297 ref += 2./((i+1)*(j+1))
298 else:
299 ref += 1./((i+1)*(j+1))
300 error=abs(res-ref)/abs(ref)
301 self.failUnless(error<=self.TOL,"surface integration for order (%s,%s,%s) failed. True value = %s, calculated = %s"%(i,j,k,ref,res))
302
303 def __test_3DT(self,dom,order,raise_tol=1.):
304 x=Function(dom).getX()
305 x_bound=FunctionOnBoundary(dom).getX()
306 for i in xrange(order+1):
307 for j in xrange(order+1):
308 for k in xrange(order+1):
309 if i+j+k<=order:
310 res=integrate(x[0]**i*x[1]**j*x[2]**k)
311 ref=1./((i+1)*(j+1)*(k+1))
312 error=abs(res-ref)/abs(ref)
313 self.failUnless(error<=self.TOL*raise_tol,"integration for order (%s,%s,%s) failed. True value = %s, calculated = %s"%(i,j,k,ref,res))
314
315 res=integrate(x_bound[0]**i*x_bound[1]**j*x_bound[2]**k)
316 ref=0
317 if i == 0:
318 ref += 2./((j+1)*(k+1))
319 else:
320 ref += 1./((j+1)*(k+1))
321 if j == 0:
322 ref += 2./((i+1)*(k+1))
323 else:
324 ref += 1./((i+1)*(k+1))
325 if k == 0:
326 ref += 2./((i+1)*(j+1))
327 else:
328 ref += 1./((i+1)*(j+1))
329 error=abs(res-ref)/abs(ref)
330 self.failUnless(error<=self.TOL*raise_tol,"surface integration for order (%s,%s,%s) failed. True value = %s, calculated = %s"%(i,j,k,ref,res))
331
332 def test_hex2D_order1(self):
333 my_dom=Rectangle(1,1,integrationOrder=1)
334 self.__test_2DQ(my_dom,1)
335 def test_hex2D_order2(self):
336 my_dom=Rectangle(1,1,integrationOrder=2)
337 self.__test_2DQ(my_dom,2)
338 def test_hex2D_order3(self):
339 my_dom=Rectangle(1,1,integrationOrder=3)
340 self.__test_2DQ(my_dom,3)
341 def test_hex2D_order4(self):
342 my_dom=Rectangle(1,1,integrationOrder=4)
343 self.__test_2DQ(my_dom,4)
344 def test_hex2D_order5(self):
345 my_dom=Rectangle(1,1,integrationOrder=5)
346 self.__test_2DQ(my_dom,5)
347 def test_hex2D_order6(self):
348 my_dom=Rectangle(1,1,integrationOrder=6)
349 self.__test_2DQ(my_dom,6)
350 def test_hex2D_order7(self):
351 my_dom=Rectangle(1,1,integrationOrder=7)
352 self.__test_2DQ(my_dom,7)
353 def test_hex2D_order8(self):
354 my_dom=Rectangle(1,1,integrationOrder=8)
355 self.__test_2DQ(my_dom,8)
356 def test_hex2D_order9(self):
357 my_dom=Rectangle(1,1,integrationOrder=9)
358 self.__test_2DQ(my_dom,9)
359 def test_hex2D_order10(self):
360 my_dom=Rectangle(1,1,integrationOrder=10)
361 self.__test_2DQ(my_dom,10)
362
363 def test_Tet2D_order1(self):
364 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=1)
365 self.__test_2DT(my_dom,1)
366 def test_Tet2D_order2(self):
367 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=2)
368 self.__test_2DT(my_dom,2)
369 def test_Tet2D_order3(self):
370 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=3)
371 self.__test_2DT(my_dom,3)
372 def test_Tet2D_order4(self):
373 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=4)
374 self.__test_2DT(my_dom,4)
375 def test_Tet2D_order5(self):
376 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=5)
377 self.__test_2DT(my_dom,5)
378 def test_Tet2D_order6(self):
379 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=6)
380 self.__test_2DT(my_dom,6)
381 def test_Tet2D_order7(self):
382 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=7)
383 self.__test_2DT(my_dom,7)
384 def test_Tet2D_order8(self):
385 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=8)
386 self.__test_2DT(my_dom,8,1./sqrt(EPSILON))
387 def test_Tet2D_order9(self):
388 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=9)
389 self.__test_2DT(my_dom,9,1./sqrt(EPSILON))
390 def test_Tet2D_order10(self):
391 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tri3.fly",optimize=False,integrationOrder=10)
392 self.__test_2DT(my_dom,10)
393
394 def test_hex3D_order1(self):
395 my_dom=Brick(1,1,1,integrationOrder=1)
396 self.__test_3DQ(my_dom,1)
397 def test_hex3D_order2(self):
398 my_dom=Brick(1,1,1,integrationOrder=2)
399 self.__test_3DQ(my_dom,2)
400 def test_hex3D_order3(self):
401 my_dom=Brick(1,1,1,integrationOrder=3)
402 self.__test_3DQ(my_dom,3)
403 def test_hex3D_order4(self):
404 my_dom=Brick(1,1,1,integrationOrder=4)
405 self.__test_3DQ(my_dom,4)
406 def test_hex3D_order5(self):
407 my_dom=Brick(1,1,1,integrationOrder=5)
408 self.__test_3DQ(my_dom,5)
409 def test_hex3D_order6(self):
410 my_dom=Brick(1,1,1,integrationOrder=6)
411 self.__test_3DQ(my_dom,6)
412 def test_hex3D_order7(self):
413 my_dom=Brick(1,1,1,integrationOrder=7)
414 self.__test_3DQ(my_dom,7)
415 def test_hex3D_order8(self):
416 my_dom=Brick(1,1,1,integrationOrder=8)
417 self.__test_3DQ(my_dom,8)
418 def test_hex3D_order9(self):
419 my_dom=Brick(1,1,1,integrationOrder=9)
420 self.__test_3DQ(my_dom,9)
421 def test_hex3D_order10(self):
422 my_dom=Brick(1,1,1,integrationOrder=10)
423 self.__test_3DQ(my_dom,10)
424
425 def test_Tet3D_order1(self):
426 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=1)
427 self.__test_3DT(my_dom,1)
428 def test_Tet3D_order2(self):
429 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=2)
430 self.__test_3DT(my_dom,2)
431 def test_Tet3D_order3(self):
432 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=3)
433 self.__test_3DT(my_dom,3)
434 def test_Tet3D_order4(self):
435 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=4)
436 self.__test_3DT(my_dom,4)
437 def test_Tet3D_order5(self):
438 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=5)
439 self.__test_3DT(my_dom,5)
440 def test_Tet3D_order6(self):
441 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=6)
442 self.__test_3DT(my_dom,6,1./sqrt(EPSILON))
443 def test_Tet3D_order7(self):
444 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=7)
445 self.__test_3DT(my_dom,7,1./sqrt(EPSILON))
446 def test_Tet3D_order8(self):
447 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=8)
448 self.__test_3DT(my_dom,8,1./sqrt(EPSILON))
449 def test_Tet3D_order9(self):
450 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=9)
451 self.__test_3DT(my_dom,9,1./sqrt(EPSILON))
452 def test_Tet3D_order10(self):
453 my_dom = ReadMesh(FINLEY_TEST_MESH_PATH+os.sep+"tet4.fly",optimize=False,integrationOrder=10)
454 self.__test_3DT(my_dom,10)
455
456 if __name__ == '__main__':
457 suite = unittest.TestSuite()
458 suite.addTest(unittest.makeSuite(Test_Generators))
459 suite.addTest(unittest.makeSuite(Test_GMSHReader))
460 suite.addTest(unittest.makeSuite(Test_Reader))
461 suite.addTest(unittest.makeSuite(Test_Integration))
462 s=unittest.TextTestRunner(verbosity=2).run(suite)

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26