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

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

internal tags nolonger default to 'internal', they now default to not being a tag.

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