source: trunk/anuga_core/source/anuga/pmesh/mesh.py @ 8279

Last change on this file since 8279 was 8010, checked in by steve, 14 years ago

Changed default boundary tag for holes to interior. Check out buildings.py in
user_manua/demos
File size: 85.0 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
24
25import numpy as num
26
27
28from anuga.coordinate_transforms.geo_reference import Geo_reference, \
29     DEFAULT_ZONE
30from  anuga.load_mesh.loadASCII import NOMAXAREA, export_mesh_file, \
31     import_mesh_file
32import anuga.alpha_shape.alpha_shape
33from anuga.geospatial_data.geospatial_data import Geospatial_data, \
34     ensure_geospatial, ensure_absolute, ensure_numeric
35from anuga.mesh_engine.mesh_engine import generate_mesh
36import anuga.utilities.log as log
37
38from anuga.file.ungenerate import load_ungenerate
39
40try: 
41    import kinds 
42except ImportError: 
43    # Hand-built mockup of the things we need from the kinds package, since it
44    # was recently removed from the standard numeric distro.  Some users may 
45    # not have it by default. 
46    class _bunch: 
47        pass 
48         
49    class _kinds(_bunch): 
50        default_float_kind = _bunch() 
51        default_float_kind.MIN = 2.2250738585072014e-308  #smallest +ve number
52        default_float_kind.MAX = 1.7976931348623157e+308 
53     
54    kinds = _kinds()
55   
56
57# 1st and third values must be the same
58# FIXME: maybe make this a switch that the user can change? - DSG
59initialconversions = ['', 'exterior', '']
60SEG_COLOUR = 'blue'
61
62
63class MeshObject:
64    """
65    An abstract superclass for the basic mesh objects, eg vertex, segment,
66    triangle.
67    """
68    def __init__(self):
69        pass
70   
71class Point(MeshObject): 
72    """
73    Define a point in a 2D space.
74    """
75    def __init__(self,X,Y):
76        __slots__ = ['x','y']
77        self.x=X
78        self.y=Y
79       
80    def DistanceToPoint(self, OtherPoint):
81        """
82        Returns the distance from this point to another
83        """
84        SumOfSquares = ((self.x - OtherPoint.x)**2) + \
85                       ((self.y - OtherPoint.y)**2)
86        return math.sqrt(SumOfSquares)
87       
88    def IsInsideCircle(self, Center, Radius):
89        """
90        Return 1 if this point is inside the circle,
91        0 otherwise
92        """
93       
94        if (self.DistanceToPoint(Center)<Radius):
95            return 1
96        else:
97            return 0
98       
99    def __repr__(self):
100        return "(%f,%f)" % (self.x,self.y) 
101
102    def cmp_xy(self, point):
103        if self.x < point.x:
104            return -1
105        elif self.x > point.x:
106            return 1
107        else:           
108            if self.y < point.y:
109                return -1
110            elif self.y > point.y:
111                return 1
112            else:
113                return 0
114       
115    def same_x_y(self, point):
116        if self.x == point.x and self.y == point.y:
117            return True
118        else:
119            return False
120       
121           
122
123class Vertex(Point):
124    """A point on the mesh.
125    Object attributes based on the Triangle program
126    """
127
128    VERTEXSQUARESIDELENGTH = 6
129
130    def __init__(self, X, Y, attributes=None):
131        __slots__ = ['x','y','attributes']
132
133        # we don't care what we get, as long as we can get float *value*
134        self.x = float(X)
135        self.y = float(Y)
136
137        self.attributes = [] 
138        if not attributes is None:
139            self.attributes = attributes
140   
141    def setAttributes(self,attributes):
142        """attributes is a list. """
143
144        self.attributes = attributes
145       
146    def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0,
147             yoffset =0, ):
148        x =  scale*(self.x + xoffset)
149        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
150        cornerOffset= self.VERTEXSQUARESIDELENGTH/2
151
152        # A hack to see the vert tags
153        # note: there will be many tags, since tags will not be removed
154        #when zooming
155        #canvas.create_text(x+ 2*cornerOffset,
156        #                   y+ 2*cornerOffset,
157        #                        text=tags)
158       
159        # This gives points info.  It is a mess though, since numbers are
160        # not deleted when zaooming in.
161        #canvas.create_text(x+ 2*cornerOffset,
162        #                   y+ 2*cornerOffset,
163        #                        text=str(x)+','+str(y))
164       
165        return canvas.create_rectangle(x-cornerOffset,
166                                       y-cornerOffset,
167                                       x+cornerOffset,
168                                       y+cornerOffset,
169                                       tags = tags,
170                                       outline=colour,
171                                       fill = 'white')
172   
173        #return tags
174     
175    def __repr__(self):
176        return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
177
178
179class Hole(Point):
180    """A region of the mesh were no triangles are generated.
181    Defined by a point in the hole enclosed by segments.
182    """
183
184    HOLECORNERLENGTH = 6
185   
186    def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0,
187             yoffset =0, ):
188        x =  scale*(self.x + xoffset)
189        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
190        cornerOffset= self.HOLECORNERLENGTH/2
191        return canvas.create_oval(x-cornerOffset,
192                                       y-cornerOffset,
193                                       x+cornerOffset,
194                                       y+cornerOffset,
195                                       tags = tags,
196                                       outline=colour,
197                                       fill = 'white')
198   
199class Region(Point):
200    """A region of the mesh.
201    Defined by a point in the region enclosed by segments. Used to tag areas.
202    """
203
204    CROSSLENGTH = 6
205    TAG = 0
206    MAXAREA = 1
207   
208    def __init__(self, X, Y, tag=None, maxArea=None):
209        """Precondition: tag is a string and maxArea is a real"""
210
211        self.x=X
212        self.y=Y   
213        self.attributes = [] # index 0 is the tag string
214                             # optional index 1 is the max triangle area
215                             # NOTE the size of this attribute is assumed
216                             # to be 1 or 2 in regionstrings2int
217        if tag is None:
218            self.attributes.append("")
219        else:
220            self.attributes.append(tag) #this is a string
221           
222        if maxArea is not None:
223            self.setMaxArea(maxArea) # maxArea is a number
224           
225    def getTag(self,):
226        return self.attributes[self.TAG]
227   
228    def setTag(self,tag):
229        self.attributes[self.TAG] = tag
230       
231    def getMaxArea(self):
232        """ Returns the Max Area of a Triangle or
233        None, if the Max Area has not been set.
234        """
235        if self.isMaxArea():
236            return self.attributes[1]
237        else:
238            return None
239   
240    def setMaxArea(self,MaxArea):
241        if MaxArea is not None:
242            if self.isMaxArea(): 
243                self.attributes[self.MAXAREA] = float(MaxArea)
244            else:
245                self.attributes.append( float(MaxArea) )
246   
247    def deleteMaxArea(self):
248        if self.isMaxArea():
249            self.attributes.pop(self.MAXAREA)
250           
251    def isMaxArea(self):
252        return len(self.attributes)> 1
253   
254    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0,
255             colour = "black"):
256        """
257        Draw a black cross, returning the objectID
258        """
259        x =  scale*(self.x + xoffset)
260        y = -1*scale*(self.y + yoffset) 
261        cornerOffset= self.CROSSLENGTH/2
262        return canvas.create_polygon(x,
263                                     y-cornerOffset,
264                                     x,
265                                     y,
266                                     x+cornerOffset,
267                                     y,
268                                     x,
269                                     y,
270                                     x,
271                                     y+cornerOffset,
272                                     x,
273                                     y,
274                                     x-cornerOffset,
275                                     y,
276                                     x,
277                                     y,
278                                     tags = tags,
279                                     outline = colour,fill = '')
280   
281    def __repr__(self):
282        if self.isMaxArea():
283            area = self.getMaxArea() 
284            return "(%f,%f,%s,%f)" % (self.x,self.y,
285                                      self.getTag(), self.getMaxArea())
286        else:
287            return "(%f,%f,%s)" % (self.x,self.y,
288                                   self.getTag())
289       
290class Segment(MeshObject):
291    """
292    Segments are edges whose presence in the triangulation is enforced.
293   
294    """
295    def __init__(self, vertex1, vertex2, tag = None ):
296        """
297        Each segment is specified by listing the vertices of its endpoints
298        The vertices are Vertex objects
299        """
300        assert(vertex1 != vertex2)
301        self.vertices = [vertex1,vertex2 ]
302       
303        if tag is None:
304            self.tag = self.__class__.default
305        else:
306            self.tag = tag #this is a string
307       
308    def __repr__(self):
309        return "[%s,%s]" % (self.vertices,self.tag)
310           
311       
312    def draw(self, canvas, tags,scale=1, xoffset=0,
313             yoffset=0, colour=SEG_COLOUR):
314        x=[]
315        y=[]
316        for end in self.vertices:
317            #end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
318            x.append(scale*(end.x + xoffset))
319            y.append(-1*scale*(end.y + yoffset)) # - since for a canvas - is up
320       
321        return canvas.create_line(x[0], y[0], x[1], y[1],
322                                  tags = tags,fill=colour)
323    def set_tag(self,tag):
324        self.tag = tag
325       
326    # Class methods
327    def set_default_tag(cls, default):
328        cls.default = default
329   
330    def get_default_tag(cls):
331        return cls.default
332   
333    set_default_tag = classmethod(set_default_tag) 
334    get_default_tag = classmethod(get_default_tag)
335
336Segment.set_default_tag("")       
337
338class Rigid_triangulation:
339    """
340    This is a triangulation that can't have triangles added or taken away.
341
342    It just represents the triangulation, not the mesh outline needed to
343    build the triangulation.
344
345      This is the guts of the data structure;
346        generated vertex list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
347        generated segment list: [(point1,point2),(p3,p4),...]
348            (Tuples of integers)
349        generated segment tag list: [tag,tag,...] list of strings
350
351        generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
352
353        generated triangle attribute list: [s1,s2,...] list of strings
354
355        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
356            tuple of triangles
357
358            Should I do the basic structure like triangle, general
359            mesh or something else? How about like triangle, since
360            that's where the info is from, and it'll fit easier into
361            this file..
362
363             Removing the loners is difficult, since all the vert's
364             after it must be removed.
365
366             This happens in set_triangulation.
367           
368    """
369
370    def __init__(self,
371                 triangles,
372                 segments,
373                 vertices,
374                 triangle_tags,
375                 triangle_neighbors,
376                 segment_tags,
377                 vertex_attributes,
378                 vertex_attribute_titles=None
379                 ):
380       
381        self.triangles = ensure_numeric(triangles) 
382        self.triangle_neighbors = ensure_numeric(triangle_neighbors)
383        self.triangle_tags = triangle_tags # list of strings
384        self.segments = ensure_numeric(segments) 
385       
386        self.segment_tags = segment_tags # list of strings
387        self.vertices = ensure_numeric(vertices)
388        # This is needed for __cmp__
389        if vertex_attributes is None:
390            self.vertex_attributes = []
391        else:
392            self.vertex_attributes = ensure_numeric(vertex_attributes)
393        if vertex_attribute_titles is None:
394            self.vertex_attribute_titles = []
395        else:
396            self.vertex_attribute_titles = vertex_attribute_titles
397       
398    def draw_triangulation(self, canvas, scale=1, xoffset=0, yoffset=0,
399             colour="green"):
400        """
401        Draw a triangle, returning the objectID
402        """
403       
404        # FIXME(DSG-DSG) This could be a data structure that is
405        # remembered and doesn't have any duplicates.
406        # but I wouldn't be able to use create_polygon.
407        # regard it as drawing a heap of segments
408       
409        for tri in self.triangles:
410            vertices = []
411            for v_index in range(3):
412                vertices.append(self.vertices[tri[v_index]])
413           
414            objectID = canvas.create_polygon(
415                scale*(vertices[1][0] + xoffset),
416                scale*-1*(vertices[1][1] + yoffset),
417                scale*(vertices[0][0] + xoffset),
418                scale*-1*(vertices[0][1] + yoffset),
419                scale*(vertices[2][0] + xoffset),
420                scale*-1*(vertices[2][1] + yoffset),
421                outline = colour,fill = ''
422                )
423           
424    def calc_mesh_area(self):
425        area = 0
426        for tri in self.triangles:
427            vertices = []
428            for v_index in range(3):
429                vertices.append(self.vertices[tri[v_index]])
430            ax = vertices[0][0]
431            ay = vertices[0][1]
432           
433            bx = vertices[1][0]
434            by = vertices[1][1]
435           
436            cx = vertices[2][0]
437            cy = vertices[2][1]
438           
439            area += abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
440        return area           
441       
442class Mesh:
443    """
444    Representation of a 2D triangular mesh.
445    User attributes describe the mesh region/segments/vertices/attributes
446
447    mesh attributes describe the mesh that is produced eg triangles and
448    vertices.
449    All point information is relative to the geo_reference passed in
450   
451   
452    """
453
454    def __repr__(self):
455        return """
456        mesh Triangles: %s 
457        mesh Attribute Titles: %s 
458        mesh Segments: %s 
459        mesh Vertices: %s 
460        user Segments: %s 
461        user Vertices: %s 
462        holes: %s 
463        regions: %s""" % (self.meshTriangles,
464                                self.attributeTitles,
465                                self.meshSegments,
466                                self.meshVertices,
467                                self.getUserSegments(),
468                                self.userVertices,
469                                self.holes,
470                                self.regions) 
471   
472    def __init__(self,
473                 userSegments=None,
474                 userVertices=None,
475                 holes=None,
476                 regions=None,
477                 geo_reference=None):
478        self.meshTriangles=[] 
479        self.attributeTitles=[] 
480        self.meshSegments=[]
481        self.meshVertices=[]
482
483        # Rigid
484        self.tri_mesh=None
485       
486       
487        self.visualise_graph = True
488
489        if userSegments is None:
490            self.userSegments=[]
491        else:
492            self.userSegments=userSegments
493        self.alphaUserSegments=[]
494           
495        if userVertices is None:
496            self.userVertices=[]
497        else:
498            self.userVertices=userVertices
499           
500        if holes is None:
501            self.holes=[]
502        else:
503            self.holes=holes
504           
505        if regions is None:
506            self.regions=[]
507        else:
508            self.regions=regions
509
510        if geo_reference is None:
511            self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0) 
512        else:
513            self.geo_reference = geo_reference
514
515        self.shape = None
516           
517    def __cmp__(self,other):
518       
519        # A dic for the initial m
520        dic = self.Mesh2triangList()
521        dic_mesh = self.Mesh2MeshList()
522        for element in dic_mesh.keys():
523            dic[element] = dic_mesh[element]
524        for element in dic.keys():
525            dic[element].sort()
526           
527        # A dic for the exported/imported m
528        dic_other = other.Mesh2triangList()
529        dic_mesh = other.Mesh2MeshList()
530        for element in dic_mesh.keys():
531            dic_other[element] = dic_mesh[element]
532        for element in dic.keys():
533            dic_other[element].sort()
534
535        #print "dsg************************8"
536        #print "dic ",dic
537        #print "*******8"
538        #print "mesh",dic_other
539        #print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
540        #print "dsg************************8"
541       
542        return (dic.__cmp__(dic_other))
543   
544    def addUserPoint(self, pointType, x,y):
545        if pointType == Vertex:
546            point = self.addUserVertex(x,y)
547        if pointType == Hole:
548            point = self._addHole(x,y)
549        if pointType == Region:
550            point = self._addRegion(x,y)
551        return point
552   
553    def addUserVertex(self, x,y):
554        v=Vertex(x, y)
555        self.userVertices.append(v)
556        return v
557
558    def _addHole(self, x,y):
559        h=Hole(x, y)
560        self.holes.append(h)
561        return h
562   
563    def add_hole(self, x,y, geo_reference=None):
564        """
565        adds a point, which represents a hole.
566
567        The point data can have it's own geo_refernece.
568        If geo_reference is None the data is asumed to be absolute
569        """
570        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
571                                                 points_geo_ref=geo_reference)
572        return self._addHole(x, y)
573
574    def _addRegion(self, x,y):
575        h=Region(x, y)
576        self.regions.append(h)
577        return h
578   
579    def add_region(self, x,y, geo_reference=None, tag=None):
580        """
581        adds a point, which represents a region.
582
583        The point data can have it's own geo_refernece.
584        If geo_reference is None the data is asumed to be absolute
585        """
586        #FIXME: have the user set the tag and resolution here,
587        # but still return the instance, just in case.
588        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
589                                                 points_geo_ref=geo_reference)
590        region =  self._addRegion(x, y)
591        if tag is not None:
592            region.setTag(tag)
593        return region
594
595    def build_grid(self,  vert_rows, vert_columns):
596        """
597        Build a grid with vert_rows number of vertex rows and
598        vert_columns number if vertex columns
599
600        Grid spacing of 1, the origin is the lower left hand corner.
601
602        FIXME(DSG-DSG) no test.
603        """
604
605        for i in range(vert_rows):
606            for j in range(vert_columns):
607                self.addUserVertex(j,i)
608        self.auto_segment()
609        self.generateMesh(mode = "Q", minAngle=20.0)
610       
611   
612    def add_vertices(self, point_data):
613        """
614        Add user vertices.
615
616        The point_data can be a list of (x,y) values, a numeric
617        array or a geospatial_data instance.
618        """
619        point_data = ensure_geospatial(point_data)
620        #print "point_data",point_data
621        # get points relative to the mesh geo_ref
622        points = point_data.get_data_points(geo_reference=self.geo_reference)
623   
624        for point in points:
625            v=Vertex(point[0], point[1])
626            self.userVertices.append(v)
627           
628    def add_hole_from_polygon(self, polygon, segment_tags=None,
629                              geo_reference=None):
630        """
631        Add a polygon with tags to the current mesh, as a region.
632        The maxArea of the region can be specified.
633
634        If a geo_reference of the polygon points is given, this is used.
635        If not;
636        The x,y info is assumed to be Easting and Northing, absolute,
637        for the meshes zone.
638
639        polygon a list of points, in meters that describe the polygon
640             (e.g. [[x1,y1],[x2,y2],...]
641        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
642
643        This returns the region instance, so if the user whats to modify
644        it they can.       
645        """
646        return self._add_area_from_polygon(polygon,
647                                           segment_tags=segment_tags,
648                                           hole=True,
649                                           geo_reference=geo_reference
650                                           )
651
652       
653    def add_region_from_polygon(self, polygon, segment_tags=None,
654                                max_triangle_area=None, geo_reference=None,
655                                region_tag=None):
656        """
657        Add a polygon with tags to the current mesh, as a region.
658        The maxArea of the region can be specified.
659
660        If a geo_reference of the polygon points is given, this is used.
661        If not;
662        The x,y info is assumed to be Easting and Northing, absolute,
663        for the meshes zone.
664
665        polygon a list of points, in meters that describe the polygon
666             (e.g. [[x1,y1],[x2,y2],...]
667        segment_tags (e.g.{'wall':[0,1,3],'ocean':[2]}) add tags to the
668        segements of the region
669        region_tags  - add a tag to all of the triangles in the region.
670
671        This returns the region instance (if a max_triangle_area is given),
672        so if the user whats to modify it they can.     
673        """
674        if max_triangle_area is None and region_tag is None:
675            create_region = False
676        else:
677            create_region = True
678           
679        region = self._add_area_from_polygon(polygon,
680                                             segment_tags=segment_tags,
681                                             geo_reference=geo_reference,
682                                             region=create_region)
683        if max_triangle_area is not None:
684            region.setMaxArea(max_triangle_area)
685        if region_tag is not None:
686            region.setTag(region_tag)
687       
688       
689        return region
690   
691       
692    def _add_area_from_polygon(self, polygon, segment_tags=None,
693                               geo_reference=None,
694                               hole=False,
695                               region=False):
696        """
697        Add a polygon with tags to the current mesh, as a region.
698        The maxArea of the region can be specified.
699
700        If a geo_reference of the polygon points is given, this is used.
701        If not;
702        The x,y info is assumed to be Easting and Northing, absolute,
703        for the meshes zone.
704
705        polygon a list of points, in meters that describe the polygon
706             (e.g. [[x1,y1],[x2,y2],...]
707        segment_tags (e.g.{'wall':[0,1,3],'ocean':[2]})add tags to the
708        segements of the region.
709
710        This returns the region instance, so if the user whats to modify
711        it they can.
712       
713        """
714        # Only import this if necessary.
715        # Trying to get pmesh working in an uncompiled environment
716        from anuga.geometry.polygon import point_in_polygon
717       
718        #get absolute values
719        if geo_reference is not None:
720            polygon = geo_reference.get_absolute(polygon)
721        # polygon is now absolute
722        #print "polygon  should be absolute",polygon
723       
724        #create points, segs and tags
725        region_dict = {}
726        region_dict['points'] = polygon
727       
728        #Create segments
729        #E.g. [[0,1], [1,2], [2,3], [3,0]]
730        #from polygon
731        #[0,1,2,3]
732        segments = []
733        N = len(polygon)
734        for i in range(N):
735            lo = i
736            hi = (lo + 1) % N
737            segments.append( [lo, hi] ) 
738        region_dict['segments'] = segments
739        region_dict['segment_tags'] = self._tag_dict2list(segment_tags, N,
740                                                          hole=hole)
741       
742   
743        self.addVertsSegs(region_dict) #this is passing absolute values
744
745        if region is True:
746            #get inner point - absolute values
747            inner_point = point_in_polygon(polygon)
748            inner = self.add_region(inner_point[0], inner_point[1],
749                                    geo_reference=None) 
750        elif hole is True:
751            #get inner point - absolute values
752            inner_point = point_in_polygon(polygon)
753            inner = self.add_hole(inner_point[0], inner_point[1],
754                                    geo_reference=None)
755        else:
756            inner = None
757           
758        return inner
759
760    def _tag_dict2list(self, tags, number_of_segs, hole=False):
761        """
762        Convert a tag dictionary from this sort of format;
763        #{'wall':[0,3],'ocean':[2]}
764
765        To a list format
766        # ['wall', '', 'ocean', 'wall']
767
768        Note: '' is a default value of nothing
769        """
770        # FIXME (DSG-DSG): Using '' as a default isn't good.
771        #Try None.
772        # Due to this default this method is too connected to
773        # _add_area_from_polygon
774
775        if hole:
776            default_tag ='interior'
777        else:
778            default_tag =''
779        segment_tags = [default_tag]*number_of_segs
780        if tags is not None:
781            for key in tags:
782                indices = tags[key]
783                for i in indices:
784                    segment_tags[i] = key
785        return segment_tags
786       
787    def add_circle(self, center, radius, segment_count=100,
788                   center_geo_reference=None, tag = None,
789                   region=False, hole=False):
790        """
791        center is a point, in absulute or relative to center_geo_ref
792        radius is the radius of the circle
793        segment_count is the number of segments used to represent the circle
794        tag is a string name, given to the segments.
795        If region is True a region object is returned.
796        If hole is True a hole object is returned.
797           (Don't have them both True.)
798
799           
800        """
801        # convert center and radius to a polygon
802        cuts = []
803        factor = 2* math.pi/segment_count
804        for cut in range(segment_count):
805             cuts.append(cut*factor)
806
807        polygon = []
808        for cut in cuts:
809           
810            x = center[0] + radius * math.cos(cut)
811            y = center[1] + radius * math.sin(cut)
812            polygon.append([x,y])
813        # build the tags
814        tags = {tag:range(segment_count)}
815       
816        return self._add_area_from_polygon(polygon, segment_tags=tags,
817                                           region=region, hole=hole,
818                                           geo_reference=center_geo_reference)
819
820    def auto_set_geo_reference(self):
821        """
822        Automatically set the georeference, based in the minimum x and y
823        user vertex values.
824
825        Not to be used with the graphical interface
826       
827        Not implemented.
828        Don't implement now.  using the new georeferenced points class
829        will change this?
830        """
831        #to do
832        # find the lower left hand corner
833        [xmin, ymin, xmax, ymax] = self.boxsize()
834
835        # set all points to that lower left hand corner.
836        # use change_geo_reference
837        new_geo = Geo_reference(self.geo_reference.get_zone(), xmin, ymin)
838        self.change_geo_reference(new_geo)
839       
840    def change_geo_reference(self, new_geo_reference):
841        """
842        Change from one geo_reference to another.
843        Not to be used with the graphical interface
844        """
845        # FIXME
846        # change georeference of;
847        #self.userVertices = self.geo_reference.change_points_geo_ref( \
848        #self.userVertices)
849        #self.holes = self.geo_reference.change_points_geo_ref(self.holes)
850        #self.regions = self.geo_reference.change_points_geo_ref(self.regions)
851        # The above will not work.
852        # since userVertices (etc) is a list of point objects,
853        #not a list of lists.
854        # add a method to the points class to fix this up.
855     
856    def add_segment(self, v1, v2, tag):
857        """
858        Don't do this function.
859        what will the v's be objects?  How is the user suppost to get them?
860
861        Indexes?  If add vertstosegs or add_region_from_polygon is called
862        more than once then the actual index is not obvious.  Making this
863        function confusing.
864        """
865        pass
866
867
868    def add_points_and_segments(self, points,
869                                  segments=None, segment_tags = None):
870        """
871        Add an outline of the mesh.
872        Points is a list of points a standard representation of points.
873        Segments is a list of tuples of integers.  Each tuple defines the
874           start and end of the segment by it's point index.
875        segment_tags is an optional dictionary which is used to add tags to
876           the segments.  The key is the tag name, value is the list of segment
877           indexes the tag will apply to.
878           eg. {'wall':[0,3],'ocean':[2]}
879           
880        """
881        # make sure the points are absolute
882        # Since addVertsSegs will deal with georeferencing.
883        points = ensure_absolute(points)
884
885        if segments is None:
886            segments = []
887            for i in range(len(points)-1):
888                segments.append([i, i+1])
889           
890        #create points, segs and tags
891        region_dict = {}
892        region_dict['points'] = points
893        region_dict['segments'] = segments
894        region_dict['segment_tags'] = self._tag_dict2list(segment_tags,
895                                                          len(segments))
896        self.addVertsSegs(region_dict)
897       
898    def addVertsSegs(self, outlineDict):
899        """
900        Add out-line (user Mesh) attributes given a dictionary of the lists
901        points: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
902        segments: [(point1,point2),(p3,p4),...] (Tuples of integers)
903        segment_tags: [S1Tag, S2Tag, ...] (list of strings)
904
905        Assume the values are in Eastings and Northings, with no reference
906        point. eg absolute
907        """
908        if not outlineDict.has_key('segment_tags'):
909            outlineDict['segment_tags'] = []
910            for i in range(len(outlineDict['segments'])):
911                outlineDict['segment_tags'].append('')
912        #print "outlineDict['segment_tags']",outlineDict['segment_tags']
913        #print "outlineDict['points']",outlineDict['points']
914        #print "outlineDict['segments']",outlineDict['segments']
915       
916        i_offset = len(self.userVertices)
917        #print "self.userVertices",self.userVertices
918        #print "index_offset",index_offset
919        for point in outlineDict['points']:
920            v=Vertex(point[0]-self.geo_reference.xllcorner,
921                     point[1]-self.geo_reference.yllcorner)
922            self.userVertices.append(v)
923           
924        for seg,seg_tag in map(None,outlineDict['segments'],
925                       outlineDict['segment_tags']):
926            segObject = Segment(self.userVertices[int(seg[0])+i_offset],
927                           self.userVertices[int(seg[1])+i_offset] )
928            if not seg_tag == '':
929                segObject.set_tag(seg_tag)
930            self.userSegments.append(segObject)
931           
932       
933    def get_triangle_count(self):
934        return len(self.getTriangulation())
935       
936    def getUserVertices(self):
937        """
938        Note: The x,y values will be relative to the mesh geo_reference
939        This returns a list of vertex objects
940        """
941        return self.userVertices
942
943    def get_user_vertices(self, absolute=True):
944        """
945        This returns a list of (x,y) values
946        (maybe it should return a geospatical object?
947        It makes it a bit confusing though.)
948        """
949        pointlist=[]
950        for vertex in self.userVertices:
951            pointlist.append([vertex.x,vertex.y])
952        spat = Geospatial_data(pointlist, geo_reference=self.geo_reference)
953        return spat.get_data_points(absolute=absolute)
954   
955    def getUserSegments(self):
956        allSegments = self.userSegments + self.alphaUserSegments
957        #print "self.userSegments",self.userSegments
958        #print "self.alphaUserSegments",self.alphaUserSegments
959        #print "allSegments",allSegments
960        return allSegments
961   
962    def deleteUserSegments(self,seg):
963        if self.userSegments.count(seg) == 0:
964            self.alphaUserSegments.remove(seg)
965            pass
966        else:
967            self.userSegments.remove(seg)
968           
969    def clearUserSegments(self):
970        self.userSegments = []
971        self.alphaUserSegments = []
972
973       #FIXME see where this is used. return an array instead
974    def getTriangulation(self):
975        #return self.meshTriangles
976        return self.tri_mesh.triangles.tolist()
977   
978    def getMeshVertices(self):
979        #return self.meshVertices
980        return self.tri_mesh.vertices
981 
982    def getMeshVerticeAttributes(self):
983        #return self.meshVertices
984        return self.tri_mesh.vertex_attributes
985   
986    def getMeshSegments(self):
987        #return self.meshSegments
988        return self.tri_mesh.segments
989   
990    def getMeshSegmentTags(self):
991        #return self.meshSegments
992        return self.tri_mesh.segment_tags
993   
994    def getHoles(self):
995        return self.holes
996   
997    def getRegions(self):
998        return self.regions
999   
1000    def isTriangulation(self):
1001        if self.meshVertices == []:
1002            return False 
1003        else:
1004            return True
1005   
1006    def addUserSegment(self, v1,v2):
1007        """
1008        PRECON: A segment between the two vertices is not already present.
1009        Check by calling isUserSegmentNew before calling this function.
1010       
1011        """
1012        s=Segment( v1,v2)
1013        self.userSegments.append(s)
1014        return s
1015       
1016    def generate_mesh(self,
1017                      maximum_triangle_area="",
1018                      minimum_triangle_angle=28.0,
1019                      verbose=True):
1020        if verbose is True:
1021            silent = ''
1022        else:
1023            silent = 'Q'
1024        self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle)
1025                                  +"a"+str(maximum_triangle_area)
1026                                  +"a")
1027        #The last a is so areas for regions will be used
1028       
1029    def generateMesh(self, mode = None, maxArea = None, minAngle= None,
1030                     isRegionalMaxAreas = True):
1031        """
1032        Based on the current user vaules, holes and regions
1033        generate a new mesh
1034        mode is a string that sets conditions on the mesh generations
1035        see triangle_instructions.txt for a definition of the commands
1036       
1037        PreCondition: maxArea is a double between 1e-20 and 1e30 or is a
1038        string.
1039        """
1040        #print "mode ",mode
1041        if mode == None:
1042            self.mode = ""
1043        else:
1044            self.mode = mode
1045       
1046        if self.mode.find('p') < 0:
1047            self.mode += 'p' #p - read a planar straight line graph.
1048            #there must be segments to use this switch
1049            # TODO throw an aception if there aren't seg's
1050            # it's more comlex than this.  eg holes
1051        if self.mode.find('z') < 0:
1052            self.mode += 'z' # z - Number all items starting from zero
1053                             # (rather than one)
1054        if self.mode.find('n'):
1055            self.mode += 'n' # n - output a list of neighboring triangles
1056        if self.mode.find('A') < 0:
1057            self.mode += 'A' # A - output region attribute list for triangles
1058
1059        if not self.mode.find('V') and not self.mode.find('Q'):
1060            self.mode += 'V' # V - output info about what Triangle is doing
1061       
1062        if self.mode.find('q') < 0 and minAngle is not None:
1063            #   print "**********8minAngle******** ",minAngle
1064            min_angle = 'q' + str(minAngle)
1065            self.mode += min_angle # z - Number all items starting from zero
1066                             # (rather than one)
1067        if maxArea != None:
1068            self.mode += 'a' + str(maxArea)
1069            try:
1070                self.mode += 'a' + '%20.20f' %maxArea
1071            except TypeError:
1072                self.mode += 'a' + str(maxArea)
1073            #print "self.mode", self.mode
1074        #FIXME (DSG-DSG) This isn't explained.
1075        if isRegionalMaxAreas:
1076            self.mode += 'a'
1077        #print "mesh#generateMesh# self.mode",self.mode 
1078        meshDict = self.Mesh2triangList()
1079
1080        #FIXME (DSG-DSG)  move below section into generate_mesh.py
1081        #                  & 4 functions eg segment_strings2ints
1082        # Actually, because of region_list.append((1.0,2.0,""))
1083        # don't move it, without careful thought
1084        #print "*************************!@!@ This is going to triangle   !@!@"
1085        #print meshDict
1086        #print "************************!@!@ This is going to triangle   !@!@"
1087
1088        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist']
1089        [meshDict['segmenttaglist'],
1090         segconverter] =  segment_strings2ints(meshDict['segmenttaglist'],
1091                                             initialconversions)
1092        #print "regionlist",meshDict['regionlist']
1093        [meshDict['regionlist'],
1094         regionconverter] =  region_strings2ints(meshDict['regionlist'])
1095        #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%regionlist",meshDict['regionlist']
1096        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist'
1097        #print "self.mode", self.mode
1098        generatedMesh = generate_mesh(
1099                              meshDict['pointlist'],
1100                              meshDict['segmentlist'],
1101                              meshDict['holelist'],
1102                              meshDict['regionlist'],
1103                              meshDict['pointattributelist'],
1104                              meshDict['segmenttaglist'],
1105                              self.mode,
1106                              meshDict['pointlist'])
1107        #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%generated",generatedMesh
1108        generatedMesh['qaa'] = 1
1109        if generatedMesh['generatedsegmentmarkerlist'] is not None:
1110            generatedMesh['generatedsegmentmarkerlist'] = \
1111              segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
1112                                   segconverter)
1113        #print "processed gen",generatedMesh['generatedsegmentmarkerlist']
1114        #print "pmesh mesh generatedMesh['generatedtriangleattributelist']", generatedMesh['generatedtriangleattributelist']
1115        if generatedMesh['generatedtriangleattributelist'] is not None:
1116            generatedMesh['generatedtriangleattributelist'] = \
1117            region_ints2strings(generatedMesh['generatedtriangleattributelist'],
1118                                  regionconverter)
1119
1120        #print "pmesh mesh generatedMesh['generatedtriangleattributelist']", generatedMesh['generatedtriangleattributelist']
1121        #FIXME (DSG-DSG)  move above section into generate_mesh.py
1122       
1123        if generatedMesh['generatedpointattributelist'] is None or \
1124               generatedMesh['generatedpointattributelist'].shape[1] ==0:
1125            self.attributeTitles = []
1126        generatedMesh['generatedpointattributetitlelist']= \
1127                                            self.attributeTitles
1128        #print "################  FROM TRIANGLE"
1129        #print "generatedMesh",generatedMesh
1130        #print "##################"
1131        self.setTriangulation(generatedMesh)
1132   
1133    def clearTriangulation(self):
1134
1135        #Clear the current generated mesh values
1136        self.meshTriangles=[] 
1137        self.meshSegments=[]
1138        self.meshVertices=[]
1139
1140    def removeDuplicatedUserVertices(self):
1141        """Pre-condition: There are no user segments
1142        This function will keep the first duplicate
1143        """
1144        assert self.getUserSegments() == []
1145        self.userVertices, counter =  self.removeDuplicatedVertices(
1146            self.userVertices)
1147        return counter
1148   
1149    def removeDuplicatedVertices(self, Vertices):
1150        """
1151        This function will keep the first duplicate, remove all others
1152        Precondition: Each vertex has a dupindex, which is the list
1153        index.
1154
1155        Note: this removes vertices that have the same x,y values,
1156        not duplicate instances in the Vertices list.
1157        """
1158        remove = []
1159        index = 0
1160        for v in Vertices:
1161            v.dupindex = index
1162            index += 1
1163        t = list(Vertices)
1164        t.sort(Point.cmp_xy)
1165   
1166        length = len(t)
1167        behind = 0
1168        ahead  = 1
1169        counter = 0
1170        while ahead < length:
1171            b = t[behind]
1172            ah = t[ahead]
1173            if (b.y == ah.y and b.x == ah.x):
1174                remove.append(ah.dupindex) 
1175            behind += 1
1176            ahead += 1
1177
1178        # remove the duplicate vertices
1179        remove.sort()
1180        remove.reverse()
1181        for i in remove:
1182            Vertices.pop(i)
1183            pass
1184
1185        #Remove the attribute that this function added
1186        for v in Vertices:
1187            del v.dupindex
1188        return Vertices,counter
1189
1190    # FIXME (DSG-DSG) Move this to geospatial
1191    def thinoutVertices(self, delta):
1192        """Pre-condition: There are no user segments
1193        This function will keep the first duplicate
1194        """
1195        assert self.getUserSegments() == []
1196        #t = self.userVertices
1197        #self.userVertices =[]
1198        boxedVertices = {}
1199        thinnedUserVertices =[]
1200        delta = round(delta,1)
1201       
1202        for v in self.userVertices :
1203            # tag is the center of the boxes
1204            tag = (round(v.x/delta,0)*delta,round(v.y/delta,0)*delta)
1205            #this creates a dict of lists of faces, indexed by tag
1206            boxedVertices.setdefault(tag,[]).append(v)
1207
1208        for [tag,verts] in boxedVertices.items():
1209            min = delta
1210            bestVert = None
1211            tagVert = Vertex(tag[0],tag[1])
1212            for v in verts:
1213                dist = v.DistanceToPoint(tagVert)
1214                if (dist<min):
1215                    min = dist
1216                    bestVert = v
1217            thinnedUserVertices.append(bestVert)
1218        self.userVertices =thinnedUserVertices
1219       
1220           
1221    def isUserSegmentNew(self, v1,v2):
1222        identicalSegs= [x for x in self.getUserSegments() \
1223                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1224        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1225       
1226        return len(identicalSegs) == 0
1227
1228       
1229    def deleteSegsOfVertex(self, delVertex):
1230        """
1231        Delete this vertex and any segments that connect to it.
1232        """
1233        #Find segments that connect to delVertex
1234        deletedSegments = []
1235        for seg in self.getUserSegments():
1236            if (delVertex in seg.vertices):
1237                deletedSegments.append(seg)
1238        # Delete segments that connect to delVertex
1239        for seg in deletedSegments:
1240            self.deleteUserSegments(seg)
1241        return deletedSegments
1242
1243   
1244    def deleteMeshObject(self, MeshObject):
1245        """
1246        Returns a list of all objects that were removed
1247        """
1248        deletedObs = []
1249        if isinstance(MeshObject, Vertex ):
1250            deletedObs = self.deleteSegsOfVertex(MeshObject)
1251            deletedObs.append(MeshObject)
1252            self.userVertices.remove(MeshObject)
1253        elif isinstance(MeshObject, Segment):
1254            deletedObs.append(MeshObject)
1255            self.deleteUserSegments(MeshObject)
1256        elif isinstance(MeshObject, Hole):
1257            deletedObs.append(MeshObject)
1258            self.holes.remove(MeshObject)
1259        elif isinstance(MeshObject, Region):
1260            deletedObs.append(MeshObject)
1261            self.regions.remove(MeshObject)         
1262        return deletedObs
1263                                                 
1264    def Mesh2triangList(self, userVertices=None,
1265                        userSegments=None,
1266                        holes=None,
1267                        regions=None):
1268        """
1269        Convert the Mesh to a dictionary of the lists needed for the
1270        triang module
1271        points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1272        pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
1273        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1274        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole)
1275        regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
1276       
1277        Note, this adds an index attribute to the user Vertex objects.
1278
1279        Used to produce output to triangle
1280        """
1281        if userVertices is None:
1282            userVertices = self.getUserVertices()
1283        if userSegments is None:
1284            userSegments = self.getUserSegments()
1285        if holes is None:
1286            holes = self.getHoles()
1287        if regions is None:
1288            regions = self.getRegions()
1289           
1290        meshDict = {}
1291       
1292        pointlist=[]
1293        pointattributelist=[]
1294        index = 0
1295        for vertex in userVertices:
1296            vertex.index = index
1297            pointlist.append((vertex.x,vertex.y))
1298            pointattributelist.append((vertex.attributes))
1299           
1300            index += 1
1301        meshDict['pointlist'] = pointlist
1302        meshDict['pointattributelist'] = pointattributelist
1303
1304        segmentlist=[]
1305        segmenttaglist=[]
1306        for seg in userSegments:
1307            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1308            segmenttaglist.append(seg.tag)
1309        meshDict['segmentlist'] =segmentlist
1310        meshDict['segmenttaglist'] =segmenttaglist
1311       
1312        holelist=[]
1313        for hole in holes:
1314            holelist.append((hole.x,hole.y)) 
1315        meshDict['holelist'] = holelist
1316       
1317        regionlist=[]
1318        for region in regions:
1319            if (region.getMaxArea() != None): 
1320                regionlist.append((region.x,region.y,region.getTag(),
1321                               region.getMaxArea()))
1322            else:
1323                regionlist.append((region.x,region.y,region.getTag()))
1324        meshDict['regionlist'] = regionlist
1325        #print "*(*("
1326        #print meshDict
1327        #print meshDict['regionlist']
1328        #print "*(*("
1329        return meshDict
1330                                               
1331    def Mesh2MeshList(self):
1332        """
1333        Convert the Mesh to a dictionary of lists describing the
1334        triangulation variables;
1335
1336        This is used by __cmp__
1337        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1338        generated point attribute list: [(a11,a12,...),(a21,a22),...]
1339            (Tuples of doubles)
1340        generated point attribute title list:[A1Title, A2Title ...]
1341            (list of strings)
1342        generated segment list: [(point1,point2),(p3,p4),...]
1343            (Tuples of integers)
1344        generated segment tag list: [tag,tag,...] list of strings
1345
1346        generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
1347
1348        generated triangle attribute list: [s1,s2,...] list of strings
1349
1350        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
1351            tuple of triangles
1352       
1353        Used to produce .tsh file
1354        """
1355       
1356        meshDict = {}   
1357        #print "old meshDict['generatedpointattributetitlelist']",meshDict['generatedpointattributetitlelist']
1358        #print "self.tri_mesh", self.tri_mesh
1359        if self.tri_mesh is not None:
1360            #print "self.tri_mesh.triangles.tolist()", self.tri_mesh.triangles.tolist()
1361            meshDict['generatedtrianglelist'] = self.tri_mesh.triangles.tolist()
1362            meshDict['generatedtriangleattributelist'] = self.tri_mesh.triangle_tags
1363            meshDict['generatedtriangleneighborlist'] = self.tri_mesh.triangle_neighbors.tolist()
1364            meshDict['generatedpointlist'] = self.tri_mesh.vertices.tolist()
1365            if  self.tri_mesh.vertex_attributes == []:
1366                meshDict['generatedpointattributelist'] = []
1367            #meshDict['generatedpointattributelist'] = self.tri_mesh.vertex_attributes
1368            meshDict['generatedpointattributetitlelist'] = \
1369                       self.tri_mesh.vertex_attribute_titles
1370            meshDict['generatedsegmentlist'] = self.tri_mesh.segments.tolist()
1371            meshDict['generatedsegmenttaglist'] = self.tri_mesh.segment_tags
1372        else:
1373            meshDict['generatedtrianglelist'] = []
1374            meshDict['generatedtriangleattributelist'] = []
1375            meshDict['generatedtriangleneighborlist'] = []
1376            meshDict['generatedpointlist'] = []
1377            meshDict['generatedpointattributelist'] = []
1378            meshDict['generatedpointattributetitlelist'] = []
1379            meshDict['generatedsegmentlist'] = []
1380            meshDict['generatedsegmenttaglist'] = []
1381           
1382        #print "new meshDict['generatedpointattributetitlelist']",meshDict['generatedpointattributetitlelist']
1383        #print "mesh.Mesh2MeshList*)*)"
1384        #print meshDict
1385        #print "mesh.Mesh2MeshList*)*)"
1386
1387        return meshDict
1388
1389                               
1390    def Mesh2MeshDic(self):
1391        """
1392        Convert the user and generated info of a mesh to a dictionary
1393        structure
1394        """
1395        dic = self.Mesh2triangList()
1396        dic_mesh = self.Mesh2MeshList()
1397        for element in dic_mesh.keys():
1398            dic[element] = dic_mesh[element]
1399        return dic
1400   
1401    def setTriangulation(self, genDict):
1402        """
1403        Set the mesh attributes given a dictionary of the lists
1404        returned from the triang module       
1405        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1406        generated point attribute list:[(P1att1,P1attt2, ...),
1407            (P2att1,P2attt2,...),...]- not implemented
1408        generated point attribute title list:[A1Title, A2Title ...]
1409            (list of strings) - not implemented
1410        generated segment list: [(point1,point2),(p3,p4),...]
1411            (Tuples of integers)
1412        generated segment marker list: [S1Tag, S2Tag, ...] (list of strings)
1413        triangle list:  [(point1,point2, point3),(p5,p4, p1),...]
1414            (Tuples of integers)
1415        triangle neighbor list: [(triangle1,triangle2, triangle3),
1416            (t5,t4, t1),...] (Tuples of integers)
1417            -1 means there's no triangle neighbor
1418        triangle attribute list: [(T1att), (T2att), ...](list of strings)
1419
1420        """
1421        # Setting up the rigid triangulation
1422        self.tri_mesh = Rigid_triangulation(
1423            genDict['generatedtrianglelist']
1424            ,genDict['generatedsegmentlist']
1425            ,genDict['generatedpointlist']
1426            ,genDict['generatedtriangleattributelist']
1427            ,genDict['generatedtriangleneighborlist']
1428            ,genDict['generatedsegmentmarkerlist']
1429            ,genDict['generatedpointattributelist']
1430            ,genDict['generatedpointattributetitlelist']
1431            )
1432           
1433    def setMesh(self, genDict):
1434        """
1435        Set the user Mesh attributes given a dictionary of the lists
1436        point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1437        point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1438        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1439        segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1440        region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1441        region attribute list: ["","reservoir",""] list of strings
1442        region max area list:[real, None, Real,...] list of None and reals
1443       
1444        mesh is an instance of a mesh object
1445        """
1446        #Clear the current user mesh values
1447        self.clearUserSegments()
1448        self.userVertices=[]
1449        self.Holes=[]
1450        self.Regions=[]
1451
1452        #print "mesh.setMesh@#@#@#"
1453        #print genDict
1454        #print "@#@#@#"
1455       
1456        #index = 0
1457        for point in genDict['pointlist']:
1458            v=Vertex(point[0], point[1])
1459            #v.index = index
1460            #index +=1
1461            self.userVertices.append(v)
1462
1463        #index = 0
1464        for seg,tag in map(None,genDict['segmentlist'],
1465                           genDict['segmenttaglist']):
1466            segObject = Segment( self.userVertices[int(seg[0])],
1467                           self.userVertices[int(seg[1])], tag = tag )
1468            #segObject.index = index
1469            #index +=1
1470            self.userSegments.append(segObject)
1471
1472# Remove the loading of attribute info.
1473# Have attribute info added using least_squares in pyvolution
1474#         index = 0
1475#         for att in genDict['pointattributelist']:
1476#             if att == None:
1477#                 self.userVertices[index].setAttributes([])
1478#             else:
1479#                 self.userVertices[index].setAttributes(att)
1480#            index += 1
1481       
1482        #index = 0
1483        for point in genDict['holelist']:
1484            h=Hole(point[0], point[1])
1485            #h.index = index
1486            #index +=1
1487            self.holes.append(h)
1488
1489        #index = 0
1490        for reg,att,maxArea in map(None,
1491                                   genDict['regionlist'],
1492                                   genDict['regionattributelist'],
1493                                   genDict['regionmaxarealist']):
1494            Object = Region( reg[0],
1495                             reg[1],
1496                             tag = att,
1497                             maxArea = maxArea)
1498            #Object.index = index
1499            #index +=1
1500            self.regions.append(Object)
1501           
1502    def Testauto_segment(self):
1503        newsegs = []
1504        s1 = Segment(self.userVertices[0],
1505                               self.userVertices[1])
1506        s2 = Segment(self.userVertices[0],
1507                               self.userVertices[2])
1508        s3 = Segment(self.userVertices[2],
1509                               self.userVertices[1])
1510        if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1511            newsegs.append(s1)
1512        if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1513            newsegs.append(s2)
1514        if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1515            newsegs.append(s3)
1516        #DSG!!!
1517        self.userSegments.extend(newsegs)
1518        return newsegs
1519
1520   
1521    def savePickle(self, currentName):
1522        fd = open(currentName, 'w')
1523        pickle.dump(self,fd)
1524        fd.close()
1525
1526    def auto_segmentFilter(self,raw_boundary=True,
1527                    remove_holes=False,
1528                    smooth_indents=False,
1529                    expand_pinch=False):
1530        """
1531        Change the filters applied on the alpha shape boundary
1532        """
1533        if self.shape is None:
1534            return [],[],0.0
1535        return self._boundary2mesh(raw_boundary=raw_boundary,
1536                    remove_holes=remove_holes,
1537                    smooth_indents=smooth_indents,
1538                    expand_pinch=expand_pinch)
1539   
1540       
1541   
1542    def auto_segment(self, alpha = None,
1543                    raw_boundary=True,
1544                    remove_holes=False,
1545                    smooth_indents=False,
1546                    expand_pinch=False): 
1547        """
1548        Precon: There must be 3 or more vertices in the userVertices structure
1549       
1550        This returns alpha_segs_no_user_segs, segs2delete, optimum_alpha
1551        """
1552        self._createBoundary(alpha=alpha)
1553        return self._boundary2mesh(raw_boundary=raw_boundary,
1554                    remove_holes=remove_holes,
1555                    smooth_indents=smooth_indents,
1556                    expand_pinch=expand_pinch)
1557
1558    def _createBoundary(self,alpha=None):
1559        """
1560        """
1561        points=[]
1562        for vertex in self.getUserVertices():
1563            points.append((vertex.x,vertex.y))
1564        self.shape = anuga.alpha_shape.alpha_shape.Alpha_Shape(points,
1565                                                               alpha=alpha)
1566
1567
1568    def _boundary2mesh(self, raw_boundary=True,
1569                    remove_holes=False,
1570                    smooth_indents=False,
1571                    expand_pinch=False):
1572        """
1573        Precon there must be a shape object.
1574        """
1575        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1576                                 remove_holes=remove_holes,
1577                                 smooth_indents=smooth_indents,
1578                                 expand_pinch=expand_pinch)
1579        boundary_segs = self.shape.get_boundary()
1580        #print "boundary_segs",boundary_segs
1581        segs2delete = self.alphaUserSegments
1582        #FIXME(DSG-DSG) this algorithm needs comments
1583        new_segs = {}
1584        #alpha_segs = []
1585        #user_segs = []
1586        for seg in boundary_segs:
1587            v1 = self.userVertices[int(seg[0])]
1588            v2 = self.userVertices[int(seg[1])]
1589            boundary_seg = Segment(v1, v2)
1590            new_segs[(v1,v2)] = boundary_seg
1591
1592        for user_seg in self.userSegments:
1593            if new_segs.has_key((user_seg.vertices[0],
1594                                user_seg.vertices[1])):
1595                del new_segs[user_seg.vertices[0],
1596                                user_seg.vertices[1]]
1597            elif new_segs.has_key((user_seg.vertices[1],
1598                                user_seg.vertices[0])):
1599                del new_segs[user_seg.vertices[1],
1600                                user_seg.vertices[0]]
1601               
1602        optimum_alpha = self.shape.get_alpha()
1603        alpha_segs_no_user_segs  = new_segs.values()
1604        self.alphaUserSegments = alpha_segs_no_user_segs
1605        return alpha_segs_no_user_segs, segs2delete, optimum_alpha
1606   
1607   
1608    def representedAlphaUserSegment(self, v1,v2):
1609        identicalSegs= [x for x in self.alphaUserSegments \
1610                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1611        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1612
1613        if identicalSegs == []:
1614            return None
1615        else:
1616            # Only return the first one.
1617            return identicalSegs[0]
1618   
1619    def representedUserSegment(self, v1,v2):
1620        identicalSegs= [x for x in self.userSegments \
1621                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1622        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1623
1624        if identicalSegs == []:
1625            return None
1626        else:
1627            # Only return the first one.
1628            return identicalSegs[0]
1629       
1630    def joinVertices(self):
1631        """
1632        Return list of segments connecting all userVertices
1633        in the order they were given
1634       
1635        Precon: There must be 3 or more vertices in the userVertices structure
1636        """
1637
1638        newsegs = []
1639       
1640        v1 = self.userVertices[0]
1641        for v2 in self.userVertices[1:]:
1642            if self.isUserSegmentNew(v1,v2):           
1643                newseg = Segment(v1, v2)
1644                newsegs.append(newseg)
1645            v1 = v2
1646
1647        #Connect last point to the first
1648        v2 = self.userVertices[0]       
1649        if self.isUserSegmentNew(v1,v2):           
1650                newseg = Segment(v1, v2)
1651                newsegs.append(newseg)
1652           
1653
1654        #Update list of user segments       
1655        #DSG!!!
1656        self.userSegments.extend(newsegs)               
1657        return newsegs
1658   
1659    def normaliseMesh(self,scale, offset, height_scale):
1660        [xmin, ymin, xmax, ymax] = self.boxsize()
1661        [attmin0, attmax0] = self.maxMinVertAtt(0)
1662        #print "[attmin0, attmax0]" ,[attmin0, attmax0]
1663        [attmin1, attmax1] = self.maxMinVertAtt(1)
1664        #print [xmin, ymin, xmax, ymax]
1665        xrange = xmax - xmin
1666        yrange = ymax - ymin
1667        if xrange > yrange:
1668            min,max = xmin, xmax
1669        else:
1670            min,max = ymin, ymax
1671           
1672        for obj in self.getUserVertices():
1673            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1674            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1675            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1676                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1677                                    (attmax0-attmin0)*height_scale
1678            if len(obj.attributes) > 1 and attmin1 != attmax1:
1679                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1680                                    (attmax1-attmin1)*height_scale
1681           
1682        for obj in self.getMeshVertices():
1683            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1684            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1685            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1686                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1687                                    (attmax0-attmin0)*height_scale
1688            if len(obj.attributes) > 1 and attmin1 != attmax1:
1689                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1690                                    (attmax1-attmin1)*height_scale
1691               
1692        for obj in self.getHoles():
1693            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1694            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1695        for obj in self.getRegions():
1696            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1697            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1698        [xmin, ymin, xmax, ymax] = self.boxsize()
1699        #print [xmin, ymin, xmax, ymax]
1700     
1701    def boxsizeVerts(self):
1702        """
1703        Returns a list of verts denoting a box or triangle that contains
1704        verts on the xmin, ymin, xmax and ymax axis.
1705        Structure: list of verts
1706        """
1707       
1708        large = kinds.default_float_kind.MAX
1709        xmin= large
1710        xmax=-large
1711        ymin= large
1712        ymax=-large
1713        for vertex in self.userVertices:
1714            if vertex.x < xmin:
1715                xmin = vertex.x
1716                xminVert = vertex
1717            if vertex.x > xmax:
1718                xmax = vertex.x
1719                xmaxVert = vertex
1720               
1721            if vertex.y < ymin:
1722                ymin = vertex.y
1723                yminVert = vertex
1724            if vertex.y > ymax:
1725                ymax = vertex.y
1726                ymaxVert = vertex
1727        verts, count = self.removeDuplicatedVertices([xminVert,
1728                                                      xmaxVert,
1729                                                      yminVert,
1730                                                      ymaxVert])
1731         
1732        return verts
1733   
1734    def boxsize(self):
1735        """
1736        Returns a list denoting a box that contains the entire
1737        structure of vertices
1738        Structure: [xmin, ymin, xmax, ymax]
1739        """
1740     
1741        large = kinds.default_float_kind.MAX
1742        xmin= large
1743        xmax=-large
1744        ymin= large
1745        ymax=-large
1746        for vertex in self.userVertices:
1747            if vertex.x < xmin:
1748                xmin = vertex.x
1749            if vertex.x > xmax:
1750                xmax = vertex.x
1751               
1752            if vertex.y < ymin:
1753                ymin = vertex.y
1754            if vertex.y > ymax:
1755                ymax = vertex.y
1756        return [xmin, ymin, xmax, ymax]
1757 
1758    def maxMinVertAtt(self, iatt):
1759        """
1760        Returns a list denoting a box that contains the entire structure
1761        of vertices
1762        Structure: [xmin, ymin, xmax, ymax]
1763        """
1764       
1765        large = kinds.default_float_kind.MAX
1766        min= large
1767        max=-large
1768        for vertex in self.userVertices:
1769            if len(vertex.attributes) > iatt:
1770                if vertex.attributes[iatt] < min:
1771                    min = vertex.attributes[iatt]
1772                if vertex.attributes[iatt] > max:
1773                    max = vertex.attributes[iatt] 
1774        for vertex in self.meshVertices:
1775            if len(vertex.attributes) > iatt:
1776                if vertex.attributes[iatt] < min:
1777                    min = vertex.attributes[iatt]
1778                if vertex.attributes[iatt] > max:
1779                    max = vertex.attributes[iatt] 
1780        return [min, max]
1781   
1782    def scaleoffset(self, WIDTH, HEIGHT):
1783        """
1784        Returns a list denoting the scale and offset terms that need to be
1785        applied when converting  mesh co-ordinates onto grid co-ordinates
1786        Structure: [scale, xoffset, yoffset]
1787        """   
1788        OFFSET = 0.05*min([WIDTH, HEIGHT])
1789        [xmin, ymin, xmax, ymax] = self.boxsize()
1790        SCALE = min([0.9*WIDTH, 0.9*HEIGHT])/max([xmax-xmin, ymax-ymin])
1791       
1792        if SCALE*xmin < OFFSET:
1793            xoffset = abs(SCALE*xmin) + OFFSET
1794        if SCALE*xmax > WIDTH - OFFSET:
1795            xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1796        if SCALE*ymin < OFFSET:
1797            b = abs(SCALE*ymin)+OFFSET
1798        if SCALE*ymax > HEIGHT-OFFSET:
1799            b = -(SCALE*ymax - HEIGHT + OFFSET)
1800        yoffset = HEIGHT - b
1801        return [SCALE, xoffset, yoffset]
1802
1803         
1804    def exportASCIIobj(self,ofile):
1805        """
1806        export a file, ofile, with the format
1807         lines:  v <x> <y> <first attribute>
1808        f <vertex #>  <vertex #> <vertex #> (of the triangles)
1809        """
1810        fd = open(ofile,'w')
1811        self.writeASCIIobj(fd)   
1812        fd.close()
1813
1814
1815    def writeASCIIobj(self,fd):
1816        fd.write(" # Triangulation as an obj file\n")
1817        numVert = str(len(self.meshVertices))
1818       
1819        index1 = 1 
1820        for vert in self.meshVertices:
1821            vert.index1 = index1
1822            index1 += 1
1823           
1824            fd.write("v "
1825                     + str(vert.x) + " "
1826                     + str(vert.y) + " "
1827                     + str(vert.attributes[0]) + "\n")
1828           
1829        for tri in self.meshTriangles:
1830            fd.write("f "
1831                     + str(tri.vertices[0].index1) + " " 
1832                     + str(tri.vertices[1].index1) + " " 
1833                     + str(tri.vertices[2].index1) + "\n")
1834           
1835    def exportASCIIsegmentoutlinefile(self,ofile):
1836        """Write the boundary user mesh info, eg
1837         vertices that are connected to segments,
1838         segments
1839        """
1840       
1841        verts = {}
1842        for seg in self.getUserSegments():
1843            verts[seg.vertices[0]] = seg.vertices[0]
1844            verts[seg.vertices[1]] = seg.vertices[1]
1845        meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values())
1846        export_mesh_file(ofile,meshDict)
1847       
1848        # exportASCIImeshfile   - this function is used
1849    def export_mesh_file(self,ofile):
1850        """
1851        export a file, ofile, with the format
1852        """
1853       
1854        dict = self.Mesh2IODict()
1855        export_mesh_file(ofile,dict)
1856
1857    # FIXME(DSG-DSG):Break this into two functions.
1858    #One for the outline points.
1859    #One for the mesh points.
1860    # Note: this function is not in the gui
1861    def exportPointsFile(self,ofile):
1862        """
1863        export a points file, ofile.
1864       
1865        """
1866       
1867        mesh_dict = self.Mesh2IODict()
1868        #point_dict = {}
1869        #point_dict['attributelist'] = {} #this will need to be expanded..
1870                                         # if attributes are brought back in.
1871        #point_dict['geo_reference'] = self.geo_reference
1872        if mesh_dict['vertices'] == []:
1873            #point_dict['pointlist'] = mesh_dict['points']
1874            geo = Geospatial_data(mesh_dict['points'],
1875                                  geo_reference=self.geo_reference)
1876        else:
1877            #point_dict['pointlist'] = mesh_dict['vertices']
1878            geo = Geospatial_data(mesh_dict['vertices'],
1879                                  geo_reference=self.geo_reference)
1880
1881        geo.export_points_file(ofile, absolute=True)
1882       
1883
1884
1885    def import_ungenerate_file(self,ofile, tag=None, region_tag=None):
1886        """
1887        Imports an ungenerate file, from arcGIS into mesh.
1888        This file describes many polygons.
1889
1890        ofile is the name of the ungenerated file.
1891        Tag is a string name to be taggged on each segment.
1892       
1893        region_tag is the tag applied to the polygon regions.
1894        if it is a string the one value will be assigned to all regions
1895        if it is a list the first value in the list will be applied to the first polygon etc.
1896        WARNING: size of list and number of polygons isn't checked
1897
1898        WARNING values are assumed to be absolute.
1899        geo-refs are not taken into account..
1900        """
1901   
1902        dict = load_ungenerate(ofile)
1903        default_tag = Segment.get_default_tag()
1904        if tag is not None:
1905            Segment.set_default_tag(str(tag))
1906       
1907        if region_tag is None:
1908            self.addVertsSegs(dict)
1909        else:           
1910            if not isinstance(region_tag, list):
1911                region_tag = [region_tag]*len(dict['polygons'])
1912            for a_tag,polygon in map(None, region_tag, dict['polygons']): 
1913                segment_tags = {tag:range(len(polygon))}
1914                self.add_region_from_polygon(polygon,segment_tags=segment_tags,
1915                                region_tag=a_tag)
1916           
1917       
1918        Segment.set_default_tag(default_tag)
1919
1920        # change the tag back to it's default
1921   
1922       
1923########### IO CONVERTERS ##################
1924        """
1925        The dict fromat for IO with .tsh files is;
1926        (the triangulation)
1927        vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
1928        vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1929        vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
1930        segments: [[v1,v2],[v3,v4],...] (lists of integers)
1931        segment_tags : [tag,tag,...] list of strings
1932        triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
1933        triangle tags: [s1,s2,...] A list of strings
1934        triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
1935       
1936        (the outline)   
1937        points: [[x1,y1],[x2,y2],...] (lists of doubles)
1938        point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1939        outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
1940        outline_segment_tags : [tag1,tag2,...] list of strings
1941        holes : [[x1,y1],...](List of doubles, one inside each hole region)
1942        regions : [ [x1,y1],...] (List of 4 doubles)
1943        region_tags : [tag1,tag2,...] (list of strings)
1944        region_max_areas: [ma1,ma2,...] (A list of doubles)
1945        {Convension: A -ve max area means no max area}
1946       
1947        """
1948     
1949
1950                               
1951    def Mesh2IODict(self):
1952        """
1953        Convert the triangulation and outline info of a mesh to a dictionary
1954        structure
1955        """
1956        dict = self.Mesh2IOTriangulationDict()
1957        dict_mesh = self.Mesh2IOOutlineDict()
1958        for element in dict_mesh.keys():
1959            dict[element] = dict_mesh[element]
1960
1961        # add the geo reference
1962        dict['geo_reference'] = self.geo_reference
1963        return dict
1964   
1965    def Mesh2IOTriangulationDict(self):
1966        """
1967        Convert the Mesh to a dictionary of lists describing the
1968        triangulation variables;
1969       
1970        Used to produce .tsh file
1971        """
1972        meshDict = {}
1973        if self.tri_mesh is not None:
1974            meshDict['triangles'] = self.tri_mesh.triangles
1975            meshDict['triangle_tags'] = self.tri_mesh.triangle_tags
1976            #print "mesh meshDict['triangle_tags']", meshDict['triangle_tags']
1977            meshDict['triangle_neighbors'] = self.tri_mesh.triangle_neighbors
1978            meshDict['vertices'] = self.tri_mesh.vertices
1979            meshDict['vertex_attributes'] = self.tri_mesh.vertex_attributes
1980            meshDict['vertex_attribute_titles'] = \
1981                               self.tri_mesh.vertex_attribute_titles
1982            meshDict['segments'] = self.tri_mesh.segments
1983            meshDict['segment_tags'] = self.tri_mesh.segment_tags
1984        else:
1985            meshDict['triangles'] = []
1986            meshDict['triangle_tags'] = []
1987            meshDict['triangle_neighbors'] = []
1988            meshDict['vertices'] = []
1989            meshDict['vertex_attributes'] = []
1990            meshDict['vertex_attribute_titles'] = []
1991            meshDict['segments'] = []
1992            meshDict['segment_tags'] = []
1993        #print "mesh.Mesh2IOTriangulationDict*)*)"
1994        #print meshDict
1995        #print "mesh.Mesh2IOTriangulationDict*)*)"
1996
1997        return meshDict
1998
1999                                                     
2000    def Mesh2IOOutlineDict(self, userVertices=None,
2001                        userSegments=None,
2002                        holes=None,
2003                        regions=None):
2004        """
2005        Convert the mesh outline to a dictionary of the lists needed for the
2006        triang module;
2007       
2008        Note, this adds an index attribute to the user Vertex objects.
2009
2010        Used to produce .tsh file and output to triangle
2011        """
2012       
2013        if userVertices is None:
2014            userVertices = self.getUserVertices()
2015        if userSegments is None:
2016            userSegments = self.getUserSegments()
2017        if holes is None:
2018            holes = self.getHoles()
2019        if regions is None:
2020            regions = self.getRegions()
2021           
2022        meshDict = {}
2023        #print "userVertices",userVertices
2024        #print "userSegments",userSegments
2025        pointlist=[]
2026        pointattributelist=[]
2027        index = 0
2028        for vertex in userVertices:
2029            vertex.index = index
2030            pointlist.append([vertex.x,vertex.y])
2031            pointattributelist.append(vertex.attributes)
2032           
2033            index += 1
2034        meshDict['points'] = pointlist
2035        meshDict['point_attributes'] = pointattributelist
2036
2037        segmentlist=[]
2038        segmenttaglist=[]
2039        for seg in userSegments:
2040            segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
2041            segmenttaglist.append(seg.tag)
2042        meshDict['outline_segments'] =segmentlist
2043        meshDict['outline_segment_tags'] =segmenttaglist
2044       
2045        holelist=[]
2046        for hole in holes:
2047            holelist.append([hole.x,hole.y]) 
2048        meshDict['holes'] = holelist
2049       
2050        regionlist=[]
2051        regiontaglist = []
2052        regionmaxarealist = []
2053        for region in regions:
2054            regionlist.append([region.x,region.y])
2055            regiontaglist.append(region.getTag())
2056           
2057            if (region.getMaxArea() != None): 
2058                regionmaxarealist.append(region.getMaxArea())
2059            else:
2060                regionmaxarealist.append(NOMAXAREA)
2061        meshDict['regions'] = regionlist
2062        meshDict['region_tags'] = regiontaglist
2063        meshDict['region_max_areas'] = regionmaxarealist
2064        #print "*(*("
2065        #print meshDict
2066        #print meshDict['regionlist']
2067        #print "*(*("
2068        return meshDict
2069
2070    def IOTriangulation2Mesh(self, genDict):
2071        """
2072        Set the mesh attributes given an tsh IO dictionary
2073        """
2074        #Clear the current generated mesh values
2075        self.tri_mesh = None
2076       
2077        self.tri_mesh = Rigid_triangulation(
2078            genDict['triangles']
2079            ,genDict['segments']
2080            ,genDict['vertices']
2081            ,genDict['triangle_tags']
2082            ,genDict['triangle_neighbors']
2083            ,genDict['segment_tags']
2084            ,genDict['vertex_attributes']
2085            ,genDict['vertex_attribute_titles']
2086            )
2087        self.attributeTitles = genDict['vertex_attribute_titles']
2088        self.maxVertexIndex = len(genDict['vertices'])
2089        #print "self.maxVertexIndex ", self.maxVertexIndex
2090       
2091    def IOOutline2Mesh(self, genDict):
2092        """
2093        Set the outline (user Mesh attributes) given a IO tsh dictionary
2094       
2095        mesh is an instance of a mesh object
2096        """
2097        #Clear the current user mesh values
2098        self.clearUserSegments()
2099        self.userVertices=[]
2100        self.Holes=[]
2101        self.Regions=[]
2102
2103        #print "mesh.IOOutline2Mesh@#@#@#"
2104        #print "genDict",genDict
2105        #print "@#@#@#"
2106       
2107        #index = 0
2108        for point in genDict['points']:
2109            v=Vertex(point[0], point[1])
2110            #v.index = index
2111            #index +=1
2112            self.userVertices.append(v)
2113
2114        #index = 0
2115        for seg,tag in map(None,genDict['outline_segments'],
2116                           genDict['outline_segment_tags']):
2117
2118            segObject = Segment( self.userVertices[int(seg[0])],
2119                           self.userVertices[int(seg[1])], tag = tag )
2120            #segObject.index = index
2121            #index +=1
2122            self.userSegments.append(segObject)
2123
2124# Remove the loading of attribute info.
2125# Have attribute info added using least_squares in pyvolution
2126#         index = 0
2127#         for att in genDict['point_attributes']:
2128#             if att == None:
2129#                 self.userVertices[index].setAttributes([])
2130#             else:
2131#                 self.userVertices[index].setAttributes(att)
2132#            index += 1
2133       
2134        #index = 0
2135        for point in genDict['holes']:
2136            h=Hole(point[0], point[1])
2137            #h.index = index
2138            #index +=1
2139            self.holes.append(h)
2140
2141        #index = 0
2142        for reg,att,maxArea in map(None,
2143                                   genDict['regions'],
2144                                   genDict['region_tags'],
2145                                   genDict['region_max_areas']):
2146            if maxArea > 0:  # maybe I should ref NOMAXAREA? Prob' not though
2147                Object = Region( reg[0],
2148                                 reg[1],
2149                                 tag = att,
2150                                 maxArea = maxArea)
2151            else:
2152                Object = Region( reg[0],
2153                                 reg[1],
2154                                 tag = att)
2155               
2156            #Object.index = index
2157            #index +=1
2158            self.regions.append(Object)
2159
2160
2161def importMeshFromFile(ofile):
2162    """returns a mesh object, made from a points file or .tsh/.msh file
2163    Often raises IOError,RuntimeError
2164    """
2165    newmesh = None
2166    if (ofile[-4:]== ".pts" or ofile[-4:]== ".txt" or \
2167        ofile[-4:]== ".csv"):
2168        geospatial = Geospatial_data(ofile)
2169        dict = {}
2170        dict['points'] = geospatial.get_data_points(absolute=False)
2171        dict['outline_segments'] = []
2172        dict['outline_segment_tags'] = []
2173        dict['regions'] = []
2174        dict['region_tags'] = []
2175        dict['region_max_areas'] = []
2176        dict['holes'] = [] 
2177        newmesh= Mesh(geo_reference = geospatial.geo_reference)
2178        newmesh.IOOutline2Mesh(dict)
2179        counter = newmesh.removeDuplicatedUserVertices()
2180        if (counter >0):
2181            log.critical("%i duplicate vertices removed from dataset" % counter)
2182    elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
2183        dict = import_mesh_file(ofile)
2184        #print "********"
2185        #print "zq mesh.dict",dict
2186        #print "********"
2187        newmesh= Mesh()
2188        newmesh.IOOutline2Mesh(dict)
2189        newmesh.IOTriangulation2Mesh(dict)
2190    else:
2191        raise RuntimeError
2192   
2193    if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
2194        newmesh.geo_reference = dict['geo_reference']
2195    return newmesh
2196
2197def loadPickle(currentName):
2198    fd = open(currentName)
2199    mesh = pickle.load(fd)
2200    fd.close()
2201    return mesh
2202   
2203def square_outline(side_length = 1,up = "top", left = "left", right = "right",
2204                   down = "bottom", regions = False):
2205   
2206        a = Vertex (0,0)
2207        b = Vertex (0,side_length)
2208        c = Vertex (side_length,0)
2209        d = Vertex (side_length,side_length)
2210     
2211        s2 = Segment(b,d, tag = up)
2212        s3 = Segment(b,a, tag = left)
2213        s4 = Segment(d,c, tag = right)
2214        s5 = Segment(a,c, tag = down)
2215
2216        if regions:
2217            e =  Vertex (side_length/2,side_length/2)
2218            s6 = Segment(a,e, tag = down + left)
2219            s7 = Segment(b,e, tag = up + left)
2220            s8 = Segment(c,e, tag = down + right)
2221            s9 = Segment(d,e, tag = up + right)
2222            r1 = Region(side_length/2,3.*side_length/4, tag = up)
2223            r2 = Region(1.*side_length/4,side_length/2, tag = left)
2224            r3 = Region(3.*side_length/4,side_length/2, tag = right)
2225            r4 = Region(side_length/2,1.*side_length/4, tag = down)
2226            mesh = Mesh(userVertices=[a,b,c,d,e],
2227                        userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
2228                        regions = [r1,r2,r3,r4])
2229        else:
2230            mesh = Mesh(userVertices=[a,b,c,d],
2231                 userSegments=[s2,s3,s4,s5])
2232     
2233        return mesh
2234
2235   
2236
2237def region_strings2ints(region_list):
2238    """Given a list of (x_int,y_int,tag_string) lists it returns a list of
2239    (x_int,y_int,tag_int) and a list to convert the tag_int's back to
2240    the tag_strings 
2241    """
2242    # Make sure "" has an index of 0
2243    region_list.reverse()
2244    region_list.append((1.0,2.0,""))
2245    region_list.reverse()
2246    convertint2string = []
2247    for i in xrange(len(region_list)):
2248        convertint2string.append(region_list[i][2])
2249        if len(region_list[i]) == 4: # there's an area value
2250            region_list[i] = (region_list[i][0],
2251                              region_list[i][1],i,region_list[i][3])
2252        elif len(region_list[i]) == 3: # no area value
2253            region_list[i] = (region_list[i][0],region_list[i][1],i)
2254        else:
2255            log.critical("The region list has a bad size")
2256            # raise an error ..
2257            raise Error
2258
2259    #remove "" from the region_list
2260    region_list.pop(0)
2261   
2262    return [region_list, convertint2string]
2263
2264def region_ints2strings(region_list,convertint2string):
2265    """Reverses the transformation of region_strings2ints
2266    """
2267    #print 'region_ints2strings region_list', region_list
2268   
2269    returned_region_list = []
2270    # may not need (not region_list[0] == [])
2271    # or region_list[0] == [0.0]
2272    if (not region_list[0] == []): # or region_list[0] == [0.0]:
2273        #print "in loop"
2274        for i in xrange(len(region_list)):
2275            temp = region_list[i]
2276            returned_region_list.append(convertint2string[int(temp[0])])
2277    return returned_region_list
2278
2279def segment_ints2strings(intlist, convertint2string):
2280    """Reverses the transformation of segment_strings2ints """
2281    stringlist = []
2282    for x in intlist:
2283        stringlist.append(convertint2string[int(x)])
2284    return stringlist
2285
2286def segment_strings2ints(stringlist, preset):
2287    """Given a list of strings return a list of 0 to n ints which represent
2288    the strings and a converting list of the strings, indexed by 0 to n ints.
2289    Also, input an initial converting list of the strings
2290    Note, the converting list starts off with
2291    ["internal boundary", "external boundary", "internal boundary"]
2292    example input and output
2293    input ["a","b","a","c"],["c"]
2294    output [[2, 1, 2, 0], ['c', 'b', 'a']]
2295
2296    the first element in the converting list will be
2297    overwritten with "".
2298    ?This will become the third item in the converting list?
2299   
2300    # note, order the initial converting list is important,
2301     since the index = the int tag
2302    """
2303    nodups = unique(stringlist)
2304    # note, order is important, the index = the int tag
2305    #preset = ["internal boundary", "external boundary"]
2306    #Remove the preset tags from the list with no duplicates
2307    nodups = [x for x in nodups if x not in preset]
2308
2309    try:
2310        nodups.remove("") # this has to go to zero
2311    except ValueError:
2312        pass
2313   
2314    # Add the preset list at the beginning of no duplicates
2315    preset.reverse()
2316    nodups.extend(preset)
2317    nodups.reverse()
2318
2319       
2320    convertstring2int = {}
2321    convertint2string = []
2322    index = 0
2323    for x in nodups:
2324        convertstring2int[x] = index
2325        convertint2string.append(x)
2326        index += 1
2327    convertstring2int[""] = 0
2328
2329    intlist = []
2330    for x in stringlist:
2331        intlist.append(convertstring2int[x])
2332    return [intlist, convertint2string]
2333
2334
2335def unique(s):
2336    """Return a list of the elements in s, but without duplicates.
2337
2338    For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
2339    unique("abcabc") some permutation of ["a", "b", "c"], and
2340    unique(([1, 2], [2, 3], [1, 2])) some permutation of
2341    [[2, 3], [1, 2]].
2342
2343    For best speed, all sequence elements should be hashable.  Then
2344    unique() will usually work in linear time.
2345
2346    If not possible, the sequence elements should enjoy a total
2347    ordering, and if list(s).sort() doesn't raise TypeError it's
2348    assumed that they do enjoy a total ordering.  Then unique() will
2349    usually work in O(N*log2(N)) time.
2350
2351    If that's not possible either, the sequence elements must support
2352    equality-testing.  Then unique() will usually work in quadratic
2353    time.
2354    """
2355
2356    n = len(s)
2357    if n == 0:
2358        return []
2359
2360    # Try using a dict first, as that's the fastest and will usually
2361    # work.  If it doesn't work, it will usually fail quickly, so it
2362    # usually doesn't cost much to *try* it.  It requires that all the
2363    # sequence elements be hashable, and support equality comparison.
2364    u = {}
2365    try:
2366        for x in s:
2367            u[x] = 1
2368    except TypeError:
2369        del u  # move on to the next method
2370    else:
2371        return u.keys()
2372
2373    # We can't hash all the elements.  Second fastest is to sort,
2374    # which brings the equal elements together; then duplicates are
2375    # easy to weed out in a single pass.
2376    # NOTE:  Python's list.sort() was designed to be efficient in the
2377    # presence of many duplicate elements.  This isn't true of all
2378    # sort functions in all languages or libraries, so this approach
2379    # is more effective in Python than it may be elsewhere.
2380    try:
2381        t = list(s)
2382        t.sort()
2383    except TypeError:
2384        del t  # move on to the next method
2385    else:
2386        assert n > 0
2387        last = t[0]
2388        lasti = i = 1
2389        while i < n:
2390            if t[i] != last:
2391                t[lasti] = last = t[i]
2392                lasti += 1
2393            i += 1
2394        return t[:lasti]
2395
2396    # Brute force is all that's left.
2397    u = []
2398    for x in s:
2399        if x not in u:
2400            u.append(x)
2401    return u
2402
2403
2404# FIXME (DSG-DSG)
2405# To do-
2406# Create a clear interface. eg
2407# have the interface methods more at the top of this file and add comments
2408# for the interface functions/methods, use function_name (not functionName),
2409
2410#Currently
2411#function_name methods assume absolute values.  Geo-refs can be passed in.
2412#
2413
2414# instead of functionName
2415if __name__ == "__main__":
2416    pass
Note: See TracBrowser for help on using the repository browser.