Context Navigation

source: inundation/ga/storm_surge/pmesh/mesh.py @ 1599

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

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

Download in other formats: