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

Last change on this file since 4888 was 4888, checked in by duncan, 16 years ago

tweaking the tagging of regions

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