source: inundation/ga/storm_surge/pmesh/mesh.py @ 364

Last change on this file since 364 was 364, checked in by duncan, 20 years ago

adding vert attribute titles
crippling loading vert attributes
File size: 60.1 KB
Line 
1#!/usr/bin/env python
2#
3"""General 2D triangular classes for triangular mesh generation.
4
5   Copyright 2003/2004
6   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
7   Geoscience Australia
8"""
9
10import load_mesh.loadASCII
11
12import sys
13import math
14import triang
15import re
16import os
17import pickle
18
19import types
20import exceptions
21
22initialconversions = ['internal', 'external','internal']
23
24#from os import sep
25#sys.path.append('..'+sep+'pmesh')
26#print "sys.path",sys.path
27
28class MeshObject:
29    """
30    An abstract superclass for the basic mesh objects, eg vertex, segment,
31    triangle.
32    """
33    def __init__(self):
34        pass
35   
36class Point(MeshObject): 
37    """
38    Define a point in a 2D space.
39    """
40    def __init__(self,X,Y):
41        __slots__ = ['x','y']
42        self.x=X
43        self.y=Y
44       
45    def DistanceToPoint(self, OtherPoint):
46        """
47        Returns the distance from this point to another
48        """
49        SumOfSquares = ((self.x - OtherPoint.x)**2) + ((self.y - OtherPoint.y)**2)
50        return math.sqrt(SumOfSquares)
51       
52    def IsInsideCircle(self, Center, Radius):
53        """
54        Return 1 if this point is inside the circle,
55        0 otherwise
56        """
57       
58        if (self.DistanceToPoint(Center)<Radius):
59            return 1
60        else:
61            return 0
62       
63    def __repr__(self):
64        return "(%f,%f)" % (self.x,self.y) 
65
66    def cmp_xy(self, point):
67        if self.x < point.x:
68            return -1
69        elif self.x > point.x:
70            return 1
71        else:           
72            if self.y < point.y:
73                return -1
74            elif self.y > point.y:
75                return 1
76            else:
77                return 0
78       
79    def same_x_y(self, point):
80        if self.x == point.x and self.y == point.y:
81            return True
82        else:
83            return False
84       
85           
86
87class Vertex(Point):
88    """
89    A point on the mesh.
90    Object attributes based on the Triangle program
91    """
92    def __init__(self,X,Y, attributes = None):
93        __slots__ = ['x','y','attributes']
94        self.x=X
95        self.y=Y       
96        self.attributes=[] 
97       
98        if attributes is None:
99            self.attributes=[]
100        else:
101            self.attributes=attributes
102   
103
104    def setAttributes(self,attributes):
105        """
106        attributes is a list.
107        """
108        self.attributes = attributes
109       
110    VERTEXSQUARESIDELENGTH = 6
111    def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, yoffset =0, ):
112        x =  scale*(self.x + xoffset)
113        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
114        #print "draw x:", x
115        #print "draw y:", y
116        cornerOffset= self.VERTEXSQUARESIDELENGTH/2
117        return canvas.create_rectangle(x-cornerOffset,
118                                       y-cornerOffset,
119                                       x+cornerOffset,
120                                       y+cornerOffset,
121                                       tags = tags,
122                                       outline=colour,
123                                       fill = 'white')
124     
125    def __repr__(self):
126        return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
127   
128class Hole(Point):
129    """
130    A region of the mesh were no triangles are generated.
131    Defined by a point in the hole enclosed by segments.
132    """
133    HOLECORNERLENGTH = 6
134    def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ):
135        x =  scale*(self.x + xoffset)
136        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
137        #print "draw x:", x
138        #print "draw y:", y
139        cornerOffset= self.HOLECORNERLENGTH/2
140        return canvas.create_oval(x-cornerOffset,
141                                       y-cornerOffset,
142                                       x+cornerOffset,
143                                       y+cornerOffset,
144                                       tags = tags,
145                                       outline=colour,
146                                       fill = 'white')
147   
148class Region(Point):
149    """
150    A region of the mesh, defined by a point in the region
151    enclosed by segments. Used to tag areas.
152    """
153    CROSSLENGTH = 6
154    TAG = 0
155    MAXAREA = 1
156   
157    def __init__(self,X,Y, tag = None, maxArea = None):
158        """Precondition: tag is a string and maxArea is a real
159        """
160        # This didn't work. 
161        #super(Region,self)._init_(self,X,Y)
162        self.x=X
163        self.y=Y   
164        self.attributes =[] # index 0 is the tag string
165                            #optoinal index 1 is the max triangle area
166                            #NOTE the size of this attribute is assumed
167                            # to be 1 or 2 in regionstrings2int
168        if tag is None:
169            self.attributes.append("")
170        else:
171            self.attributes.append(tag) #this is a string
172           
173        if maxArea is not None:
174            self.setMaxArea(maxArea) # maxArea is a number
175           
176    def getTag(self,):
177        return self.attributes[self.TAG]
178   
179    def setTag(self,tag):
180        self.attributes[self.TAG] = tag
181       
182    def getMaxArea(self):
183        """ Returns the Max Area of a Triangle or
184        None, if the Max Area has not been set.
185        """
186        if self.isMaxArea():
187            return self.attributes[1]
188        else:
189            return None
190   
191    def setMaxArea(self,MaxArea):
192        if self.isMaxArea(): 
193            self.attributes[self.MAXAREA] = float(MaxArea)
194        else:
195            self.attributes.append( float(MaxArea) )
196   
197    def deleteMaxArea(self):
198        if self.isMaxArea():
199            self.attributes.pop(self.MAXAREA)
200           
201    def isMaxArea(self):
202        return len(self.attributes)> 1
203   
204    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "black"):
205        """
206        Draw a black cross, returning the objectID
207        """
208        x =  scale*(self.x + xoffset)
209        y = -1*scale*(self.y + yoffset) 
210        cornerOffset= self.CROSSLENGTH/2
211        return canvas.create_polygon(x,
212                                     y-cornerOffset,
213                                     x,
214                                     y,
215                                     x+cornerOffset,
216                                     y,
217                                     x,
218                                     y,
219                                     x,
220                                     y+cornerOffset,
221                                     x,
222                                     y,
223                                     x-cornerOffset,
224                                     y,
225                                     x,
226                                     y,
227                                     tags = tags,
228                                     outline = colour,fill = '')
229   
230    def __repr__(self):
231        if self.isMaxArea():
232            area = self.getMaxArea() 
233            return "(%f,%f,%s,%f)" % (self.x,self.y,
234                                      self.getTag(), self.getMaxArea())
235        else:
236            return "(%f,%f,%s)" % (self.x,self.y,
237                                   self.getTag())
238       
239class Triangle(MeshObject):
240    """
241    A triangle element, defined by 3 vertices.
242    Attributes based on the Triangle program.
243    """
244
245    def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ):
246        """
247        Vertices, the initial arguments, are listed in counterclockwise order.
248        """
249        self.vertices= [vertex1,vertex2, vertex3 ]
250       
251        if attribute is None:
252            self.attribute =""
253        else:
254            self.attribute = attribute #this is a string
255           
256        if neighbors is None:
257            self.neighbors=[]
258        else:
259            self.neighbors=neighbors
260
261    def getVertices(self):
262        return self.vertices
263   
264    def calcArea(self):
265        ax = self.vertices[0].x
266        ay = self.vertices[0].y
267       
268        bx = self.vertices[1].x
269        by = self.vertices[1].y
270       
271        cx = self.vertices[2].x
272        cy = self.vertices[2].y
273       
274        return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
275   
276    def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
277        """
278        neighbor1 is the triangle opposite vertex1 and so on.
279        Null represents no neighbor
280        """
281        self.neighbors = [neighbor1, neighbor2, neighbor3]
282
283    def setAttribute(self,attribute):
284        """
285        neighbor1 is the triangle opposite vertex1 and so on.
286        Null represents no neighbor
287        """
288        self.attribute = attribute
289       
290    def __repr__(self):
291        return "[%s,%s]" % (self.vertices,self.attribute)
292       
293
294    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"):
295        """
296        Draw a red triagle, returning the objectID
297        """
298        return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
299                                     scale*-1*(self.vertices[1].y + yoffset),
300                                     scale*(self.vertices[0].x + xoffset),
301                                     scale*-1*(self.vertices[0].y + yoffset),
302                                     scale*(self.vertices[2].x + xoffset),
303                                     scale*-1*(self.vertices[2].y + yoffset),
304                                     tags = tags,
305                                     outline = colour,fill = '')
306
307class Segment(MeshObject):
308    """
309    Segments are edges whose presence in the triangulation is enforced.
310   
311    """
312    def __init__(self, vertex1, vertex2, marker = None ):
313        """
314        Each segment is specified by listing the vertices of its endpoints
315        """
316
317        assert(vertex1 != vertex2)
318        self.vertices = [vertex1,vertex2 ]
319       
320        if marker is None:
321            self.marker = self.__class__.default
322        else:
323            self.marker = marker #this is a string
324       
325    def __repr__(self):
326        return "[%s,%s]" % (self.vertices,self.marker)
327           
328       
329    def draw(self, canvas, tags,scale=1 , xoffset=0 , yoffset=0,colour='blue' ):
330        x=[]
331        y=[]
332        for end in self.vertices:
333            #end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
334            x.append(scale*(end.x + xoffset))
335            y.append(-1*scale*(end.y + yoffset))  # - since for a canvas - is up
336        return canvas.create_line(x[0], y[0], x[1], y[1], tags = tags,fill=colour)
337    # Class methods
338    def set_default_tag(cls, default):
339        cls.default = default
340   
341    def get_default_tag(cls):
342        return cls.default
343   
344    set_default_tag = classmethod(set_default_tag) 
345    get_default_tag = classmethod(get_default_tag)
346
347Segment.set_default_tag("")       
348class Mesh:
349    """
350    Representation of a 2D triangular mesh.
351    User attributes describe the mesh region/segments/vertices/attributes
352
353    mesh attributes describe the mesh that is produced eg triangles and vertices.
354   
355    The Mesh holds user information to define the
356    """
357
358    def __repr__(self):
359        return """
360        mesh Triangles: %s 
361        mesh Segments: %s 
362        mesh Vertices: %s 
363        user Segments: %s 
364        user Vertices: %s 
365        holes: %s 
366        regions: %s""" % (self.meshTriangles,
367                                self.meshSegments,
368                                self.meshVertices,
369                                self.userSegments,
370                                self.userVertices,
371                                self.holes,
372                                self.regions) 
373   
374    def __init__(self, userSegments=None, userVertices=None, holes=None, regions=None):
375        self.meshTriangles=[] 
376        self.meshSegments=[]
377        self.meshVertices=[]
378       
379        if userSegments is None:
380            self.userSegments=[]
381        else:
382            self.userSegments=userSegments
383           
384        if userVertices is None:
385            self.userVertices=[]
386        else:
387            self.userVertices=userVertices
388           
389        if holes is None:
390            self.holes=[]
391        else:
392            self.holes=holes
393           
394        if regions is None:
395            self.regions=[]
396        else:
397            self.regions=regions
398    def __cmp__(self,other):
399       
400        # A dic for the initial m
401        dic = self.Mesh2triangList()
402        dic_mesh = self.Mesh2MeshList()
403        for element in dic_mesh.keys():
404            dic[element] = dic_mesh[element]
405       
406        # A dic for the exported/imported m
407        dic_other = other.Mesh2triangList()
408        dic_mesh = other.Mesh2MeshList()
409        for element in dic_mesh.keys():
410            dic_other[element] = dic_mesh[element]
411
412        #print "dsg************************8"
413        #print "dic ",dic
414        #print "*******8"
415        #print "mesh",dic_other
416        #print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
417        #print "dsg************************8"
418       
419        return (dic.__cmp__(dic_other))
420       
421    def generateMesh(self, mode = None, maxArea = None, isRegionalMaxAreas = True):
422        """
423        Based on the current user vaules, holes and regions
424        generate a new mesh
425        mode is a string that sets conditions on the mesh generations
426        see triangle_instructions.txt for a definition of the commands
427        PreCondition: maxArea is a double
428        """
429        if mode == None:
430            self.mode = ""
431        else:
432            self.mode = mode
433       
434        if not re.match('p',self.mode):
435            self.mode += 'p' #p - read a planar straight line graph.
436            #there must be segments to use this switch
437            # TODO throw an aception if there aren't seg's
438            # it's more comlex than this.  eg holes
439        if not re.match('z',self.mode):
440            self.mode += 'z' # z - Number all items starting from zero (rather than one)
441        if not re.match('n',self.mode):
442            self.mode += 'n' # n - output a list of neighboring triangles
443        if not re.match('A',self.mode):
444            self.mode += 'A' # A - output region attribute list for triangles
445        if not re.match('V',self.mode):
446            self.mode += 'V' # V - output info about what Triangle is doing
447       
448        if maxArea != None:
449            self.mode += 'a' + str(maxArea)
450
451        if isRegionalMaxAreas:
452            self.mode += 'a'
453           
454        meshDict = self.Mesh2triangList()
455        #print "!@!@ This is going to triangle   !@!@"
456        #print meshDict
457        #print "!@!@ This is going to triangle   !@!@"
458
459        #print "meshDict['segmentmarkerlist']", meshDict['segmentmarkerlist']
460        #initialconversions = ['internal boundary', 'external boundary','internal boundary']
461        [meshDict['segmentmarkerlist'],
462         segconverter] =  segment_strings2ints(meshDict['segmentmarkerlist'],
463                                             initialconversions)
464        #print "regionlist",meshDict['regionlist']
465        [meshDict['regionlist'],
466         regionconverter] =  region_strings2ints(meshDict['regionlist'])
467        #print "regionlist",meshDict['regionlist']
468       
469        #print "meshDict['segmentmarkerlist']", meshDict['segmentmarkerlist'
470        generatedMesh = triang.genMesh(
471                              meshDict['pointlist'],
472                              meshDict['segmentlist'],
473                              meshDict['holelist'],
474                              meshDict['regionlist'],
475                              meshDict['pointattributelist'],
476                              meshDict['segmentmarkerlist'],
477                              [],  # since the trianglelist isn't used
478                              self.mode)
479        #print "generated",generatedMesh['generatedsegmentmarkerlist']
480        generatedMesh['generatedsegmentmarkerlist'] = \
481                     segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
482                                  segconverter)
483        #print "processed gen",generatedMesh['generatedsegmentmarkerlist']
484        generatedMesh['generatedtriangleattributelist'] = \
485         region_ints2strings(generatedMesh['generatedtriangleattributelist'],
486                                  regionconverter)
487       
488        self.setTriangulation(generatedMesh)
489   
490    def addUserPoint(self, pointType, x,y):
491        if pointType == Vertex:
492            point = self.addUserVertex(x,y)
493        if pointType == Hole:
494            point = self.addHole(x,y)
495        if pointType == Region:
496            point = self.addRegion(x,y)
497        return point
498   
499    def addUserVertex(self, x,y):
500        v=Vertex(x, y)
501        self.userVertices.append(v)
502        return v
503
504    def addHole(self, x,y):
505        h=Hole(x, y)
506        self.holes.append(h)
507        return h
508   
509    def addRegion(self, x,y):
510        h=Region(x, y)
511        self.regions.append(h)
512        return h
513   
514    def getUserVertices(self):
515        return self.userVertices
516   
517    def getUserSegments(self):
518        return self.userSegments
519   
520    def getTriangulation(self):
521        return self.meshTriangles
522   
523    def getMeshVertices(self):
524        return self.meshVertices
525 
526    def getMeshSegments(self):
527        return self.meshSegments
528   
529    def getHoles(self):
530        return self.holes
531   
532    def getRegions(self):
533        return self.regions
534   
535    def isTriangulation(self):
536        if self.meshVertices == []:
537            return False 
538        else:
539            return True
540   
541    def addUserSegment(self, v1,v2):
542        """
543        PRECON: A segment between the two vertices is not already present.
544        Check by calling isUserSegmentNew before calling this function.
545       
546        """
547        s=Segment( v1,v2)
548        self.userSegments.append(s)
549        return s
550
551   
552    def clearTriangulation(self):
553
554        #Clear the current generated mesh values
555        self.meshTriangles=[] 
556        self.meshSegments=[]
557        self.meshVertices=[]
558
559    def removeDuplicatedUserVertices(self):
560        """Pre-condition: There are no user segments
561        This function will keep the first duplicate
562        """
563        assert self.userSegments == []
564        self.userVertices, counter =  self.removeDuplicatedVertices(self.userVertices)
565        return counter
566   
567    def removeDuplicatedVertices(self, Vertices):
568        """
569        This function will keep the first duplicate, remove all others
570        Precondition: Each vertex has a dupindex, which is the list
571        index.
572        """
573        remove = []
574        index = 0
575        for v in Vertices:
576            v.dupindex = index
577            index += 1
578        t = list(Vertices)
579        t.sort(Point.cmp_xy)
580   
581        length = len(t)
582        behind = 0
583        ahead  = 1
584        counter = 0
585        while ahead < length:
586            b = t[behind]
587            ah = t[ahead]
588            if (b.y == ah.y and b.x == ah.x):
589                remove.append(ah.dupindex) 
590            behind += 1
591            ahead += 1
592
593        # remove the duplicate vertices
594        remove.sort()
595        remove.reverse()
596        for i in remove:
597            Vertices.pop(i)
598            pass
599
600        #Remove the attribute that this function added
601        for v in Vertices:
602            del v.dupindex
603        return Vertices,counter
604   
605    def thinoutVertices(self, delta):
606        """Pre-condition: There are no user segments
607        This function will keep the first duplicate
608        """
609        assert self.userSegments == []
610        #t = self.userVertices
611        #self.userVertices =[]
612        boxedVertices = {}
613        thinnedUserVertices =[]
614        delta = round(delta,1)
615       
616        for v in self.userVertices :
617            # marker is the center of the boxes
618            marker = (round(v.x/delta,0)*delta,round(v.y/delta,0)*delta)
619            #this creates a dict of lists of faces, indexed by marker
620            boxedVertices.setdefault(marker,[]).append(v)
621
622        for [marker,verts] in boxedVertices.items():
623            min = delta
624            bestVert = None
625            markerVert = Vertex(marker[0],marker[1])
626            for v in verts:
627                dist = v.DistanceToPoint(markerVert)
628                if (dist<min):
629                    min = dist
630                    bestVert = v
631            thinnedUserVertices.append(bestVert)
632        self.userVertices =thinnedUserVertices
633       
634           
635    def isUserSegmentNew(self, v1,v2):
636        identicalSegs= [x for x in self.userSegments if (x.vertices[0] == v1 and x.vertices[1] == v2) or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
637       
638        return len(identicalSegs) == 0
639
640    def deleteSegsOfVertex(self, delVertex):
641        """
642        Delete this vertex and any segments that connect to it.
643        """
644        #Find segments that connect to delVertex
645        deletedSegments = []
646        for seg in self.userSegments:
647            if (delVertex in seg.vertices):
648                deletedSegments.append(seg)
649        # Delete segments that connect to delVertex
650        for seg in deletedSegments:
651            self.userSegments.remove(seg)
652        # Delete delVertex
653        # self.userVertices.remove(delVertex)
654        return deletedSegments
655
656   
657    def deleteMeshObject(self, MeshObject):
658        """
659        Returns a list of all objects that were removed
660        """
661        deletedObs = []
662        if isinstance(MeshObject, Vertex ):
663            deletedObs = self.deleteSegsOfVertex(MeshObject)
664            deletedObs.append(MeshObject)
665            self.userVertices.remove(MeshObject)
666        elif isinstance(MeshObject, Segment):
667            deletedObs.append(MeshObject)
668            self.userSegments.remove(MeshObject)
669        elif isinstance(MeshObject, Hole):
670            deletedObs.append(MeshObject)
671            self.holes.remove(MeshObject)
672        elif isinstance(MeshObject, Region):
673            deletedObs.append(MeshObject)
674            self.regions.remove(MeshObject)         
675        return deletedObs
676                                                 
677    def Mesh2triangList(self):
678        """
679        Convert the Mesh to a dictionary of the lists needed for the triang modul;
680        points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
681        pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
682        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
683        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole region)
684        regionlist: [ (x1,y1,index),...] (Tuple of 3 doubles)
685        Note, this adds an index attribute to the user Vertex objects.
686        """
687       
688        meshDict = {}
689       
690        pointlist=[]
691        pointattributelist=[]
692        index = 0
693
694        for vertex in self.userVertices:
695            vertex.index = index
696            pointlist.append((vertex.x,vertex.y))
697            pointattributelist.append((vertex.attributes))
698           
699            index += 1
700        meshDict['pointlist'] = pointlist
701        meshDict['pointattributelist'] = pointattributelist
702
703        segmentlist=[]
704        segmentmarkerlist=[]
705        for seg in self.userSegments:
706            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
707            segmentmarkerlist.append(seg.marker)
708        meshDict['segmentlist'] =segmentlist
709        meshDict['segmentmarkerlist'] =segmentmarkerlist
710
711        holelist=[]
712        for hole in self.holes:
713            holelist.append((hole.x,hole.y)) 
714        meshDict['holelist'] = holelist
715       
716        regionlist=[]
717        for region in self.regions:
718            if (region.getMaxArea() != None): 
719                regionlist.append((region.x,region.y,region.getTag(),
720                               region.getMaxArea()))
721            else:
722                regionlist.append((region.x,region.y,region.getTag()))
723        meshDict['regionlist'] = regionlist
724        #print "*(*("
725        #print meshDict
726        #print meshDict['regionlist']
727        #print "*(*("
728        return meshDict
729                                               
730    def Mesh2MeshList(self):
731        """
732        Convert the Mesh to a dictionary of lists describing the triangulation variables;
733        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
734        generated point attribute list: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
735        generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
736        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole region)
737        regionlist: [ (x1,y1,index),...] (Tuple of 3 doubles)
738        Note, this adds an index attribute to the user Vertex objects.
739        """
740       
741        meshDict = {}       
742        pointlist=[]
743        pointattributelist=[]
744       
745        index = 0
746        for vertex in self.meshVertices:
747            vertex.index = index
748            pointlist.append((vertex.x,vertex.y))
749            pointattributelist.append((vertex.attributes))           
750            index += 1
751           
752        meshDict['generatedpointlist'] = pointlist
753        meshDict['generatedpointattributelist'] = pointattributelist
754
755        #segments
756        segmentlist=[]
757        segmentmarkerlist=[]
758        for seg in self.meshSegments:
759            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
760            segmentmarkerlist.append(seg.marker)
761        meshDict['generatedsegmentlist'] =segmentlist
762        meshDict['generatedsegmentmarkerlist'] =segmentmarkerlist
763
764        # Make sure that the indexation is correct
765        index = 0
766        for tri in self.meshTriangles:
767            tri.index = index           
768            index += 1
769
770        trianglelist = []
771        triangleattributelist = []
772        triangleneighborlist = []
773        for tri in self.meshTriangles: 
774            trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,tri.vertices[2].index)) 
775            triangleattributelist.append(tri.attribute)
776            neighborlist = [-1,-1,-1]
777            for neighbor,index in map(None,tri.neighbors,
778                                      range(len(tri.neighbors))):
779                if neighbor:   
780                    neighborlist[index] = neighbor.index
781            triangleneighborlist.append(neighborlist)
782       
783        meshDict['generatedtrianglelist'] = trianglelist
784        meshDict['generatedtriangleattributelist'] = triangleattributelist
785        meshDict['generatedtriangleneighborlist'] = triangleneighborlist
786       
787        #print "*(*("
788        #print meshDict
789        #print "*(*("
790        return meshDict
791
792                               
793    def Mesh2MeshDic(self):
794        """
795        Convert the user and generated info of a mesh to a dictionary
796        structure
797        """
798        dic = self.Mesh2triangList()
799        dic_mesh = self.Mesh2MeshList()
800        for element in dic_mesh.keys():
801            dic[element] = dic_mesh[element]
802        return dic
803   
804    def setTriangulation(self, genDict):
805        """
806        Set the mesh attributes given a dictionary of the lists
807        returned from the triang module       
808        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
809        generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
810        generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
811        generated segment marker list: [S1Marker, S2Marker, ...] (list of ints)
812        triangle list:  [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
813        triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
814        triangle attribute list: [(T1att), (T2att), ...] (list of a list of strings)
815        """
816
817        #Clear the current generated mesh values
818        self.meshTriangles=[] 
819        self.meshSegments=[]
820        self.meshVertices=[]
821
822        #print "@#@#@#"
823        #print genDict
824        #print "@#@#@#"
825       
826        index = 0
827        for point in genDict['generatedpointlist']:
828            v=Vertex(point[0], point[1])
829            v.index = index
830            index +=1
831            self.meshVertices.append(v)
832
833        index = 0
834        for seg,marker in map(None,genDict['generatedsegmentlist'],genDict['generatedsegmentmarkerlist']):
835            segObject = Segment( self.meshVertices[seg[0]],
836                           self.meshVertices[seg[1]], marker = marker )
837            segObject.index = index
838            index +=1
839            self.meshSegments.append(segObject)
840
841        index = 0
842        for triangle in genDict['generatedtrianglelist']:
843            tObject =Triangle( self.meshVertices[triangle[0]],
844                        self.meshVertices[triangle[1]],
845                        self.meshVertices[triangle[2]] )
846            tObject.index = index
847            index +=1
848            self.meshTriangles.append(tObject)
849
850        index = 0
851        for att in genDict['generatedtriangleattributelist']:
852            if att == []:
853                self.meshTriangles[index].setAttribute("")
854            else:
855                # Note, is the first attribute always the region att?
856                # haven't confirmed this
857                self.meshTriangles[index].setAttribute(att[0])
858            index += 1
859           
860        index = 0
861        for att in genDict['generatedpointattributelist']:
862            if att == None:
863                self.meshVertices[index].setAttributes([])
864            else:
865                self.meshVertices[index].setAttributes(att)
866            index += 1
867   
868        index = 0
869        for triangle in genDict['generatedtriangleneighborlist']:
870            # Build a list of triangle object neighbors
871            ObjectNeighbor = []
872            for neighbor in triangle:
873                if ( neighbor != -1):
874                    ObjectNeighbor.append(self.meshTriangles[neighbor])
875                else:
876                    ObjectNeighbor.append(None)
877            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],ObjectNeighbor[1],ObjectNeighbor[2])
878            index += 1
879
880
881    def setMesh(self, genDict):
882        """
883        Set the user Mesh attributes given a dictionary of the lists
884        point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
885        point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
886        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
887        segment marker list: [S1Marker, S2Marker, ...] (list of ints)
888        region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
889        region attribute list: ["","reservoir",""] list of strings
890        region max area list:[real, None, Real,...] list of None and reals
891       
892        mesh is an instance of a mesh object
893        """
894
895        #Clear the current user mesh values
896        self.userSegments=[]
897        self.userVertices=[]
898        self.Holes=[]
899        self.Regions=[]
900
901        #print "@#@#@#"
902        #print genDict
903        #print "@#@#@#"
904       
905        index = 0
906        for point in genDict['pointlist']:
907            v=Vertex(point[0], point[1])
908            v.index = index
909            index +=1
910            self.userVertices.append(v)
911
912        index = 0
913        for seg,marker in map(None,genDict['segmentlist'],genDict['segmentmarkerlist']):
914            segObject = Segment( self.userVertices[seg[0]],
915                           self.userVertices[seg[1]], marker = marker )
916            segObject.index = index
917            index +=1
918            self.userSegments.append(segObject)
919
920# Remove the loading of attribute info.
921# Have attribute info added using least_squares in pyvolution
922#         index = 0
923#         for att in genDict['pointattributelist']:
924#             if att == None:
925#                 self.userVertices[index].setAttributes([])
926#             else:
927#                 self.userVertices[index].setAttributes(att)
928#            index += 1
929       
930        index = 0
931        for point in genDict['holelist']:
932            h=Hole(point[0], point[1])
933            h.index = index
934            index +=1
935            self.holes.append(h)
936
937        index = 0
938        for reg,att,maxArea in map(None,
939                                   genDict['regionlist'],
940                                   genDict['regionattributelist'],
941                                   genDict['regionmaxarealist']):
942            Object = Region( reg[0],
943                             reg[1],
944                             tag = att,
945                             maxArea = maxArea)
946            Object.index = index
947            index +=1
948            self.regions.append(Object)
949
950       
951    def TestautoSegment(self):
952        newsegs = []
953        s1 = Segment(self.userVertices[0],
954                               self.userVertices[1])
955        s2 = Segment(self.userVertices[0],
956                               self.userVertices[2])
957        s3 = Segment(self.userVertices[2],
958                               self.userVertices[1])
959        if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
960            newsegs.append(s1)
961        if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
962            newsegs.append(s2)
963        if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
964            newsegs.append(s3)
965        self.userSegments.extend(newsegs)
966        return newsegs
967
968   
969    def savePickle(self, currentName):
970        fd = open(currentName, 'w')
971        pickle.dump(self,fd)
972        fd.close()
973
974    def autoSegment(self):
975        """
976        initially work by running an executable
977        Later compile the c code with a python wrapper.
978
979        Precon: There must be 3 or more vertices in the userVertices structure
980        """
981        newsegs = []
982        inputfile = 'hull_in.txt'
983        outputfile = inputfile + '-alf'
984        #write vertices to file
985        fd = open(inputfile,'w')
986        for v in self.userVertices:
987            fd.write(str(v.x))
988            fd.write(' ')
989            fd.write(str(v.y))
990            fd.write('\n')
991        fd.close()
992
993        # os.system('SET PATH=..\pmesh;%PATH%') # this doesn't work
994        #print "os.getcwd()",os.getcwd()    # this could be usefull
995        #command = "k:\\EQRMsens\\" + s
996        #os.chdir(command)
997        #command = 'env'
998        #os.system(command)
999       
1000        #run hull executable
1001        #warning need to compile hull for the current operating system
1002        command = 'hull.exe -A -i ' + inputfile
1003        os.system(command)
1004       
1005        #read results into this object
1006        fd = open(outputfile)
1007        lines = fd.readlines()
1008        fd.close()
1009        #print "(*(*(*("
1010        #print lines
1011        #print "(*(*(*("
1012        lines.pop(0) #remove the first (title) line
1013        for line in lines:
1014            vertindexs = line.split()
1015            #print 'int(vertindexs[0])', int(vertindexs[0])
1016            #print 'int(vertindexs[1])', int(vertindexs[1])
1017            #print 'self.userVertices[int(vertindexs[0])]' ,self.userVertices[int(vertindexs[0])]
1018            #print 'self.userVertices[int(vertindexs[1])]' ,self.userVertices[int(vertindexs[1])]
1019            v1 = self.userVertices[int(vertindexs[0])]
1020            v2 = self.userVertices[int(vertindexs[1])]
1021           
1022            if self.isUserSegmentNew(v1,v2):
1023                newseg = Segment(v1, v2)
1024                newsegs.append(newseg)
1025        self.userSegments.extend(newsegs)
1026        return newsegs
1027     
1028    def joinVertices(self):
1029        """
1030        Return list of segments connecting all userVertices
1031        in the order they were given
1032       
1033        Precon: There must be 3 or more vertices in the userVertices structure
1034        """
1035
1036        newsegs = []
1037       
1038        v1 = self.userVertices[0]
1039        for v2 in self.userVertices[1:]:
1040            if self.isUserSegmentNew(v1,v2):           
1041                newseg = Segment(v1, v2)
1042                newsegs.append(newseg)
1043            v1 = v2
1044
1045        #Connect last point to the first
1046        v2 = self.userVertices[0]       
1047        if self.isUserSegmentNew(v1,v2):           
1048                newseg = Segment(v1, v2)
1049                newsegs.append(newseg)
1050           
1051
1052        #Update list of user segments       
1053        self.userSegments.extend(newsegs)               
1054        return newsegs
1055   
1056    def normaliseMesh(self,scale, offset, height_scale):
1057        [xmin, ymin, xmax, ymax] = self.boxsize()
1058        [attmin0, attmax0] = self.maxMinVertAtt(0)
1059        #print "[attmin0, attmax0]" ,[attmin0, attmax0]
1060        [attmin1, attmax1] = self.maxMinVertAtt(1)
1061        #print [xmin, ymin, xmax, ymax]
1062        xrange = xmax - xmin
1063        yrange = ymax - ymin
1064        if xrange > yrange:
1065            min,max = xmin, xmax
1066        else:
1067            min,max = ymin, ymax
1068           
1069        for obj in self.getUserVertices():
1070            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1071            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1072            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1073                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1074                                    (attmax0-attmin0)*height_scale
1075            if len(obj.attributes) > 1 and attmin1 != attmax1:
1076                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1077                                    (attmax1-attmin1)*height_scale
1078           
1079        for obj in self.getMeshVertices():
1080            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1081            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1082            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1083                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1084                                    (attmax0-attmin0)*height_scale
1085            if len(obj.attributes) > 1 and attmin1 != attmax1:
1086                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1087                                    (attmax1-attmin1)*height_scale
1088               
1089        for obj in self.getHoles():
1090            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1091            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1092        for obj in self.getRegions():
1093            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1094            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1095        [xmin, ymin, xmax, ymax] = self.boxsize()
1096        #print [xmin, ymin, xmax, ymax]
1097     
1098    def boxsize(self):
1099        """
1100        Returns a list denoting a box that contains the entire structure of vertices
1101        Structure: [xmin, ymin, xmax, ymax]
1102        """
1103        # FIXME dsg!!! large is a hack
1104        #You want the kinds package, part of Numeric:
1105        #In [2]: import kinds
1106       
1107        #In [3]: kinds.default_float_kind.M
1108        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1109    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
1110        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1111
1112        #In [3]: kinds.default_float_kind.MIN
1113        #Out[3]: 2.2250738585072014e-308
1114
1115        large = 1e100
1116        xmin= large
1117        xmax=-large
1118        ymin= large
1119        ymax=-large
1120        for vertex in self.userVertices:
1121            if vertex.x < xmin:
1122                xmin = vertex.x
1123            if vertex.x > xmax:
1124                xmax = vertex.x
1125               
1126            if vertex.y < ymin:
1127                ymin = vertex.y
1128            if vertex.y > ymax:
1129                ymax = vertex.y
1130        return [xmin, ymin, xmax, ymax]
1131 
1132    def maxMinVertAtt(self, iatt):
1133        """
1134        Returns a list denoting a box that contains the entire structure of vertices
1135        Structure: [xmin, ymin, xmax, ymax]
1136        """
1137        # FIXME dsg!!! large is a hack
1138        #You want the kinds package, part of Numeric:
1139        #In [2]: import kinds
1140       
1141        #In [3]: kinds.default_float_kind.M
1142        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1143    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
1144        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1145
1146        #In [3]: kinds.default_float_kind.MIN
1147        #Out[3]: 2.2250738585072014e-308
1148
1149        large = 1e100
1150        min= large
1151        max=-large
1152        for vertex in self.userVertices:
1153            if len(vertex.attributes) > iatt:
1154                if vertex.attributes[iatt] < min:
1155                    min = vertex.attributes[iatt]
1156                if vertex.attributes[iatt] > max:
1157                    max = vertex.attributes[iatt] 
1158        for vertex in self.meshVertices:
1159            if len(vertex.attributes) > iatt:
1160                if vertex.attributes[iatt] < min:
1161                    min = vertex.attributes[iatt]
1162                if vertex.attributes[iatt] > max:
1163                    max = vertex.attributes[iatt] 
1164        return [min, max]
1165   
1166    def scaleoffset(self, WIDTH, HEIGHT):
1167        """
1168        Returns a list denoting the scale and offset terms that need to be
1169        applied when converting  mesh co-ordinates onto grid co-ordinates
1170        Structure: [scale, xoffset, yoffset]
1171        """   
1172        OFFSET = 0.05*min([WIDTH, HEIGHT])
1173        [xmin, ymin, xmax, ymax] = self.boxsize()
1174        SCALE = min([0.9*WIDTH, 0.9*HEIGHT])/max([xmax-xmin, ymax-ymin])
1175       
1176        if SCALE*xmin < OFFSET:
1177            xoffset = abs(SCALE*xmin) + OFFSET
1178        if SCALE*xmax > WIDTH - OFFSET:
1179            xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1180        if SCALE*ymin < OFFSET:
1181            b = abs(SCALE*ymin)+OFFSET
1182        if SCALE*ymax > HEIGHT-OFFSET:
1183            b = -(SCALE*ymax - HEIGHT + OFFSET)
1184        yoffset = HEIGHT - b
1185        return [SCALE, xoffset, yoffset]
1186           
1187    def plotMeshTriangle(self,tag = 0,WIDTH = 400,HEIGHT = 400):
1188        """
1189        Plots all node connections.
1190        tag = 0 (no node numbers), tag = 1 (node numbers)
1191        """
1192
1193        try:
1194            from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
1195
1196            [SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
1197           
1198            root = Tk()
1199            frame = Frame(root)
1200            frame.pack()
1201            button = Button(frame, text="OK", fg="red", command=frame.quit)
1202            button.pack(side=BOTTOM)
1203            canvas = Canvas(frame,bg="white", width=WIDTH, height=HEIGHT)
1204            canvas.pack()
1205            text = Label(frame, width=20, height=10, text='triangle mesh')
1206            text.pack()
1207
1208            #print self.meshTriangles
1209            for triangle in self.meshTriangles:
1210                triangle.draw(canvas,1,
1211                              scale = SCALE,
1212                              xoffset = xoffset,
1213                              yoffset = yoffset )
1214               
1215            root.mainloop()
1216
1217        except:
1218            print "Unexpected error:", sys.exc_info()[0]
1219            raise
1220
1221            #print """
1222            #node::plot Failed.
1223            #Most probably, the Tkinter module is not available.
1224            #"""
1225
1226    def plotUserSegments(self,tag = 0,WIDTH = 400,HEIGHT = 400):
1227        """
1228        Plots all node connections.
1229        tag = 0 (no node numbers), tag = 1 (node numbers)
1230        """
1231
1232        try:
1233            from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
1234           
1235            [SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
1236
1237            root = Tk()
1238            frame = Frame(root)
1239            frame.pack()
1240            button = Button(frame, text="OK", fg="red", command=frame.quit)
1241            button.pack(side=BOTTOM)
1242            canvas = Canvas(frame, bg="white", width=WIDTH, height=HEIGHT)
1243            canvas.pack()
1244            text = Label(frame, width=20, height=10, text='user segments')
1245            text.pack()
1246           
1247            for segment in self.userSegments:
1248                segment.draw(canvas,SCALE, xoffset, yoffset )
1249               
1250            root.mainloop()
1251
1252        except:
1253            print "Unexpected error:", sys.exc_info()[0]
1254            raise
1255
1256            #print """
1257            #node::plot Failed.
1258            #Most probably, the Tkinter module is not available.
1259            #"""
1260     
1261    def exportASCIItrianglulationfile(self,ofile):
1262        """
1263        export a file, ofile, with the format
1264       
1265        First line:  <# of vertices> <# of attributes>
1266        Following lines:  <vertex #> <x> <y> [attributes]
1267        One line:  <# of triangles>
1268        Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
1269        One line:  <# of segments>
1270        Following lines:  <segment #> <vertex #>  <vertex #> [boundary marker]
1271        """
1272        fd = open(ofile,'w')
1273        gen_dict = self.Mesh2MeshList()
1274        load_mesh.loadASCII.write_ASCII_trianglulation(fd,gen_dict)
1275        self.writeASCIImesh(fd,
1276                            self.userVertices,
1277                            self.userSegments,
1278                            self.holes,
1279                            self.regions)   
1280        fd.close()
1281
1282    def exportASCIIsegmentoutlinefile(self,ofile):
1283        """
1284        export a file, ofile, with no triangulation and only vertices connected to segments.
1285        """
1286        fd = open(ofile,'w')
1287        meshDict = {}
1288       
1289        meshDict['generatedpointlist'] = []
1290        meshDict['generatedpointattributelist'] = []
1291        meshDict['generatedsegmentlist'] = []
1292        meshDict['generatedsegmentmarkerlist'] = []
1293
1294        meshDict['generatedtrianglelist'] = [] 
1295        meshDict['generatedtriangleattributelist'] = []
1296        meshDict['generatedtriangleneighborlist'] = []
1297       
1298        load_mesh.loadASCII.write_ASCII_trianglulation(fd,meshDict)
1299        self.writeASCIIsegmentoutline(fd,
1300                            self.userVertices,
1301                            self.userSegments,
1302                            self.holes,
1303                            self.regions)   
1304        fd.close()
1305
1306    def exportASCIIobj(self,ofile):
1307        """
1308        export a file, ofile, with the format
1309         lines:  v <x> <y> <first attribute>
1310        f <vertex #>  <vertex #> <vertex #> (of the triangles)
1311        """
1312        fd = open(ofile,'w')
1313        self.writeASCIIobj(fd)   
1314        fd.close()
1315
1316
1317    def writeASCIIobj(self,fd):
1318        fd.write(" # Triangulation as an obj file\n")
1319        numVert = str(len(self.meshVertices))
1320       
1321        index1 = 1 
1322        for vert in self.meshVertices:
1323            vert.index1 = index1
1324            index1 += 1
1325           
1326            fd.write("v "
1327                     + str(vert.x) + " "
1328                     + str(vert.y) + " "
1329                     + str(vert.attributes[0]) + "\n")
1330           
1331        for tri in self.meshTriangles:
1332            fd.write("f "
1333                     + str(tri.vertices[0].index1) + " " 
1334                     + str(tri.vertices[1].index1) + " " 
1335                     + str(tri.vertices[2].index1) + "\n")
1336               
1337    def writeASCIIsegmentoutline(self,
1338                       fd,
1339                       userVertices,
1340                       userSegments,
1341                       holes,
1342                       regions):
1343        """Write the user mesh info, only with vertices that are connected to segs
1344        """
1345        verts = []
1346        #dupindex = 0
1347        for seg in self.userSegments:
1348            verts.append(seg.vertices[0])
1349            verts.append(seg.vertices[1])
1350        #print 'verts>',verts
1351
1352        verts, count = self.removeDuplicatedVertices(verts)
1353        #print 'verts no dups>',verts
1354        self.writeASCIImesh(fd,
1355                            verts,
1356                            self.userSegments,
1357                            self.holes,
1358                            self.regions)   
1359           
1360    def exportASCIImeshfile(self,ofile):
1361        """
1362        export a file, ofile, with the format
1363       
1364        First line:  <# of user vertices> <# of attributes>
1365        Following lines:  <x> <y> [attributes]
1366        One line:  <# of segments>
1367        Following lines:  <segment #> <vertex #>  <vertex #> [boundary marker]
1368        """
1369       
1370        fd = open(ofile,'w')
1371        self.writeASCIImesh(fd,
1372                            self.userVertices,
1373                            self.userSegments,
1374                            self.holes,
1375                            self.regions)       
1376        fd.close()
1377
1378    def writeASCIImesh(self,
1379                       fd,
1380                       userVertices,
1381                       userSegments,
1382                       holes,
1383                       regions):
1384        numVert = str(len(userVertices))
1385        if (numVert == "0"):
1386            numVertAttrib = "0"
1387        else:
1388            numVertAttrib = str(len(userVertices[0].attributes))
1389        fd.write(numVert + " " + numVertAttrib + " # <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
1390       
1391        # <x> <y> [attributes]
1392        index = 0 
1393        for vert in userVertices:
1394            vert.index = index
1395            index += 1
1396            attlist = ""
1397            for att in vert.attributes:
1398                attlist = attlist + str(att)+" "
1399            attlist.strip()
1400            fd.write(str(vert.index) + " "
1401                     +str(vert.x) + " "
1402                     + str(vert.y) + " "
1403                     + attlist + "\n")
1404           
1405        #One line:  <# of segments>
1406        fd.write(str(len(userSegments)) + 
1407                 " # <segment #> <vertex #>  <vertex #> [boundary marker] ...Mesh Segments..." + "\n")
1408       
1409        #Following lines:  <vertex #>  <vertex #> [boundary marker]
1410        index = 0
1411        for seg in userSegments:
1412            fd.write(str(index) + " "
1413                     + str(seg.vertices[0].index) + " " 
1414                     + str(seg.vertices[1].index) + " " 
1415                     + str(seg.marker) + "\n")
1416            index += 1
1417
1418        #One line:  <# of holes>
1419        fd.write(str(len(holes)) + 
1420                 " # <Hole #> <x> <y> ...Mesh Holes..." + "\n")   
1421        # <x> <y>
1422        index = 0
1423        for h in holes:
1424            fd.write(str(index) + " "
1425                     + str(h.x) + " "
1426                     + str(h.y) + "\n")
1427            index += 1
1428       
1429        #One line:  <# of regions>
1430        fd.write(str(len(regions)) + 
1431                 " # <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")   
1432        # <index> <x> <y> <tag>
1433        index = 0
1434        for r in regions:
1435            fd.write(str(index) + " " 
1436                     + str(r.x) + " "
1437                     + str(r.y)+ " "
1438                     + str(r.getTag()) + "\n")
1439            index += 1
1440        index = 0
1441        # <index> [<MaxArea>|""]
1442        for r in regions:
1443            area = r.getMaxArea()
1444            if area == None:
1445                area = ""
1446            else:
1447                area = str(area)
1448           
1449            fd.write(str(index) + " " + area + "\n")
1450            index += 1
1451               
1452    def exportxyafile(self,ofile):
1453        """
1454        export a file, ofile, with the format
1455       
1456        First line:  <# of vertices> <# of attributes>
1457        Following lines:  <vertex #> <x> <y> [attributes]
1458        """
1459        #load_mesh.loadASCII
1460
1461        if self.meshVertices == []:
1462            Vertices = self.userVertices
1463        else:
1464            Vertices = self.meshVertices
1465           
1466        numVert = str(len(Vertices))
1467       
1468        if Vertices == []:
1469            raise RuntimeError
1470        numVertAttrib = str(len(Vertices[0].attributes))
1471        title = numVert + " " + numVertAttrib + " # <vertex #> <x> <y> [attributes]"
1472
1473        #Convert the Vertices to pointlist and pointattributelist
1474        xya_dict = {}
1475        pointattributes = []
1476        points = []
1477       
1478        for vert in Vertices:
1479            points.append([vert.x,vert.y])
1480            pointattributes.append(vert.attributes)
1481           
1482        xya_dict['pointlist'] = points
1483        xya_dict['pointattributelist'] = pointattributes
1484
1485        load_mesh.loadASCII.export_xya_file(ofile, xya_dict, title, delimiter = " ")
1486
1487
1488def importMeshFromFile(ofile):
1489    """returns a mesh object, made from a .xya or .txh file
1490    Often raises SyntaxError, IOError
1491    """
1492    newmesh = None
1493    if ofile[-4:]== ".xya":
1494        print "loading " + ofile
1495        try:
1496            dict = load_mesh.loadASCII.load_xya_file(ofile, delimiter = ',')
1497        except SyntaxError:
1498            dict = load_mesh.loadASCII.load_xya_file(ofile, delimiter = ' ')
1499        newmesh= Mesh()
1500        newmesh.setMesh(dict)
1501        counter = newmesh.removeDuplicatedUserVertices()
1502        if (counter >0):
1503            print "%i duplicate vertices removed from dataset" % (counter) 
1504    elif ofile[-4:]== ".tsh":
1505        dict = load_mesh.loadASCII.import_trianglulation(ofile)
1506        newmesh= Mesh()
1507        newmesh.setMesh(dict)
1508        newmesh.setTriangulation(dict)
1509    else:
1510        raise RuntimeError
1511    return newmesh
1512
1513def loadPickle(currentName):
1514    fd = open(currentName)
1515    mesh = pickle.load(fd)
1516    fd.close()
1517    return mesh
1518   
1519def square_outline(side_length = 1,up = "top", left = "left", right = "right",
1520                   down = "bottom"):
1521   
1522        a = Vertex (0,0)
1523        b = Vertex (0,side_length)
1524        c = Vertex (side_length,0)
1525        d = Vertex (side_length,side_length)
1526     
1527        s2 = Segment(b,d, marker = up)
1528        s3 = Segment(b,a, marker = left)
1529        s4 = Segment(d,c, marker = right)
1530        s5 = Segment(a,c, marker = down)
1531     
1532        return Mesh(userVertices=[a,b,c,d],
1533                 userSegments=[s2,s3,s4,s5]) 
1534
1535def region_strings2ints(region_list):
1536    """Given a list of (x_int,y_int,marker_string) lists it returns a list of
1537    (x_int,y_int,marker_int) and a list to convert the marker_int's back to
1538    the marker_strings 
1539    """
1540    # Make sure "" has an index of 0
1541    region_list.reverse()
1542    region_list.append((1.0,2.0,""))
1543    region_list.reverse()
1544    convertint2string = []
1545    for i in xrange(len(region_list)):
1546        convertint2string.append(region_list[i][2])
1547        if len(region_list[i]) == 4: # there's an area value
1548            region_list[i] = (region_list[i][0],
1549                              region_list[i][1],i,region_list[i][3])
1550        elif len(region_list[i]) == 3: # no area value
1551            region_list[i] = (region_list[i][0],region_list[i][1],i)
1552        else:
1553            print "The region list has a bad size"
1554            # raise an error ..
1555            raise Error
1556
1557    #remove "" from the region_list
1558    region_list.pop(0)
1559   
1560    return [region_list, convertint2string]
1561
1562def region_ints2strings(region_list,convertint2string):
1563    """Reverses the transformation of region_strings2ints
1564    """
1565    if region_list[0] != []:
1566        for i in xrange(len(region_list)):
1567            region_list[i] = [convertint2string[int(region_list[i][0])]]
1568    return region_list
1569
1570def segment_ints2strings(intlist, convertint2string):
1571    """Reverses the transformation of segment_strings2ints """
1572    stringlist = []
1573    for x in intlist:
1574        stringlist.append(convertint2string[x])
1575    return stringlist
1576
1577def segment_strings2ints(stringlist, preset):
1578    """Given a list of strings return a list of 0 to n ints which represent
1579    the strings and a converting list of the strings, indexed by 0 to n ints.
1580    Also, input an initial converting list of the strings
1581    Note, the converting list starts off with
1582    ["internal boundary", "external boundary"]
1583    example input and output
1584    input ["a","b","a","c"],["c"]
1585    output [[2, 1, 2, 0], ['c', 'b', 'a']]
1586
1587    the first element in the converting list will not be used.  it is
1588    overwritten with ""
1589   
1590    # note, order the initial converting list is important,
1591     since the index = the int marker
1592    """
1593    nodups = unique(stringlist)
1594    # note, order is important, the index = the int marker
1595    #preset = ["internal boundary", "external boundary"]
1596    #Remove the preset tags from the list with no duplicates
1597    nodups = [x for x in nodups if x not in preset]
1598
1599    try:
1600        nodups.remove("") # this has to go to zero
1601    except ValueError:
1602        pass
1603   
1604    # Add the preset list at the beginning of no duplicates
1605    preset.reverse()
1606    nodups.extend(preset)
1607    nodups.reverse()
1608
1609       
1610    convertstring2int = {}
1611    convertint2string = []
1612    index = 0
1613    for x in nodups:
1614        convertstring2int[x] = index
1615        convertint2string.append(x)
1616        index += 1
1617    convertstring2int[""] = 0
1618
1619    intlist = []
1620    for x in stringlist:
1621        intlist.append(convertstring2int[x])
1622    return [intlist, convertint2string]
1623
1624
1625def unique(s):
1626    """Return a list of the elements in s, but without duplicates.
1627
1628    For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
1629    unique("abcabc") some permutation of ["a", "b", "c"], and
1630    unique(([1, 2], [2, 3], [1, 2])) some permutation of
1631    [[2, 3], [1, 2]].
1632
1633    For best speed, all sequence elements should be hashable.  Then
1634    unique() will usually work in linear time.
1635
1636    If not possible, the sequence elements should enjoy a total
1637    ordering, and if list(s).sort() doesn't raise TypeError it's
1638    assumed that they do enjoy a total ordering.  Then unique() will
1639    usually work in O(N*log2(N)) time.
1640
1641    If that's not possible either, the sequence elements must support
1642    equality-testing.  Then unique() will usually work in quadratic
1643    time.
1644    """
1645
1646    n = len(s)
1647    if n == 0:
1648        return []
1649
1650    # Try using a dict first, as that's the fastest and will usually
1651    # work.  If it doesn't work, it will usually fail quickly, so it
1652    # usually doesn't cost much to *try* it.  It requires that all the
1653    # sequence elements be hashable, and support equality comparison.
1654    u = {}
1655    try:
1656        for x in s:
1657            u[x] = 1
1658    except TypeError:
1659        del u  # move on to the next method
1660    else:
1661        return u.keys()
1662
1663    # We can't hash all the elements.  Second fastest is to sort,
1664    # which brings the equal elements together; then duplicates are
1665    # easy to weed out in a single pass.
1666    # NOTE:  Python's list.sort() was designed to be efficient in the
1667    # presence of many duplicate elements.  This isn't true of all
1668    # sort functions in all languages or libraries, so this approach
1669    # is more effective in Python than it may be elsewhere.
1670    try:
1671        t = list(s)
1672        t.sort()
1673    except TypeError:
1674        del t  # move on to the next method
1675    else:
1676        assert n > 0
1677        last = t[0]
1678        lasti = i = 1
1679        while i < n:
1680            if t[i] != last:
1681                t[lasti] = last = t[i]
1682                lasti += 1
1683            i += 1
1684        return t[:lasti]
1685
1686    # Brute force is all that's left.
1687    u = []
1688    for x in s:
1689        if x not in u:
1690            u.append(x)
1691    return u
1692
1693         
1694
1695if __name__ == "__main__":
1696# THIS CAN BE DELETED
1697    list = ["internal boundary","sea","river inlet",
1698            "","sea","","moat","internal boundary"]
1699    #list = ["sea","river inlet","","sea","","moat"]
1700    list = []
1701    print list
1702    preset = [""]
1703    print preset
1704    [intlist, converter] = segment_strings2ints(list, preset)
1705    print "intlist",intlist
1706    print "converter",converter
1707    newlist = segment_ints2strings(intlist, converter)
1708    print newlist
1709    print "***************************"
1710    list = [(4,5,"lo"),(3,2,"mo"),(3,2,"ro")]
1711    #list = []
1712    print "list",list
1713    [newlist, converter] = region_strings2ints(list)
1714    print "newlist",newlist
1715    trilist = []
1716    for x in newlist:
1717        trilist.append([x[2]])
1718        print "trilist",trilist
1719    oldlist = region_ints2strings(trilist, converter) 
1720    print "oldlist",oldlist   
Note: See TracBrowser for help on using the repository browser.