Context Navigation

source: inundation/pmesh/mesh.py @ 2299

Last change on this file since 2299 was 2297, checked in by duncan, 19 years ago
new pmesh method, add_hole_from_polygon
File size: 125.0 KB

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

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

Download in other formats: