source: anuga_core/source/anuga/pmesh/mesh.py @ 4150

Last change on this file since 4150 was 4144, checked in by duncan, 18 years ago

Changing the mesh generator and triangle interface to work with Numeric arrys, as input. Fixed to work on Linux
File size: 127.8 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 persistant 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 sys
17import math
18import re
19import os
20import pickle
21
22import types
23import exceptions
24from Numeric import array, Float, Int
25
26
27#import load_mesh
28from anuga.coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
29from anuga.utilities.polygon import point_in_polygon
30import load_mesh.loadASCII
31import anuga.alpha_shape.alpha_shape
32from anuga.geospatial_data.geospatial_data import Geospatial_data, ensure_geospatial, ensure_absolute, ensure_numeric
33from anuga.mesh_engine.mesh_engine import generate_mesh
34
35#import anuga.mesh_engine_b.mesh_engine as triang
36
37try: 
38    import kinds 
39except ImportError: 
40    # Hand-built mockup of the things we need from the kinds package, since it 
41    # was recently removed from the standard Numeric distro.  Some users may 
42    # not have it by default. 
43    class _bunch: 
44        pass 
45         
46    class _kinds(_bunch): 
47        default_float_kind = _bunch() 
48        default_float_kind.MIN = 2.2250738585072014e-308  #smallest +ve number
49        default_float_kind.MAX = 1.7976931348623157e+308 
50     
51    kinds = _kinds()
52   
53SET_COLOUR='red'
54
55#FIXME: this is not tested.
56from anuga.utilities.numerical_tools import gradient
57
58
59
60# 1st and third values must be the same
61# FIXME: maybe make this a switch that the user can change? - DSG
62initialconversions = ['', 'exterior', '']
63
64#from os import sep
65#sys.path.append('..'+sep+'pmesh')
66#print "sys.path",sys.path
67
68class MeshObject:
69    """
70    An abstract superclass for the basic mesh objects, eg vertex, segment,
71    triangle.
72    """
73    def __init__(self):
74        pass
75   
76class Point(MeshObject): 
77    """
78    Define a point in a 2D space.
79    """
80    def __init__(self,X,Y):
81        __slots__ = ['x','y']
82        self.x=X
83        self.y=Y
84       
85    def DistanceToPoint(self, OtherPoint):
86        """
87        Returns the distance from this point to another
88        """
89        SumOfSquares = ((self.x - OtherPoint.x)**2) + ((self.y - OtherPoint.y)**2)
90        return math.sqrt(SumOfSquares)
91       
92    def IsInsideCircle(self, Center, Radius):
93        """
94        Return 1 if this point is inside the circle,
95        0 otherwise
96        """
97       
98        if (self.DistanceToPoint(Center)<Radius):
99            return 1
100        else:
101            return 0
102       
103    def __repr__(self):
104        return "(%f,%f)" % (self.x,self.y) 
105
106    def cmp_xy(self, point):
107        if self.x < point.x:
108            return -1
109        elif self.x > point.x:
110            return 1
111        else:           
112            if self.y < point.y:
113                return -1
114            elif self.y > point.y:
115                return 1
116            else:
117                return 0
118       
119    def same_x_y(self, point):
120        if self.x == point.x and self.y == point.y:
121            return True
122        else:
123            return False
124       
125           
126
127class Vertex(Point):
128    """
129    A point on the mesh.
130    Object attributes based on the Triangle program
131    """
132    def __init__(self,X,Y, attributes = None):
133        __slots__ = ['x','y','attributes']
134       
135        assert (type(X) == types.FloatType or type(X) == types.IntType)
136        assert (type(Y) == types.FloatType or type(Y) == types.IntType)
137        self.x=X
138        self.y=Y       
139        self.attributes=[] 
140       
141        if attributes is None:
142            self.attributes=[]
143        else:
144            self.attributes=attributes
145   
146
147    def setAttributes(self,attributes):
148        """
149        attributes is a list.
150        """
151        self.attributes = attributes
152       
153    VERTEXSQUARESIDELENGTH = 6
154    def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, yoffset =0, ):
155        x =  scale*(self.x + xoffset)
156        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
157        #print "draw x:", x
158        #print "draw y:", y
159        cornerOffset= self.VERTEXSQUARESIDELENGTH/2
160
161        # A hack to see the vert tags
162        # note: there will be many tags, since tags will not be removed
163        #when zooming
164        #canvas.create_text(x+ 2*cornerOffset,
165        #                   y+ 2*cornerOffset,
166        #                        text=tags)
167       
168        return canvas.create_rectangle(x-cornerOffset,
169                                       y-cornerOffset,
170                                       x+cornerOffset,
171                                       y+cornerOffset,
172                                       tags = tags,
173                                       outline=colour,
174                                       fill = 'white')
175   
176        #return tags
177     
178    def __repr__(self):
179        return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
180   
181class Hole(Point):
182    """
183    A region of the mesh were no triangles are generated.
184    Defined by a point in the hole enclosed by segments.
185    """
186
187    HOLECORNERLENGTH = 6
188   
189    def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ):
190        x =  scale*(self.x + xoffset)
191        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
192        #print "draw x:", x
193        #print "draw y:", y
194        cornerOffset= self.HOLECORNERLENGTH/2
195        return canvas.create_oval(x-cornerOffset,
196                                       y-cornerOffset,
197                                       x+cornerOffset,
198                                       y+cornerOffset,
199                                       tags = tags,
200                                       outline=colour,
201                                       fill = 'white')
202   
203class Region(Point):
204    """
205    A region of the mesh, defined by a point in the region
206    enclosed by segments. Used to tag areas.
207    """
208    CROSSLENGTH = 6
209    TAG = 0
210    MAXAREA = 1
211   
212    def __init__(self,X,Y, tag = None, maxArea = None):
213        """Precondition: tag is a string and maxArea is a real
214        """
215        # This didn't work. 
216        #super(Region,self)._init_(self,X,Y)
217        self.x=X
218        self.y=Y   
219        self.attributes =[] # index 0 is the tag string
220                            #optoinal index 1 is the max triangle area
221                            #NOTE the size of this attribute is assumed
222                            # to be 1 or 2 in regionstrings2int
223        if tag is None:
224            self.attributes.append("")
225        else:
226            self.attributes.append(tag) #this is a string
227           
228        if maxArea is not None:
229            self.setMaxArea(maxArea) # maxArea is a number
230           
231    def getTag(self,):
232        return self.attributes[self.TAG]
233   
234    def setTag(self,tag):
235        self.attributes[self.TAG] = tag
236       
237    def getMaxArea(self):
238        """ Returns the Max Area of a Triangle or
239        None, if the Max Area has not been set.
240        """
241        if self.isMaxArea():
242            return self.attributes[1]
243        else:
244            return None
245   
246    def setMaxArea(self,MaxArea):
247        if MaxArea is not None:
248            if self.isMaxArea(): 
249                self.attributes[self.MAXAREA] = float(MaxArea)
250            else:
251                self.attributes.append( float(MaxArea) )
252   
253    def deleteMaxArea(self):
254        if self.isMaxArea():
255            self.attributes.pop(self.MAXAREA)
256           
257    def isMaxArea(self):
258        return len(self.attributes)> 1
259   
260    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "black"):
261        """
262        Draw a black cross, returning the objectID
263        """
264        x =  scale*(self.x + xoffset)
265        y = -1*scale*(self.y + yoffset) 
266        cornerOffset= self.CROSSLENGTH/2
267        return canvas.create_polygon(x,
268                                     y-cornerOffset,
269                                     x,
270                                     y,
271                                     x+cornerOffset,
272                                     y,
273                                     x,
274                                     y,
275                                     x,
276                                     y+cornerOffset,
277                                     x,
278                                     y,
279                                     x-cornerOffset,
280                                     y,
281                                     x,
282                                     y,
283                                     tags = tags,
284                                     outline = colour,fill = '')
285   
286    def __repr__(self):
287        if self.isMaxArea():
288            area = self.getMaxArea() 
289            return "(%f,%f,%s,%f)" % (self.x,self.y,
290                                      self.getTag(), self.getMaxArea())
291        else:
292            return "(%f,%f,%s)" % (self.x,self.y,
293                                   self.getTag())
294       
295class Triangle(MeshObject):
296    """
297    A triangle element, defined by 3 vertices.
298    Attributes based on the Triangle program.
299    """
300
301    def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ):
302        """
303        Vertices, the initial arguments, are listed in counterclockwise order.
304        """
305        self.vertices= [vertex1,vertex2, vertex3 ]
306       
307        if attribute is None:
308            self.attribute =""
309        else:
310            self.attribute = attribute #this is a string
311           
312        if neighbors is None:
313            self.neighbors=[]
314        else:
315            self.neighbors=neighbors
316
317    def replace(self,new_triangle):
318        self = new_triangle
319
320    def longestSideID(self):
321        ax = self.vertices[0].x
322        ay = self.vertices[0].y
323       
324        bx = self.vertices[1].x
325        by = self.vertices[1].y
326       
327        cx = self.vertices[2].x
328        cy = self.vertices[2].y
329
330        lenA = ((cx-bx)**2+(cy-by)**2)**0.5
331        lenB = ((ax-cx)**2+(ay-cy)**2)**0.5
332        lenC = ((bx-ax)**2+(by-ay)**2)**0.5
333 
334        len = [lenA,lenB,lenC]
335        return len.index(max(len))
336
337    def rotate(self,offset):
338        """
339        permute the order of the sides of the triangle
340        offset must be 0,1 or 2
341        """
342
343        if offset == 0:
344            pass
345        else:
346            if offset == 1:
347                self.vertices = [self.vertices[1],self.vertices[2],self.vertices[0]]
348                self.neighbors = [self.neighbors[1],self.neighbors[2],self.neighbors[0]]
349            if offset == 2:
350                self.vertices = [self.vertices[2],self.vertices[0],self.vertices[1]]
351                self.neighbors = [self.neighbors[2],self.neighbors[0],self.neighbors[1]]
352
353    def rotate_longest_side(self):
354        self.rotate(self.longestSideID())
355
356    def getVertices(self):
357        return self.vertices
358   
359    def get_vertices(self):
360        """
361        Return a list of the vertices.  The x and y values will be relative
362        Easting and Northings for the zone of the current geo_ref.
363        """
364        return self.vertices
365   
366    def calcArea(self):
367        ax = self.vertices[0].x
368        ay = self.vertices[0].y
369       
370        bx = self.vertices[1].x
371        by = self.vertices[1].y
372       
373        cx = self.vertices[2].x
374        cy = self.vertices[2].y
375       
376        return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
377   
378    def calcP(self):
379        #calculate the perimeter
380        ax = self.vertices[0].x
381        ay = self.vertices[0].y
382       
383        bx = self.vertices[1].x
384        by = self.vertices[1].y
385       
386        cx = self.vertices[2].x
387        cy = self.vertices[2].y
388
389        a =  ((cx-bx)**2+(cy-by)**2)**0.5 
390        b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
391        c =  ((bx-ax)**2+(by-ay)**2)**0.5 
392
393        return a+b+c
394           
395    def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
396        """
397        neighbor1 is the triangle opposite vertex1 and so on.
398        Null represents no neighbor
399        """
400        self.neighbors = [neighbor1, neighbor2, neighbor3]
401
402    def setAttribute(self,attribute):
403        """
404        neighbor1 is the triangle opposite vertex1 and so on.
405        Null represents no neighbor
406        """
407        self.attribute = attribute #this is a string
408       
409    def __repr__(self):
410        return "[%s,%s]" % (self.vertices,self.attribute)
411       
412
413    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"):
414        """
415        Draw a triangle, returning the objectID
416        """
417        return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
418                                     scale*-1*(self.vertices[1].y + yoffset),
419                                     scale*(self.vertices[0].x + xoffset),
420                                     scale*-1*(self.vertices[0].y + yoffset),
421                                     scale*(self.vertices[2].x + xoffset),
422                                     scale*-1*(self.vertices[2].y + yoffset),
423                                     tags = tags,
424                                     outline = colour,fill = '')
425
426class Segment(MeshObject):
427    """
428    Segments are edges whose presence in the triangulation is enforced.
429   
430    """
431    def __init__(self, vertex1, vertex2, tag = None ):
432        """
433        Each segment is specified by listing the vertices of its endpoints
434        The vertices are Vertex objects
435        """
436        assert(vertex1 != vertex2)
437        self.vertices = [vertex1,vertex2 ]
438       
439        if tag is None:
440            self.tag = self.__class__.default
441        else:
442            self.tag = tag #this is a string
443       
444    def __repr__(self):
445        return "[%s,%s]" % (self.vertices,self.tag)
446           
447       
448    def draw(self, canvas, tags,scale=1 , xoffset=0 , yoffset=0,colour='blue' ):
449        x=[]
450        y=[]
451        for end in self.vertices:
452            #end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
453            x.append(scale*(end.x + xoffset))
454            y.append(-1*scale*(end.y + yoffset)) # - since for a canvas - is up
455       
456        return canvas.create_line(x[0], y[0], x[1], y[1],
457                                  tags = tags,fill=colour)
458    def set_tag(self,tag):
459        self.tag = tag
460       
461    # Class methods
462    def set_default_tag(cls, default):
463        cls.default = default
464   
465    def get_default_tag(cls):
466        return cls.default
467   
468    set_default_tag = classmethod(set_default_tag) 
469    get_default_tag = classmethod(get_default_tag)
470
471Segment.set_default_tag("")       
472
473class Rigid_triangulation:
474    """
475    This is a triangulation that can't have triangles added or taken away.
476
477    It just represents the triangulation, not the mesh outline needed to
478    build the triangulation.
479
480      This is the guts of the data structure;
481        generated vertex list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
482        generated segment list: [(point1,point2),(p3,p4),...]
483            (Tuples of integers)
484        generated segment tag list: [tag,tag,...] list of strings
485
486        generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
487
488        generated triangle attribute list: [s1,s2,...] list of strings
489
490        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
491            tuple of triangles
492
493            Should I do the basic structure like triangle, general
494            mesh or something else? How about like triangle, since
495            that's where the info is from, and it'll fit easier into
496            this file..
497
498             Removing the loners is difficult, since all the vert's
499             after it must be removed.
500
501             This happens in set_triangulation.
502           
503    """
504
505    def __init__(self,
506                 triangles,
507                 segments,
508                 vertices,
509                 triangle_tags,
510                 triangle_neighbor,
511                 segment_tags,
512                 ):
513       
514        self.triangles=ensure_numeric(triangles) 
515        self.triangle_neighbor=ensure_numeric(triangle_neighbor)
516        self.triangle_tags=triangle_tags
517        self.segments=ensure_numeric(segments) 
518       
519        self.segment_tags=segment_tags # string
520        self.vertices=ensure_numeric(vertices)
521
522       
523       
524class Mesh:
525    """
526    Representation of a 2D triangular mesh.
527    User attributes describe the mesh region/segments/vertices/attributes
528
529    mesh attributes describe the mesh that is produced eg triangles and vertices.
530    All point information is relative to the geo_reference passed in
531   
532   
533    """
534
535    def __repr__(self):
536        return """
537        mesh Triangles: %s 
538        mesh Attribute Titles: %s 
539        mesh Segments: %s 
540        mesh Vertices: %s 
541        user Segments: %s 
542        user Vertices: %s 
543        holes: %s 
544        regions: %s""" % (self.meshTriangles,
545                                self.attributeTitles,
546                                self.meshSegments,
547                                self.meshVertices,
548                                self.getUserSegments(),
549                                self.userVertices,
550                                self.holes,
551                                self.regions) 
552   
553    def __init__(self,
554                 userSegments=None,
555                 userVertices=None,
556                 holes=None,
557                 regions=None,
558                 geo_reference=None):
559        self.meshTriangles=[] 
560        self.attributeTitles=[] 
561        self.meshSegments=[]
562        self.meshVertices=[]
563
564        self.setID={}
565        #a dictionary of names.
566        #multiple sets are allowed, but the gui does not yet
567        #support this
568       
569        self.setID['None']=0
570        #contains the names of the sets pointing to the indexes
571        #in the list.
572       
573        self.sets=[[]]
574        #Contains the lists of triangles (triangle sets)
575
576       
577        self.visualise_graph = True
578
579        if userSegments is None:
580            self.userSegments=[]
581        else:
582            self.userSegments=userSegments
583        self.alphaUserSegments=[]
584           
585        if userVertices is None:
586            self.userVertices=[]
587        else:
588            self.userVertices=userVertices
589           
590        if holes is None:
591            self.holes=[]
592        else:
593            self.holes=holes
594           
595        if regions is None:
596            self.regions=[]
597        else:
598            self.regions=regions
599
600        if geo_reference is None:
601            self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0) 
602        else:
603            self.geo_reference = geo_reference
604
605        self.shape = None
606           
607    def __cmp__(self,other):
608       
609        # A dic for the initial m
610        dic = self.Mesh2triangList()
611        dic_mesh = self.Mesh2MeshList()
612        for element in dic_mesh.keys():
613            dic[element] = dic_mesh[element]
614        for element in dic.keys():
615            dic[element].sort()
616           
617        # A dic for the exported/imported m
618        dic_other = other.Mesh2triangList()
619        dic_mesh = other.Mesh2MeshList()
620        for element in dic_mesh.keys():
621            dic_other[element] = dic_mesh[element]
622        for element in dic.keys():
623            dic_other[element].sort()
624
625        #print "dsg************************8"
626        #print "dic ",dic
627        #print "*******8"
628        #print "mesh",dic_other
629        #print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
630        #print "dsg************************8"
631       
632        return (dic.__cmp__(dic_other))
633   
634    def addUserPoint(self, pointType, x,y):
635        if pointType == Vertex:
636            point = self.addUserVertex(x,y)
637        if pointType == Hole:
638            point = self._addHole(x,y)
639        if pointType == Region:
640            point = self._addRegion(x,y)
641        return point
642   
643    def addUserVertex(self, x,y):
644        v=Vertex(x, y)
645        self.userVertices.append(v)
646        return v
647
648    def _addHole(self, x,y):
649        h=Hole(x, y)
650        self.holes.append(h)
651        return h
652   
653    def add_hole(self, x,y, geo_reference=None):
654        """
655        adds a point, which represents a hole.
656
657        The point data can have it's own geo_refernece.
658        If geo_reference is None the data is asumed to be absolute
659        """
660        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
661                                                 points_geo_ref=geo_reference)
662        return self._addHole(x, y)
663
664    def _addRegion(self, x,y):
665        h=Region(x, y)
666        self.regions.append(h)
667        return h
668   
669    def add_region(self, x,y, geo_reference=None):
670        """
671        adds a point, which represents a region.
672
673        The point data can have it's own geo_refernece.
674        If geo_reference is None the data is asumed to be absolute
675        """
676        #FIXME: have the user set the tag and resolution here,
677        # but still return the instance, just in case.
678        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
679                                                 points_geo_ref=geo_reference)
680        return self._addRegion(x, y)
681
682    # Depreciated
683    def addRegionEN(self, x,y):
684        print "depreciated, use add_region"
685        return self.add_region(x,y)
686
687   
688    def add_vertices(self, point_data):
689        """
690        Add user vertices.
691
692        The point_data can be a list of (x,y) values, a numeric
693        array or a geospatial_data instance.
694        """
695        point_data = ensure_geospatial(point_data)
696        #print "point_data",point_data
697        # get points relative to the mesh geo_ref
698        points = point_data.get_data_points(geo_reference=self.geo_reference)
699   
700        for point in points:
701            v=Vertex(point[0], point[1])
702            self.userVertices.append(v)
703           
704    def add_hole_from_polygon(self, polygon, tags=None, geo_reference=None):
705        """
706        Add a polygon with tags to the current mesh, as a region.
707        The maxArea of the region can be specified.
708
709        If a geo_reference of the polygon points is given, this is used.
710        If not;
711        The x,y info is assumed to be Easting and Northing, absolute,
712        for the meshes zone.
713
714        polygon a list of points, in meters that describe the polygon
715             (e.g. [[x1,y1],[x2,y2],...]
716        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
717
718        This returns the region instance, so if the user whats to modify
719        it they can.       
720        """
721        return self._add_area_from_polygon(polygon, tags=tags,
722                                           hole=True,
723                                           geo_reference=geo_reference
724                                           )
725
726       
727    def add_region_from_polygon(self, polygon, tags=None,
728                                max_triangle_area=None, geo_reference=None):
729        """
730        Add a polygon with tags to the current mesh, as a region.
731        The maxArea of the region can be specified.
732
733        If a geo_reference of the polygon points is given, this is used.
734        If not;
735        The x,y info is assumed to be Easting and Northing, absolute,
736        for the meshes zone.
737
738        polygon a list of points, in meters that describe the polygon
739             (e.g. [[x1,y1],[x2,y2],...]
740        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
741
742        This returns the region instance (if a max_triangle_area is given),
743        so if the user whats to modify it they can.     
744        """
745        if max_triangle_area is None:
746            create_region = False
747        else:
748            create_region = True
749           
750        region = self._add_area_from_polygon(polygon, tags=tags,
751                                             geo_reference=geo_reference,
752                                             region=create_region)
753        if max_triangle_area is not None:
754            region.setMaxArea(max_triangle_area)
755        return region
756   
757       
758    def _add_area_from_polygon(self, polygon, tags=None,
759                               geo_reference=None,
760                               hole=False,
761                               region=False):
762        """
763        Add a polygon with tags to the current mesh, as a region.
764        The maxArea of the region can be specified.
765
766        If a geo_reference of the polygon points is given, this is used.
767        If not;
768        The x,y info is assumed to be Easting and Northing, absolute,
769        for the meshes zone.
770
771        polygon a list of points, in meters that describe the polygon
772             (e.g. [[x1,y1],[x2,y2],...]
773        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
774
775        This returns the region instance, so if the user whats to modify
776        it they can.
777       
778        """
779        #get absolute values
780        if geo_reference is not None:
781            polygon = geo_reference.get_absolute(polygon)
782        # polygon is now absolute
783        #print "polygon  should be absolute",polygon
784       
785        #create points, segs and tags
786        region_dict = {}
787        region_dict['points'] = polygon
788       
789        #Create segments
790        #E.g. [[0,1], [1,2], [2,3], [3,0]]
791        #from polygon
792        #[0,1,2,3]
793        segments = []
794        N = len(polygon)
795        for i in range(N):
796            lo = i
797            hi = (lo + 1) % N
798            segments.append( [lo, hi] ) 
799        region_dict['segments'] = segments
800        region_dict['segment_tags'] = self._tag_dict2list(tags, N)
801       
802   
803        self.addVertsSegs(region_dict) #this is passing absolute values
804
805        if region is True:
806            #get inner point - absolute values
807            inner_point = point_in_polygon(polygon)
808            inner = self.add_region(inner_point[0], inner_point[1],
809                                    geo_reference=None) 
810        elif hole is True:
811            #get inner point - absolute values
812            inner_point = point_in_polygon(polygon)
813            inner = self.add_hole(inner_point[0], inner_point[1],
814                                    geo_reference=None)
815        else:
816            inner = None
817           
818        return inner
819
820    def _tag_dict2list(self, tags, number_of_segs):
821        """
822        Convert a tag dictionary from this sort of format;
823        #{'wall':[0,3],'ocean':[2]}
824
825        To a list format
826        # ['wall', '', 'ocean', 'wall']
827
828        Note: '' is a default value of nothing
829        """
830        # FIXME (DSG-DSG): Using '' as a default isn't good.
831        #Try None.
832        # Due to this default this method is too connected to
833        # _add_area_from_polygon
834       
835        segment_tags = ['']*number_of_segs
836        if tags is not None:
837            for key in tags:
838                indices = tags[key]
839                for i in indices:
840                    segment_tags[i] = key
841        return segment_tags
842       
843    def add_circle(self, center, radius, segment_count=100,
844                   center_geo_reference=None, tag = None,
845                   region=False, hole=False):
846        """
847        center is a point, in absulute or relative to center_geo_ref
848        radius is the radius of the circle
849        segment_count is the number of segments used to represent the circle
850        tag is a string name, given to the segments.
851        If region is True a region object is returned.
852        If hole is True a hole object is returned.
853           (Don't have them both True.)
854
855           
856        """
857        # convert center and radius to a polygon
858        cuts = []
859        factor = 2* math.pi/segment_count
860        for cut in range(segment_count):
861             cuts.append(cut*factor)
862
863        polygon = []
864        for cut in cuts:
865           
866            x = center[0] + radius * math.cos(cut)
867            y = center[1] + radius * math.sin(cut)
868            polygon.append([x,y])
869        # build the tags
870        tags = {tag:range(segment_count)}
871       
872        return self._add_area_from_polygon(polygon, tags=tags,
873                                           region=region, hole=hole,
874                                           geo_reference=center_geo_reference)
875
876    def auto_set_geo_reference(self):
877        """
878        Automatically set the georeference, based in the minimum x and y
879        user vertex values.
880
881        Not to be used with the graphical interface
882       
883        Not implemented.
884        Don't implement now.  using the new georeferenced points class
885        will change this?
886        """
887        #to do
888        # find the lower left hand corner
889        [xmin, ymin, xmax, ymax] = self.boxsize()
890
891        # set all points to that lower left hand corner.
892        # use change_geo_reference
893        new_geo = Geo_reference(self.geo_reference.get_zone(), xmin, ymin)
894        self.change_geo_reference(new_geo)
895       
896    def change_geo_reference(self, new_geo_reference):
897        """
898        Change from one geo_reference to another.
899        Not to be used with the graphical interface
900        """
901        # FIXME
902        # change georeference of;
903        #self.userVertices = self.geo_reference.change_points_geo_ref( \
904        #self.userVertices)
905        #self.holes = self.geo_reference.change_points_geo_ref(self.holes)
906        #self.regions = self.geo_reference.change_points_geo_ref(self.regions)
907        # The above will not work.
908        # since userVertices (etc) is a list of point objects,
909        #not a list of lists.
910        # add a method to the points class to fix this up.
911     
912    def add_segment(self, v1, v2, tag):
913        """
914        Don't do this function.
915        what will the v's be objects?  How is the user suppost to get them?
916
917        Indexes?  If add vertstosegs or add_region_from_polygon is called
918        more than once then the actual index is not obvious.  Making this
919        function confusing.
920        """
921        pass
922
923
924    def add_points_and_segments(self, points,
925                                  segments, segment_tags = None):
926        """
927        Add an outline of the mesh.
928        Vertices is a list of points/ a standard representation of points.
929        Segments is a list of tuples of integers.  Each tuple defines the
930           start and end of the segment by it's vertex index, in relation to
931           the list of vertices.
932        segment_tags is an optional dictionary which is used to add tags to
933           the segments.  The key is the tag name, value is the list of segment
934           indexes the tag will apply to.
935           eg. {'wall':[0,3],'ocean':[2]}
936           
937        """
938        #make sure the points are absolute
939        points = ensure_absolute(points)
940       
941        #create points, segs and tags
942        region_dict = {}
943        region_dict['points'] = points
944        region_dict['segments'] = segments
945        region_dict['segment_tags'] = self._tag_dict2list(segment_tags,
946                                                          len(segments))
947        self.addVertsSegs(region_dict)
948       
949    def addVertsSegs(self, outlineDict):
950        """
951        Add out-line (user Mesh) attributes given a dictionary of the lists
952        points: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
953        segments: [(point1,point2),(p3,p4),...] (Tuples of integers)
954        segment_tags: [S1Tag, S2Tag, ...] (list of strings)
955
956        Assume the values are in Eastings and Northings, with no reference
957        point. eg absolute
958        """
959        if not outlineDict.has_key('segment_tags'):
960            outlineDict['segment_tags'] = []
961            for i in range(len(outlineDict['segments'])):
962                outlineDict['segment_tags'].append('')
963        #print "outlineDict['segment_tags']",outlineDict['segment_tags']
964        #print "outlineDict['points']",outlineDict['points']
965        #print "outlineDict['segments']",outlineDict['segments']
966       
967        i_offset = len(self.userVertices)
968        #print "self.userVertices",self.userVertices
969        #print "index_offset",index_offset
970        for point in outlineDict['points']:
971            v=Vertex(point[0]-self.geo_reference.xllcorner,
972                     point[1]-self.geo_reference.yllcorner)
973            self.userVertices.append(v)
974           
975        for seg,seg_tag in map(None,outlineDict['segments'],
976                       outlineDict['segment_tags']):
977            segObject = Segment(self.userVertices[int(seg[0])+i_offset],
978                           self.userVertices[int(seg[1])+i_offset] )
979            if not seg_tag == '':
980                segObject.set_tag(seg_tag)
981            self.userSegments.append(segObject)
982           
983       
984    def get_triangle_count(self):
985        return len(self.getTriangulation())
986       
987    def getUserVertices(self):
988        """
989        Note: The x,y values will be relative to the mesh geo_reference
990        This returns a list of vertex objects
991        """
992        return self.userVertices
993
994    def get_user_vertices(self, absolute=True):
995        """
996        This returns a list of (x,y) values
997        (maybe it should return a geospatical object?
998        It makes it a bit confusing though.)
999        """
1000        pointlist=[]
1001        for vertex in self.userVertices:
1002            pointlist.append([vertex.x,vertex.y])
1003        spat = Geospatial_data(pointlist, geo_reference=self.geo_reference)
1004        return spat.get_data_points(absolute=absolute)
1005   
1006    def getUserSegments(self):
1007        allSegments = self.userSegments + self.alphaUserSegments
1008        #print "self.userSegments",self.userSegments
1009        #print "self.alphaUserSegments",self.alphaUserSegments
1010        #print "allSegments",allSegments
1011        return allSegments
1012   
1013    def deleteUserSegments(self,seg):
1014        if self.userSegments.count(seg) == 0:
1015            self.alphaUserSegments.remove(seg)
1016            pass
1017        else:
1018            self.userSegments.remove(seg)
1019           
1020    def clearUserSegments(self):
1021        self.userSegments = []
1022        self.alphaUserSegments = []
1023               
1024    def getTriangulation(self):
1025        return self.meshTriangles
1026   
1027    def getMeshVertices(self):
1028        return self.meshVertices
1029 
1030    def getMeshSegments(self):
1031        return self.meshSegments
1032   
1033    def getHoles(self):
1034        return self.holes
1035   
1036    def getRegions(self):
1037        return self.regions
1038   
1039    def isTriangulation(self):
1040        if self.meshVertices == []:
1041            return False 
1042        else:
1043            return True
1044   
1045    def addUserSegment(self, v1,v2):
1046        """
1047        PRECON: A segment between the two vertices is not already present.
1048        Check by calling isUserSegmentNew before calling this function.
1049       
1050        """
1051        s=Segment( v1,v2)
1052        self.userSegments.append(s)
1053        return s
1054       
1055    def generate_mesh(self,
1056                      maximum_triangle_area="",
1057                      minimum_triangle_angle=28.0,
1058                      verbose=True):
1059        if verbose is True:
1060            silent = ''
1061        else:
1062            silent = 'Q'
1063        self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle)
1064                                  +"a"+str(maximum_triangle_area)
1065                                  +"a")
1066        #The last a is so areas for regions will be used
1067       
1068    def generateMesh(self, mode = None, maxArea = None, minAngle= None,
1069                     isRegionalMaxAreas = True):
1070        """
1071        Based on the current user vaules, holes and regions
1072        generate a new mesh
1073        mode is a string that sets conditions on the mesh generations
1074        see triangle_instructions.txt for a definition of the commands
1075        PreCondition: maxArea is a double
1076        """
1077        #print "mode ",mode
1078        if mode == None:
1079            self.mode = ""
1080        else:
1081            self.mode = mode
1082       
1083        if self.mode.find('p') < 0:
1084            self.mode += 'p' #p - read a planar straight line graph.
1085            #there must be segments to use this switch
1086            # TODO throw an aception if there aren't seg's
1087            # it's more comlex than this.  eg holes
1088        if self.mode.find('z') < 0:
1089            self.mode += 'z' # z - Number all items starting from zero
1090                             # (rather than one)
1091        if self.mode.find('n'):
1092            self.mode += 'n' # n - output a list of neighboring triangles
1093        if self.mode.find('A') < 0:
1094            self.mode += 'A' # A - output region attribute list for triangles
1095
1096        if not self.mode.find('V') and not self.mode.find('Q'):
1097            self.mode += 'V' # V - output info about what Triangle is doing
1098       
1099        if self.mode.find('q') < 0 and minAngle is not None:
1100            #   print "**********8minAngle******** ",minAngle
1101            min_angle = 'q' + str(minAngle)
1102            self.mode += min_angle # z - Number all items starting from zero
1103                             # (rather than one)
1104        if maxArea != None:
1105            self.mode += 'a' + str(maxArea)
1106        #FIXME (DSG-DSG) This isn't explained.
1107        if isRegionalMaxAreas:
1108            self.mode += 'a'
1109        #print "mesh#generateMesh# self.mode",self.mode 
1110        meshDict = self.Mesh2triangList()
1111
1112        #FIXME (DSG-DSG)  move below section into generate_mesh.py
1113        #                  & 4 functions eg segment_strings2ints
1114        #print "*************************!@!@ This is going to triangle   !@!@"
1115        #print meshDict
1116        #print "************************!@!@ This is going to triangle   !@!@"
1117
1118        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist']
1119        [meshDict['segmenttaglist'],
1120         segconverter] =  segment_strings2ints(meshDict['segmenttaglist'],
1121                                             initialconversions)
1122        #print "regionlist",meshDict['regionlist']
1123        [meshDict['regionlist'],
1124         regionconverter] =  region_strings2ints(meshDict['regionlist'])
1125        #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%regionlist",meshDict['regionlist']
1126        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist'
1127        generatedMesh = generate_mesh(
1128                              meshDict['pointlist'],
1129                              meshDict['segmentlist'],
1130                              meshDict['holelist'],
1131                              meshDict['regionlist'],
1132                              meshDict['pointattributelist'],
1133                              meshDict['segmenttaglist'],
1134                              self.mode,
1135                              meshDict['pointlist'])
1136        #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%generated",generatedMesh
1137        generatedMesh['generatedsegmentmarkerlist'] = \
1138             segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
1139                                  segconverter)
1140        #print "processed gen",generatedMesh['generatedsegmentmarkerlist']
1141        #print "generatedtriangleattributelist",generatedMesh['generatedtriangleattributelist']
1142        generatedMesh['generatedtriangleattributelist'] = \
1143         region_ints2strings(generatedMesh['generatedtriangleattributelist'],
1144                                  regionconverter)
1145
1146        #FIXME (DSG-DSG)  move above section into generate_mesh.py
1147
1148        if len(generatedMesh['generatedpointattributelist'][0])==0:
1149            self.attributeTitles = []
1150        generatedMesh['generatedpointattributetitlelist']= \
1151                                            self.attributeTitles
1152        #print "################  FROM TRIANGLE"
1153        #print "generatedMesh",generatedMesh
1154        #print "##################"
1155       
1156        self.setTriangulation(generatedMesh)
1157   
1158    def clearTriangulation(self):
1159
1160        #Clear the current generated mesh values
1161        self.meshTriangles=[] 
1162        self.meshSegments=[]
1163        self.meshVertices=[]
1164
1165    def removeDuplicatedUserVertices(self):
1166        """Pre-condition: There are no user segments
1167        This function will keep the first duplicate
1168        """
1169        assert self.getUserSegments() == []
1170        self.userVertices, counter =  self.removeDuplicatedVertices(
1171            self.userVertices)
1172        return counter
1173   
1174    def removeDuplicatedVertices(self, Vertices):
1175        """
1176        This function will keep the first duplicate, remove all others
1177        Precondition: Each vertex has a dupindex, which is the list
1178        index.
1179
1180        Note: this removes vertices that have the same x,y values,
1181        not duplicate instances in the Vertices list.
1182        """
1183        remove = []
1184        index = 0
1185        for v in Vertices:
1186            v.dupindex = index
1187            index += 1
1188        t = list(Vertices)
1189        t.sort(Point.cmp_xy)
1190   
1191        length = len(t)
1192        behind = 0
1193        ahead  = 1
1194        counter = 0
1195        while ahead < length:
1196            b = t[behind]
1197            ah = t[ahead]
1198            if (b.y == ah.y and b.x == ah.x):
1199                remove.append(ah.dupindex) 
1200            behind += 1
1201            ahead += 1
1202
1203        # remove the duplicate vertices
1204        remove.sort()
1205        remove.reverse()
1206        for i in remove:
1207            Vertices.pop(i)
1208            pass
1209
1210        #Remove the attribute that this function added
1211        for v in Vertices:
1212            del v.dupindex
1213        return Vertices,counter
1214   
1215    def thinoutVertices(self, delta):
1216        """Pre-condition: There are no user segments
1217        This function will keep the first duplicate
1218        """
1219        assert self.getUserSegments() == []
1220        #t = self.userVertices
1221        #self.userVertices =[]
1222        boxedVertices = {}
1223        thinnedUserVertices =[]
1224        delta = round(delta,1)
1225       
1226        for v in self.userVertices :
1227            # tag is the center of the boxes
1228            tag = (round(v.x/delta,0)*delta,round(v.y/delta,0)*delta)
1229            #this creates a dict of lists of faces, indexed by tag
1230            boxedVertices.setdefault(tag,[]).append(v)
1231
1232        for [tag,verts] in boxedVertices.items():
1233            min = delta
1234            bestVert = None
1235            tagVert = Vertex(tag[0],tag[1])
1236            for v in verts:
1237                dist = v.DistanceToPoint(tagVert)
1238                if (dist<min):
1239                    min = dist
1240                    bestVert = v
1241            thinnedUserVertices.append(bestVert)
1242        self.userVertices =thinnedUserVertices
1243       
1244           
1245    def isUserSegmentNew(self, v1,v2):
1246        identicalSegs= [x for x in self.getUserSegments() \
1247                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1248        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1249       
1250        return len(identicalSegs) == 0
1251
1252       
1253    def deleteSegsOfVertex(self, delVertex):
1254        """
1255        Delete this vertex and any segments that connect to it.
1256        """
1257        #Find segments that connect to delVertex
1258        deletedSegments = []
1259        for seg in self.getUserSegments():
1260            if (delVertex in seg.vertices):
1261                deletedSegments.append(seg)
1262        # Delete segments that connect to delVertex
1263        for seg in deletedSegments:
1264            self.deleteUserSegments(seg)
1265        return deletedSegments
1266
1267   
1268    def deleteMeshObject(self, MeshObject):
1269        """
1270        Returns a list of all objects that were removed
1271        """
1272        deletedObs = []
1273        if isinstance(MeshObject, Vertex ):
1274            deletedObs = self.deleteSegsOfVertex(MeshObject)
1275            deletedObs.append(MeshObject)
1276            self.userVertices.remove(MeshObject)
1277        elif isinstance(MeshObject, Segment):
1278            deletedObs.append(MeshObject)
1279            self.deleteUserSegments(MeshObject)
1280        elif isinstance(MeshObject, Hole):
1281            deletedObs.append(MeshObject)
1282            self.holes.remove(MeshObject)
1283        elif isinstance(MeshObject, Region):
1284            deletedObs.append(MeshObject)
1285            self.regions.remove(MeshObject)         
1286        return deletedObs
1287                                                 
1288    def Mesh2triangList(self, userVertices=None,
1289                        userSegments=None,
1290                        holes=None,
1291                        regions=None):
1292        """
1293        Convert the Mesh to a dictionary of the lists needed for the
1294        triang module
1295        points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1296        pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
1297        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1298        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole)
1299        regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
1300       
1301        Note, this adds an index attribute to the user Vertex objects.
1302
1303        Used to produce output to triangle
1304        """
1305        if userVertices is None:
1306            userVertices = self.getUserVertices()
1307        if userSegments is None:
1308            userSegments = self.getUserSegments()
1309        if holes is None:
1310            holes = self.getHoles()
1311        if regions is None:
1312            regions = self.getRegions()
1313           
1314        meshDict = {}
1315       
1316        pointlist=[]
1317        pointattributelist=[]
1318        index = 0
1319        for vertex in userVertices:
1320            vertex.index = index
1321            pointlist.append((vertex.x,vertex.y))
1322            pointattributelist.append((vertex.attributes))
1323           
1324            index += 1
1325        meshDict['pointlist'] = pointlist
1326        meshDict['pointattributelist'] = pointattributelist
1327
1328        segmentlist=[]
1329        segmenttaglist=[]
1330        for seg in userSegments:
1331            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1332            segmenttaglist.append(seg.tag)
1333        meshDict['segmentlist'] =segmentlist
1334        meshDict['segmenttaglist'] =segmenttaglist
1335       
1336        holelist=[]
1337        for hole in holes:
1338            holelist.append((hole.x,hole.y)) 
1339        meshDict['holelist'] = holelist
1340       
1341        regionlist=[]
1342        for region in regions:
1343            if (region.getMaxArea() != None): 
1344                regionlist.append((region.x,region.y,region.getTag(),
1345                               region.getMaxArea()))
1346            else:
1347                regionlist.append((region.x,region.y,region.getTag()))
1348        meshDict['regionlist'] = regionlist
1349        #print "*(*("
1350        #print meshDict
1351        #print meshDict['regionlist']
1352        #print "*(*("
1353        return meshDict
1354                                               
1355    def Mesh2MeshList(self):
1356        """
1357        Convert the Mesh to a dictionary of lists describing the
1358        triangulation variables;
1359        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1360        generated point attribute list: [(a11,a12,...),(a21,a22),...]
1361            (Tuples of doubles)
1362        generated point attribute title list:[A1Title, A2Title ...]
1363            (list of strings)
1364        generated segment list: [(point1,point2),(p3,p4),...]
1365            (Tuples of integers)
1366        generated segment tag list: [tag,tag,...] list of strings
1367
1368        generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
1369
1370        generated triangle attribute list: [s1,s2,...] list of strings
1371
1372        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
1373            tuple of triangles
1374       
1375        Used to produce .tsh file
1376        """
1377       
1378        meshDict = {}       
1379        pointlist=[]
1380        pointattributelist=[]
1381
1382
1383        self.maxVertexIndex=0
1384        for vertex in self.meshVertices:
1385            vertex.index = self.maxVertexIndex
1386            pointlist.append((vertex.x,vertex.y))
1387            pointattributelist.append((vertex.attributes))           
1388            self.maxVertexIndex += 1
1389
1390        meshDict['generatedpointlist'] = pointlist
1391        meshDict['generatedpointattributelist'] = pointattributelist
1392        meshDict['generatedpointattributetitlelist'] = self.attributeTitles
1393        #segments
1394        segmentlist=[]
1395        segmenttaglist=[]
1396        for seg in self.meshSegments:
1397            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1398            segmenttaglist.append(seg.tag)
1399        meshDict['generatedsegmentlist'] =segmentlist
1400        meshDict['generatedsegmenttaglist'] =segmenttaglist
1401
1402        # Make sure that the indexation is correct
1403        index = 0
1404        for tri in self.meshTriangles:
1405            tri.index = index           
1406            index += 1
1407
1408        trianglelist = []
1409        triangleattributelist = []
1410        triangleneighborlist = []
1411        for tri in self.meshTriangles: 
1412            trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,
1413                                 tri.vertices[2].index)) 
1414            triangleattributelist.append([tri.attribute])
1415            neighborlist = [-1,-1,-1]
1416            for neighbor,index in map(None,tri.neighbors,
1417                                      range(len(tri.neighbors))):
1418                if neighbor:   
1419                    neighborlist[index] = neighbor.index
1420            triangleneighborlist.append(neighborlist)
1421       
1422        meshDict['generatedtrianglelist'] = trianglelist
1423        meshDict['generatedtriangleattributelist'] = triangleattributelist
1424        meshDict['generatedtriangleneighborlist'] = triangleneighborlist
1425       
1426        #print "mesh.Mesh2MeshList*)*)"
1427        #print meshDict
1428        #print "mesh.Mesh2MeshList*)*)"
1429
1430        return meshDict
1431
1432                               
1433    def Mesh2MeshDic(self):
1434        """
1435        Convert the user and generated info of a mesh to a dictionary
1436        structure
1437        """
1438        dic = self.Mesh2triangList()
1439        dic_mesh = self.Mesh2MeshList()
1440        for element in dic_mesh.keys():
1441            dic[element] = dic_mesh[element]
1442        return dic
1443   
1444    def setTriangulation(self, genDict):
1445        """
1446        Set the mesh attributes given a dictionary of the lists
1447        returned from the triang module       
1448        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1449        generated point attribute list:[(P1att1,P1attt2, ...),
1450            (P2att1,P2attt2,...),...]
1451        generated point attribute title list:[A1Title, A2Title ...]
1452            (list of strings)
1453        generated segment list: [(point1,point2),(p3,p4),...]
1454            (Tuples of integers)
1455        generated segment marker list: [S1Tag, S2Tag, ...] (list of ints)
1456        triangle list:  [(point1,point2, point3),(p5,p4, p1),...]
1457            (Tuples of integers)
1458        triangle neighbor list: [(triangle1,triangle2, triangle3),
1459            (t5,t4, t1),...] (Tuples of integers)
1460            -1 means there's no triangle neighbor
1461        triangle attribute list: [(T1att), (T2att), ...]
1462
1463
1464            (list of a list of strings)
1465        lone point list:[point1, ...] (list of integers)
1466        """
1467        #Clear the current generated mesh values
1468        self.meshTriangles=[]
1469        self.attributeTitles=[]
1470        self.meshSegments=[]
1471        self.meshVertices=[]
1472
1473        #print "mesh.setTriangulation@#@#@#"
1474        #print genDict
1475        #print "@#@#@#"
1476
1477        self.maxVertexIndex = 0
1478        for point in genDict['generatedpointlist']:
1479            v=Vertex(point[0], point[1])
1480            v.index =  self.maxVertexIndex
1481            self.maxVertexIndex +=1
1482            self.meshVertices.append(v)
1483
1484        self.attributeTitles = genDict['generatedpointattributetitlelist']
1485
1486        index = 0
1487        for seg,marker in map(None,genDict['generatedsegmentlist'],
1488                              genDict['generatedsegmentmarkerlist']):
1489            segObject = Segment( self.meshVertices[int(seg[0])],
1490                           self.meshVertices[int(seg[1])], tag = marker )
1491            segObject.index = index
1492            index +=1
1493            self.meshSegments.append(segObject)
1494
1495        index = 0
1496        for triangle in genDict['generatedtrianglelist']:
1497            tObject =Triangle( self.meshVertices[int(triangle[0])],
1498                        self.meshVertices[int(triangle[1])],
1499                        self.meshVertices[int(triangle[2])] )
1500            tObject.index = index
1501            index +=1
1502            self.meshTriangles.append(tObject)
1503
1504        index = 0
1505        for att in genDict['generatedtriangleattributelist']:
1506            if att == []:
1507                self.meshTriangles[index].setAttribute("")
1508            else:
1509                self.meshTriangles[index].setAttribute(att[0])
1510            index += 1
1511           
1512        index = 0
1513        for att in genDict['generatedpointattributelist']:
1514            if att == None:
1515                self.meshVertices[index].setAttributes([])
1516            else:
1517                self.meshVertices[index].setAttributes(att)
1518            index += 1
1519   
1520        index = 0
1521        for triangle in genDict['generatedtriangleneighborlist']:
1522            # Build a list of triangle object neighbors
1523            ObjectNeighbor = []
1524            for neighbor in triangle:
1525                if ( neighbor != -1):
1526                    ObjectNeighbor.append(self.meshTriangles[int(neighbor)])
1527                else:
1528                    ObjectNeighbor.append(None)
1529            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
1530                                                   ObjectNeighbor[1],
1531                                                   ObjectNeighbor[2])
1532            index += 1
1533        genDict['lonepointlist'].sort() 
1534        genDict['lonepointlist'].reverse()       
1535        for loner in genDict['lonepointlist']:
1536            # Remove the loner vertex
1537            #print "Removing the loner", loner
1538            self.meshVertices.pop(loner)
1539
1540    def setMesh(self, genDict):
1541        """
1542        Set the user Mesh attributes given a dictionary of the lists
1543        point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1544        point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1545        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1546        segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1547        region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1548        region attribute list: ["","reservoir",""] list of strings
1549        region max area list:[real, None, Real,...] list of None and reals
1550       
1551        mesh is an instance of a mesh object
1552        """
1553        #Clear the current user mesh values
1554        self.clearUserSegments()
1555        self.userVertices=[]
1556        self.Holes=[]
1557        self.Regions=[]
1558
1559        #print "mesh.setMesh@#@#@#"
1560        #print genDict
1561        #print "@#@#@#"
1562       
1563        #index = 0
1564        for point in genDict['pointlist']:
1565            v=Vertex(point[0], point[1])
1566            #v.index = index
1567            #index +=1
1568            self.userVertices.append(v)
1569
1570        #index = 0
1571        for seg,tag in map(None,genDict['segmentlist'],
1572                           genDict['segmenttaglist']):
1573            segObject = Segment( self.userVertices[int(seg[0])],
1574                           self.userVertices[int(seg[1])], tag = tag )
1575            #segObject.index = index
1576            #index +=1
1577            self.userSegments.append(segObject)
1578
1579# Remove the loading of attribute info.
1580# Have attribute info added using least_squares in pyvolution
1581#         index = 0
1582#         for att in genDict['pointattributelist']:
1583#             if att == None:
1584#                 self.userVertices[index].setAttributes([])
1585#             else:
1586#                 self.userVertices[index].setAttributes(att)
1587#            index += 1
1588       
1589        #index = 0
1590        for point in genDict['holelist']:
1591            h=Hole(point[0], point[1])
1592            #h.index = index
1593            #index +=1
1594            self.holes.append(h)
1595
1596        #index = 0
1597        for reg,att,maxArea in map(None,
1598                                   genDict['regionlist'],
1599                                   genDict['regionattributelist'],
1600                                   genDict['regionmaxarealist']):
1601            Object = Region( reg[0],
1602                             reg[1],
1603                             tag = att,
1604                             maxArea = maxArea)
1605            #Object.index = index
1606            #index +=1
1607            self.regions.append(Object)
1608           
1609    def Testauto_segment(self):
1610        newsegs = []
1611        s1 = Segment(self.userVertices[0],
1612                               self.userVertices[1])
1613        s2 = Segment(self.userVertices[0],
1614                               self.userVertices[2])
1615        s3 = Segment(self.userVertices[2],
1616                               self.userVertices[1])
1617        if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1618            newsegs.append(s1)
1619        if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1620            newsegs.append(s2)
1621        if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1622            newsegs.append(s3)
1623        #DSG!!!
1624        self.userSegments.extend(newsegs)
1625        return newsegs
1626
1627   
1628    def savePickle(self, currentName):
1629        fd = open(currentName, 'w')
1630        pickle.dump(self,fd)
1631        fd.close()
1632
1633    def auto_segmentFilter(self,raw_boundary=True,
1634                    remove_holes=False,
1635                    smooth_indents=False,
1636                    expand_pinch=False):
1637        """
1638        Change the filters applied on the alpha shape boundary
1639        """
1640        if self.shape is None:
1641            return [],[],0.0
1642        return self._boundary2mesh(raw_boundary=raw_boundary,
1643                    remove_holes=remove_holes,
1644                    smooth_indents=smooth_indents,
1645                    expand_pinch=expand_pinch)
1646   
1647       
1648   
1649    def auto_segment(self, alpha = None,
1650                    raw_boundary=True,
1651                    remove_holes=False,
1652                    smooth_indents=False,
1653                    expand_pinch=False): 
1654        """
1655        Precon: There must be 3 or more vertices in the userVertices structure
1656        """
1657        self._createBoundary(alpha=alpha)
1658        return self._boundary2mesh(raw_boundary=raw_boundary,
1659                    remove_holes=remove_holes,
1660                    smooth_indents=smooth_indents,
1661                    expand_pinch=expand_pinch)
1662
1663    def _createBoundary(self,alpha=None):
1664        """
1665        """
1666        points=[]
1667        for vertex in self.getUserVertices():
1668            points.append((vertex.x,vertex.y))
1669        self.shape = anuga.alpha_shape.alpha_shape.Alpha_Shape(points, alpha = alpha)
1670
1671
1672    def _boundary2mesh(self, raw_boundary=True,
1673                    remove_holes=False,
1674                    smooth_indents=False,
1675                    expand_pinch=False):
1676        """
1677        Precon there must be a shape object.
1678        """
1679        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1680                                 remove_holes=remove_holes,
1681                                 smooth_indents=smooth_indents,
1682                                 expand_pinch=expand_pinch)
1683        boundary_segs = self.shape.get_boundary()
1684        #print "boundary_segs",boundary_segs
1685        segs2delete = self.alphaUserSegments
1686        #FIXME(DSG-DSG) this algorithm needs comments
1687        new_segs = {}
1688        #alpha_segs = []
1689        #user_segs = []
1690        for seg in boundary_segs:
1691            v1 = self.userVertices[int(seg[0])]
1692            v2 = self.userVertices[int(seg[1])]
1693            boundary_seg = Segment(v1, v2)
1694            new_segs[(v1,v2)] = boundary_seg
1695
1696        for user_seg in self.userSegments:
1697            if new_segs.has_key((user_seg.vertices[0],
1698                                user_seg.vertices[1])):
1699                del new_segs[user_seg.vertices[0],
1700                                user_seg.vertices[1]]
1701            elif new_segs.has_key((user_seg.vertices[1],
1702                                user_seg.vertices[0])):
1703                del new_segs[user_seg.vertices[1],
1704                                user_seg.vertices[0]]
1705               
1706        optimum_alpha = self.shape.get_alpha()
1707        alpha_segs_no_user_segs  = new_segs.values()
1708        self.alphaUserSegments = alpha_segs_no_user_segs
1709        return alpha_segs_no_user_segs, segs2delete, optimum_alpha
1710   
1711    def _boundary2mesh_old(self, raw_boundary=True,
1712                    remove_holes=False,
1713                    smooth_indents=False,
1714                    expand_pinch=False):
1715        """
1716        Precon there must be a shape object.
1717        """
1718        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1719                                 remove_holes=remove_holes,
1720                                 smooth_indents=smooth_indents,
1721                                 expand_pinch=expand_pinch)
1722        boundary_segs = self.shape.get_boundary()
1723        #print "boundary_segs",boundary_segs
1724        segs2delete = self.alphaUserSegments
1725
1726        #FIXME(DSG-DSG) this algorithm needs comments
1727        #FIXME(DSG-DSG) can it be sped up?  It's slow
1728        new_segs = []
1729        alpha_segs = []
1730        user_segs = []
1731        for seg in boundary_segs:
1732            v1 = self.userVertices[int(seg[0])]
1733            v2 = self.userVertices[int(seg[1])]
1734            alpha_seg = self.representedAlphaUserSegment(v1, v2)
1735            user_seg = self.representedUserSegment(v1, v2)
1736            #DSG!!!
1737            assert not(not (alpha_seg == None) and not (user_seg == None))
1738            if not alpha_seg == None:
1739                alpha_segs.append(alpha_seg)
1740            elif not user_seg  == None:
1741                user_segs.append(user_seg)
1742            else:
1743                unique_seg = Segment(v1, v2)
1744                new_segs.append(unique_seg)
1745               
1746            for seg in alpha_segs:
1747                try:
1748                    segs2delete.remove(seg)
1749                except:
1750                    pass
1751       
1752        self.alphaUserSegments = []
1753        self.alphaUserSegments.extend(new_segs)
1754        self.alphaUserSegments.extend(alpha_segs)
1755
1756        optimum_alpha = self.shape.get_alpha()
1757        # need to draw newsegs
1758        return new_segs, segs2delete, optimum_alpha
1759   
1760    def representedAlphaUserSegment(self, v1,v2):
1761        identicalSegs= [x for x in self.alphaUserSegments \
1762                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1763        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1764
1765        if identicalSegs == []:
1766            return None
1767        else:
1768            # Only return the first one.
1769            return identicalSegs[0]
1770   
1771    def representedUserSegment(self, v1,v2):
1772        identicalSegs= [x for x in self.userSegments \
1773                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1774        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1775
1776        if identicalSegs == []:
1777            return None
1778        else:
1779            # Only return the first one.
1780            return identicalSegs[0]
1781       
1782    def joinVertices(self):
1783        """
1784        Return list of segments connecting all userVertices
1785        in the order they were given
1786       
1787        Precon: There must be 3 or more vertices in the userVertices structure
1788        """
1789
1790        newsegs = []
1791       
1792        v1 = self.userVertices[0]
1793        for v2 in self.userVertices[1:]:
1794            if self.isUserSegmentNew(v1,v2):           
1795                newseg = Segment(v1, v2)
1796                newsegs.append(newseg)
1797            v1 = v2
1798
1799        #Connect last point to the first
1800        v2 = self.userVertices[0]       
1801        if self.isUserSegmentNew(v1,v2):           
1802                newseg = Segment(v1, v2)
1803                newsegs.append(newseg)
1804           
1805
1806        #Update list of user segments       
1807        #DSG!!!
1808        self.userSegments.extend(newsegs)               
1809        return newsegs
1810   
1811    def normaliseMesh(self,scale, offset, height_scale):
1812        [xmin, ymin, xmax, ymax] = self.boxsize()
1813        [attmin0, attmax0] = self.maxMinVertAtt(0)
1814        #print "[attmin0, attmax0]" ,[attmin0, attmax0]
1815        [attmin1, attmax1] = self.maxMinVertAtt(1)
1816        #print [xmin, ymin, xmax, ymax]
1817        xrange = xmax - xmin
1818        yrange = ymax - ymin
1819        if xrange > yrange:
1820            min,max = xmin, xmax
1821        else:
1822            min,max = ymin, ymax
1823           
1824        for obj in self.getUserVertices():
1825            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1826            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1827            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1828                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1829                                    (attmax0-attmin0)*height_scale
1830            if len(obj.attributes) > 1 and attmin1 != attmax1:
1831                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1832                                    (attmax1-attmin1)*height_scale
1833           
1834        for obj in self.getMeshVertices():
1835            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1836            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1837            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1838                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1839                                    (attmax0-attmin0)*height_scale
1840            if len(obj.attributes) > 1 and attmin1 != attmax1:
1841                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1842                                    (attmax1-attmin1)*height_scale
1843               
1844        for obj in self.getHoles():
1845            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1846            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1847        for obj in self.getRegions():
1848            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1849            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1850        [xmin, ymin, xmax, ymax] = self.boxsize()
1851        #print [xmin, ymin, xmax, ymax]
1852     
1853    def boxsizeVerts(self):
1854        """
1855        Returns a list of verts denoting a box or triangle that contains
1856        verts on the xmin, ymin, xmax and ymax axis.
1857        Structure: list of verts
1858        """
1859       
1860        large = kinds.default_float_kind.MAX
1861        xmin= large
1862        xmax=-large
1863        ymin= large
1864        ymax=-large
1865        for vertex in self.userVertices:
1866            if vertex.x < xmin:
1867                xmin = vertex.x
1868                xminVert = vertex
1869            if vertex.x > xmax:
1870                xmax = vertex.x
1871                xmaxVert = vertex
1872               
1873            if vertex.y < ymin:
1874                ymin = vertex.y
1875                yminVert = vertex
1876            if vertex.y > ymax:
1877                ymax = vertex.y
1878                ymaxVert = vertex
1879        verts, count = self.removeDuplicatedVertices([xminVert,
1880                                                      xmaxVert,
1881                                                      yminVert,
1882                                                      ymaxVert])
1883         
1884        return verts
1885   
1886    def boxsize(self):
1887        """
1888        Returns a list denoting a box that contains the entire
1889        structure of vertices
1890        Structure: [xmin, ymin, xmax, ymax]
1891        """
1892     
1893        large = kinds.default_float_kind.MAX
1894        xmin= large
1895        xmax=-large
1896        ymin= large
1897        ymax=-large
1898        for vertex in self.userVertices:
1899            if vertex.x < xmin:
1900                xmin = vertex.x
1901            if vertex.x > xmax:
1902                xmax = vertex.x
1903               
1904            if vertex.y < ymin:
1905                ymin = vertex.y
1906            if vertex.y > ymax:
1907                ymax = vertex.y
1908        return [xmin, ymin, xmax, ymax]
1909 
1910    def maxMinVertAtt(self, iatt):
1911        """
1912        Returns a list denoting a box that contains the entire structure
1913        of vertices
1914        Structure: [xmin, ymin, xmax, ymax]
1915        """
1916       
1917        large = kinds.default_float_kind.MAX
1918        min= large
1919        max=-large
1920        for vertex in self.userVertices:
1921            if len(vertex.attributes) > iatt:
1922                if vertex.attributes[iatt] < min:
1923                    min = vertex.attributes[iatt]
1924                if vertex.attributes[iatt] > max:
1925                    max = vertex.attributes[iatt] 
1926        for vertex in self.meshVertices:
1927            if len(vertex.attributes) > iatt:
1928                if vertex.attributes[iatt] < min:
1929                    min = vertex.attributes[iatt]
1930                if vertex.attributes[iatt] > max:
1931                    max = vertex.attributes[iatt] 
1932        return [min, max]
1933   
1934    def scaleoffset(self, WIDTH, HEIGHT):
1935        """
1936        Returns a list denoting the scale and offset terms that need to be
1937        applied when converting  mesh co-ordinates onto grid co-ordinates
1938        Structure: [scale, xoffset, yoffset]
1939        """   
1940        OFFSET = 0.05*min([WIDTH, HEIGHT])
1941        [xmin, ymin, xmax, ymax] = self.boxsize()
1942        SCALE = min([0.9*WIDTH, 0.9*HEIGHT])/max([xmax-xmin, ymax-ymin])
1943       
1944        if SCALE*xmin < OFFSET:
1945            xoffset = abs(SCALE*xmin) + OFFSET
1946        if SCALE*xmax > WIDTH - OFFSET:
1947            xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1948        if SCALE*ymin < OFFSET:
1949            b = abs(SCALE*ymin)+OFFSET
1950        if SCALE*ymax > HEIGHT-OFFSET:
1951            b = -(SCALE*ymax - HEIGHT + OFFSET)
1952        yoffset = HEIGHT - b
1953        return [SCALE, xoffset, yoffset]
1954
1955         
1956    def exportASCIIobj(self,ofile):
1957        """
1958        export a file, ofile, with the format
1959         lines:  v <x> <y> <first attribute>
1960        f <vertex #>  <vertex #> <vertex #> (of the triangles)
1961        """
1962        fd = open(ofile,'w')
1963        self.writeASCIIobj(fd)   
1964        fd.close()
1965
1966
1967    def writeASCIIobj(self,fd):
1968        fd.write(" # Triangulation as an obj file\n")
1969        numVert = str(len(self.meshVertices))
1970       
1971        index1 = 1 
1972        for vert in self.meshVertices:
1973            vert.index1 = index1
1974            index1 += 1
1975           
1976            fd.write("v "
1977                     + str(vert.x) + " "
1978                     + str(vert.y) + " "
1979                     + str(vert.attributes[0]) + "\n")
1980           
1981        for tri in self.meshTriangles:
1982            fd.write("f "
1983                     + str(tri.vertices[0].index1) + " " 
1984                     + str(tri.vertices[1].index1) + " " 
1985                     + str(tri.vertices[2].index1) + "\n")
1986           
1987    def exportASCIIsegmentoutlinefile(self,ofile):
1988        """Write the boundary user mesh info, eg
1989         vertices that are connected to segments,
1990         segments
1991        """
1992       
1993        verts = {}
1994        for seg in self.getUserSegments():
1995            verts[seg.vertices[0]] = seg.vertices[0]
1996            verts[seg.vertices[1]] = seg.vertices[1]
1997        meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values())
1998        load_mesh.loadASCII.export_mesh_file(ofile,meshDict)
1999       
2000        # exportASCIImeshfile   - this function is used
2001    def export_mesh_file(self,ofile):
2002        """
2003        export a file, ofile, with the format
2004        """
2005       
2006        dict = self.Mesh2IODict()
2007        load_mesh.loadASCII.export_mesh_file(ofile,dict)
2008
2009    # FIXME(DSG-DSG):Break this into two functions.
2010    #One for the outline points.
2011    #One for the mesh points.
2012    # Note: this function is not in the gui
2013    def exportPointsFile(self,ofile):
2014        """
2015        export a points (.xya or .pts)  file, ofile.
2016       
2017        """
2018       
2019        mesh_dict = self.Mesh2IODict()
2020        point_dict = {}
2021        point_dict['attributelist'] = {} #this will need to be expanded..
2022                                         # if attributes are brought back in.
2023        point_dict['geo_reference'] = self.geo_reference
2024        if mesh_dict['vertices'] == []:
2025            point_dict['pointlist'] = mesh_dict['points']
2026        else:
2027            point_dict['pointlist'] = mesh_dict['vertices']
2028
2029        load_mesh.loadASCII.export_points_file(ofile,point_dict)
2030
2031
2032    def import_ungenerate_file(self,ofile, tag=None):
2033        """
2034        Imports an ungenerate file, from arcGIS into mesh.
2035
2036        ofile is the name of the ungenerated file.
2037        Tag is a string name to be taggged on each segment.
2038
2039        WARNING values are assumed to be absolute.
2040        geo-refs are not taken into account..
2041        """
2042   
2043        dict = importUngenerateFile(ofile)
2044        default_tag = Segment.get_default_tag()
2045        if tag is not None:
2046            Segment.set_default_tag(str(tag))
2047        self.addVertsSegs(dict)
2048        Segment.set_default_tag(default_tag)
2049
2050        # change the tag back to  it's default
2051   
2052       
2053########### IO CONVERTERS ##################
2054        """
2055        The dict fromat for IO with .tsh files is;
2056        (the triangulation)
2057        vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
2058        vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
2059        vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
2060        segments: [[v1,v2],[v3,v4],...] (lists of integers)
2061        segment_tags : [tag,tag,...] list of strings
2062        triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
2063        triangle tags: [s1,s2,...] A list of strings
2064        triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
2065       
2066        (the outline)   
2067        points: [[x1,y1],[x2,y2],...] (lists of doubles)
2068        point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
2069        outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
2070        outline_segment_tags : [tag1,tag2,...] list of strings
2071        holes : [[x1,y1],...](List of doubles, one inside each hole region)
2072        regions : [ [x1,y1],...] (List of 4 doubles)
2073        region_tags : [tag1,tag2,...] (list of strings)
2074        region_max_areas: [ma1,ma2,...] (A list of doubles)
2075        {Convension: A -ve max area means no max area}
2076       
2077        """
2078     
2079
2080                               
2081    def Mesh2IODict(self):
2082        """
2083        Convert the triangulation and outline info of a mesh to a dictionary
2084        structure
2085        """
2086        dict = self.Mesh2IOTriangulationDict()
2087        dict_mesh = self.Mesh2IOOutlineDict()
2088        for element in dict_mesh.keys():
2089            dict[element] = dict_mesh[element]
2090
2091        # add the geo reference
2092        dict['geo_reference'] = self.geo_reference
2093        return dict
2094   
2095    def Mesh2IOTriangulationDict(self):
2096        """
2097        Convert the Mesh to a dictionary of lists describing the
2098        triangulation variables;
2099       
2100        Used to produce .tsh file
2101        """
2102       
2103        meshDict = {}       
2104        vertices=[]
2105        vertex_attributes=[]
2106
2107        self.maxVertexIndex=0
2108        for vertex in self.meshVertices:
2109            vertex.index = self.maxVertexIndex
2110            vertices.append([vertex.x,vertex.y])
2111            vertex_attributes.append(vertex.attributes)           
2112            self.maxVertexIndex += 1
2113
2114        meshDict['vertices'] = vertices
2115        meshDict['vertex_attributes'] = vertex_attributes
2116        meshDict['vertex_attribute_titles'] = self.attributeTitles
2117        #segments
2118        segments=[]
2119        segment_tags=[]
2120        for seg in self.meshSegments:
2121            segments.append([seg.vertices[0].index,seg.vertices[1].index])
2122            segment_tags.append(seg.tag)
2123        meshDict['segments'] =segments
2124        meshDict['segment_tags'] =segment_tags
2125
2126        # Make sure that the indexation is correct
2127        index = 0
2128        for tri in self.meshTriangles:
2129            tri.index = index           
2130            index += 1
2131
2132        triangles = []
2133        triangle_tags = []
2134        triangle_neighbors = []
2135        for tri in self.meshTriangles: 
2136            triangles.append([tri.vertices[0].index,
2137                              tri.vertices[1].index,
2138                              tri.vertices[2].index]) 
2139            triangle_tags.append(tri.attribute)
2140            neighborlist = [-1,-1,-1]
2141            for neighbor,index in map(None,tri.neighbors,
2142                                      range(len(tri.neighbors))):
2143                if neighbor:   
2144                    neighborlist[index] = neighbor.index
2145            triangle_neighbors.append(neighborlist)
2146       
2147        meshDict['triangles'] = triangles
2148        meshDict['triangle_tags'] = triangle_tags
2149        meshDict['triangle_neighbors'] = triangle_neighbors
2150       
2151        #print "mesh.Mesh2IOTriangulationDict*)*)"
2152        #print meshDict
2153        #print "mesh.Mesh2IOTriangulationDict*)*)"
2154
2155        return meshDict
2156
2157                                                     
2158    def Mesh2IOOutlineDict(self, userVertices=None,
2159                        userSegments=None,
2160                        holes=None,
2161                        regions=None):
2162        """
2163        Convert the mesh outline to a dictionary of the lists needed for the
2164        triang module;
2165       
2166        Note, this adds an index attribute to the user Vertex objects.
2167
2168        Used to produce .tsh file and output to triangle
2169        """
2170        if userVertices is None:
2171            userVertices = self.getUserVertices()
2172        if userSegments is None:
2173            userSegments = self.getUserSegments()
2174        if holes is None:
2175            holes = self.getHoles()
2176        if regions is None:
2177            regions = self.getRegions()
2178           
2179        meshDict = {}
2180        #print "userVertices",userVertices
2181        #print "userSegments",userSegments
2182        pointlist=[]
2183        pointattributelist=[]
2184        index = 0
2185        for vertex in userVertices:
2186            vertex.index = index
2187            pointlist.append([vertex.x,vertex.y])
2188            pointattributelist.append(vertex.attributes)
2189           
2190            index += 1
2191        meshDict['points'] = pointlist
2192        meshDict['point_attributes'] = pointattributelist
2193
2194        segmentlist=[]
2195        segmenttaglist=[]
2196        for seg in userSegments:
2197            segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
2198            segmenttaglist.append(seg.tag)
2199        meshDict['outline_segments'] =segmentlist
2200        meshDict['outline_segment_tags'] =segmenttaglist
2201       
2202        holelist=[]
2203        for hole in holes:
2204            holelist.append([hole.x,hole.y]) 
2205        meshDict['holes'] = holelist
2206       
2207        regionlist=[]
2208        regiontaglist = []
2209        regionmaxarealist = []
2210        for region in regions:
2211            regionlist.append([region.x,region.y])
2212            regiontaglist.append(region.getTag())
2213           
2214            if (region.getMaxArea() != None): 
2215                regionmaxarealist.append(region.getMaxArea())
2216            else:
2217                regionmaxarealist.append( load_mesh.loadASCII.NOMAXAREA)
2218        meshDict['regions'] = regionlist
2219        meshDict['region_tags'] = regiontaglist
2220        meshDict['region_max_areas'] = regionmaxarealist
2221        #print "*(*("
2222        #print meshDict
2223        #print meshDict['regionlist']
2224        #print "*(*("
2225        return meshDict
2226
2227    def IOTriangulation2Mesh(self, genDict):
2228        """
2229        Set the mesh attributes given an tsh IO dictionary
2230        """
2231        #Clear the current generated mesh values
2232        self.meshTriangles=[]
2233        self.attributeTitles=[]
2234        self.meshSegments=[]
2235        self.meshVertices=[]
2236
2237        #print "mesh.setTriangulation@#@#@#"
2238        #print genDict
2239        #print "@#@#@#"
2240
2241        self.maxVertexIndex = 0
2242        for point in genDict['vertices']:
2243            v=Vertex(point[0], point[1])
2244            v.index =  self.maxVertexIndex
2245            self.maxVertexIndex +=1
2246            self.meshVertices.append(v)
2247
2248        self.attributeTitles = genDict['vertex_attribute_titles']
2249
2250        index = 0
2251        for seg,tag in map(None,genDict['segments'],genDict['segment_tags']):
2252            segObject = Segment( self.meshVertices[int(seg[0])],
2253                           self.meshVertices[int(seg[1])], tag = tag )
2254            segObject.index = index
2255            index +=1
2256            self.meshSegments.append(segObject)
2257
2258        index = 0
2259        for triangle in genDict['triangles']:
2260            tObject =Triangle( self.meshVertices[int(triangle[0])],
2261                        self.meshVertices[int(triangle[1])],
2262                        self.meshVertices[int(triangle[2])] )
2263            tObject.index = index
2264            index +=1
2265            self.meshTriangles.append(tObject)
2266
2267        index = 0
2268        for att in genDict['triangle_tags']:
2269            if att == []:
2270                self.meshTriangles[index].setAttribute("")
2271            else:
2272                self.meshTriangles[index].setAttribute(att)
2273            index += 1
2274           
2275        index = 0
2276        for att in genDict['vertex_attributes']:
2277            if att == None:
2278                self.meshVertices[index].setAttributes([])
2279            else:
2280                self.meshVertices[index].setAttributes(att)
2281            index += 1
2282   
2283        index = 0
2284        for triangle in genDict['triangle_neighbors']:
2285            # Build a list of triangle object neighbors
2286            ObjectNeighbor = []
2287            for neighbor in triangle:
2288                if ( neighbor != -1):
2289                    ObjectNeighbor.append(self.meshTriangles[int(neighbor)])
2290                else:
2291                    ObjectNeighbor.append(None)
2292            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
2293                                                   ObjectNeighbor[1],
2294                                                   ObjectNeighbor[2])
2295            index += 1
2296
2297
2298    def IOOutline2Mesh(self, genDict):
2299        """
2300        Set the outline (user Mesh attributes) given a IO tsh dictionary
2301       
2302        mesh is an instance of a mesh object
2303        """
2304        #Clear the current user mesh values
2305        self.clearUserSegments()
2306        self.userVertices=[]
2307        self.Holes=[]
2308        self.Regions=[]
2309
2310        #print "mesh.IOOutline2Mesh@#@#@#"
2311        #print "genDict",genDict
2312        #print "@#@#@#"
2313       
2314        #index = 0
2315        for point in genDict['points']:
2316            v=Vertex(point[0], point[1])
2317            #v.index = index
2318            #index +=1
2319            self.userVertices.append(v)
2320
2321        #index = 0
2322        for seg,tag in map(None,genDict['outline_segments'],
2323                           genDict['outline_segment_tags']):
2324
2325            segObject = Segment( self.userVertices[int(seg[0])],
2326                           self.userVertices[int(seg[1])], tag = tag )
2327            #segObject.index = index
2328            #index +=1
2329            self.userSegments.append(segObject)
2330
2331# Remove the loading of attribute info.
2332# Have attribute info added using least_squares in pyvolution
2333#         index = 0
2334#         for att in genDict['point_attributes']:
2335#             if att == None:
2336#                 self.userVertices[index].setAttributes([])
2337#             else:
2338#                 self.userVertices[index].setAttributes(att)
2339#            index += 1
2340       
2341        #index = 0
2342        for point in genDict['holes']:
2343            h=Hole(point[0], point[1])
2344            #h.index = index
2345            #index +=1
2346            self.holes.append(h)
2347
2348        #index = 0
2349        for reg,att,maxArea in map(None,
2350                                   genDict['regions'],
2351                                   genDict['region_tags'],
2352                                   genDict['region_max_areas']):
2353            if maxArea > 0:  # maybe I should ref NOMAXAREA? Prob' not though
2354                Object = Region( reg[0],
2355                                 reg[1],
2356                                 tag = att,
2357                                 maxArea = maxArea)
2358            else:
2359                Object = Region( reg[0],
2360                                 reg[1],
2361                                 tag = att)
2362               
2363            #Object.index = index
2364            #index +=1
2365            self.regions.append(Object)
2366 
2367############################################
2368
2369       
2370    def refineSet(self,setName):
2371        Triangles = self.sets[self.setID[setName]]
2372        Refine(self,Triangles)
2373
2374    def selectAllTriangles(self):
2375        A=[]
2376        A.extend(self.meshTriangles)
2377        if not('All' in self.setID.keys()):
2378            self.setID['All']=len(self.sets)
2379            self.sets.append(A)
2380        else:
2381            self.sets[self.setID['All']]=A
2382        return 'All'
2383        # and objectIDs
2384
2385
2386    def clearSelection(self):
2387        A = []
2388        if not('None' in self.setID.keys()):
2389            self.setID['None']=len(self.sets)
2390            self.sets.append(A)
2391        return 'None'
2392
2393    def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
2394    #FIXME Draws over previous triangles - may bloat canvas
2395        Triangles = self.sets[self.setID[setName]]
2396        for triangle in Triangles:
2397            triangle.draw(canvas,1,
2398                          scale = SCALE,
2399                          colour = colour)
2400           
2401    def undrawSet(self,canvas,setName,SCALE,colour='green'):
2402    #FIXME Draws over previous lines - may bloat canvas
2403        Triangles = self.sets[self.setID[setName]]
2404        for triangle in Triangles:
2405            triangle.draw(canvas,1,
2406                          scale = SCALE,
2407                          colour = colour)
2408
2409    def weed(self,Vertices,Segments):
2410        #Depreciated
2411        #weed out existing duplicates
2412        print 'len(self.getUserSegments())'
2413        print len(self.getUserSegments())
2414        print 'len(self.getUserVertices())'
2415        print len(self.getUserVertices())
2416
2417        point_keys = {}
2418        for vertex in Vertices:
2419            point = (vertex.x,vertex.y)
2420            point_keys[point]=vertex
2421        #inlined would looks very ugly
2422
2423        line_keys = {}
2424        for segment in Segments:
2425            vertex1 = segment.vertices[0]
2426            vertex2 = segment.vertices[1]
2427            point1 = (vertex1.x,vertex1.y)
2428            point2 = (vertex2.x,vertex2.y)
2429            segment.vertices[0]=point_keys[point1]
2430            segment.vertices[1]=point_keys[point2]
2431            vertex1 = segment.vertices[0]
2432            vertex2 = segment.vertices[1]
2433            point1 = (vertex1.x,vertex1.y)
2434            point2 = (vertex2.x,vertex2.y)
2435            line1 = (point1,point2)
2436            line2 = (point2,point1)
2437            if not (line_keys.has_key(line1) \
2438                 or line_keys.has_key(line2)):
2439                 line_keys[line1]=segment
2440        Vertices=point_keys.values()
2441        Segments=line_keys.values()
2442        return Vertices,Segments
2443
2444    def segs_to_dict(self,segments):
2445        dict={}
2446        for segment in segments:
2447            vertex1 = segment.vertices[0]
2448            vertex2 = segment.vertices[1]
2449            point1 = (vertex1.x,vertex1.y)
2450            point2 = (vertex2.x,vertex2.y)
2451            line = (point1,point2)
2452            dict[line]=segment
2453        return dict
2454
2455    def seg2line(self,s):
2456        return ((s.vertices[0].x,s.vertices[0].y,)\
2457                (s.vertices[1].x,s.vertices[1].y))
2458
2459    def line2seg(self,line,tag=None):
2460        point0 = self.point2ver(line[0])
2461        point1 = self.point2ver(line[1])
2462        return Segment(point0,point1,tag=tag)
2463
2464    def ver2point(self,vertex):
2465        return (vertex.x,vertex.y)
2466
2467    def point2ver(self,point):
2468        return Vertex(point[0],point[1])
2469
2470    def smooth_polySet(self,min_radius=0.05):
2471        #for all pairs of connecting segments:
2472        #    propose a new segment that replaces the 2
2473
2474        #    If the difference between the new segment
2475        #    and the old lines is small: replace the
2476        #    old lines.
2477
2478        seg2line = self.seg2line
2479        ver2point= self.ver2point
2480        line2seg = self.line2seg
2481        point2ver= self.point2ver
2482
2483        #create dictionaries of lines -> segments
2484        userSegments = self.segs_to_dict(self.userSegments)
2485        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2486
2487        #lump user and alpha segments
2488        for key in alphaSegments.keys():
2489            userSegments[key]=alphaSegments[key]
2490
2491        #point_keys = tuple -> vertex
2492        #userVertices = vertex -> [line,line] - lines from that node
2493        point_keys = {}
2494        userVertices={}
2495        for vertex in self.getUserVertices():
2496            point = ver2point(vertex)
2497            if not point_keys.has_key(point):
2498                point_keys[point]=vertex
2499                userVertices[vertex]=[]
2500        for key in userSegments.keys():
2501            line = key
2502            point_0 = key[0]
2503            point_1 = key[1]
2504            userVertices[point_keys[point_0]].append(line)
2505            userVertices[point_keys[point_1]].append(line)
2506
2507        for point in point_keys.keys():
2508            try:
2509            #removed keys can cause keyerrors
2510                vertex = point_keys[point]
2511                lines = userVertices[vertex]
2512   
2513                #if there are 2 lines on the node
2514                if len(lines)==2:
2515                    line_0 = lines[0]
2516                    line_1 = lines[1]
2517   
2518                    #if the tags are the the same on the 2 lines
2519                    if userSegments[line_0].tag == userSegments[line_1].tag:
2520                        tag = userSegments[line_0].tag
2521   
2522                        #point_a is one of the next nodes, point_b is the other
2523                        if point==line_0[0]:
2524                            point_a = line_0[1]
2525                        if point==line_0[1]:
2526                            point_a = line_0[0]
2527                        if point==line_1[0]:
2528                            point_b = line_1[1]
2529                        if point==line_1[1]:
2530                            point_b = line_1[0]
2531   
2532   
2533                        #line_2 is proposed
2534                        line_2 = (point_a,point_b)
2535
2536                        #calculate the area of the triangle between
2537                        #the two existing segments and the proposed
2538                        #new segment
2539                        ax = point_a[0]
2540                        ay = point_a[1]
2541                        bx = point_b[0]
2542                        by = point_b[1]
2543                        cx = point[0]
2544                        cy = point[1]
2545                        area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
2546
2547                        #calculate the perimeter
2548                        len_a =  ((cx-bx)**2+(cy-by)**2)**0.5 
2549                        len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
2550                        len_c =  ((bx-ax)**2+(by-ay)**2)**0.5 
2551                        perimeter = len_a+len_b+len_c
2552
2553                        #calculate the radius
2554                        r = area/(2*perimeter)
2555
2556                        #if the radius is small: then replace the existing
2557                        #segments with the new one
2558                        if r < min_radius:
2559                            if len_c < min_radius: append = False
2560                            else: append = True
2561                            #if the new seg is also time, don't add it
2562                            if append:
2563                                segment = self.line2seg(line_2,tag=tag)
2564
2565                            list_a=userVertices[point_keys[point_a]]
2566                            list_b=userVertices[point_keys[point_b]]
2567
2568                            if line_0 in list_a:
2569                                list_a.remove(line_0)
2570                            else:
2571                                list_a.remove(line_1)
2572
2573                            if line_0 in list_b:
2574                                list_b.remove(line_0)
2575                            else:
2576                                list_b.remove(line_1)
2577
2578                            if append:
2579                                list_a.append(line_2)
2580                                list_b.append(line_2)
2581                            else:
2582                                if len(list_a)==0:
2583                                    userVertices.pop(point_keys[point_a])
2584                                    point_keys.pop(point_a)
2585                                if len(list_b)==0:
2586                                    userVertices.pop(point_keys[point_b])
2587                                    point_keys.pop(point_b)
2588
2589                            userVertices.pop(point_keys[point])
2590                            point_keys.pop(point)
2591                            userSegments.pop(line_0)
2592                            userSegments.pop(line_1)
2593
2594                            if append:
2595                                userSegments[line_2]=segment
2596            except:
2597                pass
2598
2599        #self.userVerticies = userVertices.keys()
2600        #self.userSegments = []
2601        #for key in userSegments.keys():
2602        #    self.userSegments.append(userSegments[key])
2603        #self.alphaUserSegments = []
2604
2605        self.userVerticies = []
2606        self.userSegments = []
2607        self.alphaUserSegments = []
2608
2609        return userVertices,userSegments,alphaSegments
2610
2611    def triangles_to_polySet(self,setName):
2612        #self.smooth_polySet()
2613
2614        seg2line = self.seg2line
2615        ver2point= self.ver2point
2616        line2seg = self.line2seg
2617        point2ver= self.point2ver
2618
2619        from Numeric import array,allclose
2620        #turn the triangles into a set
2621        Triangles = self.sets[self.setID[setName]]
2622        Triangles_dict = {}
2623        for triangle in Triangles:
2624            Triangles_dict[triangle]=None
2625 
2626
2627        #create a dict of points to vertexes (tuple -> object)
2628        #also create a set of vertexes (object -> True)
2629        point_keys = {}
2630        userVertices={}
2631        for vertex in self.getUserVertices():
2632            point = ver2point(vertex)
2633            if not point_keys.has_key(point):
2634                point_keys[point]=vertex
2635                userVertices[vertex]=True
2636
2637        #create a dict of lines to segments (tuple -> object)
2638        userSegments = self.segs_to_dict(self.userSegments)
2639        #append the userlines in an affine linespace
2640        affine_lines = Affine_Linespace()
2641        for line in userSegments.keys():
2642            affine_lines.append(line)
2643        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2644        for line in alphaSegments.keys():
2645            affine_lines.append(line)
2646       
2647        for triangle in Triangles:
2648            for i in (0,1,2):
2649                #for every triangles neighbour:
2650                if not Triangles_dict.has_key(triangle.neighbors[i]):
2651                #if the neighbour is not in the set:
2652                    a = triangle.vertices[i-1]
2653                    b = triangle.vertices[i-2]
2654                    #Get possible matches:
2655                    point_a = ver2point(a)
2656                    point_b = ver2point(b)
2657                    midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
2658                    line = (point_a,point_b)
2659                    tag = None
2660
2661
2662                    #this bit checks for matching lines
2663                    possible_lines = affine_lines[line] 
2664                    possible_lines = unique(possible_lines)
2665                    found = 0                           
2666                    for user_line in possible_lines:
2667                        if self.point_on_line(midpoint,user_line):
2668                            found+=1
2669                            assert found<2
2670                            if userSegments.has_key(user_line):
2671                                parent_segment = userSegments.pop(user_line)
2672                            if alphaSegments.has_key(user_line):
2673                                parent_segment = alphaSegments.pop(user_line)
2674                            tag = parent_segment.tag
2675                            offspring = [line]
2676                            offspring.extend(self.subtract_line(user_line,line))
2677                            affine_lines.remove(user_line)
2678                            for newline in offspring:
2679                                line_vertices = []
2680                                for point in newline:
2681                                    if point_keys.has_key(point):
2682                                        vert = point_keys[point]
2683                                    else:
2684                                        vert = Vertex(point[0],point[1])
2685                                        userVertices[vert]=True
2686                                        point_keys[point]=vert
2687                                    line_vertices.append(vert)
2688                                segment = Segment(line_vertices[0],line_vertices[1],tag)
2689                                userSegments[newline]=segment
2690                                affine_lines.append(newline)
2691                            #break
2692                    assert found<2
2693
2694
2695
2696                    #if no matching lines
2697                    if not found:
2698                        line_vertices = []
2699                        for point in line:
2700                            if point_keys.has_key(point):
2701                                vert = point_keys[point]
2702                            else:
2703                                vert = Vertex(point[0],point[1])
2704                                userVertices[vert]=True
2705                                point_keys[point]=vert
2706                            line_vertices.append(vert)
2707                        segment = Segment(line_vertices[0],line_vertices[1],tag)
2708                        userSegments[line]=segment
2709                        affine_lines.append(line)
2710       
2711        self.userVerticies = []
2712        self.userSegments = []
2713        self.alphaUserSegments = []
2714
2715        return userVertices,userSegments,alphaSegments
2716
2717    def subtract_line(self,parent,child):
2718        #Subtracts child from parent
2719        #Requires that the child is a
2720        #subline of parent to work.
2721
2722        from Numeric import allclose,dot,array
2723        A= parent[0]
2724        B= parent[1]
2725        a = child[0]
2726        b = child[1]
2727
2728        A_array = array(parent[0])
2729        B_array = array(parent[1])
2730        a_array   = array(child[0])
2731        b_array   = array(child[1])
2732
2733        assert not A == B
2734        assert not a == b
2735
2736        answer = []
2737
2738        #if the new line does not share a
2739        #vertex with the old one
2740        if not (allclose(A_array,a_array)\
2741             or allclose(B_array,b_array)\
2742             or allclose(A_array,b_array)\
2743             or allclose(a_array,B_array)):
2744            if dot(A_array-a_array,A_array-a_array) \
2745            < dot(A_array-b_array,A_array-b_array):
2746                sibling1 = (A,a)
2747                sibling2 = (B,b)
2748                return [sibling1,sibling2]
2749            else:
2750                sibling1 = (A,b)
2751                sibling2 = (B,a)
2752                return [sibling1,sibling2]
2753
2754        elif allclose(A_array,a_array):
2755            if allclose(B_array,b_array):
2756                return []
2757            else:
2758                sibling = (b,B)
2759                return [sibling]
2760        elif allclose(B_array,b_array):
2761            sibling = (a,A)
2762            return [sibling]
2763
2764        elif allclose(A_array,b_array):
2765            if allclose(B,a):
2766                return []
2767            else:
2768                sibling = (a,B)
2769                return [sibling]
2770        elif allclose(a_array,B_array):
2771            sibling = (b,A)
2772            return [sibling]
2773
2774    def point_on_line(self,point,line):
2775        #returns true within a tolerance of 3 degrees
2776        x=point[0]
2777        y=point[1]
2778        x0=line[0][0]
2779        x1=line[1][0]
2780        y0=line[0][1]
2781        y1=line[1][1]
2782        from Numeric import array, dot, allclose
2783        from math import sqrt
2784        tol = 3. #DEGREES
2785        tol = tol*3.1415/180
2786
2787        a = array([x - x0, y - y0]) 
2788        a_normal = array([a[1], -a[0]])     
2789        len_a_normal = sqrt(sum(a_normal**2)) 
2790
2791        b = array([x1 - x0, y1 - y0])         
2792        len_b = sqrt(sum(b**2)) 
2793   
2794        if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
2795            #Point is somewhere on the infinite extension of the line
2796
2797            len_a = sqrt(sum(a**2))                                     
2798            if dot(a, b) >= 0 and len_a <= len_b:
2799               return True
2800            else:   
2801               return False
2802        else:
2803          return False
2804
2805    def line_length(self,line):
2806        x0=line[0][0]
2807        x1=line[1][0]
2808        y0=line[0][1]
2809        y1=line[1][1]
2810        return ((x1-x0)**2-(y1-y0)**2)**0.5     
2811
2812    def threshold(self,setName,min=None,max=None,attribute_name = 'elevation'):
2813        """
2814        threshold using  d
2815        """
2816        triangles = self.sets[self.setID[setName]]
2817        A = []
2818
2819        if attribute_name in self.attributeTitles:
2820            i = self.attributeTitles.index(attribute_name)
2821        else: i = -1#no attribute
2822        if not max == None:
2823            for t in triangles:
2824                if (min<self.av_att(t,i)<max):
2825                    A.append(t)
2826        else:
2827            for t in triangles:
2828                if (min<self.av_att(t,i)):
2829                    A.append(t)
2830        self.sets[self.setID[setName]] = A
2831
2832    def general_threshold(self,setName,min=None,max=None\
2833              ,attribute_name = 'elevation',function=None):
2834        """
2835        Thresholds the triangles
2836        """
2837        from visual.graph import arange,ghistogram,color as colour
2838        triangles = self.sets[self.setID[setName]]
2839        A = []
2840        data=[]
2841        #data is for the graph
2842
2843        if attribute_name in self.attributeTitles:
2844            i = self.attributeTitles.index(attribute_name)
2845        else: i = -1
2846        if not max == None:
2847            for t in triangles:
2848                value=function(t,i)
2849                if (min<value<max):
2850                    A.append(t)
2851                data.append(value)
2852        else:
2853            for t in triangles:
2854                value=function(t,i)
2855                if (min<value):
2856                    A.append(t)
2857                data.append(value)
2858        self.sets[self.setID[setName]] = A
2859
2860        if self.visualise_graph:
2861            if len(data)>0:
2862                max=data[0]
2863                min=data[0]
2864                for value in data:
2865                    if value > max:
2866                        max = value
2867                    if value < min:
2868                        min = value
2869
2870                inc = (max-min)/100
2871
2872                histogram = ghistogram(bins=arange(min,max,inc),\
2873                             color = colour.red)
2874                histogram.plot(data=data)
2875       
2876    def av_att(self,triangle,i):
2877        if i==-1: return 1
2878        else:
2879            #evaluates the average attribute of the vertices of a triangle.
2880            V = triangle.getVertices()
2881            a0 = (V[0].attributes[i])
2882            a1 = (V[1].attributes[i])
2883            a2 = (V[2].attributes[i])
2884            return (a0+a1+a2)/3
2885
2886    def Courant_ratio(self,triangle,index):
2887        """
2888        Uses the courant threshold
2889        """
2890        e = self.av_att(triangle,index)
2891        A = triangle.calcArea()
2892        P = triangle.calcP()
2893        r = A/(2*P)
2894        e = max(0.1,abs(e))
2895        return r/e**0.5
2896
2897    def Gradient(self,triangle,index):
2898        V = triangle.vertices
2899        x0, y0, x1, y1, x2, y2, q0, q1, q2 = V[0].x,V[0].y,V[1].x,V[1].y,V[2].x,V[2].y,V[0].attributes[index],V[1].attributes[index],V[2].attributes[index]
2900        grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
2901        if ((grad_x**2)+(grad_y**2))**(0.5)<0:
2902            print ((grad_x**2)+(grad_y**2))**(0.5)
2903        return ((grad_x**2)+(grad_y**2))**(0.5)
2904   
2905
2906    def append_triangle(self,triangle):
2907        self.meshTriangles.append(triangle)
2908
2909    def replace_triangle(self,triangle,replacement):
2910        i = self.meshTriangles.index(triangle)
2911        self.meshTriangles[i]=replacement
2912        assert replacement in self.meshTriangles
2913
2914def importUngenerateFile(ofile):
2915    """
2916    import a file, ofile, with the format
2917    [poly]
2918    poly format:
2919    First line:  <# of vertices> <x centroid> <y centroid>
2920    Following lines: <x> <y>
2921    last line:  "END"
2922
2923    Note: These are clockwise.
2924    """
2925    fd = open(ofile,'r')
2926    Dict = readUngenerateFile(fd)
2927    fd.close()
2928    return Dict
2929
2930def readUngenerateFile(fd):
2931    """
2932    import a file, ofile, with the format
2933    [poly]
2934    poly format:
2935    First line:  <# of polynomial> <x centroid> <y centroid>
2936    Following lines: <x> <y>
2937    last line:  "END"
2938    """
2939    END_DELIMITER = 'END'
2940   
2941    points = []
2942    segments = []
2943   
2944    isEnd = False
2945    line = fd.readline() #not used <# of polynomial> <x> <y>
2946    while not isEnd:
2947        line = fd.readline()
2948        fragments = line.split()
2949        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2950        points.append(vert)
2951        PreviousVertIndex = len(points)-1
2952        firstVertIndex = PreviousVertIndex
2953       
2954        line = fd.readline() #Read the next line
2955        while not line.startswith(END_DELIMITER): 
2956            #print "line >" + line + "<"
2957            fragments = line.split()
2958            vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2959            points.append(vert)
2960            thisVertIndex = len(points)-1
2961            segment = [PreviousVertIndex,thisVertIndex]
2962            segments.append(segment)
2963            PreviousVertIndex = thisVertIndex
2964            line = fd.readline() #Read the next line
2965            i =+ 1
2966        # If the last and first segments are the same,
2967        # Remove the last segment and the last vertex
2968        # then add a segment from the second last vert to the 1st vert
2969        thisVertIndex = len(points)-1
2970        firstVert = points[firstVertIndex]
2971        thisVert = points[thisVertIndex]
2972        #print "firstVert",firstVert
2973        #print "thisVert",thisVert
2974        if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
2975            points.pop()
2976            segments.pop()
2977            thisVertIndex = len(points)-1
2978            segments.append([thisVertIndex, firstVertIndex])
2979       
2980        line = fd.readline() # read <# of polynomial> <x> <y> OR END
2981        #print "line >>" + line + "<<"
2982        if line.startswith(END_DELIMITER):
2983            isEnd = True
2984   
2985    #print "points", points       
2986    #print "segments", segments
2987    ungenerated_dict = {}
2988    ungenerated_dict['points'] = points
2989    ungenerated_dict['segments'] = segments
2990    return ungenerated_dict
2991
2992def importMeshFromFile(ofile):
2993    """returns a mesh object, made from a .xya/.pts or .tsh/.msh file
2994    Often raises IOError,RuntimeError
2995    """
2996    newmesh = None
2997    if (ofile[-4:]== ".xya" or ofile[-4:]== ".pts"):
2998        dict = load_mesh.loadASCII.import_points_file(ofile)
2999        dict['points'] = dict['pointlist']
3000        dict['outline_segments'] = []
3001        dict['outline_segment_tags'] = []
3002        dict['regions'] = []
3003        dict['region_tags'] = []
3004        dict['region_max_areas'] = []
3005        dict['holes'] = [] 
3006        newmesh= Mesh(geo_reference = dict['geo_reference'])
3007        newmesh.IOOutline2Mesh(dict)
3008        counter = newmesh.removeDuplicatedUserVertices()
3009        if (counter >0):
3010            print "%i duplicate vertices removed from dataset" % (counter)
3011    elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
3012        dict = load_mesh.loadASCII.import_mesh_file(ofile)
3013        #print "********"
3014        #print "zq mesh.dict",dict
3015        #print "********"
3016        newmesh= Mesh()
3017        newmesh.IOOutline2Mesh(dict)
3018        newmesh.IOTriangulation2Mesh(dict)
3019    else:
3020        raise RuntimeError
3021   
3022    if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
3023        newmesh.geo_reference = dict['geo_reference']
3024    return newmesh
3025
3026def loadPickle(currentName):
3027    fd = open(currentName)
3028    mesh = pickle.load(fd)
3029    fd.close()
3030    return mesh
3031   
3032def square_outline(side_length = 1,up = "top", left = "left", right = "right",
3033                   down = "bottom", regions = False):
3034   
3035        a = Vertex (0,0)
3036        b = Vertex (0,side_length)
3037        c = Vertex (side_length,0)
3038        d = Vertex (side_length,side_length)
3039     
3040        s2 = Segment(b,d, tag = up)
3041        s3 = Segment(b,a, tag = left)
3042        s4 = Segment(d,c, tag = right)
3043        s5 = Segment(a,c, tag = down)
3044
3045        if regions:
3046            e =  Vertex (side_length/2,side_length/2)
3047            s6 = Segment(a,e, tag = down + left)
3048            s7 = Segment(b,e, tag = up + left)
3049            s8 = Segment(c,e, tag = down + right)
3050            s9 = Segment(d,e, tag = up + right)
3051            r1 = Region(side_length/2,3.*side_length/4, tag = up)
3052            r2 = Region(1.*side_length/4,side_length/2, tag = left)
3053            r3 = Region(3.*side_length/4,side_length/2, tag = right)
3054            r4 = Region(side_length/2,1.*side_length/4, tag = down)
3055            mesh = Mesh(userVertices=[a,b,c,d,e],
3056                        userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
3057                        regions = [r1,r2,r3,r4])
3058        else:
3059            mesh = Mesh(userVertices=[a,b,c,d],
3060                 userSegments=[s2,s3,s4,s5])
3061     
3062        return mesh
3063
3064   
3065
3066def region_strings2ints(region_list):
3067    """Given a list of (x_int,y_int,tag_string) lists it returns a list of
3068    (x_int,y_int,tag_int) and a list to convert the tag_int's back to
3069    the tag_strings 
3070    """
3071    # Make sure "" has an index of 0
3072    region_list.reverse()
3073    region_list.append((1.0,2.0,""))
3074    region_list.reverse()
3075    convertint2string = []
3076    for i in xrange(len(region_list)):
3077        convertint2string.append(region_list[i][2])
3078        if len(region_list[i]) == 4: # there's an area value
3079            region_list[i] = (region_list[i][0],
3080                              region_list[i][1],i,region_list[i][3])
3081        elif len(region_list[i]) == 3: # no area value
3082            region_list[i] = (region_list[i][0],region_list[i][1],i)
3083        else:
3084            print "The region list has a bad size"
3085            # raise an error ..
3086            raise Error
3087
3088    #remove "" from the region_list
3089    region_list.pop(0)
3090   
3091    return [region_list, convertint2string]
3092
3093def region_ints2strings(region_list,convertint2string):
3094    """Reverses the transformation of region_strings2ints
3095    """
3096    if region_list[0] != []:
3097        for i in xrange(len(region_list)):
3098            region_list[i] = [convertint2string[int(region_list[i][0])]]
3099    return region_list
3100
3101def segment_ints2strings(intlist, convertint2string):
3102    """Reverses the transformation of segment_strings2ints """
3103    stringlist = []
3104    for x in intlist:
3105        stringlist.append(convertint2string[x])
3106    return stringlist
3107
3108def segment_strings2ints(stringlist, preset):
3109    """Given a list of strings return a list of 0 to n ints which represent
3110    the strings and a converting list of the strings, indexed by 0 to n ints.
3111    Also, input an initial converting list of the strings
3112    Note, the converting list starts off with
3113    ["internal boundary", "external boundary", "internal boundary"]
3114    example input and output
3115    input ["a","b","a","c"],["c"]
3116    output [[2, 1, 2, 0], ['c', 'b', 'a']]
3117
3118    the first element in the converting list will be
3119    overwritten with "".
3120    ?This will become the third item in the converting list?
3121   
3122    # note, order the initial converting list is important,
3123     since the index = the int tag
3124    """
3125    nodups = unique(stringlist)
3126    # note, order is important, the index = the int tag
3127    #preset = ["internal boundary", "external boundary"]
3128    #Remove the preset tags from the list with no duplicates
3129    nodups = [x for x in nodups if x not in preset]
3130
3131    try:
3132        nodups.remove("") # this has to go to zero
3133    except ValueError:
3134        pass
3135   
3136    # Add the preset list at the beginning of no duplicates
3137    preset.reverse()
3138    nodups.extend(preset)
3139    nodups.reverse()
3140
3141       
3142    convertstring2int = {}
3143    convertint2string = []
3144    index = 0
3145    for x in nodups:
3146        convertstring2int[x] = index
3147        convertint2string.append(x)
3148        index += 1
3149    convertstring2int[""] = 0
3150
3151    intlist = []
3152    for x in stringlist:
3153        intlist.append(convertstring2int[x])
3154    return [intlist, convertint2string]
3155
3156
3157
3158
3159def unique(s):
3160    """Return a list of the elements in s, but without duplicates.
3161
3162    For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
3163    unique("abcabc") some permutation of ["a", "b", "c"], and
3164    unique(([1, 2], [2, 3], [1, 2])) some permutation of
3165    [[2, 3], [1, 2]].
3166
3167    For best speed, all sequence elements should be hashable.  Then
3168    unique() will usually work in linear time.
3169
3170    If not possible, the sequence elements should enjoy a total
3171    ordering, and if list(s).sort() doesn't raise TypeError it's
3172    assumed that they do enjoy a total ordering.  Then unique() will
3173    usually work in O(N*log2(N)) time.
3174
3175    If that's not possible either, the sequence elements must support
3176    equality-testing.  Then unique() will usually work in quadratic
3177    time.
3178    """
3179
3180    n = len(s)
3181    if n == 0:
3182        return []
3183
3184    # Try using a dict first, as that's the fastest and will usually
3185    # work.  If it doesn't work, it will usually fail quickly, so it
3186    # usually doesn't cost much to *try* it.  It requires that all the
3187    # sequence elements be hashable, and support equality comparison.
3188    u = {}
3189    try:
3190        for x in s:
3191            u[x] = 1
3192    except TypeError:
3193        del u  # move on to the next method
3194    else:
3195        return u.keys()
3196
3197    # We can't hash all the elements.  Second fastest is to sort,
3198    # which brings the equal elements together; then duplicates are
3199    # easy to weed out in a single pass.
3200    # NOTE:  Python's list.sort() was designed to be efficient in the
3201    # presence of many duplicate elements.  This isn't true of all
3202    # sort functions in all languages or libraries, so this approach
3203    # is more effective in Python than it may be elsewhere.
3204    try:
3205        t = list(s)
3206        t.sort()
3207    except TypeError:
3208        del t  # move on to the next method
3209    else:
3210        assert n > 0
3211        last = t[0]
3212        lasti = i = 1
3213        while i < n:
3214            if t[i] != last:
3215                t[lasti] = last = t[i]
3216                lasti += 1
3217            i += 1
3218        return t[:lasti]
3219
3220    # Brute force is all that's left.
3221    u = []
3222    for x in s:
3223        if x not in u:
3224            u.append(x)
3225    return u
3226
3227"""Refines triangles
3228
3229   Implements the #triangular bisection?# algorithm.
3230 
3231   
3232"""
3233
3234def Refine(mesh, triangles):
3235    """
3236    Given a general_mesh, and a triangle number, split
3237    that triangle in the mesh in half. Then to prevent
3238    vertices and edges from meeting, keep refining
3239    neighbouring triangles until the mesh is clean.
3240    """
3241    state = BisectionState(mesh)
3242    for triangle in triangles:
3243        if not state.refined_triangles.has_key(triangle):
3244            triangle.rotate_longest_side()
3245            state.start(triangle)
3246            Refine_mesh(mesh, state)
3247
3248def Refine_mesh(mesh, state):
3249    """
3250    """
3251    state.getState(mesh)
3252    refine_triangle(mesh,state)
3253    state.evolve()
3254    if not state.end:
3255        Refine_mesh(mesh,state)
3256
3257def refine_triangle(mesh,state):
3258    split(mesh,state.current_triangle,state.new_point)
3259    if state.case == 'one':
3260        state.r[3]=state.current_triangle#triangle 2
3261
3262        new_triangle_id = len(mesh.meshTriangles)-1
3263        new_triangle = mesh.meshTriangles[new_triangle_id]
3264
3265        split(mesh,new_triangle,state.old_point)
3266        state.r[2]=new_triangle#triangle 1.2
3267        state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
3268        r = state.r
3269        state.repairCaseOne()
3270
3271    if state.case == 'two':
3272        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3273
3274        new_triangle = state.current_triangle
3275
3276        split(mesh,new_triangle,state.old_point)
3277
3278        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
3279        state.r[4]=new_triangle#triangle 2.2
3280        r = state.r
3281
3282        state.repairCaseTwo()
3283
3284    if state.case == 'vertex':
3285        state.r[2]=state.current_triangle#triangle 2
3286        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3287        r = state.r
3288        state.repairCaseVertex()
3289       
3290    if state.case == 'start':
3291        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3292        state.r[3]=state.current_triangle#triangle 2
3293
3294    if state.next_case == 'boundary':
3295        state.repairCaseBoundary()
3296
3297
3298def split(mesh, triangle, new_point):
3299        """
3300        Given a mesh, triangle_id and a new point,
3301        split the corrosponding triangle into two
3302        new triangles and update the mesh.
3303        """
3304
3305        new_triangle1 = Triangle(new_point,triangle.vertices[0],triangle.vertices[1],attribute = triangle.attribute, neighbors = None)
3306        new_triangle2 = Triangle(new_point,triangle.vertices[2],triangle.vertices[0],attribute = triangle.attribute, neighbors = None)
3307
3308        new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
3309        new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
3310
3311        mesh.meshTriangles.append(new_triangle1)
3312
3313        triangle.vertices = new_triangle2.vertices
3314        triangle.neighbors = new_triangle2.neighbors
3315
3316
3317class State:
3318
3319    def __init__(self):
3320        pass
3321
3322class BisectionState(State):
3323
3324
3325    def __init__(self,mesh):
3326        self.len = len(mesh.meshTriangles)
3327        self.refined_triangles = {}
3328        self.mesh = mesh
3329        self.current_triangle = None
3330        self.case = 'start'
3331        self.end = False
3332        self.r = [None,None,None,None,None]
3333
3334    def start(self, triangle):
3335        self.current_triangle = triangle
3336        self.case = 'start'
3337        self.end = False
3338        self.r = [None,None,None,None,None]
3339
3340    def getState(self,mesh):
3341        if not self.case == 'vertex':
3342            self.new_point=self.getNewVertex(mesh, self.current_triangle)
3343            #self.neighbour=self.getNeighbour(mesh, self.current_triangle)
3344            self.neighbour = self.current_triangle.neighbors[0]
3345            if not self.neighbour is None:
3346                self.neighbour.rotate_longest_side()
3347            self.next_case = self.get_next_case(mesh,self.neighbour,self.current_triangle)
3348        if self.case == 'vertex':
3349            self.new_point=self.old_point
3350
3351
3352    def evolve(self):
3353        if self.case == 'vertex':
3354            self.end = True
3355
3356        self.last_case = self.case
3357        self.case = self.next_case
3358
3359        self.old_point = self.new_point
3360        self.current_triangle = self.neighbour
3361
3362        if self.case == 'boundary':
3363            self.end = True
3364        self.refined_triangles[self.r[2]]=1
3365        self.refined_triangles[self.r[3]]=1
3366        if not self.r[4] is None:
3367            self.refined_triangles[self.r[4]]=1
3368        self.r[0]=self.r[2]
3369        self.r[1]=self.r[3]
3370
3371
3372    def getNewVertex(self,mesh,triangle):
3373        coordinate1 = triangle.vertices[1]
3374        coordinate2 = triangle.vertices[2]
3375        a = ([coordinate1.x*1.,coordinate1.y*1.])
3376        b = ([coordinate2.x*1.,coordinate2.y*1.])
3377        attributes = []
3378        for i in range(len(coordinate1.attributes)):
3379            att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
3380            attributes.append(att)
3381        new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
3382        newVertex = Vertex(new_coordinate[0],new_coordinate[1], attributes = attributes)
3383        mesh.maxVertexIndex+=1
3384        newVertex.index = mesh.maxVertexIndex
3385        mesh.meshVertices.append(newVertex)
3386        return newVertex
3387
3388    def get_next_case(self, mesh,neighbour,triangle):
3389        """
3390        Given the locations of two neighbouring triangles,
3391        examine the interior indices of their vertices (i.e.
3392        0,1 or 2) to determine what how the neighbour needs
3393        to be refined.
3394        """
3395        if (neighbour is None):
3396                next_case = 'boundary'
3397        else:
3398                if triangle.vertices[1].x==neighbour.vertices[2].x:
3399                    if triangle.vertices[1].y==neighbour.vertices[2].y:
3400                        next_case = 'vertex'
3401                if triangle.vertices[1].x==neighbour.vertices[0].x:
3402                    if triangle.vertices[1].y==neighbour.vertices[0].y:
3403                        next_case = 'two'
3404                if triangle.vertices[1].x==neighbour.vertices[1].x:
3405                    if triangle.vertices[1].y==neighbour.vertices[1].y:
3406                        next_case = 'one'
3407        return next_case
3408
3409
3410
3411    def repairCaseVertex(self):
3412
3413        r = self.r
3414
3415
3416        self.link(r[0],r[2])
3417        self.repair(r[0])
3418
3419        self.link(r[1],r[3])
3420        self.repair(r[1])
3421
3422        self.repair(r[2])
3423
3424        self.repair(r[3])
3425
3426
3427    def repairCaseOne(self):
3428        r = self.rkey
3429
3430
3431        self.link(r[0],r[2])
3432        self.repair(r[0])
3433
3434        self.link(r[1],r[4])
3435        self.repair(r[1])
3436
3437        self.repair(r[4])
3438
3439    def repairCaseTwo(self):
3440        r = self.r
3441
3442        self.link(r[0],r[4])
3443        self.repair(r[0])
3444
3445        self.link(r[1],r[3])
3446        self.repair(r[1])
3447
3448        self.repair(r[4])
3449
3450    def repairCaseBoundary(self):
3451        r = self.r
3452        self.repair(r[2])
3453        self.repair(r[3])
3454
3455
3456
3457    def repair(self,triangle):
3458        """
3459        Given a triangle that knows its neighbours, this will
3460        force the neighbours to comply.
3461
3462        However, it needs to compare the vertices of triangles
3463        for this implementation
3464
3465        But it doesn't work for invalid neighbour structures
3466        """
3467        n=triangle.neighbors
3468        for i in (0,1,2):
3469            if not n[i] is None:
3470                for j in (0,1,2):#to find which side of the list is broken
3471                    if not (n[i].vertices[j] in triangle.vertices):
3472                    #ie if j is the side of n that needs fixing
3473                        k = j
3474                n[i].neighbors[k]=triangle
3475
3476
3477
3478    def link(self,triangle1,triangle2):
3479        """
3480        make triangle1 neighbors point to t
3481                #count = 0riangle2
3482        """
3483        count = 0
3484        for i in (0,1,2):#to find which side of the list is broken
3485            if not (triangle1.vertices[i] in triangle2.vertices):
3486                j = i
3487                count+=1
3488        assert count == 1
3489        triangle1.neighbors[j]=triangle2
3490
3491class Discretised_Tuple_Set:
3492    """
3493    if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
3494    a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
3495    a[(10000)]=[] #NOT KEYERROR
3496
3497    a.append[(0.01)]
3498    => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
3499
3500    #NOT IMPLIMENTED
3501    a.remove[(0.01)]
3502    => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
3503
3504    a.remove[(0.17)]
3505    => {(0.0):[(0.02),(0.01)],0.2:[]}
3506    #NOT IMPLIMENTED
3507    at a.precision = 2:
3508        a.round_up_rel[0.0]=
3509        a.round_flat[0.0]=
3510        a.round_down_rel[0.0]=
3511
3512        a.up((0.1,2.04))=
3513
3514    If t_rel = 0, nothing gets rounded into
3515    two bins. If t_rel = 0.5, everything does.
3516
3517    Ideally, precision can be set high, so that multiple
3518    entries are rarely in the same bin. And t_rel should
3519    be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
3520    so that it is rare to put items in mutiple bins.
3521
3522    Ex bins per entry = product(a,b,c...,n)
3523    a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
3524    b = 1 or 2 ...
3525
3526    BUT!!! to avoid missing the right bin:
3527    (-10)**(precision+1)*t_rel must be greater than the
3528    greatest possible variation that an identical element
3529    can display.
3530
3531
3532    Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
3533    but not .6 - this looks wrong, but note that *everything* will round,
3534    so .6 wont be missed as everything close to it will check in .7 and .5.
3535    """
3536    def __init__(self,p_rel = 6,t_rel = 0.01):
3537        self.__p_rel__ = p_rel
3538        self.__t_rel__ = t_rel
3539
3540        self.__p_abs__ = p_rel+1
3541        self.__t_abs__ = t_rel
3542
3543        assert t_rel <= 0.5
3544        self.__items__ = {}
3545        from math import frexp
3546        self.frexp = frexp
3547        roundings = [self.round_up_rel,\
3548        self.round_down_rel,self.round_flat_rel,\
3549        self.round_down_abs,self.round_up_abs,\
3550        self.round_flat_abs]#
3551
3552        self.roundings = roundings
3553
3554    def __repr__(self):
3555        return '%s'%self.__items__
3556
3557    def rounded_keys(self,key):
3558        key = tuple(key)
3559        keys = [key]
3560        keys = self.__rounded_keys__(key)
3561        return (keys)
3562
3563    def __rounded_keys__(self,key):
3564        keys = []
3565        rounded_key=list(key)
3566        rounded_values=list(key)
3567
3568        roundings = list(self.roundings)
3569
3570        #initialise rounded_values
3571        round = roundings.pop(0)
3572        for i in range(len(rounded_values)):
3573            rounded_key[i]=round(key[i])
3574            rounded_values[i]={}
3575            rounded_values[i][rounded_key[i]]=None
3576        keys.append(tuple(rounded_key))
3577       
3578        for round in roundings:
3579            for i in range(len(rounded_key)):
3580                rounded_value=round(key[i])
3581                if not rounded_values[i].has_key(rounded_value):
3582                    #ie unless round_up_rel = round_down_rel
3583                    #so the keys stay unique
3584                    for j in range(len(keys)):
3585                        rounded_key = list(keys[j])
3586                        rounded_key[i]=rounded_value
3587                        keys.append(tuple(rounded_key))
3588        return keys
3589
3590    def append(self,item):
3591        keys = self.rounded_keys(item)
3592        for key in keys:
3593            if self.__items__.has_key(key):
3594                self.__items__[key].append(item)
3595            else:
3596                self.__items__[key]=[item]
3597
3598    def __getitem__(self,key):
3599        answer = []
3600        keys = self.rounded_keys(key)
3601        for key in keys:
3602            if self.__items__.has_key(key):
3603                answer.extend(self.__items__[key])
3604        #if len(answer)==0:
3605        #    raise KeyError#FIXME or return KeyError
3606        #                  #FIXME or just return []?
3607        else:
3608            return answer #FIXME or unique(answer)?
3609
3610    def __delete__(self,item):
3611        keys = self.rounded_keys(item)
3612        answer = False
3613        #if any of the possible keys contains
3614        #a list, return true
3615        for key in keys:       
3616            if self.__items__.has_key(key):
3617                if item in self.__items__[key]:
3618                    self.__items__[key].remove(item)
3619
3620    def remove(self,item):
3621        self.__delete__(item)
3622
3623    def __contains__(self,item):
3624
3625        keys = self.rounded_keys(item)
3626        answer = False
3627        #if any of the possible keys contains
3628        #a list, return true
3629        for key in keys:       
3630            if self.__items__.has_key(key):
3631                if item in self.__items__[key]:
3632                    answer = True
3633        return answer
3634
3635
3636    def has_item(self,item):
3637        return self.__contains__(item)
3638
3639    def round_up_rel2(self,value):
3640         t_rel=self.__t_rel__
3641         #Rounding up the value
3642         m,e = self.frexp(value)
3643         m = m/2
3644         e = e + 1
3645         #m is the mantissa, e the exponent
3646         # 0.5 < |m| < 1.0
3647         m = m+t_rel*(10**-(self.__p_rel__))
3648         #bump m up
3649         m = round(m,self.__p_rel__)
3650         return m*(2.**e)
3651
3652    def round_down_rel2(self,value):
3653         t_rel=self.__t_rel__
3654         #Rounding down the value
3655         m,e = self.frexp(value)
3656         m = m/2
3657         e = e + 1
3658         #m is the mantissa, e the exponent
3659         # 0.5 < m < 1.0
3660         m = m-t_rel*(10**-(self.__p_rel__))
3661         #bump the |m| down, by 5% or whatever
3662         #self.p_rel dictates
3663         m = round(m,self.__p_rel__)
3664         return m*(2.**e)
3665
3666    def round_flat_rel2(self,value):
3667    #redundant
3668         m,e = self.frexp(value)
3669         m = m/2
3670         e = e + 1
3671         m = round(m,self.__p_rel__)
3672         return m*(2.**e)
3673
3674    def round_up_rel(self,value):
3675         t_rel=self.__t_rel__
3676         #Rounding up the value
3677         m,e = self.frexp(value)
3678         #m is the mantissa, e the exponent
3679         # 0.5 < |m| < 1.0
3680         m = m+t_rel*(10**-(self.__p_rel__))
3681         #bump m up
3682         m = round(m,self.__p_rel__)
3683         return m*(2.**e)
3684
3685    def round_down_rel(self,value):
3686         t_rel=self.__t_rel__
3687         #Rounding down the value
3688         m,e = self.frexp(value)
3689         #m is the mantissa, e the exponent
3690         # 0.5 < m < 1.0
3691         m = m-t_rel*(10**-(self.__p_rel__))
3692         #bump the |m| down, by 5% or whatever
3693         #self.p_rel dictates
3694         m = round(m,self.__p_rel__)
3695         return m*(2.**e)
3696
3697    def round_flat_rel(self,value):
3698    #redundant
3699         m,e = self.frexp(value)
3700         m = round(m,self.__p_rel__)
3701         return m*(2.**e)
3702
3703    def round_up_abs(self,value):
3704         t_abs=self.__t_abs__
3705         #Rounding up the value
3706         m = value+t_abs*(10**-(self.__p_abs__))
3707         #bump m up
3708         m = round(m,self.__p_abs__)
3709         return m
3710
3711    def round_down_abs(self,value):
3712         t_abs=self.__t_abs__
3713         #Rounding down the value
3714         m = value-t_abs*(10**-(self.__p_abs__))
3715         #bump the |m| down, by 5% or whatever
3716         #self.p_rel dictates
3717         m = round(m,self.__p_abs__)
3718         return m
3719
3720    def round_flat_abs(self,value):
3721    #redundant?
3722         m = round(value,self.__p_abs__)
3723         return m
3724
3725    def keys(self):
3726        return self.__items__.keys()
3727
3728
3729class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
3730    """
3731    This is a discretised tuple set, but
3732    based on a mapping. The mapping MUST
3733    return a sequence.
3734
3735    example:
3736    def weight(animal):
3737        return [animal.weight]
3738   
3739    a = Mapped_Discretised_Tuple_Set(weight)
3740    a.append[cow]
3741    a.append[fox]
3742    a.append[horse]
3743
3744    a[horse]    -> [cow,horse]
3745    a[dog]      -> [fox]
3746    a[elephant] -> []
3747    """
3748    def __init__(self,mapping,p_rel = 6, t_rel=0.01):
3749        Discretised_Tuple_Set.__init__\
3750         (self,p_rel,t_rel = t_rel)
3751        self.mapping = mapping
3752
3753    def rounded_keys(self,key):
3754        mapped_key = tuple(self.mapping(key))
3755        keys = self.__rounded_keys__(mapped_key)
3756        return keys
3757
3758class Affine_Linespace(Mapped_Discretised_Tuple_Set):
3759    """
3760    The affine linespace creates a way to record and compare lines.
3761    Precision is a bit of a hack, but it creates a way to avoid
3762    misses caused by round offs (between two lines of different
3763    lenghts, the short one gets rounded off more).
3764    I am starting to think that a quadratic search would be faster.
3765    Nearly.
3766    """
3767    def __init__(self,p_rel=4,t_rel=0.2):
3768        Mapped_Discretised_Tuple_Set.__init__\
3769            (self,self.affine_line,\
3770            p_rel=p_rel,t_rel=t_rel)
3771
3772        roundings = \
3773        [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
3774        self.roundings = roundings
3775        #roundings = \
3776        #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
3777        #self.roundings = roundings
3778
3779    def affine_line(self,line):
3780        point_1 = line[0]
3781        point_2 = line[1]
3782        #returns the equation of a line
3783        #between two points, in the from
3784        #(a,b,-c), as in ax+by-c=0
3785        #or line *dot* (x,y,1) = (0,0,0)
3786
3787        #Note that it normalises the line
3788        #(a,b,-c) so Line*Line = 1.
3789        #This does not change the mathematical
3790        #properties, but it makes comparism
3791        #easier.
3792
3793        #There are probably better algorithms.
3794        x1 = point_1[0]
3795        x2 = point_2[0]
3796        y1 = point_1[1]
3797        y2 = point_2[1]
3798        dif_x = x1-x2
3799        dif_y = y1-y2
3800
3801        if dif_x == dif_y == 0:
3802            msg = 'points are the same'
3803            raise msg
3804        elif abs(dif_x)>=abs(dif_y):
3805            alpha = (-dif_y)/dif_x
3806            #a = alpha * b
3807            b = -1.
3808            c = (x1*alpha+x2*alpha+y1+y2)/2.
3809            a = alpha*b
3810        else:
3811            beta = dif_x/(-dif_y)
3812            #b = beta * a
3813            a = 1.
3814            c = (x1+x2+y1*beta+y2*beta)/2.
3815            b = beta*a
3816        mag = abs(a)+abs(b)
3817        #This does not change the mathematical
3818        #properties, but it makes comparism possible.
3819
3820        #note that the gradient is b/a, or (a/b)**-1.
3821        #so
3822
3823        #if a == 0:
3824        #    sign_a = 1.
3825        #else:
3826        #    sign_a = a/((a**2)**0.5)
3827        #if b == 0:
3828        #    sign_b = 1.
3829        #else:
3830        #    sign_b = b/((b**2)**0.5)
3831        #if c == 0:
3832        #    sign_c = 1.
3833        #else:
3834        #    sign_c = c/((c**2)**0.5)
3835        #a = a/mag*sign_a
3836        #b = b/mag*sign_b
3837        #c = c/mag*sign_c
3838        a = a/mag
3839        b = b/mag
3840        c = c/mag
3841        return a,b,c
3842
3843
3844# FIXME (DSG-DSG)
3845# To do-
3846# Create a clear interface. eg
3847# have the interface methods more at the top of this file and add comments
3848# for the interface functions/methods, use function_name (not functionName),
3849
3850#Currently
3851#function_name methods assume absolute values.  Geo-refs can be passed in.
3852#
3853
3854# instead of functionName
3855if __name__ == "__main__":
3856    #from mesh import *
3857    # THIS CAN BE DELETED
3858    m = Mesh()
3859    dict = importUngenerateFile("ungen_test.txt")
3860    m.addVertsSegs(dict)
3861    print m3
Note: See TracBrowser for help on using the repository browser.