Context Navigation

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

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

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

Download in other formats: