Context Navigation

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

Last change on this file since 4382 was 4382, checked in by duncan, 18 years ago
add option of sellecting a hole in urs_ungridded2sww
File size: 128.2 KB

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

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

Download in other formats: