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

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

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