source: anuga_core/source/anuga/pmesh/test_mesh.py @ 4203

Last change on this file since 4203 was 4165, checked in by duncan, 18 years ago

removing read/write of .xya files from loadASCII.py

File size: 85.8 KB
Line 
1#!/usr/bin/env python
2#
3import tempfile
4import unittest
5
6#try:
7from anuga.pmesh.mesh import *
8#except ImportError: 
9#    from mesh import *
10
11
12from load_mesh.loadASCII import *
13from anuga.coordinate_transforms.geo_reference import Geo_reference
14from anuga.geospatial_data.geospatial_data import Geospatial_data
15from anuga.utilities.polygon import  is_inside_polygon ### inside_polygon
16
17class meshTestCase(unittest.TestCase):
18    def setUp(self):
19        pass
20   
21    def tearDown(self):
22        pass
23
24    def testPointDistance(self):
25        a = Point(0.0, 0.0)
26        b = Point(0.0, 10.0)
27       
28        self.failUnless( a.DistanceToPoint(b) == 10.0,
29                        'Point DistanceToPoint is wrong!')
30   
31    def testVertexDistance(self):
32        a = Vertex (0.0, 0.0)
33        b = Vertex (0.0, 10.0)
34       
35        self.failUnless( a.DistanceToPoint(b) == 10.0,
36                        'Point DistanceToPoint is wrong!')
37       
38    def testTriangle(self):
39        a = Vertex (0.0, 0.0)
40        b = Vertex (0.0, 2.0)
41        c = Vertex (2.0,0.0)
42        d = Vertex (0.0, 4.0)
43        e = Vertex (2.0, 2.0)
44        f = Vertex (4.0,0.0)
45       
46        t1 = Triangle(b,a,c)       
47        t2 = Triangle(b,c,e)     
48        t3 = Triangle(e,c,f)     
49        t4 = Triangle(d,b,e)
50        t2.setNeighbors(t3,t4,t1)
51       
52        self.failUnless( t2.neighbors[2].vertices[0] == b, 'Triangle initialisation is wrong!')
53   
54       
55    def testSegment(self):
56        a = Vertex (0.0, 0.0)
57        b = Vertex (0.0, 10.0)
58        s = Segment(a,b, tag = 20)     
59       
60        self.failUnless( s.vertices[0].DistanceToPoint(s.vertices[1]) == 10.0,
61                        'vertices in a segment are wrong')
62       
63        self.failUnless( s.tag == 20.0,
64                        'tag in a segment are wrong')
65
66    def testdeleteUserVertex(self):
67
68       
69        mesh = Mesh()
70        a = mesh.addUserVertex(0.0, 0.0)
71        b = mesh.addUserVertex (0.0, 2.0)
72        c = mesh.addUserVertex (2.0,0.0)
73       
74        s1 = mesh.addUserSegment(a,b)
75        s2 = mesh.addUserSegment(a,c)
76        s3 = mesh.addUserSegment(c,b)
77
78        mesh.deleteMeshObject (a) 
79        self.failUnless(mesh.userSegments[0] == s3,
80                        'Bad segment. ')       
81        self.failUnless(len(mesh.userSegments) ==1,
82                        'Segments not deleted.')
83        self.failUnless(len(mesh.userVertices) == 2,
84                        'Vertex not deleted.')
85       
86 
87    # FIXME add test for minAngle   
88    def testgenerateMesh(self):
89        a = Vertex (0.0, 0.0)
90        d = Vertex (0.0, 4.0)
91        f = Vertex (4.0,0.0)
92
93        s1 = Segment(a,d)
94        s2 = Segment(d,f)
95        s3 = Segment(a,f)
96
97        r1 = Region(0.3, 0.3,tag = 1.3,maxArea = .6)
98        #print r1
99        m = Mesh(userVertices=[a,d,f], userSegments=[s1,s2,s3], regions=[r1] )
100       
101        m.generateMesh("Q", maxArea = 2.1 )         
102
103        #print m
104
105        #m.plotMeshTriangle()
106
107        result = 1.414214
108        delta  = 0.00001
109       
110        self.failUnless((m.meshTriangles[1].vertices[0].x < result + delta) or
111                        (m.meshTriangles[1].vertices[0].x > result - delta),
112                        'generated mesh is wrong!')
113
114    def test_regionalMaxArea(self):
115        v0 = Vertex (0.0, 0.0)
116        v1 = Vertex (6.0, 0.0)
117        v2 = Vertex (6.0,6.0)
118        v3 = Vertex (0.0,6.0)
119
120        s1 = Segment(v0,v1)
121        s2 = Segment(v1,v2)
122        s3 = Segment(v3,v2)
123        s4 = Segment(v3,v0)
124        s5 = Segment(v2,v0)
125
126        r1 = Region(3, 1,tag = 1.3)
127        #print r1
128        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5], regions=[r1] )
129       
130        m.generateMesh("Q", maxArea = 36 )         
131
132        #m.plotMeshTriangle()
133        #print "len(m.meshTriangles)",len(m.meshTriangles)
134
135        self.failUnless(len(m.meshTriangles) == 2, 
136                        'test_regionalMaxArea 1:generated mesh is wrong!')
137       
138        ## Another test case
139        r1 = Region(3, 1,tag = 1.3)
140        r2 = Region(1, 3,tag = 1.3, maxArea = 8)
141        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
142                 regions=[r1,r2] )
143        m.generateMesh("Q", maxArea = 36 )
144       
145        self.failUnless(len(m.meshTriangles) >= 6,
146                        'testregion_with_maxarea 2: # of tris is wrong!')   
147       
148               
149        ## Another test case
150        r1 = Region(3, 1, tag = 1.3, maxArea = 8)
151        r2 = Region(1, 3, tag = 1.3, maxArea = 8)
152        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
153                 regions=[r1,r2] )
154        m.generateMesh("Q", maxArea = 36 ) 
155        #print "len(m.meshTriangles)",len(m.meshTriangles)
156       
157        self.failUnless(len(m.meshTriangles) >= 8,
158                        'testregion_with_maxarea 3: # of tris is wrong!')
159               
160               
161        ## Another test case
162        r1 = Region(3, 1, tag = 1.3 )
163        r2 = Region(1, 3, tag = 1.3, maxArea = 8)
164        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
165                 regions=[r1,r2] )
166        m.generateMesh("Q", maxArea = 8 ) 
167        self.failUnless(len(m.meshTriangles) >= 8,
168                        'testregion_with_maxarea 4: # of tris is wrong!')   
169
170       
171        ## Another test case
172        r1 = Region(3, 1,tag = 1.3, maxArea = 8)
173        r2 = Region(1, 3,tag = 1.3, maxArea = 8)
174        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
175                 regions=[r1,r2] )
176        m.generateMesh("Q", maxArea = 36,isRegionalMaxAreas = False )     
177        self.failUnless(len(m.meshTriangles) == 2, 
178                        'test_regionalMaxArea 5:generated mesh is wrong!')
179       
180        ## Another test case
181        r1 = Region(3, 1,tag = 1.3, maxArea = 8)
182        r2 = Region(1, 3,tag = 1.3, maxArea = 8)
183        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
184                 regions=[r1,r2] )
185        m.generateMesh("Q",isRegionalMaxAreas = False )
186        self.failUnless(len(m.meshTriangles) == 2, 
187                        'test_regionalMaxArea 5:generated mesh is wrong!')
188       
189    def test_generate_mesh(self):
190        v0 = Vertex (0.0, 0.0)
191        v1 = Vertex (6.0, 0.0)
192        v2 = Vertex (6.0,6.0)
193        v3 = Vertex (0.0,6.0)
194
195        s1 = Segment(v0,v1)
196        s2 = Segment(v1,v2)
197        s3 = Segment(v3,v2)
198        s4 = Segment(v3,v0)
199        s5 = Segment(v2,v0)
200
201        r1 = Region(3, 1,tag = 1.3)
202        #print r1
203        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
204                 regions=[r1] )
205       
206        m.generate_mesh(maximum_triangle_area=36,verbose=False)         
207
208        self.failUnless(len(m.meshTriangles) == 2, 
209                        'test_regionalMaxArea 1:generated mesh is wrong!')
210       
211        ## Another test case
212        r1 = Region(3, 1,tag = 1.3)
213        r2 = Region(1, 3,tag = 1.3, maxArea = 8)
214        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
215                 regions=[r1,r2] )
216        m.generate_mesh(maximum_triangle_area=36,verbose=False) 
217       
218        self.failUnless(len(m.meshTriangles) >= 6,
219                        'testregion_with_maxarea 2: # of tris is wrong!')   
220               
221        ## Another test case
222        r1 = Region(3, 1, tag = 1.3, maxArea = 8)
223        r2 = Region(1, 3, tag = 1.3, maxArea = 8)
224        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
225                 regions=[r1,r2] )
226        m.generate_mesh(maximum_triangle_area=36,verbose=False)         
227        #print "len(m.meshTriangles)",len(m.meshTriangles)
228       
229        self.failUnless(len(m.meshTriangles) >= 8,
230                        'testregion_with_maxarea 3: # of tris is wrong!')
231                         
232        ## Another test case
233        r1 = Region(3, 1, tag = 1.3 )
234        r2 = Region(1, 3, tag = 1.3, maxArea = 8)
235        m = Mesh(userVertices=[v0,v1,v2,v3], userSegments=[s1,s2,s3,s4,s5],
236                 regions=[r1,r2] )
237        m.generate_mesh(maximum_triangle_area=8,verbose=False)   
238        self.failUnless(len(m.meshTriangles) >= 8,
239                        'testregion_with_maxarea 4: # of tris is wrong!')   
240
241        ## Another test case r1 = Region(3, 1,tag = 1.3, maxArea = 8)
242        r2 = Region(1, 3,tag = 1.3, maxArea = 8)
243        m = Mesh(userVertices=[v0,v1,v2,v3],
244        userSegments=[s1,s2,s3,s4,s5], regions=[r1,r2] )
245        m.generate_mesh(verbose=False)
246        #print "en(m.meshTriangles)", len(m.meshTriangles)
247        self.failUnless(len(m.meshTriangles) >= 8,
248        'You have issues!')
249       
250    def testdeleteUserVertex(self):
251        mesh = Mesh()
252        a = mesh.addUserVertex(0.0, 0.0)
253        b = mesh.addUserVertex (0.0, 2.0)
254        c = mesh.addUserVertex (2.0,0.0)
255       
256        s1 = mesh.addUserSegment(a,b)
257        s2 = mesh.addUserSegment(a,c)
258        s3 = mesh.addUserSegment(c,b)
259
260        mesh.deleteMeshObject (s2)
261       
262        #print ",s2 in mesh.userSegments" ,s2 in mesh.userSegments
263        self.failUnless(not(s2 in mesh.userSegments),
264                        'Bad segment. ')       
265        self.failUnless(len(mesh.userSegments) ==2,
266                        'Segments not deleted.')
267        self.failUnless(len(mesh.userVertices) == 3,
268                        'Vertex deleted, instead of segment.')
269
270    def testTriangleArea(self):
271        a = Vertex (10.0, 10.0)
272        b = Vertex (10.0, 20.0)
273        c = Vertex (20.0,10.0)
274       
275        d = Vertex (-20.0, 0.0)
276        e = Vertex (-20.0, -20.0)
277        f = Vertex (20.0,-20.0)
278       
279        t1 = Triangle(b,a,c)     
280        t2 = Triangle(e,d,f)
281       
282#         print "t1", t1
283#         print "t1 area ", t1.calcArea()
284#         print "t2", t2
285#         print "t2 area ", t2.calcArea()
286        self.failUnless( t1.calcArea() == 50 and t2.calcArea() == 400, 'Triangle area is wrong!')
287    def testisUserSegmentNew (self):
288        mesh = Mesh()
289        a = mesh.addUserVertex(0.0, 0.0)
290        b = mesh.addUserVertex (0.0, 2.0)
291        c = mesh.addUserVertex (2.0,0.0)
292        d = mesh.addUserVertex (2.0,3.0)
293       
294        s1 = mesh.addUserSegment(a,b)
295        s2 = mesh.addUserSegment(a,c)
296        s3 = mesh.addUserSegment(c,b)
297
298        self.failUnless(mesh.isUserSegmentNew(a,d) ,
299                        'Segment should be new. ')
300        self.failUnless(not(mesh.isUserSegmentNew(a,b)) ,
301                        'Segment should not be new. ')
302
303
304    def testisUserSegmentNewII (self):
305        mesh = Mesh()
306        a = mesh.addUserVertex(0.0, 0.0)
307        b = mesh.addUserVertex (0.0, 2.0)
308        c = mesh.addUserVertex (2.0,0.0)
309        d = mesh.addUserVertex (2.0,3.0)
310       
311        s1 = mesh.addUserSegment(a,b)
312        s2 = mesh.addUserSegment(a,c)
313        s3 = mesh.addUserSegment(c,b)
314
315        self.failUnless(mesh.representedUserSegment(a,d) == None,
316                        'Segment should be new. ')
317        self.failUnless(mesh.representedUserSegment(a,b) == s1 ,
318                        'Segment should not be new. ')
319       
320    def testauto_segment(self):
321        p0 = Vertex (0.0, 0.0)
322        p1 = Vertex (0.0, 4.0)
323        p2 = Vertex (4.0,4.0)
324        p3 = Vertex (4.0,0.0)
325
326        s1 = Segment(p0,p1)
327       
328        m = Mesh(userVertices=[p0, p1, p2, p3], userSegments=[s1] ) 
329        m.auto_segment()
330       
331        #print 'Len', len(m.userSegments)
332        self.failUnless(len(m.getUserSegments()) == 4 ,
333                        'userSegments is wrong!')
334     
335        m.auto_segment()
336        self.failUnless(len(m.getUserSegments()) == 4 ,
337                        'userSegments is wrong!')
338     
339    def testauto_segmentII(self):
340        p1 = Vertex (3.0, 4.0)
341        p2 = Vertex (3.0,2.0)
342        p3 = Vertex (3.0,0.0)
343        p4 = Vertex (6.0, 4.0)
344        p5 = Vertex (6.0,2.0)
345        p0 = Vertex (6.0,0.0)
346
347
348        s1 = Segment(p2,p3)
349        s2 = Segment(p4,p5)
350       
351        m = Mesh(userVertices=[p0, p1, p2, p3, p4, p5],
352                 userSegments=[s1, s2])     
353
354        m.auto_segment()
355       
356        s3 = m.representedAlphaUserSegment(p3,p0)
357        self.failUnless(not (s3 == None) ,
358                        'userSegments is wrong!')
359
360       
361        s6 = m.representedAlphaUserSegment(p1,p4)       
362        self.failUnless(not (s6 == None) ,
363                        'userSegments is wrong!')
364       
365        # remove a segment, add a point, auto_segment
366        m.alphaUserSegments.remove(s3)
367        p6 = Vertex (1.0, 2.0)
368        m.userVertices.append(p6)
369       
370        m.auto_segment()
371       
372        s1_now = m.representedUserSegment(p3,p2)
373        self.failUnless(s1_now == s1 ,
374                        'userSegments is wrong!')
375       
376        s2_now = m.representedUserSegment(p5,p4)       
377        self.failUnless(s2_now == s2 ,
378                        'userSegments is wrong!')
379       
380        s3 = m.representedAlphaUserSegment(p3,p6)       
381        self.failUnless(not (s3 == None) ,
382                        'userSegments is wrong!')
383       
384        s4 = m.representedAlphaUserSegment(p3,p6)       
385        self.failUnless(not (s4 == None) ,
386                        'userSegments is wrong!')
387       
388        s5 = m.representedAlphaUserSegment(p4,p6)       
389        self.failUnless(s5 == None ,
390                        'userSegments is wrong!')
391        #print m
392       
393    def testRegions(self):
394        a = Vertex (0.0, 0.0)
395        b = Vertex (0.0, 2.0)
396        c = Vertex (2.0,0.0)
397        d = Vertex (0.0, 4.0)
398        e = Vertex (2.0, 2.0)
399        f = Vertex (4.0,0.0)
400        g = Vertex (0.0,-2.0)
401       
402        Segment.set_default_tag("")
403        s1 = Segment(a,b)
404        s2 = Segment(b,e)
405        s3 = Segment(e,c)
406        s4 = Segment(c,a)
407        s5 = Segment(b,d)
408        s6 = Segment(e,d)
409        s7 = Segment(e,f)
410        s8 = Segment(c,f)
411        s9 = Segment(g,c)
412        s10 = Segment(g,a)
413
414        r1 = Region(0.1,0.1,tag="22")
415        r2 = Region(0.1,2.1,tag="11")
416        r3 = Region(2.1,0.1)
417       
418        m = Mesh(userVertices=[a,b,c,d,e,f,g], userSegments=[s1,s2,s3,s4,s5,s6,s7,s8,s9,s10], regions=[r1,r2,r3] )
419        m.generateMesh("Q", maxArea = 2.1 )
420        #print m
421        Triangulation =  m.getTriangulation()
422        #print Triangulation[0].attribute
423        #print Triangulation[1].attribute
424        #print Triangulation[2].attribute
425        #print Triangulation[3].attribute
426        #print Triangulation[4].attribute
427       
428        self.failUnless(Triangulation[0].attribute == "" and
429                        Triangulation[1].attribute == "22" and
430                        Triangulation[2].attribute == "" and
431                        Triangulation[3].attribute == "11" and
432                        Triangulation[4].attribute == "22" ,
433                        'region attributes are wrong!')   
434
435    def test_vertexAttribs(self):
436        a = Vertex (0.0, 0.0, attributes = [12.0,2.0])
437        d = Vertex (0.0, 4.0, attributes = [9.0,7.0])
438        f = Vertex (4.0,0.0, attributes = [14.0,3.0])
439   
440        Segment.set_default_tag("")
441        s1 = Segment(a,d)
442        s2 = Segment(d,f)
443        s3 = Segment(a,f)
444     
445        r1 = Region(0.3, 0.3, tag = 88.9)
446     
447        m = Mesh(userVertices=[a,d,f], userSegments=[s1,s2,s3], regions=[r1])
448
449        m.generateMesh("Q", maxArea = 2.1)
450
451        vert = m.getMeshVertices()
452       
453        self.failUnless(vert[0].attributes == [12.0, 2.0] and
454                        vert[1].attributes == [9.0, 7.0] and
455                        vert[2].attributes == [14.0,3.0] and
456                        vert[3].attributes == [12.232233047033631, 4.4142135623730949] and
457                        vert[4].attributes == [13.0, 2.5] ,
458                        'vertex attributes are wrong!')
459
460       
461    def test_vertexAttribs2(self):
462   
463        a = Vertex (0.0, 0.0)
464        d = Vertex (0.0, 4.0)
465        f = Vertex (4.0,0.0)
466   
467        Segment.set_default_tag("")
468        s1 = Segment(a,d)
469        s2 = Segment(d,f)
470        s3 = Segment(a,f)
471     
472        r1 = Region(0.3, 0.3, tag = 88.9)
473     
474        m = Mesh(userVertices=[a,d,f], userSegments=[s1,s2,s3], regions=[r1])
475
476        m.generateMesh("Q", maxArea = 2.1 )
477
478        vert = m.getMeshVertices()
479        self.failUnless(vert[0].attributes == [] and
480                        vert[1].attributes == [] and
481                        vert[2].attributes == [] and
482                        vert[3].attributes == [] and
483                        vert[4].attributes == [],
484                        'vertex attributes are wrong!')
485
486    def test_segtag(self):
487   
488        a = Vertex (0.0, 0.0)
489        d = Vertex (0.0, 4.0)
490        f = Vertex (4.0,0.0)
491   
492        s1 = Segment(a,d,tag = 5)
493        s2 = Segment(d,f,tag = 7)
494        s3 = Segment(a,f,tag = 9)
495     
496        m = Mesh(userVertices=[a,d,f], userSegments=[s1,s2,s3])
497
498        m.generateMesh("Q", maxArea = 2.1 )
499
500        #m.export_mesh_file("from_test_mesh.tsh")
501        seg = m.getMeshSegments()
502        #print "seg",seg
503        #print "seg[0].tag"
504        #print seg[0].tag
505        #print "seg[0].tag"
506       
507        self.failUnless(seg[0].tag == 5 and
508                        seg[1].tag == 7 and
509                        seg[2].tag == 9 and
510                        seg[3].tag == 7 and
511                        seg[4].tag == 9,
512                        'seg tags are wrong')
513           
514
515    def test_segtag2(self):
516   
517        a = Vertex (0.0, 0.0)
518        d = Vertex (0.0, 4.0)
519        f = Vertex (4.0,0.0)
520        e = Vertex (1.0,1.0)
521   
522        s1 = Segment(a,d)
523        s2 = Segment(d,f)
524        s3 = Segment(a,f)
525        s4 = Segment(a,e)
526     
527        m = Mesh(userVertices=[a,d,f,e], userSegments=[s1,s2,s3,s4])
528
529        m.generateMesh("Q", maxArea = 2.1)
530
531        seg = m.getMeshSegments()
532        self.failUnless(seg[0].tag == "exterior" and
533                        seg[1].tag == "exterior" and
534                        seg[2].tag == "exterior" and
535                        seg[3].tag == "" and
536                        seg[4].tag == "exterior",
537                        '2nd seg tags are wrong')
538
539    def test_asciiFile(self):
540   
541        a = Vertex (0.0, 0.0)  #, attributes = [1.1])
542        d = Vertex (0.0, 4.0)  #, attributes = [1.2])
543        f = Vertex (4.0,0.0)  #, attributes = [1.3])
544        e = Vertex (1.0,1.0)  #, attributes = [1.4])
545   
546        s1 = Segment(a,d)
547        s2 = Segment(d,f)
548        s3 = Segment(a,f)
549        s4 = Segment(a,e)
550     
551        m = Mesh(userVertices=[a,d,f,e], userSegments=[s1,s2,s3,s4])
552
553        m.generateMesh("Q", maxArea = 2.1 )
554        seg = m.getMeshSegments()
555       
556        fileName = tempfile.mktemp(".tsh")
557        m.export_mesh_file(fileName)
558        file = open(fileName)
559        lFile = file.read().split('\n')
560        file.close()
561        os.remove(fileName)
562       
563        #print "@^@^"
564        #for l in lFile:
565        #    print l,"<"
566        #print "@^@^"
567
568        # no need to check the title again
569        #self.failUnless(lFile[0] == "5 0 # <vertex #> <x> <y> [attributes] ...Triangulation Vertices..."
570          #              ,'Ascii file is wrong, vertex title')
571        self.failUnless(lFile[1] == "0 0.0 0.0 " and #1.1 " and
572                        lFile[2] == "1 0.0 4.0 " and #1.2 " and
573                        lFile[3] == "2 4.0 0.0 " and #1.3 " and
574                        lFile[4] == "3 1.0 1.0 " and #1.4 " and
575                        lFile[5] == "4 2.0 2.0 "  #1.25 "
576                        ,
577                        'Ascii file is wrong, vertex')
578       
579        #self.failUnless(lFile[6] == "# attribute column titles ...Triangulation Vertex Titles..."
580          #              ,'Ascii file is wrong, attribute column title')
581        self.failUnless(lFile[8] == "0 3 2 4 -1 2 3  " and
582                        lFile[9] == "1 1 0 3 3 2 -1  " and
583                        lFile[10] == "2 3 4 1 -1 1 0  " and
584                        lFile[11] == "3 2 3 0 1 -1 0  "
585                        ,
586                        'Ascii file is wrong, triangle') 
587
588        self.failUnless( lFile[13] == "0 0 1 exterior" and
589                        lFile[14] == "1 1 4 exterior" and
590                        lFile[15] == "2 2 0 exterior" and
591                        lFile[16] == "3 0 3 " and
592                        lFile[17] == "4 4 2 exterior" ,
593                        'Ascii file is wrong, segment')
594       
595       # self.failUnless(lFile[18] == '4 0 # <vertex #> <x> <y> [attributes] ...Mesh Vertices...',
596        #                'Ascii file is wrong, Mesh Vertices Title')
597       
598        self.failUnless(lFile[19] == '0 0.0 0.0 ' and #1.1 ' and
599                        lFile[20] == '1 0.0 4.0 ' and #1.2 ' and
600                        lFile[21] == '2 4.0 0.0 ' and #1.3 ' and
601                        lFile[22] == '3 1.0 1.0 ' #1.4 '
602                        ,
603                        'Ascii file is wrong, Mesh Vertices II')
604       
605        self.failUnless(lFile[24] == '0 0 1 ' and
606                        lFile[25] == '1 1 2 ' and
607                        lFile[26] == '2 0 2 ' and
608                        lFile[27] == '3 0 3 '
609                        ,'Ascii file is wrong, Mesh Segments')       
610
611 
612    def test_ascii_file(self):
613   
614        a = Vertex (0.0, 0.0) #, attributes = [1.1])
615        d = Vertex (0.0, 4.0) #, attributes = [1.2])
616        f = Vertex (4.0,0.0) #, attributes = [1.3])
617        e = Vertex (1.0,1.0) #, attributes = [1.4])
618   
619        s1 = Segment(a,d)
620        s2 = Segment(d,f)
621        s3 = Segment(a,f)
622        s4 = Segment(a,e)
623     
624        m = Mesh(userVertices=[a,d,f,e], userSegments=[s1,s2,s3,s4])
625
626        m.generateMesh("Q", maxArea = 2.1 )
627
628        seg = m.getMeshSegments()
629       
630        fileName = tempfile.mktemp(".tsh")
631        m.export_mesh_file(fileName)
632        file = open(fileName)
633        lFile = file.read().split('\n')
634        file.close()
635        os.remove(fileName)
636       
637        #print "@^@^"
638        #for l in lFile:
639        #    print l,"<"
640        #print "@^@^"
641        self.failUnless(lFile[0] == "5 0 # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices..."
642                        ,
643                        'Ascii file is wrong, vertex title')
644        self.failUnless(lFile[1] == "0 0.0 0.0 " and #1.1 " and
645                        lFile[2] == "1 0.0 4.0 " and #1.2 " and
646                        lFile[3] == "2 4.0 0.0 " and #1.3 " and
647                        lFile[4] == "3 1.0 1.0 " and #1.4 " and
648                        lFile[5] == "4 2.0 2.0 "  #1.25 "
649                        ,
650                        'Ascii file is wrong, vertex')
651       
652        self.failUnless(lFile[6] == "# attribute column titles ...Triangulation Vertex Titles..."
653                        ,
654                        'Ascii file is wrong, attribute column title')
655        self.failUnless(lFile[7] == "4 # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." and
656                        lFile[8] == "0 3 2 4 -1 2 3  " and
657                        lFile[9] == "1 1 0 3 3 2 -1  " and
658                        lFile[10] == "2 3 4 1 -1 1 0  " and
659                        lFile[11] == "3 2 3 0 1 -1 0  "
660                        ,
661                        'Ascii file is wrong, triangle') 
662
663        self.failUnless(lFile[12] == "5 # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments..." and
664                        lFile[13] == "0 0 1 exterior" and
665                        lFile[14] == "1 1 4 exterior" and
666                        lFile[15] == "2 2 0 exterior" and
667                        lFile[16] == "3 0 3 " and
668                        lFile[17] == "4 4 2 exterior" ,
669                        'Ascii file is wrong, segment')
670       
671        self.failUnless(lFile[18] == '4 0 # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices...',
672                        'Ascii file is wrong, Mesh Vertices Title')
673       
674        self.failUnless(lFile[19] == '0 0.0 0.0 ' and #1.1 ' and
675                        lFile[20] == '1 0.0 4.0 ' and #1.2 ' and
676                        lFile[21] == '2 4.0 0.0 ' and #1.3 ' and
677                        lFile[22] == '3 1.0 1.0 ' #1.4 '
678                        ,
679                        'Ascii file is wrong, Mesh Vertices II')
680       
681        self.failUnless(lFile[23] == '4 # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments...' and
682                        lFile[24] == '0 0 1 ' and
683                        lFile[25] == '1 1 2 ' and
684                        lFile[26] == '2 0 2 ' and
685                        lFile[27] == '3 0 3 ' and
686                        lFile[28] == '0 # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes...' and
687                        lFile[29] == '0 # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions...'
688                        ,
689                        'Ascii file is wrong, Mesh Segments')       
690 
691
692    def test_thinoutVertices(self):
693
694        v1 = Vertex(-20,-20)
695        v2 = Vertex(-11,-11)
696        v3 = Vertex(-10,-10)
697        v4 = Vertex(-9,-1)
698        v5 = Vertex(-8,2)
699        v6 = Vertex(6,3)
700        v7 = Vertex(12,9)
701        v8 = Vertex(15,3)
702        v9 = Vertex(24,3)
703        m = Mesh(userVertices = [v1,v2,v3,v4,v5,v6,v7,v8,v9])
704        m.thinoutVertices(10)
705         
706        self.failUnless(v1 in m.userVertices,
707                        'test_thinoutVertices, test 1 failed')
708        self.failUnless(v3 in m.userVertices,
709                        'test_thinoutVertices, test 2 failed')
710        self.failUnless(v4 in m.userVertices,
711                        'test_thinoutVertices, test 3 failed')
712        self.failUnless(v6 in m.userVertices,
713                        'test_thinoutVertices, test 4 failed')
714        self.failUnless(v7 in m.userVertices,
715                        'test_thinoutVertices, test 5 failed')
716        self.failUnless(v9 in m.userVertices,
717                        'test_thinoutVertices, test 6 failed')
718        self.failUnless(v5 not in m.userVertices,
719                        'test_thinoutVertices, test 7 failed')
720        self.failUnless(v2 not in m.userVertices,
721                        'test_thinoutVertices, test 8 failed')
722        self.failUnless(v8 not in m.userVertices,
723                        'test_thinoutVertices, test 9 failed')
724
725    def test_same_x_y(self):
726        v = Point(7,8)
727        f = Point(7,8)
728        f.same_x_y(v)
729
730        self.failUnless(f.same_x_y(v),
731                        'same_x_y True failed')
732        e = Point(7,9)
733        self.failUnless(not f.same_x_y(e),
734                        'same_x_y False failed')
735
736    def test_import_tsh(self):
737       
738        a_att = [5.0,2.0]
739        d_att =[4.0,2.0]
740        f_att = [3.0,2.0]
741        e_att = [2.0,2.0]
742        a_xy = [0.0, 0.0]
743        a = Vertex ( a_xy[0],a_xy[1]) #, attributes =a_att)
744        d = Vertex (0.0, 4.0) #, attributes =d_att)
745        f = Vertex (4.0,0.0) #, attributes =f_att)
746        e = Vertex (1.0,1.0) #, attributes =e_att)
747   
748        s1 = Segment(a,d, tag = "50")
749        s2 = Segment(d,f, tag = "40")
750        s3 = Segment(a,f, tag = "30")
751        s4 = Segment(a,e, tag = "20")
752     
753        r1 = Region(0.3, 0.3,tag = "1.3")
754        geo = Geo_reference(8.9,8.9,65)
755        m = Mesh(userVertices=[a,d,f,e],
756                 userSegments=[s1,s2,s3,s4],
757                 regions=[r1],
758                 geo_reference=geo)
759
760        m.generateMesh("Q", maxArea = 2.1)
761        fileName = tempfile.mktemp(".tsh")
762        #print "dgs!!!"
763        #print "****************** fileName", fileName
764        m.export_mesh_file(fileName)
765        #print "******************"
766        #print "m", m
767        #print "******************"
768        m_returned = importMeshFromFile(fileName)
769        #print "m_returned",m_returned
770        #print "******************"
771        #print "****************** fileName", fileName
772        os.remove(fileName)
773        self.failUnless(0 == m.__cmp__(m_returned),
774                        'loading and saving of a mesh failed')
775        # do this when .msh supports geo_refs
776        #self.failUnless(m.geo_reference == m_returned.geo_reference,
777        #                'loading and saving of a mesh geo refs failed')
778
779    def test_import_mesh(self):
780       
781        a_att = [5.0,2.0]
782        d_att =[4.0,2.0]
783        f_att = [3.0,2.0]
784        e_att = [2.0,2.0]
785        a_xy = [0.0, 0.0]
786        a = Vertex ( a_xy[0],a_xy[1]) #, attributes =a_att)
787        d = Vertex (0.0, 4.0) #, attributes =d_att)
788        f = Vertex (4.0,0.0) #, attributes =f_att)
789        e = Vertex (1.0,1.0) #, attributes =e_att)
790   
791        s1 = Segment(a,d, tag = "50")
792        s2 = Segment(d,f, tag = "40")
793        s3 = Segment(a,f, tag = "30")
794        s4 = Segment(a,e, tag = "20")
795     
796        r1 = Region(0.3, 0.3,tag = "1.3")
797        geo = Geo_reference(65,8.9,8.9)
798        m = Mesh(userVertices=[a,d,f,e],
799                 userSegments=[s1,s2,s3,s4],
800                 regions=[r1],
801                 geo_reference=geo)
802
803        m.generateMesh("Q", maxArea = 2.1)
804        fileName = tempfile.mktemp(".msh")
805        #print "dgs!!!"
806        #print "****************** fileName", fileName
807        m.export_mesh_file(fileName)
808        #print "******************"
809        #print "m", m
810        #print "******************"
811        m_returned = importMeshFromFile(fileName)
812        #print "m_returned",m_returned
813        #print "******************"
814        #print "****************** fileName", fileName
815        os.remove(fileName)
816        #print "m.geo_reference",m.geo_reference
817        #print "m_returned.geo_reference,",m_returned.geo_reference
818        self.failUnless(0 == m.__cmp__(m_returned),
819                        'loading and saving of a mesh failed')
820        self.failUnless(m.geo_reference == m_returned.geo_reference,
821                        'loading and saving of a mesh geo refs failed')
822
823    def test_normaliseMesh(self):
824       
825        a_att = [5.0,2.0]
826        d_att =[4.0,2.0]
827        f_att = [3.0,2.0]
828        e_att = [2.0,2.0]
829        a_xy = [10.0, 10.0]
830        a = Vertex ( a_xy[0],a_xy[1], attributes =a_att)
831        d = Vertex (15.0, 10.0, attributes =d_att)
832        f = Vertex (10.0,20.0, attributes =f_att)
833        e = Vertex (15.0,20.0, attributes =e_att)
834   
835        s1 = Segment(a,d, tag = 50)
836        s2 = Segment(d,e, tag = 40)
837        s3 = Segment(e,f, tag = 30)
838        s4 = Segment(f,a, tag = 20)
839     
840        r1 = Region(0.3, 0.3,tag = 1.3)
841        m = Mesh(userVertices=[a,d,f,e],
842                 userSegments=[s1,s2,s3,s4],
843                 regions=[r1])
844        m.normaliseMesh(1,0,1)
845        [xmin, ymin, xmax, ymax] = m.boxsize()
846        [attmin, attmax] = m.maxMinVertAtt(0)
847        self.failUnless(attmin == 0.0 and attmax == 1.0,
848                        'normalise failed')
849        self.failUnless(xmin == 0.0 and ymin == 0.0 and xmax == 0.5 and ymax == 1.0,
850                        'normalise failed')
851        m.normaliseMesh(200,-100,5)
852        [xmin, ymin, xmax, ymax] = m.boxsize()
853        [attmin, attmax] = m.maxMinVertAtt(0)
854        self.failUnless(attmin == 0.0 and attmax == 5.0,
855                        'normalise failed')
856        self.failUnless(xmin == -100.0 and ymin == -100.0 and xmax == 0.0 and ymax == 100.0,
857                        'normalise failed')
858       
859    def test_exportASCIIsegmentoutlinefile(self):
860        a = Vertex (0,0)
861        b = Vertex (0,3)
862        c = Vertex (3,3)
863        d = Vertex (1,2)
864        e = Vertex (3,1)
865     
866        s1 = Segment(a,b, tag = "50")
867        s2 = Segment(b,c, tag = "40")
868        s3 = Segment(c,a, tag = "30")
869     
870        r1 = Region(2, 1,tag = "1.3")
871        h1 = Hole(1,4)
872        m = Mesh(userVertices=[a,b,c,d,e],
873                 userSegments=[s1,s2,s3],
874                 regions=[r1],
875                 holes = [h1])     
876
877        # vertex e is outside of the outline, so
878        # it is a loner and it is removed.
879        m.generateMesh("Q", maxArea = 2.1)
880        #print "mesh ***************dsg*", m
881
882        fileName = tempfile.mktemp(".tsh")
883        m.exportASCIIsegmentoutlinefile(fileName)
884       
885        m_returned = importMeshFromFile(fileName)
886
887        #print "m_returned ****",m_returned
888        #print "****************** fileName", fileName
889        os.remove(fileName)
890
891        #Trim mesh, so it should like like m_returned
892        m.meshVertices = []
893        m.meshTriangles = []
894        m.meshSegments = []
895        m.userVertices=[a,b,c]
896        #print "mesh ***************dsg*", m
897        #print "(m.__cmp__(m_returned)", m.__cmp__(m_returned)
898        self.failUnless(0 == m.__cmp__(m),
899                        'test_exportASCIIsegmentoutlinefile:loading and saving of a mesh failed')
900        # Having problems with this on linux.
901        #The ordering of the dictionary values wasn't the same as the windows
902        #returned value (verts.values())
903        #self.failUnless(0 == m.__cmp__(m_returned),
904        #                'test_exportASCIIsegmentoutlinefile:loading and saving of a mesh failed')
905       
906        self.failUnless(3 == len(m_returned.userVertices),
907                        'segmentoutlinefile:IO of a mesh failed')
908        self.failUnless(len(m.userSegments) == len(m_returned.userSegments),
909                        'segmentoutlinefile:IO of a mesh failed')
910        for i in range(len(m.userSegments)):
911            self.failUnless(m.userSegments[i].vertices[0].x ==
912                            m_returned.userSegments[i].vertices[0].x,
913                        'loading and saving of a mesh outline fialed')
914            self.failUnless(m.userSegments[i].vertices[0].y ==
915                            m_returned.userSegments[i].vertices[0].y,
916                        'loading and saving of a mesh outline fialed')
917            self.failUnless(m.userSegments[i].vertices[1].x ==
918                            m_returned.userSegments[i].vertices[1].x,
919                        'loading and saving of a mesh outline fialed')
920            self.failUnless(m.userSegments[i].vertices[1].y ==
921                            m_returned.userSegments[i].vertices[1].y,
922                        'loading and saving of a mesh outline fialed')
923
924 
925    def test_exportASCIIsegmentoutlinefile2(self):
926        a = Vertex (0,0)
927        b = Vertex (0,1)
928        c = Vertex (1,0)
929        d = Vertex (1,1)
930        e = Vertex (0.5,0.5)
931        f  = Vertex (0.6,0.6)
932     
933        s1 = Segment(a,e, tag = "50")
934        s2 = Segment(b,e, tag = "40")
935        s3 = Segment(c,e, tag = "30")
936        s4 = Segment(d,e, tag = "30")
937     
938        r1 = Region(2, 1,tag = "1.3")
939        h1 = Hole(1,4)
940        m = Mesh(userVertices=[a,b,c,d,e],
941                 userSegments=[s1,s2,s3,s4],
942                 regions=[r1],
943                 holes = [h1])     
944       
945        fileName = tempfile.mktemp(".tsh")
946        m.exportASCIIsegmentoutlinefile(fileName)
947       
948        m_returned = importMeshFromFile(fileName)
949        #print "****************** fileName", fileName
950        os.remove(fileName)
951
952        #Trim mesh, so it should look like m_returned
953        m.meshVertices = []
954        m.meshTriangles = []
955        m.meshSegments = []
956        m.userVertices=[a,e,d,b,c]
957        #print "mesh ***************dsg*", m
958        #print "(m.__cmp__(m_returned)", m.__cmp__(m_returned)
959        self.failUnless(0 == m.__cmp__(m),
960                        'loading and saving of a mesh failed')
961
962        self.failUnless(5 == len(m_returned.userVertices),
963                        'segmentoutlinefile:IO of a mesh failed')
964        self.failUnless(len(m.userSegments) == len(m_returned.userSegments),
965                        'segmentoutlinefile:IO of a mesh failed')
966        for i in range(len(m.userSegments)):
967            self.failUnless(m.userSegments[i].vertices[0].x ==
968                            m_returned.userSegments[i].vertices[0].x,
969                        'loading and saving of a mesh outline fialed')
970            self.failUnless(m.userSegments[i].vertices[0].y ==
971                            m_returned.userSegments[i].vertices[0].y,
972                        'loading and saving of a mesh outline fialed')
973            self.failUnless(m.userSegments[i].vertices[1].x ==
974                            m_returned.userSegments[i].vertices[1].x,
975                        'loading and saving of a mesh outline fialed')
976            self.failUnless(m.userSegments[i].vertices[1].y ==
977                            m_returned.userSegments[i].vertices[1].y,
978                        'loading and saving of a mesh outline fialed')
979
980
981    def test_loadxy(self):
982        """
983        To test the mesh side of loading xya files.
984        Not the loading of xya files
985        """
986        import os
987        import tempfile
988       
989        fileName = tempfile.mktemp(".csv")
990        file = open(fileName,"w")
991        file.write("x,y,elevation, speed \n\
9921.0, 0.0, 10.0, 0.0\n\
9930.0, 1.0, 0.0, 10.0\n\
9941.0, 0.0, 10.4, 40.0\n")
995        file.close()
996        #print fileName
997        m = importMeshFromFile(fileName)
998        os.remove(fileName)
999        self.failUnless(m.userVertices[0].x == 1.0,
1000                        'loadxy, test 1 failed')
1001        self.failUnless(m.userVertices[0].y == 0.0,
1002                        'loadxy, test 2 failed')
1003        #self.failUnless(m.userVertices[0].attributes == [10.0,0.0],
1004        #                'loadxy, test 2.2 failed')
1005        self.failUnless(m.userVertices[1].x == 0.0,
1006                        'loadxy, test 3 failed')
1007        self.failUnless(m.userVertices[1].y == 1.0,
1008                        'loadxy, test 4 failed')
1009        #self.failUnless(m.userVertices[1].attributes == [0.0,10.0],
1010        #                'loadxy, test 5 failed')
1011       
1012    def exportPointsFile(self):
1013        a = Vertex (0,0)
1014        b = Vertex (0,3)
1015        c = Vertex (3,3)
1016        d = Vertex (1,2)
1017        e = Vertex (3,1)
1018        #f = Vertex (3,1)
1019     
1020        s1 = Segment(a,b, tag = 50)
1021        s2 = Segment(b,c, tag = 40)
1022        s3 = Segment(c,a, tag = 30)
1023     
1024        r1 = Region(2, 1,tag = 1.3)
1025        h1 = Hole(1,4)
1026        # Warning mesh can't produce this type of data structure its self
1027        m = Mesh(userVertices=[a,b,c,d,e],
1028                 userSegments=[s1,s2,s3],
1029                 regions=[r1],
1030                 holes = [h1])
1031       
1032        fileName = tempfile.mktemp(".txt")
1033        #fileName = 't.xya'
1034        #os.remove(fileName)
1035        m.exportPointsFile(fileName)
1036        file = open(fileName)
1037        lFile = file.read().split('\n')
1038        file.close()
1039
1040        os.remove(fileName)
1041        self.failUnless(lFile[0] == "x,y" and
1042                        lFile[1] == "0,0" and
1043                        lFile[2] == "0,3" and
1044                        lFile[3] == "3,3" 
1045                        ,
1046                        'exported Ascii xya file is wrong')
1047        self.failUnless(lFile[4] == "1,2" and
1048                        lFile[5] == "3,1" 
1049                        ,
1050                        'exported Ascii xya file is wrong')
1051       
1052        # vertex e is outside of the outline, so
1053        # it is a loner and it is removed.
1054        m.generateMesh("Q", maxArea = 2.1)
1055        fileName = tempfile.mktemp(".xya")
1056        #fileName = 't.xya'
1057        #m.export_mesh_file('m.tsh')
1058        m.exportPointsFile(fileName)
1059        file = open(fileName)
1060        lFile = file.read().split('\n')
1061        file.close()
1062        os.remove(fileName)
1063       
1064        self.failUnless(lFile[0] == "" and
1065                        lFile[1] == "0.0,0.0" and
1066                        lFile[2] == "0.0,3.0" and
1067                        lFile[3] == "3.0,3.0" and
1068                        lFile[4] == "1.0,2.0"
1069                        ,
1070                        'exported Ascii xya file is wrong')
1071     
1072    def to_be_test_lone_vert_in_mesh_gen_c_layer(self):
1073        # currently just a copy of the above test
1074        a = Vertex (0,0)
1075        b = Vertex (0,3)
1076        c = Vertex (3,3)
1077        d = Vertex (1,2)
1078        e = Vertex (3,1)
1079        #f = Vertex (3,1)
1080     
1081        s1 = Segment(a,b, tag = 50)
1082        s2 = Segment(b,c, tag = 40)
1083        s3 = Segment(c,a, tag = 30)
1084     
1085        r1 = Region(2, 1,tag = 1.3)
1086        h1 = Hole(1,4)
1087        # Warning mesh can't produce this type of data structure its self
1088        m = Mesh(userVertices=[a,b,c,d,e],
1089                 userSegments=[s1,s2,s3],
1090                 regions=[r1],
1091                 holes = [h1])
1092       
1093        fileName = tempfile.mktemp(".csv")
1094        #fileName = 't.xya'
1095        #os.remove(fileName)
1096        m.exportPointsFile(fileName)
1097        file = open(fileName)
1098        lFile = file.read().split('\n')
1099        file.close()
1100
1101        os.remove(fileName)
1102        self.failUnless(lFile[0] == "x,y" and
1103                        lFile[1] == "0,0" and
1104                        lFile[2] == "0,3" and
1105                        lFile[3] == "3,3" 
1106                        ,
1107                        'exported Ascii xya file is wrong')
1108        self.failUnless(lFile[4] == "1,2" and
1109                        lFile[5] == "3,1" 
1110                        ,
1111                        'exported Ascii xya file is wrong')
1112       
1113        # vertex e is outside of the outline, so
1114        # it is a loner and it is removed.
1115        m.generateMesh("Q", maxArea = 2.1)
1116        fileName = tempfile.mktemp(".csv")
1117        #fileName = 't.xya'
1118        #m.export_mesh_file('m.tsh')
1119        m.exportPointsFile(fileName)
1120        file = open(fileName)
1121        lFile = file.read().split('\n')
1122        file.close()
1123        os.remove(fileName)
1124       
1125        self.failUnless(lFile[0] == "x,y" and
1126                        lFile[1] == "0.0,0.0" and
1127                        lFile[2] == "0.0,3.0" and
1128                        lFile[3] == "3.0,3.0" and
1129                        lFile[4] == "1.0,2.0"
1130                        ,
1131                        'exported Ascii xya file is wrong')
1132       
1133    def NOT_test_exportPointsFilefile2(self):
1134        #geospatial needs at least one point
1135        m = Mesh()
1136       
1137        fileName = tempfile.mktemp(".csv")
1138        m.exportPointsFile(fileName)
1139        file = open(fileName)
1140        lFile = file.read().split('\n')
1141        file.close()
1142
1143        os.remove(fileName)
1144        #print "************* test_mesh exportPointsFilefile"
1145        #print "lFile",lFile
1146        #print "************* test_mesh exportPointsFilefile"
1147        self.failUnless(lFile[0] == "" 
1148                        ,
1149                        'exported Ascii xya file is wrong')
1150       
1151    def test_strings2ints(self):
1152        list = ["sea","river inlet","","sea","","moat"]
1153        preset = ["moat", "internal boundary"]
1154        [intlist, converter] = segment_strings2ints(list,preset )
1155        self.failUnless(intlist == [2,3 ,0 ,2 ,0 ,0 ]
1156                        ,
1157                        'test_strings2ints produces bad intlist')
1158        self.failUnless(converter == ['moat', 'internal boundary',
1159                                      'sea', 'river inlet']
1160                        ,
1161                        'test_strings2ints produces bad converter')
1162       
1163    def test_ints2strings(self):
1164        list = ["internal boundary","sea","river inlet",
1165            "","sea","","moat","internal boundary"]
1166        outlist = ['internal boundary', 'sea', 'river inlet', 'moat',
1167                   'sea', 'moat', 'moat', 'internal boundary']
1168        preset = ["moat", "internal boundary"]
1169        [intlist, converter] = segment_strings2ints(list,preset )
1170        newlist = segment_ints2strings(intlist, converter)
1171        self.failUnless(outlist == newlist
1172                        ,
1173                        'test_strings2ints produces bad intlist')
1174        self.failUnless(converter == ['moat', 'internal boundary',
1175                                      'sea', 'river inlet']
1176                        ,
1177                        'test_strings2ints produces bad converter')
1178       
1179    def test_ints2strings2(self):
1180        list = ["","",""]
1181        preset = ["moat", "internal boundary"]
1182        [intlist, converter] = segment_strings2ints(list,preset )
1183        newlist = segment_ints2strings(intlist, converter)
1184        outlist = ['moat', 'moat', 'moat']
1185        self.failUnless(outlist == newlist
1186                        ,
1187                        'test_strings2ints produces bad intlist')
1188        self.failUnless(converter == ['moat', 'internal boundary']
1189                        ,
1190                        'test_strings2ints produces bad converter')
1191
1192       
1193    def test_removeDuplicatedVertices(self):
1194        a = Vertex (0,0)
1195        a.index = 0
1196        b = Vertex (0,3)
1197        b.index = 1
1198        c = Vertex (3,3)
1199        c.index = 2
1200        d = Vertex (1,1)
1201        d.index = 3
1202        e = Vertex (3,1)
1203        e.index = 4
1204        f = Vertex (1,1)
1205        f.index = 5
1206        g = Vertex (1,1)
1207        g.index = 6
1208        inputVerts_noDups = [a,b,c,d,e]
1209       
1210        m = Mesh(userVertices=[a,b,c,d,e,f,g])
1211        counter = m.removeDuplicatedUserVertices()
1212        UserVerts = m.getUserVertices()
1213       
1214         
1215        self.failUnless(UserVerts == inputVerts_noDups,
1216                            'duplicate verts not removed')
1217        #for userVert, inputVert in map(None, UserVerts, inputVerts_noDups):
1218        #    self.failUnless(userVert.x == inputVert.x,
1219        #                    'x duplicate verts not removed')
1220        #    self.failUnless(userVert.y == inputVert.y,
1221        #                    'y duplicate verts not removed')
1222
1223       
1224    def test_ungenerateFileLoading(self):
1225       
1226        fileName = tempfile.mktemp(".txt")
1227        file = open(fileName,"w")
1228        file.write("         1       ??      ??\n\
1229       0.0       0.0\n\
1230       1.0       0.0\n\
1231       1.0       1.0\n\
1232       0.0       1.0\n\
1233       0.0       0.0\n\
1234END\n\
1235         2      ?? ??\n\
1236       10.0       10.0\n\
1237       10.0       20.0\n\
1238       20.0       20.0\n\
1239       10.0       10.0\n\
1240END\n\
1241END\n")
1242        file.close()
1243       
1244       
1245        a = Vertex (0.0, 0.0) #, attributes = [1.1])
1246        b = Vertex (0.0, 40.0) #, attributes = [1.2])
1247        c = Vertex (40.0,40.0) #, attributes = [1.3])
1248        d = Vertex (40.0,0.0) #, attributes = [1.4])
1249   
1250        s1 = Segment(a,b)
1251        s2 = Segment(b,c)
1252        s3 = Segment(c,d)
1253        s4 = Segment(d,a)
1254     
1255        m = Mesh(userVertices=[a,b,c,d], userSegments=[s1,s2,s3,s4])
1256        dict = importUngenerateFile(fileName)
1257        #os.remove(fileName)
1258
1259        tag = "DSG"
1260        Segment.set_default_tag(tag)
1261        m.addVertsSegs(dict)
1262
1263        # have to reset this , since it's a class attribute
1264        Segment.set_default_tag("")
1265           
1266        self.failUnless(len(m.userSegments) ==11,
1267                        'Wrong segment list length.')
1268        self.failUnless(len(m.userVertices) == 11,
1269                        'Wrong vertex list length.')
1270        self.failUnless(m.userSegments[10].vertices[0] == m.userVertices[10],
1271                        'bad vertex on segment.')
1272        self.failUnless(m.userSegments[10].vertices[1] == m.userVertices[8],
1273                        'Bad segment.')
1274        self.failUnless(m.userSegments[10].tag == tag,
1275                        'wrong tag.')
1276
1277        ## let's test the method
1278        a = Vertex (0.0, 0.0) #, attributes = [1.1])
1279        b = Vertex (0.0, 40.0) #, attributes = [1.2])
1280        c = Vertex (40.0,40.0) #, attributes = [1.3])
1281        d = Vertex (40.0,0.0) #, attributes = [1.4])
1282   
1283        s1 = Segment(a,b)
1284        s2 = Segment(b,c)
1285        s3 = Segment(c,d)
1286        s4 = Segment(d,a)
1287     
1288        m = Mesh(userVertices=[a,b,c,d], userSegments=[s1,s2,s3,s4])
1289
1290        tag = "DSG"       
1291        initial_tag = "PIG"
1292        Segment.set_default_tag(initial_tag)
1293        m.import_ungenerate_file(fileName, tag=tag)
1294
1295        os.remove(fileName)
1296
1297        self.failUnless(Segment.get_default_tag() == initial_tag,
1298                        'Wrong segment list length.')
1299       
1300
1301        # have to reset this , since it's a class attribute
1302        Segment.set_default_tag("")
1303           
1304        self.failUnless(len(m.userSegments) ==11,
1305                        'Wrong segment list length.')
1306        self.failUnless(len(m.userVertices) == 11,
1307                        'Wrong vertex list length.')
1308        self.failUnless(m.userSegments[10].vertices[0] == m.userVertices[10],
1309                        'bad vertex on segment.')
1310        self.failUnless(m.userSegments[10].vertices[1] == m.userVertices[8],
1311                        'Bad segment.')
1312        self.failUnless(m.userSegments[10].tag == tag,
1313                        'wrong tag.')
1314       
1315    def test_ungenerateFileLoadingII(self):
1316       
1317        fileName = tempfile.mktemp(".txt")
1318        file = open(fileName,"w")
1319        file.write("         1       ??      ??\n\
1320       0.0       0.0\n\
1321       1.0       0.0\n\
1322       1.0       1.0\n\
1323       0.0       1.0\n\
1324       0.0       0.0\n\
1325END\n\
1326         2      ?? ??\n\
1327       10.0       10.0\n\
1328       10.0       20.0\n\
1329       20.0       20.0\n\
1330END\n\
1331END\n")
1332        file.close()
1333       
1334       
1335        a = Vertex (0.0, 0.0) #, attributes = [1.1])
1336        b = Vertex (0.0, 40.0) #, attributes = [1.2])
1337        c = Vertex (40.0,40.0) #, attributes = [1.3])
1338        d = Vertex (40.0,0.0) #, attributes = [1.4])
1339   
1340        s1 = Segment(a,b)
1341        s2 = Segment(b,c)
1342        s3 = Segment(c,d)
1343        s4 = Segment(d,a)
1344     
1345        m = Mesh(userVertices=[a,b,c,d], userSegments=[s1,s2,s3,s4])
1346        dict = importUngenerateFile(fileName)
1347        #os.remove(fileName)
1348
1349        tag = "DSG"
1350        Segment.set_default_tag(tag)
1351        m.addVertsSegs(dict)
1352
1353        self.failUnless(len(m.userSegments) ==10,
1354                        'Wrong segment list length.')
1355        self.failUnless(len(m.userVertices) == 11,
1356                        'Wrong vertex list length.')
1357
1358        # Test the method
1359        a = Vertex (0.0, 0.0) #, attributes = [1.1])
1360        b = Vertex (0.0, 40.0) #, attributes = [1.2])
1361        c = Vertex (40.0,40.0) #, attributes = [1.3])
1362        d = Vertex (40.0,0.0) #, attributes = [1.4])
1363   
1364        s1 = Segment(a,b)
1365        s2 = Segment(b,c)
1366        s3 = Segment(c,d)
1367        s4 = Segment(d,a)
1368     
1369        m = Mesh(userVertices=[a,b,c,d], userSegments=[s1,s2,s3,s4])
1370        tag = "DSG"       
1371        initial_tag = "PIG"
1372        Segment.set_default_tag(initial_tag)
1373        m.import_ungenerate_file(fileName, tag=tag)
1374
1375        os.remove(fileName)
1376
1377        self.failUnless(Segment.get_default_tag() == initial_tag,
1378                        'Wrong segment list length.')
1379       
1380       
1381        self.failUnless(len(m.userSegments) ==10,
1382                        'Wrong segment list length.')
1383        self.failUnless(len(m.userVertices) == 11,
1384                        'Wrong vertex list length.')
1385       
1386        # have to reset this , since it's a class attribute
1387        Segment.set_default_tag("")
1388       
1389    def test_addVertsSegs(self):
1390        m = Mesh()
1391        Segment.set_default_tag("food")
1392        dict = {}
1393        dict['points'] = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]]
1394        dict['segments'] = [[0, 1], [1, 2]]
1395        dict['segment_tags'] = ['','do-op']
1396        m.addVertsSegs(dict)
1397        # have to reset this , since it's a class attribute
1398        Segment.set_default_tag("")
1399
1400       
1401        self.failUnless(len(m.userSegments) ==2,
1402                        'Wrong segment list length.')
1403        self.failUnless(len(m.userVertices) == 3,
1404                        'Wrong vertex list length.')
1405        self.failUnless(m.userSegments[0].tag =='food',
1406                        'Wrong segment tag length.')
1407        self.failUnless(m.userSegments[1].tag =='do-op',
1408                        'Wrong segment tag.')
1409       
1410    def test_addVertsSegs2(self):
1411        geo = Geo_reference(56,5,10)
1412        m = Mesh(geo_reference=geo)
1413        dict = {}
1414        dict['points'] = [[2.0, 1.0], [3.0, 1.0], [2.0, 2.0]]
1415        dict['segments'] = [[0, 1], [1, 2], [2,0]]
1416        dict['segment_tags'] = ['','do-op','']
1417        m.addVertsSegs(dict)
1418
1419    def test_addVertsSegs_done_twice(self):
1420        m = Mesh()
1421        dict = {}
1422        dict['points'] = [[0.0, 0.0], [5.0, 0.0], [5.0, 5.0]]
1423        dict['segments'] = [[0, 1], [1, 2], [2,0]]
1424        dict['segment_tags'] = ['0','1','2']
1425        m.addVertsSegs(dict)
1426       
1427        dict = {}
1428        dict['points'] = [[2.0, 1.0], [4.0, 1.0], [4.0, 3.0]]
1429        dict['segments'] = [[0, 1], [1, 2], [2,0]]
1430        dict['segment_tags'] = ['3','4','5']
1431        m.addVertsSegs(dict)
1432
1433       
1434        self.failUnless(m.userSegments[5].vertices[0].y == 3,
1435                        'Wrong vertex connected.')
1436        self.failUnless(m.userSegments[5].vertices[1].y == 1,
1437                        'Wrong vertex connected.')
1438           
1439    def test_add_points_and_segments(self):
1440        m = Mesh()
1441        Segment.set_default_tag("food")
1442        dict = {}
1443        points =  [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]]
1444        segments = [[0, 1], [1, 2]]
1445        segment_tags = {'do-op':[1]}
1446        m.add_points_and_segments(points,
1447                                    segments, segment_tags)
1448        # have to reset this , since it's a class attribute
1449        Segment.set_default_tag("")
1450
1451       
1452        self.failUnless(len(m.userSegments) ==2,
1453                        'Wrong segment list length.')
1454        self.failUnless(len(m.userVertices) == 3,
1455                        'Wrong vertex list length.')
1456        self.failUnless(m.userSegments[0].tag =='food',
1457                        'Wrong segment tag length.')
1458        self.failUnless(m.userSegments[1].tag =='do-op',
1459                        'Wrong segment tag.')
1460       
1461    def test_exportASCIImeshfile(self):
1462   
1463        #a_att = [5,2]
1464        #d_att =[4,2]
1465        #f_att = [3,2]
1466        #e_att = [2,2]
1467        a_xy = [0.0, 0.0]
1468        a = Vertex ( a_xy[0],a_xy[1]) #, attributes =a_att)
1469        d = Vertex (0.0, 4.0) #, attributes =d_att)
1470        f = Vertex (4.0,0.0) #, attributes =f_att)
1471        e = Vertex (1.0,1.0) #, attributes =e_att)
1472   
1473        s1 = Segment(a,d, tag = "50")
1474        s2 = Segment(d,f, tag = "40")
1475        s3 = Segment(a,f, tag = "30")
1476        s4 = Segment(a,e, tag = "20")
1477     
1478        r1 = Region(0.3, 0.3,tag = "1.3", maxArea = 36)
1479
1480
1481        h1 = Hole(0.2,0.6)
1482       
1483        m = Mesh(userVertices=[a,d,f,e],
1484                 userSegments=[s1,s2,s3,s4],
1485                 regions=[r1],
1486                 holes=[h1])
1487
1488        seg = m.getUserSegments()
1489        points = m.getUserVertices()
1490        holes = m.getHoles()
1491        regions = m.getRegions()
1492        fileName = tempfile.mktemp(".tsh")
1493        m.export_mesh_file(fileName)
1494        #print "***************************fileName", fileName
1495        new_m = importMeshFromFile(fileName)
1496        os.remove(fileName)
1497       
1498
1499        #print '**@@@@@******'
1500        #print "new_m",new_m
1501        #print '**@@@@@******'
1502        #print "m",m
1503        #print '**@@@@@******'
1504       
1505        self.failUnless( new_m == m,
1506                         'loadASCIITestCase failed. test new 1')
1507           
1508    def test_Mesh2MeshList(self):
1509
1510        a_att = [5,2]
1511        d_att =[4,2]
1512        f_att = [3,2]
1513        e_att = [2,2]
1514        a_xy = [0.0, 0.0]
1515        a = Vertex ( a_xy[0],a_xy[1]) #, attributes =a_att)
1516        d = Vertex (0.0, 4.0) #, attributes =d_att)
1517        f = Vertex (4.0,0.0) #, attributes =f_att)
1518        e = Vertex (1.0,1.0) #, attributes =e_att)
1519   
1520        s1 = Segment(a,d, tag = "50")
1521        s2 = Segment(d,f, tag = "40")
1522        s3 = Segment(a,f, tag = "30")
1523        s4 = Segment(a,e, tag = "20")
1524     
1525        r1 = Region(0.3, 0.3,tag = "1.3", maxArea = 45)
1526        m = Mesh(userVertices=[a,d,f,e],
1527                 userSegments=[s1,s2,s3,s4],
1528                 regions=[r1])
1529
1530        m.generateMesh("Qa2.1")
1531
1532        seg = m.getMeshSegments()
1533        points = m.getMeshVertices()
1534        dict = m.Mesh2MeshList()
1535        #print "dict",dict
1536        # test not finished...
1537 
1538    def test_Mesh2IOTriangulationDict(self):
1539
1540        a_att = [5,2]
1541        d_att =[4,2]
1542        f_att = [3,2]
1543        e_att = [2,2]
1544        a_xy = [0.0, 0.0]
1545        a = Vertex ( a_xy[0],a_xy[1] , attributes =a_att)
1546        d = Vertex (0.0, 4.0 , attributes =d_att)
1547        f = Vertex (4.0,0.0 , attributes =f_att)
1548        e = Vertex (1.0,1.0 , attributes =e_att)
1549   
1550        s1 = Segment(a,d, tag = "50")
1551        s2 = Segment(d,f, tag = "40")
1552        s3 = Segment(a,f, tag = "30")
1553        s4 = Segment(a,e, tag = "20")
1554     
1555        r1 = Region(0.3, 0.3,tag = "1.3", maxArea = 45)
1556        m = Mesh(userVertices=[a,d,f,e],
1557                 userSegments=[s1,s2,s3,s4],
1558                 regions=[r1])
1559        titles = ['ele','friction']
1560        m.attributeTitles = titles
1561        m.generateMesh("Qa2.1")
1562
1563        seg = m.getMeshSegments()
1564        verts = m.getMeshVertices()
1565        dict = m.Mesh2IOTriangulationDict()
1566        #print "dict",dict
1567       
1568        self.failUnless( dict['vertex_attribute_titles'] == titles,
1569                         'test_Mesh2IOTriangulationDict failed. test 1')
1570        answer = [a_xy,[0.0, 4.0],[4.0,0.0], [1.0,1.0], [2.0,2.0]]
1571        #print "answer",answer
1572        #print "dict['vertices']",dict['vertices']
1573       
1574        self.failUnless( dict['vertices'] == answer,
1575                         'test_Mesh2IOTriangulationDict failed. test 2')
1576
1577       
1578        for pimport,pactual,pimpatt in map(None,dict['vertices'],
1579                                           verts,dict['vertex_attributes']):
1580            #assert all_close( pimport, (pactual.x,pactual.y))
1581            self.failUnless( pimport == [pactual.x,pactual.y],
1582                        'test_Mesh2IOTriangulationDict failed. test 2.1')
1583            self.failUnless( pimpatt == pactual.attributes,
1584                        'test_Mesh2IOTriangulationDict failed. test 2.2')
1585        self.failUnless( dict['segments'][0] == [0,1],
1586                        'test_Mesh2IOTriangulationDict failed. test 3')
1587        for segimp,segactual in map(None,dict['segment_tags'],seg):
1588            self.failUnless( segimp == segactual.tag,
1589                        'test_Mesh2IOTriangulationDict failed. test 4')
1590        #print " dict['triangles'][0]", dict['triangles'][0]
1591        self.failUnless( dict['triangles'][0] == [3,2,4],
1592                        'test_Mesh2IOTriangulationDict failed. test 5')
1593        self.failUnless( dict['triangle_neighbors'][0] == [-1,2,3],
1594                        'test_Mesh2IOTriangulationDict failed. test 6')
1595        self.failUnless( dict['triangle_tags'][0] == "1.3",
1596                         'test_Mesh2IOTriangulationDict failed. test 7')
1597
1598 
1599    def test_Mesh2IODict(self):
1600
1601        a_att = [5,2]
1602        d_att =[4,2]
1603        f_att = [3,2]
1604        e_att = [2,2]
1605        a_xy = [0.0, 0.0]
1606        a = Vertex ( a_xy[0],a_xy[1] , attributes =a_att)
1607        d = Vertex (0.0, 4.0 , attributes =d_att)
1608        f = Vertex (4.0,0.0 , attributes =f_att)
1609        e = Vertex (1.0,1.0 , attributes =e_att)
1610   
1611        s1 = Segment(a,d, tag = "50")
1612        s2 = Segment(d,f, tag = "40")
1613        s3 = Segment(a,f, tag = "30")
1614        s4 = Segment(a,e, tag = "20")
1615     
1616        r1 = Region(0.3, 0.3,tag = "1.3", maxArea = 45)
1617        m = Mesh(userVertices=[a,d,f,e],
1618                 userSegments=[s1,s2,s3,s4],
1619                 regions=[r1])
1620        titles = ['ele','friction']
1621        m.attributeTitles = titles
1622        m.generateMesh("Qa2.1")
1623
1624        seg = m.getMeshSegments()
1625        verts = m.getMeshVertices()
1626        dict = m.Mesh2IODict()
1627        #print "dict",dict
1628       
1629        self.failUnless( dict['vertex_attribute_titles'] == titles,
1630                         'test_Mesh2IOTriangulationDict failed. test 1')
1631        answer = [a_xy,[0.0, 4.0],[4.0,0.0], [1.0,1.0], [2.0,2.0]]
1632        #print "answer",answer
1633        #print "dict['vertices']",dict['vertices']
1634       
1635        self.failUnless( dict['vertices'] == answer,
1636                         'test_Mesh2IOTriangulationDict failed. test 2')
1637
1638       
1639        for pimport,pactual,pimpatt in map(None,dict['vertices'],
1640                                           verts,dict['vertex_attributes']):
1641            #assert all_close( pimport, (pactual.x,pactual.y))
1642            self.failUnless( pimport == [pactual.x,pactual.y],
1643                        'test_Mesh2IODict failed. test 2.1')
1644            self.failUnless( pimpatt == pactual.attributes,
1645                        'test_Mesh2IODict failed. test 2.2')
1646        self.failUnless( dict['segments'][0] == [0,1],
1647                        'test_Mesh2IODict failed. test 3')
1648        for segimp,segactual in map(None,dict['segment_tags'],seg):
1649            self.failUnless( segimp == segactual.tag,
1650                        'test_Mesh2IODict failed. test 4')
1651        #print " dict['triangles'][0]", dict['triangles'][0]
1652        self.failUnless( dict['triangles'][0] == [3,2,4],
1653                        'test_Mesh2IODict failed. test 5')
1654        self.failUnless( dict['triangle_neighbors'][0] == [-1,2,3],
1655                        'test_Mesh2IODict failed. test 6')
1656        self.failUnless( dict['triangle_tags'][0] == "1.3",
1657                         'test_Mesh2IODict failed. test 7')
1658
1659        seg = m.getUserSegments()
1660        points = m.getUserVertices()
1661        holes = m.getHoles()
1662        regions = m.getRegions()
1663       
1664        for pimport,pactual,pimpatt in map(None,dict['points'],points,dict['point_attributes']):
1665            self.failUnless( pimport == [pactual.x,pactual.y],
1666                        'test_Mesh2IODict failed. test 1')
1667            self.failUnless( pimpatt == pactual.attributes,
1668                        'test_Mesh2IODict failed. test 1.1')
1669        self.failUnless( dict['outline_segments'][0] == [0,1],
1670                        'test_Mesh2IODict failed. test 3')
1671        for segimp,segactual in map(None,dict['outline_segment_tags'],seg):
1672            self.failUnless( segimp == segactual.tag,
1673                        'test_Mesh2IODict failed. test 4')
1674        for holeimp,holeactual in map(None,dict['holes'],holes):
1675            self.failUnless( holeimp == [holeactual.x,holeactual.y],
1676                        'test_Mesh2IODict failed. test 5')
1677       
1678        for regimp,regactual,regattimp, regmaxarea in map(None,dict['regions'],regions, dict['region_tags'], dict['region_max_areas']):
1679            self.failUnless( regimp == [regactual.x,regactual.y],
1680                        'loadASCIITestCase failed. test 6')
1681            self.failUnless( regattimp == regactual.getTag(),
1682                        'loadASCIITestCase failed. test 7')
1683            self.failUnless( regmaxarea == regactual.getMaxArea(),
1684                        'loadASCIITestCase failed. test 7')
1685   
1686           
1687       
1688    def test_Mesh2IOOutlineDict(self):
1689
1690        a_att = [5,2]
1691        d_att =[4,2]
1692        f_att = [3,2]
1693        e_att = [2,2]
1694        a_xy = [0.0, 0.0]
1695        a = Vertex ( a_xy[0],a_xy[1] , attributes =a_att)
1696        d = Vertex (0.0, 4.0 , attributes =d_att)
1697        f = Vertex (4.0,0.0 , attributes =f_att)
1698        e = Vertex (1.0,1.0 , attributes =e_att)
1699   
1700        s1 = Segment(a,d, tag = "50")
1701        s2 = Segment(d,f, tag = "40")
1702        s3 = Segment(a,f, tag = "30")
1703        s4 = Segment(a,e, tag = "20")
1704     
1705        r1 = Region(0.3, 0.3,tag = "1.3", maxArea = 45)
1706        m = Mesh(userVertices=[a,d,f,e],
1707                 userSegments=[s1,s2,s3,s4],
1708                 regions=[r1])
1709        titles = ['ele','friction']
1710        m.attributeTitles = titles
1711        m.generateMesh("Qa2.1")
1712
1713        seg = m.getMeshSegments()
1714        verts = m.getMeshVertices()
1715        dict = m.Mesh2IOOutlineDict()
1716       
1717        seg = m.getUserSegments()
1718        points = m.getUserVertices()
1719        holes = m.getHoles()
1720        regions = m.getRegions()
1721       
1722        for pimport,pactual,pimpatt in map(None,dict['points'],points,dict['point_attributes']):
1723            self.failUnless( pimport == [pactual.x,pactual.y],
1724                        'loadASCIITestCase failed. test 1')
1725            self.failUnless( pimpatt == pactual.attributes,
1726                        'loadASCIITestCase failed. test 1.1')
1727        self.failUnless( dict['outline_segments'][0] == [0,1],
1728                        'loadASCIITestCase failed. test 3')
1729        for segimp,segactual in map(None,dict['outline_segment_tags'],seg):
1730            self.failUnless( segimp == segactual.tag,
1731                        'loadASCIITestCase failed. test 4')
1732        for holeimp,holeactual in map(None,dict['holes'],holes):
1733            self.failUnless( holeimp == [holeactual.x,holeactual.y],
1734                        'loadASCIITestCase failed. test 5')
1735        #for regimp,regactual in map(None,dict['regions'],regions):
1736         #   self.failUnless( [regimp[0],regimp[1]]==[regactual.x,regactual.y],
1737          #              'loadASCIITestCase failed. test 6')
1738           # self.failUnless( regimp[2] == regactual.getTag(),
1739            #            'loadASCIITestCase failed. test 7')
1740            #self.failUnless( regimp[3] == regactual.getMaxArea(),
1741             #           'loadASCIITestCase failed. test 7')
1742
1743           
1744        for regimp,regactual,regattimp, regmaxarea in map(None,dict['regions'],regions, dict['region_tags'], dict['region_max_areas']):
1745            self.failUnless( regimp == [regactual.x,regactual.y],
1746                        'loadASCIITestCase failed. test 6')
1747            self.failUnless( regattimp == regactual.getTag(),
1748                        'loadASCIITestCase failed. test 7')
1749            self.failUnless( regmaxarea == regactual.getMaxArea(),
1750                        'loadASCIITestCase failed. test 7')
1751
1752
1753#___________beginning of Peters tests
1754
1755    def test_set_stuff(self):
1756        """
1757        Documentation
1758        """
1759        #making a test mesh
1760        p0=[0.,2.]
1761        p1=[1.,2.]
1762        p2=[0.,1.]
1763        p3=[1.,1.]
1764        p4=[0.,0.]
1765        p5=[2.,0.]
1766        p6=[-1.,1.]
1767        point_list = [p0,p1,p2,p3,p4,p5,p6]
1768
1769        a0=[0]
1770        a1=[0]
1771        a2=[100]
1772        a3=[0]
1773        a4=[0]
1774        a5=[0]
1775        a6=[0]
1776        attribute=[a0,a1,a2,a3,a4,a5,a6]
1777       
1778        t0=[0,3,1]
1779        t1=[0,2,3]
1780        t2=[2,4,3]
1781        t3=[4,5,3]
1782        t4=[1,3,5]
1783        t5=[2,6,4]
1784
1785        n0=[4,-1,2]
1786        n1=[2,0,-1]
1787        n2=[3,1,5]
1788        n3=[4,2,-1]
1789        n4=[3,-1,0]
1790        n5=[-1,2,-1]
1791
1792        tri_list = [t0,t1,t2,t3,t4,t5]
1793        n_list = [n0,n1,n2,n3,n4,n5]
1794        for i in range(6):
1795            for j in (0,1,2):
1796                a=attribute[tri_list[i][j]]
1797                tri_list[i][j]=point_list[tri_list[i][j]]
1798                tri_list[i][j]=Vertex(tri_list[i][j][0]\
1799                                      ,tri_list[i][j][1],a)
1800            neighbours=n_list[i]
1801            tri_list[i]=Triangle(tri_list[i][0],\
1802                                 tri_list[i][1],tri_list[i][2]\
1803                                ,neighbors=neighbours)
1804
1805        #testing selectAll
1806        mesh = Mesh()
1807        mesh.attributeTitles=['attrib']
1808
1809        mesh.meshTriangles=tri_list
1810
1811        mesh.selectAllTriangles()
1812        A=mesh.sets[mesh.setID['All']]
1813        assert list_comp(tri_list,A)
1814
1815       #testing threshold
1816        mesh = Mesh()
1817        mesh.attributeTitles=['attrib']
1818
1819        mesh.meshTriangles=tri_list
1820        mesh.selectAllTriangles()
1821        mesh.threshold('All',min=30,max=35,attribute_name = 'attrib')
1822        A = [tri_list[1],tri_list[2],tri_list[5]]
1823        B = mesh.sets[mesh.setID['All']]
1824        assert list_comp(A,B)
1825       
1826
1827        A = [tri_list[3],tri_list[2],tri_list[5]]
1828        assert not list_comp(A,B)
1829
1830        #testing
1831
1832    def test_Discretised_Tuple_Set_rounding(self):
1833        #This is the hardest bit of DST
1834
1835        tol = 0.1
1836        a=Discretised_Tuple_Set(p_rel=1,t_rel= tol)
1837        m = 0.541
1838        m_up = 0.6
1839        m_down = 0.5
1840        assert m_up == a.round_up_rel(m)
1841        assert m_down == a.round_down_rel(m)
1842
1843        tol = 0.1
1844        a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
1845        m = 0.539
1846        m_up = 0.5
1847        m_down = 0.5
1848        assert m_up == a.round_up_rel(m)
1849        assert m_down == a.round_down_rel(m)
1850
1851        tol = 0.5
1852        a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
1853
1854
1855        m = 0.6
1856        m_up = 0.7
1857        m_down = 0.5
1858        assert m_up == a.round_up_rel(m)
1859        assert m_down == a.round_down_rel(m)
1860
1861        m = 0.599
1862        m_up = 0.6
1863        m_down = 0.5
1864        assert m_up == a.round_up_rel(m)
1865        assert m_down == a.round_down_rel(m)
1866
1867    def test_Discretised_Tuple_Set_get(self):
1868       
1869        tol = 0.25
1870        a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
1871        b = (1.1,1.1)
1872        a.append(b)
1873        list = [(1.2,1.),(1.,1.),(1.,1.2),(1.2,1.2)]
1874        for key in list:
1875            assert a[key][0]==b
1876            assert len(a[key])==1
1877       
1878        c = (2.1,1.)
1879        a.append(c)
1880        assert a[(2.,1.)][0]==c
1881        assert a[(2.2,1.)][0]==c
1882
1883    def test_mapped_Discretised_Tuple_Set(self):
1884
1885        def map(sequence):
1886            return [len(sequence)]
1887
1888        tol = 0.5
1889        a=Mapped_Discretised_Tuple_Set(map,p_rel=1,t_rel = tol)
1890        b = range(20)
1891        a.append(b)
1892        assert b in a[range(17)] 
1893        assert b in a[range(22)]
1894
1895        tol = 0.01
1896        a=Mapped_Discretised_Tuple_Set(map,p_rel=1,t_rel = tol)
1897        b = range(20)
1898        a.append(b)
1899        assert b in a[range(20)] 
1900        assert b in a[range(19)] 
1901        assert not range(17) in a
1902
1903#___________end of Peters tests
1904
1905    def test_add_region_from_polygon(self):
1906        m=Mesh()
1907        region = m.add_region_from_polygon([[0,0],[1,0],[0,1]],
1908                                  max_triangle_area = 88)
1909        self.failUnless(len(m.regions)==1,
1910                        'FAILED!')
1911        self.failUnless(region.getMaxArea()==88,
1912                        'FAILED!')
1913        self.failUnless(len(m.getUserSegments())==3,
1914                        'FAILED!')
1915        self.failUnless(len(m.userVertices)==3,
1916                        'FAILED!')
1917       
1918    def test_add_region_from_polygon2(self):
1919        m=Mesh()
1920        m.add_region_from_polygon([[0,0],[1,0],[1,1],[0,1]],
1921                               {'tagin':[0,1],'bom':[2]},
1922                                  max_triangle_area=10)
1923        self.failUnless(len(m.regions)==1,
1924                        'FAILED!')
1925        segs = m.getUserSegments()
1926        self.failUnless(len(segs)==4,
1927                        'FAILED!')
1928        self.failUnless(len(m.userVertices)==4,
1929                        'FAILED!') 
1930        self.failUnless(segs[0].tag=='tagin',
1931                        'FAILED!') 
1932        self.failUnless(segs[1].tag=='tagin',
1933                        'FAILED!') 
1934         
1935        self.failUnless(segs[2].tag=='bom',
1936                        'FAILED!') 
1937        self.failUnless(segs[3].tag=='',
1938                        'FAILED!') 
1939       
1940    def test_add_region_from_polygon3(self):
1941        x=-500
1942        y=-1000
1943        m=Mesh(geo_reference=Geo_reference(56,x,y))
1944
1945        # These are the absolute values
1946        polygon_absolute = [[0,0],[1,0],[1,1],[0,1]]
1947       
1948        x_p = -10
1949        y_p = -40
1950        geo_ref_poly = Geo_reference(56, x_p, y_p)
1951        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
1952       
1953        poly_point = m.add_region_from_polygon(polygon,
1954                                               {'tagin':[0,1],'bom':[2]},
1955                                               geo_reference=geo_ref_poly,
1956                                               max_triangle_area=10)
1957        # poly_point values are relative to the mesh geo-ref
1958        # make them absolute
1959        self.failUnless(is_inside_polygon([poly_point.x+x,poly_point.y+y],
1960                                       polygon_absolute, closed = False),
1961                        'FAILED!')
1962               
1963        self.failUnless(len(m.regions)==1,
1964                        'FAILED!')
1965        segs = m.getUserSegments()
1966        self.failUnless(len(segs)==4,
1967                        'FAILED!')
1968        self.failUnless(len(m.userVertices)==4,
1969                        'FAILED!') 
1970        self.failUnless(segs[0].tag=='tagin',
1971                        'FAILED!') 
1972        self.failUnless(segs[1].tag=='tagin',
1973                        'FAILED!') 
1974         
1975        self.failUnless(segs[2].tag=='bom',
1976                        'FAILED!') 
1977        self.failUnless(segs[3].tag=='',
1978                        'FAILED!')
1979        verts = m.getUserVertices()
1980        #print "User verts",verts
1981        #print 'polygon',polygon
1982        #vert values are relative
1983        for point,new_point in map(None,polygon,verts):
1984            point_x = point[0] + geo_ref_poly.get_xllcorner()
1985            new_point_x = new_point.x + m.geo_reference.get_xllcorner()
1986            point_y = point[1] + geo_ref_poly.get_yllcorner()
1987            #print "new_point.y",new_point.y
1988            #print "m.geo_ref.get_yllcorner()",m.geo_reference.get_yllcorner()
1989            new_point_y = new_point.y + m.geo_reference.get_yllcorner()
1990            #print "point_y",point_y
1991            #print "new_point_y",new_point_y
1992           
1993            self.failUnless(point_x == new_point_x, ' failed')
1994            self.failUnless(point_y == new_point_y, ' failed')
1995           
1996         
1997    def test_add_region_from_polygon4(self):
1998        x=50000
1999        y=1000
2000        m=Mesh(geo_reference=Geo_reference(56,x,y))
2001        polygon = [[0,0],[1,0],[1,1],[0,1]]
2002       
2003        m.add_region_from_polygon(polygon,
2004                               {'tagin':[0,1],'bom':[2]},
2005                                  max_triangle_area=10)
2006        self.failUnless(len(m.regions)==1,
2007                        'FAILED!')
2008        segs = m.getUserSegments()
2009        self.failUnless(len(segs)==4,
2010                        'FAILED!')
2011        self.failUnless(len(m.userVertices)==4,
2012                        'FAILED!') 
2013        self.failUnless(segs[0].tag=='tagin',
2014                        'FAILED!') 
2015        self.failUnless(segs[1].tag=='tagin',
2016                        'FAILED!') 
2017         
2018        self.failUnless(segs[2].tag=='bom',
2019                        'FAILED!') 
2020        self.failUnless(segs[3].tag=='',
2021                        'FAILED!')
2022        verts = m.getUserVertices()
2023        #print "User verts",verts
2024        #print 'polygon',polygon
2025        #vert values are relative
2026        for point,new_point in map(None,polygon,verts):
2027            point_x = point[0] 
2028            new_point_x = new_point.x + m.geo_reference.get_xllcorner()
2029            #print "point_x",point_x
2030            #print "new_point_x",new_point_x
2031            point_y = point[1] 
2032            new_point_y = new_point.y + m.geo_reference.get_yllcorner()
2033           
2034            self.failUnless(point_x == new_point_x, ' failed')
2035            self.failUnless(point_y == new_point_y, ' failed')
2036
2037
2038    def test_add_hole_from_polygon(self):
2039        x=-500
2040        y=-1000
2041        m=Mesh(geo_reference=Geo_reference(56,x,y))
2042
2043        # These are the absolute values
2044        polygon_absolute = [[0,0],[1,0],[1,1],[0,1]]
2045       
2046        x_p = -10
2047        y_p = -40
2048        geo_ref_poly = Geo_reference(56, x_p, y_p)
2049        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
2050       
2051        poly_point = m.add_hole_from_polygon(polygon,
2052                                               {'tagin':[0,1],'bom':[2]},
2053                                               geo_reference=geo_ref_poly)
2054        # poly_point values are relative to the mesh geo-ref
2055        # make them absolute
2056        #print "poly_point.x+x",poly_point.x+x
2057        #print "poly_point.y+y",poly_point.y+y
2058        #print "polygon_absolute", polygon_absolute
2059        self.failUnless(is_inside_polygon([poly_point.x+x,poly_point.y+y],
2060                                       polygon_absolute, closed = False),
2061                        'FAILED!')
2062               
2063        self.failUnless(len(m.holes)==1,
2064                        'FAILED!')
2065        segs = m.getUserSegments()
2066        self.failUnless(len(segs)==4,
2067                        'FAILED!')
2068        self.failUnless(len(m.userVertices)==4,
2069                        'FAILED!') 
2070        self.failUnless(segs[0].tag=='tagin',
2071                        'FAILED!') 
2072        self.failUnless(segs[1].tag=='tagin',
2073                        'FAILED!') 
2074         
2075        self.failUnless(segs[2].tag=='bom',
2076                        'FAILED!') 
2077        self.failUnless(segs[3].tag=='',
2078                        'FAILED!')
2079        verts = m.getUserVertices()
2080        #print "User verts",verts
2081        #print 'polygon',polygon
2082        #vert values are relative
2083        for point,new_point in map(None,polygon,verts):
2084            point_x = point[0] + geo_ref_poly.get_xllcorner()
2085            new_point_x = new_point.x + m.geo_reference.get_xllcorner()
2086            point_y = point[1] + geo_ref_poly.get_yllcorner()
2087            #print "new_point.y",new_point.y
2088            #print "m.geo_ref.get_yllcorner()",m.geo_reference.get_yllcorner()
2089            new_point_y = new_point.y + m.geo_reference.get_yllcorner()
2090            #print "point_y",point_y
2091            #print "new_point_y",new_point_y
2092           
2093            self.failUnless(point_x == new_point_x, ' failed')
2094            self.failUnless(point_y == new_point_y, ' failed')
2095
2096    def test_add_circle(self):
2097        x=-500
2098        y=-1000
2099        m=Mesh(geo_reference=Geo_reference(56,x,y))
2100
2101        # These are the absolute values
2102        tag = 'hey'
2103        segment_count = 104
2104        radius = 30
2105        circle_center_absolute = [100,80]       
2106        x_p = -.666
2107        y_p = -.777
2108        geo_ref_poly = Geo_reference(56, x_p, y_p)
2109        circle_center = \
2110                geo_ref_poly.change_points_geo_ref(circle_center_absolute)
2111        circle_center = circle_center[0] #make a list of lists a list
2112        poly_point = m.add_circle(circle_center, radius, segment_count,
2113                                  tag=tag,
2114                                  region=True,
2115                                  center_geo_reference=geo_ref_poly)
2116        # poly_point values are relative to the mesh geo-ref
2117        # make them absolute
2118        #print "poly_point.x+x",poly_point.x+x
2119        #print "polygon_absolute", polygon_absolute
2120     
2121       
2122        #m.export_mesh_file("aaat.msh")
2123       
2124        self.failUnless(len(m.regions)==1,
2125                        'FAILED!')
2126        segs = m.getUserSegments()
2127        self.failUnless(len(segs)==segment_count,
2128                        'FAILED!')
2129        self.failUnless(len(m.userVertices)==segment_count,
2130                        'FAILED!') 
2131        self.failUnless(segs[0].tag==tag,
2132                        'FAILED!') 
2133        self.failUnless(segs[1].tag==tag,
2134                        'FAILED!') 
2135         
2136        verts = m.getUserVertices()
2137       
2138        #m.export_mesh_file("aaat.msh")
2139       
2140    def FIXMEtest_auto_set_geo_reference(self):
2141        x=50000
2142        y=1000
2143        m=Mesh(geo_reference=Geo_reference(56,x,y))
2144        polygon = [[0,0],[1,0],[1,1],[0,1]]
2145       
2146        m.add_region_from_polygon(polygon,
2147                               {'tagin':[0,1],'bom':[2]},
2148                                  max_triangle_area=10)
2149        m.auto_set_geo_reference()
2150       
2151 
2152    def test_duplicat_verts_are_removed(self):
2153   
2154     
2155        a = Vertex ( 0.0 ,0.0)
2156        b = Vertex (0.0, 4.0)
2157        c = Vertex (4.0,4.0)
2158        d = Vertex (4.0,0.0)
2159        e = Vertex (4.0,0.0) # duplicate point
2160   
2161        s1 = Segment(a,b, tag = "50")
2162        s2 = Segment(b,c, tag = "40")
2163        s3 = Segment(c,d, tag = "30")
2164        s4 = Segment(d,e, tag = "no where seg")
2165        s5 = Segment(e,a, tag = "20")
2166
2167       
2168        m = Mesh(userVertices=[a,b,c,d,e],
2169                 userSegments=[s1,s2,s3,s4,s5])
2170
2171        seg = m.getUserSegments()
2172        points = m.getUserVertices()
2173        holes = m.getHoles()
2174        regions = m.getRegions()
2175        #fileName = tempfile.mktemp(".tsh")
2176        #fileName = "badmesh.tsh"
2177        #m.export_mesh_file(fileName)
2178        #print "***************************fileName", fileName
2179        #new_m = importMeshFromFile(fileName)
2180        #os.remove(fileName)
2181       
2182        m.generateMesh("Q", maxArea = 2000.1 )
2183
2184        #m.export_mesh_file("from_test_mesh.tsh")
2185        seg = m.getMeshSegments()
2186        self.failUnless(4==len(seg),
2187                        'FAILED!') 
2188
2189        vert = m.getMeshVertices() 
2190        self.failUnless(4==len(vert),
2191                        'FAILED!')
2192 
2193    def test_duplicat_verts_are_removedII(self):
2194   
2195     
2196        a = Vertex ( 0.0 ,0.0)
2197        b = Vertex (0.0, 4.0)
2198        c = Vertex (4.0,4.0)
2199        d = Vertex (4.0,0.0)
2200        e = Vertex (4.0,0.0) # duplicate point
2201        f = Vertex (49.0,0.0) # unused point
2202   
2203        s1 = Segment(a,b, tag = "50")
2204        s2 = Segment(b,c, tag = "40")
2205        s3 = Segment(c,d, tag = "30")
2206        s4 = Segment(d,e, tag = "no where seg")
2207        s5 = Segment(e,a, tag = "20")
2208
2209       
2210        m = Mesh(userVertices=[a,b,c,d,e,f],
2211                 userSegments=[s1,s2,s3,s4,s5])
2212
2213        seg = m.getUserSegments()
2214        points = m.getUserVertices()
2215        holes = m.getHoles()
2216        regions = m.getRegions()
2217        #fileName = tempfile.mktemp(".tsh")
2218        #fileName = "badmesh.tsh"
2219        #m.export_mesh_file(fileName)
2220        #print "***************************fileName", fileName
2221        #new_m = importMeshFromFile(fileName)
2222        #os.remove(fileName)
2223       
2224        m.generateMesh("Q", maxArea = 2000.1 )
2225
2226        #m.export_mesh_file("from_test_mesh.tsh")
2227        seg = m.getMeshSegments()
2228        self.failUnless(4==len(seg),
2229                        'FAILED!') 
2230
2231        vert = m.getMeshVertices() 
2232        self.failUnless(4==len(vert),
2233                        'FAILED!')
2234   
2235    def test_add_vertices(self):
2236        points_ab = [[0.1,1],[0.4,.2],[7,5],[10,5]]
2237        geo =  Geo_reference(56,23,21)
2238        points = geo.change_points_geo_ref(points_ab)
2239        spat = Geospatial_data(points, geo_reference=geo)
2240       
2241        geo_mesh =  Geo_reference(56,100,200)
2242        m = Mesh(geo_reference=geo_mesh)
2243        m.add_vertices(spat)
2244
2245        vert = m.getUserVertices()
2246        #print "vert",vert
2247        self.failUnless(4==len(vert),
2248                        'FAILED!')
2249        vert= m.get_user_vertices(absolute=True)
2250        self.failUnless(vert==points_ab,
2251                        'FAILED!')
2252
2253   
2254    def test_add_vertices_more(self):
2255        points = [[0.1,1],[0.4,.2],[7,5],[10,5]]
2256        #spat = Geospatial_data(points)
2257       
2258        m = Mesh()
2259        m.add_vertices(points)
2260
2261        vert = m.getUserVertices()
2262        #print "vert",vert
2263        self.failUnless(4==len(vert),
2264                        'FAILED!')
2265        vert= m.get_user_vertices(absolute=True)
2266        self.failUnless(vert==points,
2267                        'FAILED!')
2268   
2269    def test_add_verticesII(self):
2270        points_lat_long = [[-33,152],[-35,152],[-35,150],[-33,150]]
2271       
2272        spat = Geospatial_data(data_points=points_lat_long,
2273                               points_are_lats_longs=True)
2274        points_ab = spat.get_data_points( absolute = True)
2275        m = Mesh()
2276        m.add_vertices(spat)
2277
2278        vert = m.getUserVertices()
2279        #print "vert",vert
2280        self.failUnless(4==len(vert),
2281                        'FAILED!')
2282        vert= m.get_user_vertices(absolute=True)
2283        self.failUnless(vert==points_ab,
2284                        'FAILED!')
2285
2286        spat = Geospatial_data(data_points=points_lat_long,
2287                               points_are_lats_longs=True)
2288        points_ab = spat.get_data_points( absolute = True)
2289        geo =  Geo_reference(56,400000,6000000)
2290        spat.set_geo_reference(geo)
2291        m = Mesh()
2292        m.add_vertices(spat)
2293
2294        vert = m.getUserVertices()
2295        #print "vert",vert
2296        self.failUnless(4==len(vert),
2297                        'FAILED!')
2298        vert= m.get_user_vertices(absolute=True)
2299        self.failUnless(vert==points_ab,
2300                        'FAILED!')
2301
2302        #geo =  Geo_reference(56,23,21)
2303        #points = geo.change_points_geo_ref(points_ab)
2304       
2305    def test_get_user_vertices(self):
2306        points_ab = [[0.1,1],[0.4,.2],[7,5],[10,5]]
2307        geo =  Geo_reference(56,23,21)
2308        points = geo.change_points_geo_ref(points_ab)
2309        spat = Geospatial_data(points, geo_reference=geo)
2310       
2311        geo_mesh =  Geo_reference(56,100,200)
2312        m = Mesh(geo_reference=geo_mesh)
2313        m.add_vertices(spat)
2314
2315        vert = m.getUserVertices()
2316        #print "vert",vert
2317        self.failUnless(4==len(vert),
2318                        'FAILED!')
2319        vert= m.get_user_vertices(absolute=True)
2320        self.failUnless(vert==points_ab,
2321                        'FAILED!')
2322        vert= m.get_user_vertices(absolute=False)
2323        points_new = m.geo_reference.get_absolute(vert)
2324       
2325        self.failUnless(points_ab==points_new,
2326                        'FAILED!')
2327       
2328       
2329def list_comp(A,B):
2330    yes = len(A)==len(B)
2331    for item in A:
2332        if not item in B:
2333            yes = False
2334    return yes
2335
2336
2337           
2338#-------------------------------------------------------------
2339if __name__ == "__main__":
2340    suite = unittest.makeSuite(meshTestCase,'test')
2341    #suite = unittest.makeSuite(meshTestCase,'test_exportASCIIsegmentoutlinefile')
2342    #suite = unittest.makeSuite(meshTestCase,'test_e')
2343    runner = unittest.TextTestRunner() #verbosity=2)
2344    runner.run(suite)
2345   
Note: See TracBrowser for help on using the repository browser.