Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: