Context Navigation

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

Last change on this file since 1039 was 1039, checked in by duncan, 20 years ago
new function for eastings and northings
File size: 95.6 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
30
31	SET_COLOUR='red'
32	DEFAULT_ZONE = 56
33
34	def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
35	"""
36	"""
37
38	det = (y2-y0)(x1-x0) - (y1-y0)(x2-x0)
39	a = (y2-y0)(q1-q0) - (y1-y0)(q2-q0)
40	a /= det
41
42	b = (x1-x0)(q2-q0) - (x2-x0)(q1-q0)
43	b /= det
44
45	return a, b
46
47	##############################################
48	#Initialise module
49
50	import compile
51	if compile.can_use_C_extension('util_ext.c'):
52	from util_ext import gradient, point_on_line
53	else:
54	gradient = gradient_python
55
56	# 1st and third values must be the same
57	# FIXME: maybe make this a switch that the user can change? - DSG
58	initialconversions = ['', 'exterior','']
59
60	#from os import sep
61	#sys.path.append('..'+sep+'pmesh')
62	#print "sys.path",sys.path
63
64	class MeshObject:
65	"""
66	An abstract superclass for the basic mesh objects, eg vertex, segment,
67	triangle.
68	"""
69	def __init__(self):
70	pass
71
72	class Point(MeshObject):
73	"""
74	Define a point in a 2D space.
75	"""
76	def __init__(self,X,Y):
77	__slots__ = ['x','y']
78	self.x=X
79	self.y=Y
80
81	def DistanceToPoint(self, OtherPoint):
82	"""
83	Returns the distance from this point to another
84	"""
85	SumOfSquares = ((self.x - OtherPoint.x)2) + ((self.y - OtherPoint.y)2)
86	return math.sqrt(SumOfSquares)
87
88	def IsInsideCircle(self, Center, Radius):
89	"""
90	Return 1 if this point is inside the circle,
91	0 otherwise
92	"""
93
94	if (self.DistanceToPoint(Center)<Radius):
95	return 1
96	else:
97	return 0
98
99	def __repr__(self):
100	return "(%f,%f)" % (self.x,self.y)
101
102	def cmp_xy(self, point):
103	if self.x < point.x:
104	return -1
105	elif self.x > point.x:
106	return 1
107	else:
108	if self.y < point.y:
109	return -1
110	elif self.y > point.y:
111	return 1
112	else:
113	return 0
114
115	def same_x_y(self, point):
116	if self.x == point.x and self.y == point.y:
117	return True
118	else:
119	return False
120
121
122
123	class Vertex(Point):
124	"""
125	A point on the mesh.
126	Object attributes based on the Triangle program
127	"""
128	def __init__(self,X,Y, attributes = None):
129	__slots__ = ['x','y','attributes']
130	assert (type(X) == types.FloatType)
131	self.x=X
132	self.y=Y
133	self.attributes=[]
134
135	if attributes is None:
136	self.attributes=[]
137	else:
138	self.attributes=attributes
139
140
141	def setAttributes(self,attributes):
142	"""
143	attributes is a list.
144	"""
145	self.attributes = attributes
146
147	VERTEXSQUARESIDELENGTH = 6
148	def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, yoffset =0, ):
149	x = scale*(self.x + xoffset)
150	y = -1scale(self.y + yoffset) # - since for a canvas - is up
151	#print "draw x:", x
152	#print "draw y:", y
153	cornerOffset= self.VERTEXSQUARESIDELENGTH/2
154
155	# A hack to see the vert tags
156	#canvas.create_text(x+ 2*cornerOffset,
157	# y+ 2*cornerOffset,
158	# text=tags)
159
160	return canvas.create_rectangle(x-cornerOffset,
161	y-cornerOffset,
162	x+cornerOffset,
163	y+cornerOffset,
164	tags = tags,
165	outline=colour,
166	fill = 'white')
167
168	#return tags
169
170	def __repr__(self):
171	return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
172
173	class Hole(Point):
174	"""
175	A region of the mesh were no triangles are generated.
176	Defined by a point in the hole enclosed by segments.
177	"""
178	HOLECORNERLENGTH = 6
179	def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ):
180	x = scale*(self.x + xoffset)
181	y = -1scale(self.y + yoffset) # - since for a canvas - is up
182	#print "draw x:", x
183	#print "draw y:", y
184	cornerOffset= self.HOLECORNERLENGTH/2
185	return canvas.create_oval(x-cornerOffset,
186	y-cornerOffset,
187	x+cornerOffset,
188	y+cornerOffset,
189	tags = tags,
190	outline=colour,
191	fill = 'white')
192
193	class Region(Point):
194	"""
195	A region of the mesh, defined by a point in the region
196	enclosed by segments. Used to tag areas.
197	"""
198	CROSSLENGTH = 6
199	TAG = 0
200	MAXAREA = 1
201
202	def __init__(self,X,Y, tag = None, maxArea = None):
203	"""Precondition: tag is a string and maxArea is a real
204	"""
205	# This didn't work.
206	#super(Region,self)._init_(self,X,Y)
207	self.x=X
208	self.y=Y
209	self.attributes =[] # index 0 is the tag string
210	#optoinal index 1 is the max triangle area
211	#NOTE the size of this attribute is assumed
212	# to be 1 or 2 in regionstrings2int
213	if tag is None:
214	self.attributes.append("")
215	else:
216	self.attributes.append(tag) #this is a string
217
218	if maxArea is not None:
219	self.setMaxArea(maxArea) # maxArea is a number
220
221	def getTag(self,):
222	return self.attributes[self.TAG]
223
224	def setTag(self,tag):
225	self.attributes[self.TAG] = tag
226
227	def getMaxArea(self):
228	""" Returns the Max Area of a Triangle or
229	None, if the Max Area has not been set.
230	"""
231	if self.isMaxArea():
232	return self.attributes[1]
233	else:
234	return None
235
236	def setMaxArea(self,MaxArea):
237	if self.isMaxArea():
238	self.attributes[self.MAXAREA] = float(MaxArea)
239	else:
240	self.attributes.append( float(MaxArea) )
241
242	def deleteMaxArea(self):
243	if self.isMaxArea():
244	self.attributes.pop(self.MAXAREA)
245
246	def isMaxArea(self):
247	return len(self.attributes)> 1
248
249	def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "black"):
250	"""
251	Draw a black cross, returning the objectID
252	"""
253	x = scale*(self.x + xoffset)
254	y = -1scale(self.y + yoffset)
255	cornerOffset= self.CROSSLENGTH/2
256	return canvas.create_polygon(x,
257	y-cornerOffset,
258	x,
259	y,
260	x+cornerOffset,
261	y,
262	x,
263	y,
264	x,
265	y+cornerOffset,
266	x,
267	y,
268	x-cornerOffset,
269	y,
270	x,
271	y,
272	tags = tags,
273	outline = colour,fill = '')
274
275	def __repr__(self):
276	if self.isMaxArea():
277	area = self.getMaxArea()
278	return "(%f,%f,%s,%f)" % (self.x,self.y,
279	self.getTag(), self.getMaxArea())
280	else:
281	return "(%f,%f,%s)" % (self.x,self.y,
282	self.getTag())
283
284	class Triangle(MeshObject):
285	"""
286	A triangle element, defined by 3 vertices.
287	Attributes based on the Triangle program.
288	"""
289
290	def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ):
291	"""
292	Vertices, the initial arguments, are listed in counterclockwise order.
293	"""
294	self.vertices= [vertex1,vertex2, vertex3 ]
295
296	if attribute is None:
297	self.attribute =""
298	else:
299	self.attribute = attribute #this is a string
300
301	if neighbors is None:
302	self.neighbors=[]
303	else:
304	self.neighbors=neighbors
305
306	def replace(self,new_triangle):
307	self = new_triangle
308
309	def longestSideID(self):
310	ax = self.vertices[0].x
311	ay = self.vertices[0].y
312
313	bx = self.vertices[1].x
314	by = self.vertices[1].y
315
316	cx = self.vertices[2].x
317	cy = self.vertices[2].y
318
319	lenA = ((cx-bx)2+(cy-by)2)**0.5
320	lenB = ((ax-cx)2+(ay-cy)2)**0.5
321	lenC = ((bx-ax)2+(by-ay)2)**0.5
322
323	len = [lenA,lenB,lenC]
324	return len.index(max(len))
325
326	def rotate(self,offset):
327	"""
328	offset must be 0,1 or 2
329	"""
330
331	if offset == 0:
332	pass
333	else:
334	if offset == 1:
335	self.vertices = [self.vertices[1],self.vertices[2],self.vertices[0]]
336	self.neighbors = [self.neighbors[1],self.neighbors[2],self.neighbors[0]]
337	if offset == 2:
338	self.vertices = [self.vertices[2],self.vertices[0],self.vertices[1]]
339	self.neighbors = [self.neighbors[2],self.neighbors[0],self.neighbors[1]]
340
341	def rotate_longest_side(self):
342	self.rotate(self.longestSideID())
343
344	def getVertices(self):
345	return self.vertices
346
347	def calcArea(self):
348	ax = self.vertices[0].x
349	ay = self.vertices[0].y
350
351	bx = self.vertices[1].x
352	by = self.vertices[1].y
353
354	cx = self.vertices[2].x
355	cy = self.vertices[2].y
356
357	return abs((bxay-axby)+(cxby-bxcy)+(axcy-cxay))/2
358
359	def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
360	"""
361	neighbor1 is the triangle opposite vertex1 and so on.
362	Null represents no neighbor
363	"""
364	self.neighbors = [neighbor1, neighbor2, neighbor3]
365
366	def setAttribute(self,attribute):
367	"""
368	neighbor1 is the triangle opposite vertex1 and so on.
369	Null represents no neighbor
370	"""
371	self.attribute = attribute
372
373	def __repr__(self):
374	return "[%s,%s]" % (self.vertices,self.attribute)
375
376
377	def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"):
378	"""
379	Draw a triangle, returning the objectID
380	"""
381	return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
382	scale-1(self.vertices[1].y + yoffset),
383	scale*(self.vertices[0].x + xoffset),
384	scale-1(self.vertices[0].y + yoffset),
385	scale*(self.vertices[2].x + xoffset),
386	scale-1(self.vertices[2].y + yoffset),
387	tags = tags,
388	outline = colour,fill = '')
389
390	class Segment(MeshObject):
391	"""
392	Segments are edges whose presence in the triangulation is enforced.
393
394	"""
395	def __init__(self, vertex1, vertex2, tag = None ):
396	"""
397	Each segment is specified by listing the vertices of its endpoints
398	The vertices are Vertex objects
399	"""
400	assert(vertex1 != vertex2)
401	self.vertices = [vertex1,vertex2 ]
402
403	if tag is None:
404	self.tag = self.__class__.default
405	else:
406	self.tag = tag #this is a string
407
408	def __repr__(self):
409	return "[%s,%s]" % (self.vertices,self.tag)
410
411
412	def draw(self, canvas, tags,scale=1 , xoffset=0 , yoffset=0,colour='blue' ):
413	x=[]
414	y=[]
415	for end in self.vertices:
416	#end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
417	x.append(scale*(end.x + xoffset))
418	y.append(-1scale(end.y + yoffset)) # - since for a canvas - is up
419
420	return canvas.create_line(x[0], y[0], x[1], y[1],
421	tags = tags,fill=colour)
422	def set_tag(self,tag):
423	self.tag = tag
424
425	# Class methods
426	def set_default_tag(cls, default):
427	cls.default = default
428
429	def get_default_tag(cls):
430	return cls.default
431
432	set_default_tag = classmethod(set_default_tag)
433	get_default_tag = classmethod(get_default_tag)
434
435	Segment.set_default_tag("")
436
437	class Mesh:
438	"""
439	Representation of a 2D triangular mesh.
440	User attributes describe the mesh region/segments/vertices/attributes
441
442	mesh attributes describe the mesh that is produced eg triangles and vertices.
443
444	The Mesh holds user information to define the
445	"""
446
447	def __repr__(self):
448	return """
449	mesh Triangles: %s
450	mesh Attribute Titles: %s
451	mesh Segments: %s
452	mesh Vertices: %s
453	user Segments: %s
454	user Vertices: %s
455	holes: %s
456	regions: %s""" % (self.meshTriangles,
457	self.attributeTitles,
458	self.meshSegments,
459	self.meshVertices,
460	self.getUserSegments(),
461	self.userVertices,
462	self.holes,
463	self.regions)
464
465	def __init__(self,
466	userSegments=None,
467	userVertices=None,
468	holes=None,
469	regions=None,
470	geo_reference=None):
471	self.meshTriangles=[]
472	self.attributeTitles=[]
473	self.meshSegments=[]
474	self.meshVertices=[]
475
476	#Peters stuff
477	# FIXME (DSG) Sets of what?
478	self.setID={}
479	self.setID['None']=0
480	self.sets=[[]]
481
482	if userSegments is None:
483	self.userSegments=[]
484	else:
485	self.userSegments=userSegments
486	self.alphaUserSegments=[]
487
488	if userVertices is None:
489	self.userVertices=[]
490	else:
491	self.userVertices=userVertices
492
493	if holes is None:
494	self.holes=[]
495	else:
496	self.holes=holes
497
498	if regions is None:
499	self.regions=[]
500	else:
501	self.regions=regions
502
503	if geo_reference is None:
504	self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
505	else:
506	self.geo_reference = geo_reference
507
508	def __cmp__(self,other):
509
510	# A dic for the initial m
511	dic = self.Mesh2triangList()
512	dic_mesh = self.Mesh2MeshList()
513	for element in dic_mesh.keys():
514	dic[element] = dic_mesh[element]
515
516	# A dic for the exported/imported m
517	dic_other = other.Mesh2triangList()
518	dic_mesh = other.Mesh2MeshList()
519	for element in dic_mesh.keys():
520	dic_other[element] = dic_mesh[element]
521
522	#print "dsg************************8"
523	#print "dic ",dic
524	#print "*******8"
525	#print "mesh",dic_other
526	#print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
527	#print "dsg************************8"
528
529	return (dic.__cmp__(dic_other))
530
531	def generateMesh(self, mode = None, maxArea = None, isRegionalMaxAreas = True):
532	"""
533	Based on the current user vaules, holes and regions
534	generate a new mesh
535	mode is a string that sets conditions on the mesh generations
536	see triangle_instructions.txt for a definition of the commands
537	PreCondition: maxArea is a double
538	"""
539	if mode == None:
540	self.mode = ""
541	else:
542	self.mode = mode
543
544	if not re.match('p',self.mode):
545	self.mode += 'p' #p - read a planar straight line graph.
546	#there must be segments to use this switch
547	# TODO throw an aception if there aren't seg's
548	# it's more comlex than this. eg holes
549	if not re.match('z',self.mode):
550	self.mode += 'z' # z - Number all items starting from zero (rather than one)
551	if not re.match('n',self.mode):
552	self.mode += 'n' # n - output a list of neighboring triangles
553	if not re.match('A',self.mode):
554	self.mode += 'A' # A - output region attribute list for triangles
555	if not re.match('V',self.mode) and not re.match('Q',self.mode):
556	self.mode += 'V' # V - output info about what Triangle is doing
557
558	if maxArea != None:
559	self.mode += 'a' + str(maxArea)
560
561	if isRegionalMaxAreas:
562	self.mode += 'a'
563
564	meshDict = self.Mesh2triangList()
565	#print "!@!@ This is going to triangle !@!@"
566	#print meshDict
567	#print "!@!@ This is going to triangle !@!@"
568
569	#print "meshDict['segmenttaglist']", meshDict['segmenttaglist']
570	[meshDict['segmenttaglist'],
571	segconverter] = segment_strings2ints(meshDict['segmenttaglist'],
572	initialconversions)
573	#print "regionlist",meshDict['regionlist']
574	[meshDict['regionlist'],
575	regionconverter] = region_strings2ints(meshDict['regionlist'])
576	#print "regionlist",meshDict['regionlist']
577	#print "meshDict['segmenttaglist']", meshDict['segmenttaglist'
578	generatedMesh = triang.genMesh(
579	meshDict['pointlist'],
580	meshDict['segmentlist'],
581	meshDict['holelist'],
582	meshDict['regionlist'],
583	meshDict['pointattributelist'],
584	meshDict['segmenttaglist'],
585	[], # since the trianglelist isn't used
586	self.mode)
587	#print "generated",generatedMesh['generatedsegmenttaglist']
588	generatedMesh['generatedsegmentmarkerlist'] = \
589	segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
590	segconverter)
591	#print "processed gen",generatedMesh['generatedsegmentmarkerlist']
592	generatedMesh['generatedtriangleattributelist'] = \
593	region_ints2strings(generatedMesh['generatedtriangleattributelist'],
594	regionconverter)
595
596
597	if len(generatedMesh['generatedpointattributelist'][0])==0:
598	self.attributeTitles = []
599	generatedMesh['generatedpointattributetitlelist']=self.attributeTitles
600
601	self.setTriangulation(generatedMesh)
602
603	def addUserPoint(self, pointType, x,y):
604	if pointType == Vertex:
605	point = self.addUserVertex(x,y)
606	if pointType == Hole:
607	point = self.addHole(x,y)
608	if pointType == Region:
609	point = self.addRegion(x,y)
610	return point
611
612	def addUserVertex(self, x,y):
613	v=Vertex(x, y)
614	self.userVertices.append(v)
615	return v
616
617	def addHole(self, x,y):
618	h=Hole(x, y)
619	self.holes.append(h)
620	return h
621
622	def addRegion(self, x,y):
623	h=Region(x, y)
624	self.regions.append(h)
625	return h
626
627	def addRegionEN(self, x,y):
628	h=Region(x-self.geo_reference.xllcorner,
629	y-self.geo_reference.yllcorner)
630	self.regions.append(h)
631	return h
632
633	def getUserVertices(self):
634	return self.userVertices
635
636	def getUserSegments(self):
637	allSegments = self.userSegments + self.alphaUserSegments
638	#print "self.userSegments",self.userSegments
639	#print "self.alphaUserSegments",self.alphaUserSegments
640	#print "allSegments",allSegments
641	return allSegments
642
643	def deleteUserSegments(self,seg):
644	if self.userSegments.count(seg) == 0:
645	self.alphaUserSegments.remove(seg)
646	pass
647	else:
648	self.userSegments.remove(seg)
649
650	def clearUserSegments(self):
651	self.userSegments = []
652	self.alphaUserSegments = []
653
654	def getTriangulation(self):
655	return self.meshTriangles
656
657	def getMeshVertices(self):
658	return self.meshVertices
659
660	def getMeshSegments(self):
661	return self.meshSegments
662
663	def getHoles(self):
664	return self.holes
665
666	def getRegions(self):
667	return self.regions
668
669	def isTriangulation(self):
670	if self.meshVertices == []:
671	return False
672	else:
673	return True
674
675	def addUserSegment(self, v1,v2):
676	"""
677	PRECON: A segment between the two vertices is not already present.
678	Check by calling isUserSegmentNew before calling this function.
679
680	"""
681	s=Segment( v1,v2)
682	self.userSegments.append(s)
683	return s
684
685	def clearTriangulation(self):
686
687	#Clear the current generated mesh values
688	self.meshTriangles=[]
689	self.meshSegments=[]
690	self.meshVertices=[]
691
692	def removeDuplicatedUserVertices(self):
693	"""Pre-condition: There are no user segments
694	This function will keep the first duplicate
695	"""
696	assert self.getUserSegments() == []
697	self.userVertices, counter = self.removeDuplicatedVertices(self.userVertices)
698	return counter
699
700	def removeDuplicatedVertices(self, Vertices):
701	"""
702	This function will keep the first duplicate, remove all others
703	Precondition: Each vertex has a dupindex, which is the list
704	index.
705	"""
706	remove = []
707	index = 0
708	for v in Vertices:
709	v.dupindex = index
710	index += 1
711	t = list(Vertices)
712	t.sort(Point.cmp_xy)
713
714	length = len(t)
715	behind = 0
716	ahead = 1
717	counter = 0
718	while ahead < length:
719	b = t[behind]
720	ah = t[ahead]
721	if (b.y == ah.y and b.x == ah.x):
722	remove.append(ah.dupindex)
723	behind += 1
724	ahead += 1
725
726	# remove the duplicate vertices
727	remove.sort()
728	remove.reverse()
729	for i in remove:
730	Vertices.pop(i)
731	pass
732
733	#Remove the attribute that this function added
734	for v in Vertices:
735	del v.dupindex
736	return Vertices,counter
737
738	def thinoutVertices(self, delta):
739	"""Pre-condition: There are no user segments
740	This function will keep the first duplicate
741	"""
742	assert self.getUserSegments() == []
743	#t = self.userVertices
744	#self.userVertices =[]
745	boxedVertices = {}
746	thinnedUserVertices =[]
747	delta = round(delta,1)
748
749	for v in self.userVertices :
750	# tag is the center of the boxes
751	tag = (round(v.x/delta,0)delta,round(v.y/delta,0)delta)
752	#this creates a dict of lists of faces, indexed by tag
753	boxedVertices.setdefault(tag,[]).append(v)
754
755	for [tag,verts] in boxedVertices.items():
756	min = delta
757	bestVert = None
758	tagVert = Vertex(tag[0],tag[1])
759	for v in verts:
760	dist = v.DistanceToPoint(tagVert)
761	if (dist<min):
762	min = dist
763	bestVert = v
764	thinnedUserVertices.append(bestVert)
765	self.userVertices =thinnedUserVertices
766
767
768	def isUserSegmentNew(self, v1,v2):
769	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) ]
770
771	return len(identicalSegs) == 0
772
773	def representedAlphaUserSegment(self, v1,v2):
774	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) ]
775
776	if identicalSegs == []:
777	return None
778	else:
779	# Only return the first one.
780	return identicalSegs[0]
781
782	def representedUserSegment(self, v1,v2):
783	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) ]
784
785	if identicalSegs == []:
786	return None
787	else:
788	# Only return the first one.
789	return identicalSegs[0]
790
791	def deleteSegsOfVertex(self, delVertex):
792	"""
793	Delete this vertex and any segments that connect to it.
794	"""
795	#Find segments that connect to delVertex
796	deletedSegments = []
797	for seg in self.getUserSegments():
798	if (delVertex in seg.vertices):
799	deletedSegments.append(seg)
800	# Delete segments that connect to delVertex
801	for seg in deletedSegments:
802	self.deleteUserSegments(seg)
803	return deletedSegments
804
805
806	def deleteMeshObject(self, MeshObject):
807	"""
808	Returns a list of all objects that were removed
809	"""
810	deletedObs = []
811	if isinstance(MeshObject, Vertex ):
812	deletedObs = self.deleteSegsOfVertex(MeshObject)
813	deletedObs.append(MeshObject)
814	self.userVertices.remove(MeshObject)
815	elif isinstance(MeshObject, Segment):
816	deletedObs.append(MeshObject)
817	self.deleteUserSegments(MeshObject)
818	elif isinstance(MeshObject, Hole):
819	deletedObs.append(MeshObject)
820	self.holes.remove(MeshObject)
821	elif isinstance(MeshObject, Region):
822	deletedObs.append(MeshObject)
823	self.regions.remove(MeshObject)
824	return deletedObs
825
826	def Mesh2triangList(self, userVertices=None,
827	userSegments=None,
828	holes=None,
829	regions=None):
830	"""
831	Convert the Mesh to a dictionary of the lists needed for the triang modul;
832	points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
833	pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
834	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
835	hole list: [(x1,y1),...](Tuples of doubles, one inside each hole region)
836	regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
837
838	Note, this adds an index attribute to the user Vertex objects.
839
840	Used to produce output to triangle
841	"""
842	if userVertices is None:
843	userVertices = self.getUserVertices()
844	if userSegments is None:
845	userSegments = self.getUserSegments()
846	if holes is None:
847	holes = self.getHoles()
848	if regions is None:
849	regions = self.getRegions()
850
851	meshDict = {}
852
853	pointlist=[]
854	pointattributelist=[]
855	index = 0
856	for vertex in userVertices:
857	vertex.index = index
858	pointlist.append((vertex.x,vertex.y))
859	pointattributelist.append((vertex.attributes))
860
861	index += 1
862	meshDict['pointlist'] = pointlist
863	meshDict['pointattributelist'] = pointattributelist
864
865	segmentlist=[]
866	segmenttaglist=[]
867	for seg in userSegments:
868	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
869	segmenttaglist.append(seg.tag)
870	meshDict['segmentlist'] =segmentlist
871	meshDict['segmenttaglist'] =segmenttaglist
872
873	holelist=[]
874	for hole in holes:
875	holelist.append((hole.x,hole.y))
876	meshDict['holelist'] = holelist
877
878	regionlist=[]
879	for region in regions:
880	if (region.getMaxArea() != None):
881	regionlist.append((region.x,region.y,region.getTag(),
882	region.getMaxArea()))
883	else:
884	regionlist.append((region.x,region.y,region.getTag()))
885	meshDict['regionlist'] = regionlist
886	#print "(("
887	#print meshDict
888	#print meshDict['regionlist']
889	#print "(("
890	return meshDict
891
892	def Mesh2MeshList(self):
893	"""
894	Convert the Mesh to a dictionary of lists describing the triangulation variables;
895	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
896	generated point attribute list: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
897	generated point attribute title list:[A1Title, A2Title ...] (list of strings)
898	generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
899	generated segment tag list: [tag,tag,...] list of strings
900
901	generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
902
903	generated triangle attribute list: [s1,s2,...] list of strings
904
905	generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....] tuple of triangles
906
907	Used to produce .tsh file
908	"""
909
910	meshDict = {}
911	pointlist=[]
912	pointattributelist=[]
913
914
915	self.maxVertexIndex=0
916	for vertex in self.meshVertices:
917	vertex.index = self.maxVertexIndex
918	pointlist.append((vertex.x,vertex.y))
919	pointattributelist.append((vertex.attributes))
920	self.maxVertexIndex += 1
921
922	meshDict['generatedpointlist'] = pointlist
923	meshDict['generatedpointattributelist'] = pointattributelist
924	meshDict['generatedpointattributetitlelist'] = self.attributeTitles
925	#segments
926	segmentlist=[]
927	segmenttaglist=[]
928	for seg in self.meshSegments:
929	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
930	segmenttaglist.append(seg.tag)
931	meshDict['generatedsegmentlist'] =segmentlist
932	meshDict['generatedsegmenttaglist'] =segmenttaglist
933
934	# Make sure that the indexation is correct
935	index = 0
936	for tri in self.meshTriangles:
937	tri.index = index
938	index += 1
939
940	trianglelist = []
941	triangleattributelist = []
942	triangleneighborlist = []
943	for tri in self.meshTriangles:
944	trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,tri.vertices[2].index))
945	triangleattributelist.append([tri.attribute])
946	neighborlist = [-1,-1,-1]
947	for neighbor,index in map(None,tri.neighbors,
948	range(len(tri.neighbors))):
949	if neighbor:
950	neighborlist[index] = neighbor.index
951	triangleneighborlist.append(neighborlist)
952
953	meshDict['generatedtrianglelist'] = trianglelist
954	meshDict['generatedtriangleattributelist'] = triangleattributelist
955	meshDict['generatedtriangleneighborlist'] = triangleneighborlist
956
957	#print "mesh.Mesh2MeshList))"
958	#print meshDict
959	#print "mesh.Mesh2MeshList))"
960
961	return meshDict
962
963
964	def Mesh2MeshDic(self):
965	"""
966	Convert the user and generated info of a mesh to a dictionary
967	structure
968	"""
969	dic = self.Mesh2triangList()
970	dic_mesh = self.Mesh2MeshList()
971	for element in dic_mesh.keys():
972	dic[element] = dic_mesh[element]
973	return dic
974
975	def setTriangulation(self, genDict):
976	"""
977	Set the mesh attributes given a dictionary of the lists
978	returned from the triang module
979	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
980	generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
981	generated point attribute title list:[A1Title, A2Title ...] (list of strings)
982	generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
983	generated segment marker list: [S1Tag, S2Tag, ...] (list of ints)
984	triangle list: [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
985	triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
986	triangle attribute list: [(T1att), (T2att), ...] (list of a list of strings)
987	"""
988	#Clear the current generated mesh values
989	self.meshTriangles=[]
990	self.attributeTitles=[]
991	self.meshSegments=[]
992	self.meshVertices=[]
993
994	#print "mesh.setTriangulation@#@#@#"
995	#print genDict
996	#print "@#@#@#"
997
998	self.maxVertexIndex = 0
999	for point in genDict['generatedpointlist']:
1000	v=Vertex(point[0], point[1])
1001	v.index = self.maxVertexIndex
1002	self.maxVertexIndex +=1
1003	self.meshVertices.append(v)
1004
1005	self.attributeTitles = genDict['generatedpointattributetitlelist']
1006
1007	index = 0
1008	for seg,marker in map(None,genDict['generatedsegmentlist'],genDict['generatedsegmentmarkerlist']):
1009	segObject = Segment( self.meshVertices[seg[0]],
1010	self.meshVertices[seg[1]], tag = marker )
1011	segObject.index = index
1012	index +=1
1013	self.meshSegments.append(segObject)
1014
1015	index = 0
1016	for triangle in genDict['generatedtrianglelist']:
1017	tObject =Triangle( self.meshVertices[triangle[0]],
1018	self.meshVertices[triangle[1]],
1019	self.meshVertices[triangle[2]] )
1020	tObject.index = index
1021	index +=1
1022	self.meshTriangles.append(tObject)
1023
1024	index = 0
1025	for att in genDict['generatedtriangleattributelist']:
1026	if att == []:
1027	self.meshTriangles[index].setAttribute("")
1028	else:
1029	# Note, is the first attribute always the region att?
1030	# haven't confirmed this
1031	#Peter - I think so (from manuel)
1032	#...the first such value is assumed to be a regional attribute...
1033	self.meshTriangles[index].setAttribute(att[0])
1034	index += 1
1035
1036	index = 0
1037	for att in genDict['generatedpointattributelist']:
1038	if att == None:
1039	self.meshVertices[index].setAttributes([])
1040	else:
1041	self.meshVertices[index].setAttributes(att)
1042	index += 1
1043
1044	index = 0
1045	for triangle in genDict['generatedtriangleneighborlist']:
1046	# Build a list of triangle object neighbors
1047	ObjectNeighbor = []
1048	for neighbor in triangle:
1049	if ( neighbor != -1):
1050	ObjectNeighbor.append(self.meshTriangles[neighbor])
1051	else:
1052	ObjectNeighbor.append(None)
1053	self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],ObjectNeighbor[1],ObjectNeighbor[2])
1054	index += 1
1055
1056
1057	def setMesh(self, genDict):
1058	"""
1059	Set the user Mesh attributes given a dictionary of the lists
1060	point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1061	point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1062	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1063	segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1064	region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1065	region attribute list: ["","reservoir",""] list of strings
1066	region max area list:[real, None, Real,...] list of None and reals
1067
1068	mesh is an instance of a mesh object
1069	"""
1070	#Clear the current user mesh values
1071	self.clearUserSegments()
1072	self.userVertices=[]
1073	self.Holes=[]
1074	self.Regions=[]
1075
1076	#print "mesh.setMesh@#@#@#"
1077	#print genDict
1078	#print "@#@#@#"
1079
1080	#index = 0
1081	for point in genDict['pointlist']:
1082	v=Vertex(point[0], point[1])
1083	#v.index = index
1084	#index +=1
1085	self.userVertices.append(v)
1086
1087	#index = 0
1088	for seg,tag in map(None,genDict['segmentlist'],genDict['segmenttaglist']):
1089	segObject = Segment( self.userVertices[seg[0]],
1090	self.userVertices[seg[1]], tag = tag )
1091	#segObject.index = index
1092	#index +=1
1093	self.userSegments.append(segObject)
1094
1095	# Remove the loading of attribute info.
1096	# Have attribute info added using least_squares in pyvolution
1097	# index = 0
1098	# for att in genDict['pointattributelist']:
1099	# if att == None:
1100	# self.userVertices[index].setAttributes([])
1101	# else:
1102	# self.userVertices[index].setAttributes(att)
1103	# index += 1
1104
1105	#index = 0
1106	for point in genDict['holelist']:
1107	h=Hole(point[0], point[1])
1108	#h.index = index
1109	#index +=1
1110	self.holes.append(h)
1111
1112	#index = 0
1113	for reg,att,maxArea in map(None,
1114	genDict['regionlist'],
1115	genDict['regionattributelist'],
1116	genDict['regionmaxarealist']):
1117	Object = Region( reg[0],
1118	reg[1],
1119	tag = att,
1120	maxArea = maxArea)
1121	#Object.index = index
1122	#index +=1
1123	self.regions.append(Object)
1124
1125	def addVertsSegs(self, outlineDict):
1126	"""
1127	Add out-line (user Mesh) attributes given a dictionary of the lists
1128	point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1129	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1130	#segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1131
1132	Assume the values are in Eastings and Northings, with no reference
1133	point
1134	"""
1135
1136	localUserVertices = []
1137	#index = 0
1138	for point in outlineDict['pointlist']:
1139	v=Vertex(point[0]-self.geo_reference.xllcorner,
1140	point[1]-self.geo_reference.yllcorner)
1141	#v.index = index
1142	#index +=1
1143	self.userVertices.append(v)
1144	localUserVertices.append(v)
1145
1146	#index = 0
1147	for seg in map(None,outlineDict['segmentlist']):
1148	segObject = Segment( localUserVertices[seg[0]],
1149	localUserVertices[seg[1]] )
1150	#segObject.index = index
1151	#index +=1
1152	self.userSegments.append(segObject)
1153	#DSG!!!
1154
1155	def TestautoSegment(self):
1156	newsegs = []
1157	s1 = Segment(self.userVertices[0],
1158	self.userVertices[1])
1159	s2 = Segment(self.userVertices[0],
1160	self.userVertices[2])
1161	s3 = Segment(self.userVertices[2],
1162	self.userVertices[1])
1163	if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1164	newsegs.append(s1)
1165	if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1166	newsegs.append(s2)
1167	if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1168	newsegs.append(s3)
1169	#DSG!!!
1170	self.userSegments.extend(newsegs)
1171	return newsegs
1172
1173
1174	def savePickle(self, currentName):
1175	fd = open(currentName, 'w')
1176	pickle.dump(self,fd)
1177	fd.close()
1178
1179	def autoSegmentHull(self):
1180	"""
1181	initially work by running an executable
1182	Later compile the c code with a python wrapper.
1183
1184	Precon: There must be 3 or more vertices in the userVertices structure
1185	"""
1186	newsegs = []
1187	inputfile = 'hull_in.txt'
1188	outputfile = inputfile + '-alf'
1189	#write vertices to file
1190	fd = open(inputfile,'w')
1191	for v in self.userVertices:
1192	fd.write(str(v.x))
1193	fd.write(' ')
1194	fd.write(str(v.y))
1195	fd.write('\n')
1196	fd.close()
1197
1198	#run hull executable
1199	#warning need to compile hull for the current operating system
1200	command = 'hull.exe -A -i ' + inputfile
1201	os.system(command)
1202
1203	#read results into this object
1204	fd = open(outputfile)
1205	lines = fd.readlines()
1206	fd.close()
1207	#print "(((*("
1208	#print lines
1209	#print "(((*("
1210	lines.pop(0) #remove the first (title) line
1211	for line in lines:
1212	vertindexs = line.split()
1213	#print 'int(vertindexs[0])', int(vertindexs[0])
1214	#print 'int(vertindexs[1])', int(vertindexs[1])
1215	#print 'self.userVertices[int(vertindexs[0])]' ,self.userVertices[int(vertindexs[0])]
1216	#print 'self.userVertices[int(vertindexs[1])]' ,self.userVertices[int(vertindexs[1])]
1217	v1 = self.userVertices[int(vertindexs[0])]
1218	v2 = self.userVertices[int(vertindexs[1])]
1219
1220	if self.isUserSegmentNew(v1,v2):
1221	newseg = Segment(v1, v2)
1222	newsegs.append(newseg)
1223	#DSG!!!
1224	self.userSegments.extend(newsegs)
1225	return newsegs
1226	def autoSegmentFilter(self,):
1227	pass
1228
1229	def autoSegment(self, alpha = None,
1230	raw_boundary=True,
1231	remove_holes=False,
1232	smooth_indents=False,
1233	expand_pinch=False):
1234	"""
1235	Precon: There must be 3 or more vertices in the userVertices structure
1236	"""
1237
1238	points=[]
1239	for vertex in self.getUserVertices():
1240	points.append((vertex.x,vertex.y))
1241	self.shape = alpha_shape.alpha_shape.Alpha_Shape(points, alpha = alpha)
1242	self.shape.set_boundary_type(raw_boundary=raw_boundary,
1243	remove_holes=remove_holes,
1244	smooth_indents=smooth_indents,
1245	expand_pinch=expand_pinch)
1246	boundary_segs = self.shape.get_boundary()
1247
1248	segs2delete = self.alphaUserSegments
1249
1250	new_segs = []
1251	alpha_segs = []
1252	user_segs = []
1253	for seg in boundary_segs:
1254	v1 = self.userVertices[int(seg[0])]
1255	v2 = self.userVertices[int(seg[1])]
1256	alpha_seg = self.representedAlphaUserSegment(v1, v2)
1257	user_seg = self.representedUserSegment(v1, v2)
1258	#DSG!!!
1259	assert not(not (alpha_seg == None) and not (user_seg == None))
1260	if not alpha_seg == None:
1261	alpha_segs.append(alpha_seg)
1262	elif not user_seg == None:
1263	user_segs.append(user_seg)
1264	else:
1265	unique_seg = Segment(v1, v2)
1266	new_segs.append(unique_seg)
1267
1268	for seg in alpha_segs:
1269	try:
1270	segs2delete.remove(seg)
1271	except:
1272	pass
1273
1274	self.alphaUserSegments = []
1275	self.alphaUserSegments.extend(new_segs)
1276	self.alphaUserSegments.extend(alpha_segs)
1277
1278	optimum_alpha = self.shape.get_alpha()
1279	# need to draw newsegs
1280	return new_segs, segs2delete, optimum_alpha
1281
1282	def joinVertices(self):
1283	"""
1284	Return list of segments connecting all userVertices
1285	in the order they were given
1286
1287	Precon: There must be 3 or more vertices in the userVertices structure
1288	"""
1289
1290	newsegs = []
1291
1292	v1 = self.userVertices[0]
1293	for v2 in self.userVertices[1:]:
1294	if self.isUserSegmentNew(v1,v2):
1295	newseg = Segment(v1, v2)
1296	newsegs.append(newseg)
1297	v1 = v2
1298
1299	#Connect last point to the first
1300	v2 = self.userVertices[0]
1301	if self.isUserSegmentNew(v1,v2):
1302	newseg = Segment(v1, v2)
1303	newsegs.append(newseg)
1304
1305
1306	#Update list of user segments
1307	#DSG!!!
1308	self.userSegments.extend(newsegs)
1309	return newsegs
1310
1311	def normaliseMesh(self,scale, offset, height_scale):
1312	[xmin, ymin, xmax, ymax] = self.boxsize()
1313	[attmin0, attmax0] = self.maxMinVertAtt(0)
1314	#print "[attmin0, attmax0]" ,[attmin0, attmax0]
1315	[attmin1, attmax1] = self.maxMinVertAtt(1)
1316	#print [xmin, ymin, xmax, ymax]
1317	xrange = xmax - xmin
1318	yrange = ymax - ymin
1319	if xrange > yrange:
1320	min,max = xmin, xmax
1321	else:
1322	min,max = ymin, ymax
1323
1324	for obj in self.getUserVertices():
1325	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1326	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1327	if len(obj.attributes) > 0 and attmin0 != attmax0:
1328	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1329	(attmax0-attmin0)*height_scale
1330	if len(obj.attributes) > 1 and attmin1 != attmax1:
1331	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1332	(attmax1-attmin1)*height_scale
1333
1334	for obj in self.getMeshVertices():
1335	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1336	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1337	if len(obj.attributes) > 0 and attmin0 != attmax0:
1338	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1339	(attmax0-attmin0)*height_scale
1340	if len(obj.attributes) > 1 and attmin1 != attmax1:
1341	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1342	(attmax1-attmin1)*height_scale
1343
1344	for obj in self.getHoles():
1345	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1346	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1347	for obj in self.getRegions():
1348	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1349	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1350	[xmin, ymin, xmax, ymax] = self.boxsize()
1351	#print [xmin, ymin, xmax, ymax]
1352
1353	def boxsizeVerts(self):
1354	"""
1355	Returns a list of verts denoting a box or triangle that contains verts on the xmin, ymin, xmax and ymax axis.
1356	Structure: list of verts
1357	"""
1358	# FIXME dsg!!! large is a hack
1359	#You want the kinds package, part of Numeric:
1360	#In [2]: import kinds
1361
1362	#In [3]: kinds.default_float_kind.M
1363	#kinds.default_float_kind.MAX kinds.default_float_kind.MIN
1364	#kinds.default_float_kind.MAX_10_EXP kinds.default_float_kind.MIN_10_EXP
1365	#kinds.default_float_kind.MAX_EXP kinds.default_float_kind.MIN_EXP
1366
1367	#In [3]: kinds.default_float_kind.MIN
1368	#Out[3]: 2.2250738585072014e-308
1369
1370	large = 1e100
1371	xmin= large
1372	xmax=-large
1373	ymin= large
1374	ymax=-large
1375	for vertex in self.userVertices:
1376	if vertex.x < xmin:
1377	xmin = vertex.x
1378	xminVert = vertex
1379	if vertex.x > xmax:
1380	xmax = vertex.x
1381	xmaxVert = vertex
1382
1383	if vertex.y < ymin:
1384	ymin = vertex.y
1385	yminVert = vertex
1386	if vertex.y > ymax:
1387	ymax = vertex.y
1388	ymaxVert = vertex
1389	verts, count = self.removeDuplicatedVertices([xminVert,xmaxVert,yminVert,ymaxVert])
1390
1391	return verts
1392
1393	def boxsize(self):
1394	"""
1395	Returns a list denoting a box that contains the entire structure of vertices
1396	Structure: [xmin, ymin, xmax, ymax]
1397	"""
1398	# FIXME dsg!!! large is a hack
1399	#You want the kinds package, part of Numeric:
1400	#In [2]: import kinds
1401
1402	#In [3]: kinds.default_float_kind.M
1403	#kinds.default_float_kind.MAX kinds.default_float_kind.MIN
1404	#kinds.default_float_kind.MAX_10_EXP kinds.default_float_kind.MIN_10_EXP
1405	#kinds.default_float_kind.MAX_EXP kinds.default_float_kind.MIN_EXP
1406
1407	#In [3]: kinds.default_float_kind.MIN
1408	#Out[3]: 2.2250738585072014e-308
1409
1410	large = 1e100
1411	xmin= large
1412	xmax=-large
1413	ymin= large
1414	ymax=-large
1415	for vertex in self.userVertices:
1416	if vertex.x < xmin:
1417	xmin = vertex.x
1418	if vertex.x > xmax:
1419	xmax = vertex.x
1420
1421	if vertex.y < ymin:
1422	ymin = vertex.y
1423	if vertex.y > ymax:
1424	ymax = vertex.y
1425	return [xmin, ymin, xmax, ymax]
1426
1427	def maxMinVertAtt(self, iatt):
1428	"""
1429	Returns a list denoting a box that contains the entire structure of vertices
1430	Structure: [xmin, ymin, xmax, ymax]
1431	"""
1432	# FIXME dsg!!! large is a hack
1433	#You want the kinds package, part of Numeric:
1434	#In [2]: import kinds
1435
1436	#In [3]: kinds.default_float_kind.M
1437	#kinds.default_float_kind.MAX kinds.default_float_kind.MIN
1438	#kinds.default_float_kind.MAX_10_EXP kinds.default_float_kind.MIN_10_EXP
1439	#kinds.default_float_kind.MAX_EXP kinds.default_float_kind.MIN_EXP
1440
1441	#In [3]: kinds.default_float_kind.MIN
1442	#Out[3]: 2.2250738585072014e-308
1443
1444	large = 1e100
1445	min= large
1446	max=-large
1447	for vertex in self.userVertices:
1448	if len(vertex.attributes) > iatt:
1449	if vertex.attributes[iatt] < min:
1450	min = vertex.attributes[iatt]
1451	if vertex.attributes[iatt] > max:
1452	max = vertex.attributes[iatt]
1453	for vertex in self.meshVertices:
1454	if len(vertex.attributes) > iatt:
1455	if vertex.attributes[iatt] < min:
1456	min = vertex.attributes[iatt]
1457	if vertex.attributes[iatt] > max:
1458	max = vertex.attributes[iatt]
1459	return [min, max]
1460
1461	def scaleoffset(self, WIDTH, HEIGHT):
1462	"""
1463	Returns a list denoting the scale and offset terms that need to be
1464	applied when converting mesh co-ordinates onto grid co-ordinates
1465	Structure: [scale, xoffset, yoffset]
1466	"""
1467	OFFSET = 0.05*min([WIDTH, HEIGHT])
1468	[xmin, ymin, xmax, ymax] = self.boxsize()
1469	SCALE = min([0.9WIDTH, 0.9HEIGHT])/max([xmax-xmin, ymax-ymin])
1470
1471	if SCALE*xmin < OFFSET:
1472	xoffset = abs(SCALE*xmin) + OFFSET
1473	if SCALE*xmax > WIDTH - OFFSET:
1474	xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1475	if SCALE*ymin < OFFSET:
1476	b = abs(SCALE*ymin)+OFFSET
1477	if SCALE*ymax > HEIGHT-OFFSET:
1478	b = -(SCALE*ymax - HEIGHT + OFFSET)
1479	yoffset = HEIGHT - b
1480	return [SCALE, xoffset, yoffset]
1481
1482	def plotMeshTriangle(self,tag = 0,WIDTH = 400,HEIGHT = 400):
1483	"""
1484	Plots all node connections.
1485	tag = 0 (no node numbers), tag = 1 (node numbers)
1486	"""
1487
1488	try:
1489	from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
1490
1491	[SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
1492
1493	root = Tk()
1494	frame = Frame(root)
1495	frame.pack()
1496	button = Button(frame, text="OK", fg="red", command=frame.quit)
1497	button.pack(side=BOTTOM)
1498	canvas = Canvas(frame,bg="white", width=WIDTH, height=HEIGHT)
1499	canvas.pack()
1500	text = Label(frame, width=20, height=10, text='triangle mesh')
1501	text.pack()
1502
1503	#print self.meshTriangles
1504	for triangle in self.meshTriangles:
1505	triangle.draw(canvas,1,
1506	scale = SCALE,
1507	xoffset = xoffset,
1508	yoffset = yoffset )
1509
1510	root.mainloop()
1511
1512	except:
1513	print "Unexpected error:", sys.exc_info()[0]
1514	raise
1515
1516	#print """
1517	#node::plot Failed.
1518	#Most probably, the Tkinter module is not available.
1519	#"""
1520
1521	def plotUserSegments(self,tag = 0,WIDTH = 400,HEIGHT = 400):
1522	"""
1523	Plots all node connections.
1524	tag = 0 (no node numbers), tag = 1 (node numbers)
1525	"""
1526
1527	try:
1528	from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
1529
1530	[SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
1531
1532	root = Tk()
1533	frame = Frame(root)
1534	frame.pack()
1535	button = Button(frame, text="OK", fg="red", command=frame.quit)
1536	button.pack(side=BOTTOM)
1537	canvas = Canvas(frame, bg="white", width=WIDTH, height=HEIGHT)
1538	canvas.pack()
1539	text = Label(frame, width=20, height=10, text='user segments')
1540	text.pack()
1541
1542	for segment in self.getUserSegments():
1543	segment.draw(canvas,SCALE, xoffset, yoffset )
1544
1545	root.mainloop()
1546
1547	except:
1548	print "Unexpected error:", sys.exc_info()[0]
1549	raise
1550
1551	#print """
1552	#node::plot Failed.
1553	#Most probably, the Tkinter module is not available.
1554	#"""
1555	# FIXME let's not use this function,
1556	# use export _mesh_file instead
1557	def export_triangulation_file(self,ofile):
1558	"""
1559	export a file, ofile, with the format
1560
1561	First line: <# of vertices> <# of attributes>
1562	Following lines: <vertex #> <x> <y> [attributes]
1563	One line: <vertex attribute titles>
1564	One line: <# of triangles>
1565	Following lines: <triangle #> <vertex #> <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
1566	One line: <# of segments>
1567	Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
1568	"""
1569	gen_dict = self.Mesh2IODict()
1570	if (ofile[-4:] == ".tsh"):
1571	fd = open(ofile,'w')
1572	load_mesh.loadASCII.write_ASCII_triangulation(fd,gen_dict)
1573	self.writeASCIImesh(fd,
1574	self.userVertices,
1575	self.getUserSegments(),
1576	self.holes,
1577	self.regions)
1578	fd.close()
1579	elif (ofile[-4:] == ".msh"):
1580	#print "mesh gen_dict",gen_dict
1581	load_mesh.loadASCII.write_msh_file(ofile, gen_dict)
1582
1583	# self.writeASCIImesh(fd, does anything use this?
1584
1585	#FIXME this function has a bug..
1586	def exportASCIIsegmentoutlinefile(self,ofile):
1587	"""
1588	export a file, ofile, with no triangulation and only vertices connected to segments.
1589	"""
1590	fd = open(ofile,'w')
1591	meshDict = {}
1592
1593	meshDict['vertices'] = []
1594	meshDict['vertex_attributes'] = []
1595	meshDict['segments'] = []
1596	meshDict['segment_tags'] = []
1597	meshDict['triangles'] = []
1598	meshDict['triangle_tags'] = []
1599	meshDict['triangle_neighbors'] = []
1600
1601	load_mesh.loadASCII.write_ASCII_triangulation(fd,meshDict)
1602	self.writeASCIIsegmentoutline(fd,
1603	self.userVertices,
1604	self.getUserSegments(),
1605	self.holes,
1606	self.regions)
1607	fd.close()
1608
1609	def exportASCIIobj(self,ofile):
1610	"""
1611	export a file, ofile, with the format
1612	lines: v <x> <y> <first attribute>
1613	f <vertex #> <vertex #> <vertex #> (of the triangles)
1614	"""
1615	fd = open(ofile,'w')
1616	self.writeASCIIobj(fd)
1617	fd.close()
1618
1619
1620	def writeASCIIobj(self,fd):
1621	fd.write(" # Triangulation as an obj file\n")
1622	numVert = str(len(self.meshVertices))
1623
1624	index1 = 1
1625	for vert in self.meshVertices:
1626	vert.index1 = index1
1627	index1 += 1
1628
1629	fd.write("v "
1630	+ str(vert.x) + " "
1631	+ str(vert.y) + " "
1632	+ str(vert.attributes[0]) + "\n")
1633
1634	for tri in self.meshTriangles:
1635	fd.write("f "
1636	+ str(tri.vertices[0].index1) + " "
1637	+ str(tri.vertices[1].index1) + " "
1638	+ str(tri.vertices[2].index1) + "\n")
1639
1640	#FIXME I think this has a bug...
1641	def writeASCIIsegmentoutline(self,
1642	fd,
1643	userVertices,
1644	userSegments,
1645	holes,
1646	regions):
1647	"""Write the user mesh info, only with vertices that are connected to segs
1648	"""
1649	verts = []
1650	#dupindex = 0
1651	for seg in self.userSegments:
1652	verts.append(seg.vertices[0])
1653	verts.append(seg.vertices[1])
1654	print 'verts>',verts
1655
1656	verts, count = self.removeDuplicatedVertices(verts)
1657	print 'verts no dups>',verts
1658	self.writeASCIImesh(fd,
1659	verts,
1660	self.getUserSegments(),
1661	self.holes,
1662	self.regions)
1663
1664	# exportASCIImeshfile
1665	def export_mesh_file(self,ofile):
1666	"""
1667	export a file, ofile, with the format
1668	"""
1669
1670	dict = self.Mesh2IODict()
1671	load_mesh.loadASCII.export_mesh_file(ofile,dict)
1672
1673	# not used
1674	def writeASCIImesh(self,
1675	fd,
1676	userVertices,
1677	userSegments,
1678	holes,
1679	regions):
1680
1681	#print "mesh.writeASCIImesh*********"
1682	#print "dict",dict
1683	#print "mesh*********"
1684	load_mesh.loadASCII.write_ASCII_outline(fd, dict)
1685
1686	#FIXME the title is wrong, need more comments
1687	def exportxyafile(self,ofile):
1688	"""
1689	export a file, ofile, with the format
1690
1691	First line: <# of vertices> <# of attributes>
1692	Following lines: <vertex #> <x> <y> [attributes]
1693	"""
1694	#load_mesh.loadASCII
1695	#FIXME, this should be a mesh2io method
1696	if self.meshVertices == []:
1697	Vertices = self.userVertices
1698	else:
1699	Vertices = self.meshVertices
1700
1701	numVert = str(len(Vertices))
1702
1703	if Vertices == []:
1704	raise RuntimeError
1705	numVertAttrib = str(len(Vertices[0].attributes))
1706	title = numVert + " " + numVertAttrib + " # <vertex #> <x> <y> [attributes]"
1707
1708	#Convert the Vertices to pointlist and pointattributelist
1709	xya_dict = {}
1710	pointattributes = []
1711	points = []
1712
1713	for vert in Vertices:
1714	points.append([vert.x,vert.y])
1715	pointattributes.append(vert.attributes)
1716
1717	xya_dict['pointlist'] = points
1718	xya_dict['pointattributelist'] = pointattributes
1719	xya_dict['geo_reference'] = self.geo_reference
1720
1721	load_mesh.loadASCII.export_xya_file(ofile, xya_dict, title, delimiter = " ")
1722
1723
1724	########### IO CONVERTERS ##################
1725	"""
1726	The dict fromat for IO with .tsh files is;
1727	(the triangulation)
1728	vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
1729	vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1730	vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
1731	segments: [[v1,v2],[v3,v4],...] (lists of integers)
1732	segment_tags : [tag,tag,...] list of strings
1733	triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
1734	triangle tags: [s1,s2,...] A list of strings
1735	triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
1736
1737	(the outline)
1738	points: [[x1,y1],[x2,y2],...] (lists of doubles)
1739	point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1740	outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
1741	outline_segment_tags : [tag1,tag2,...] list of strings
1742	holes : [[x1,y1],...](List of doubles, one inside each hole region)
1743	regions : [ [x1,y1],...] (List of 4 doubles)
1744	region_tags : [tag1,tag2,...] (list of strings)
1745	region_max_areas: [ma1,ma2,...] (A list of doubles)
1746	{Convension: A -ve max area means no max area}
1747
1748	"""
1749
1750
1751
1752	def Mesh2IODict(self):
1753	"""
1754	Convert the triangulation and outline info of a mesh to a dictionary
1755	structure
1756	"""
1757	dict = self.Mesh2IOTriangulationDict()
1758	dict_mesh = self.Mesh2IOOutlineDict()
1759	for element in dict_mesh.keys():
1760	dict[element] = dict_mesh[element]
1761
1762	# add the geo reference
1763	dict['geo_reference'] = self.geo_reference
1764	return dict
1765
1766	def Mesh2IOTriangulationDict(self):
1767	"""
1768	Convert the Mesh to a dictionary of lists describing the
1769	triangulation variables;
1770
1771	Used to produce .tsh file
1772	"""
1773
1774	meshDict = {}
1775	vertices=[]
1776	vertex_attributes=[]
1777
1778
1779	self.maxVertexIndex=0
1780	for vertex in self.meshVertices:
1781	vertex.index = self.maxVertexIndex
1782	vertices.append([vertex.x,vertex.y])
1783	vertex_attributes.append(vertex.attributes)
1784	self.maxVertexIndex += 1
1785
1786	meshDict['vertices'] = vertices
1787	meshDict['vertex_attributes'] = vertex_attributes
1788	meshDict['vertex_attribute_titles'] = self.attributeTitles
1789	#segments
1790	segments=[]
1791	segment_tags=[]
1792	for seg in self.meshSegments:
1793	segments.append([seg.vertices[0].index,seg.vertices[1].index])
1794	segment_tags.append(seg.tag)
1795	meshDict['segments'] =segments
1796	meshDict['segment_tags'] =segment_tags
1797
1798	# Make sure that the indexation is correct
1799	index = 0
1800	for tri in self.meshTriangles:
1801	tri.index = index
1802	index += 1
1803
1804	triangles = []
1805	triangle_tags = []
1806	triangle_neighbors = []
1807	for tri in self.meshTriangles:
1808	triangles.append([tri.vertices[0].index,tri.vertices[1].index,tri.vertices[2].index])
1809	triangle_tags.append([tri.attribute])
1810	neighborlist = [-1,-1,-1]
1811	for neighbor,index in map(None,tri.neighbors,
1812	range(len(tri.neighbors))):
1813	if neighbor:
1814	neighborlist[index] = neighbor.index
1815	triangle_neighbors.append(neighborlist)
1816
1817	meshDict['triangles'] = triangles
1818	meshDict['triangle_tags'] = triangle_tags
1819	meshDict['triangle_neighbors'] = triangle_neighbors
1820
1821	#print "mesh.Mesh2IOTriangulationDict))"
1822	#print meshDict
1823	#print "mesh.Mesh2IOTriangulationDict))"
1824
1825	return meshDict
1826
1827
1828	def Mesh2IOOutlineDict(self, userVertices=None,
1829	userSegments=None,
1830	holes=None,
1831	regions=None):
1832	"""
1833	Convert the mesh outline to a dictionary of the lists needed for the
1834	triang module;
1835
1836	Note, this adds an index attribute to the user Vertex objects.
1837
1838	Used to produce .tsh file and output to triangle
1839	"""
1840	if userVertices is None:
1841	userVertices = self.getUserVertices()
1842	if userSegments is None:
1843	userSegments = self.getUserSegments()
1844	if holes is None:
1845	holes = self.getHoles()
1846	if regions is None:
1847	regions = self.getRegions()
1848
1849	meshDict = {}
1850	#print "userVertices",userVertices
1851	#print "userSegments",userSegments
1852	pointlist=[]
1853	pointattributelist=[]
1854	index = 0
1855	for vertex in userVertices:
1856	vertex.index = index
1857	pointlist.append([vertex.x,vertex.y])
1858	pointattributelist.append(vertex.attributes)
1859
1860	index += 1
1861	meshDict['points'] = pointlist
1862	meshDict['point_attributes'] = pointattributelist
1863
1864	segmentlist=[]
1865	segmenttaglist=[]
1866	for seg in userSegments:
1867	segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
1868	segmenttaglist.append(seg.tag)
1869	meshDict['outline_segments'] =segmentlist
1870	meshDict['outline_segment_tags'] =segmenttaglist
1871
1872	holelist=[]
1873	for hole in holes:
1874	holelist.append([hole.x,hole.y])
1875	meshDict['holes'] = holelist
1876
1877	regionlist=[]
1878	regiontaglist = []
1879	regionmaxarealist = []
1880	for region in regions:
1881	regionlist.append([region.x,region.y])
1882	regiontaglist.append(region.getTag())
1883
1884	if (region.getMaxArea() != None):
1885	regionmaxarealist.append(region.getMaxArea())
1886	else:
1887	regionmaxarealist.append( load_mesh.loadASCII.NOMAXAREA)
1888	meshDict['regions'] = regionlist
1889	meshDict['region_tags'] = regiontaglist
1890	meshDict['region_max_areas'] = regionmaxarealist
1891	#print "(("
1892	#print meshDict
1893	#print meshDict['regionlist']
1894	#print "(("
1895	return meshDict
1896
1897
1898	def IOTriangulation2Mesh(self, genDict):
1899	"""
1900	Set the mesh attributes given an tsh IO dictionary
1901	"""
1902	#Clear the current generated mesh values
1903	self.meshTriangles=[]
1904	self.attributeTitles=[]
1905	self.meshSegments=[]
1906	self.meshVertices=[]
1907
1908	#print "mesh.setTriangulation@#@#@#"
1909	#print genDict
1910	#print "@#@#@#"
1911
1912	self.maxVertexIndex = 0
1913	for point in genDict['vertices']:
1914	v=Vertex(point[0], point[1])
1915	v.index = self.maxVertexIndex
1916	self.maxVertexIndex +=1
1917	self.meshVertices.append(v)
1918
1919	self.attributeTitles = genDict['vertex_attribute_titles']
1920
1921	index = 0
1922	for seg,tag in map(None,genDict['segments'],genDict['segment_tags']):
1923	segObject = Segment( self.meshVertices[seg[0]],
1924	self.meshVertices[seg[1]], tag = tag )
1925	segObject.index = index
1926	index +=1
1927	self.meshSegments.append(segObject)
1928
1929	index = 0
1930	for triangle in genDict['triangles']:
1931	tObject =Triangle( self.meshVertices[triangle[0]],
1932	self.meshVertices[triangle[1]],
1933	self.meshVertices[triangle[2]] )
1934	tObject.index = index
1935	index +=1
1936	self.meshTriangles.append(tObject)
1937
1938	index = 0
1939	for att in genDict['triangle_tags']:
1940	if att == []:
1941	self.meshTriangles[index].setAttribute("")
1942	else:
1943	# Note, is the first attribute always the region att?
1944	# haven't confirmed this
1945	#Peter - I think so (from manuel)
1946	#...the first such value is assumed to be a regional attribute...
1947	self.meshTriangles[index].setAttribute(att[0])
1948	index += 1
1949
1950	index = 0
1951	for att in genDict['vertex_attributes']:
1952	if att == None:
1953	self.meshVertices[index].setAttributes([])
1954	else:
1955	self.meshVertices[index].setAttributes(att)
1956	index += 1
1957
1958	index = 0
1959	for triangle in genDict['triangle_neighbors']:
1960	# Build a list of triangle object neighbors
1961	ObjectNeighbor = []
1962	for neighbor in triangle:
1963	if ( neighbor != -1):
1964	ObjectNeighbor.append(self.meshTriangles[neighbor])
1965	else:
1966	ObjectNeighbor.append(None)
1967	self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],ObjectNeighbor[1],ObjectNeighbor[2])
1968	index += 1
1969
1970
1971	def IOOutline2Mesh(self, genDict):
1972	"""
1973	Set the outline (user Mesh attributes) given a IO tsh dictionary
1974
1975	mesh is an instance of a mesh object
1976	"""
1977	#Clear the current user mesh values
1978	self.clearUserSegments()
1979	self.userVertices=[]
1980	self.Holes=[]
1981	self.Regions=[]
1982
1983	#print "mesh.IOOutline2Mesh@#@#@#"
1984	#print "genDict",genDict
1985	#print "@#@#@#"
1986
1987	#index = 0
1988	for point in genDict['points']:
1989	v=Vertex(point[0], point[1])
1990	#v.index = index
1991	#index +=1
1992	self.userVertices.append(v)
1993
1994	#index = 0
1995	for seg,tag in map(None,genDict['outline_segments'],genDict['outline_segment_tags']):
1996	segObject = Segment( self.userVertices[seg[0]],
1997	self.userVertices[seg[1]], tag = tag )
1998	#segObject.index = index
1999	#index +=1
2000	self.userSegments.append(segObject)
2001
2002	# Remove the loading of attribute info.
2003	# Have attribute info added using least_squares in pyvolution
2004	# index = 0
2005	# for att in genDict['point_attributes']:
2006	# if att == None:
2007	# self.userVertices[index].setAttributes([])
2008	# else:
2009	# self.userVertices[index].setAttributes(att)
2010	# index += 1
2011
2012	#index = 0
2013	for point in genDict['holes']:
2014	h=Hole(point[0], point[1])
2015	#h.index = index
2016	#index +=1
2017	self.holes.append(h)
2018
2019	#index = 0
2020	for reg,att,maxArea in map(None,
2021	genDict['regions'],
2022	genDict['region_tags'],
2023	genDict['region_max_areas']):
2024	if maxArea > 0: # maybe I should ref NOMAXAREA? Prob' not though
2025	Object = Region( reg[0],
2026	reg[1],
2027	tag = att,
2028	maxArea = maxArea)
2029	else:
2030	Object = Region( reg[0],
2031	reg[1],
2032	tag = att)
2033
2034	#Object.index = index
2035	#index +=1
2036	self.regions.append(Object)
2037
2038	############################################
2039
2040
2041	def refineSet(self,setName):
2042	Triangles = self.sets[self.setID[setName]]
2043	Refine(self,Triangles)
2044
2045	def selectAllTriangles(self):
2046	A=[]
2047	A.extend(self.meshTriangles)
2048	if not('All' in self.setID.keys()):
2049	self.setID['All']=len(self.sets)
2050	self.sets.append(A)
2051	else:
2052	self.sets[self.setID['All']]=A
2053	return 'All'
2054	# and objectIDs
2055
2056
2057	def clearSelection(self):
2058	A = []
2059	if not('None' in self.setID.keys()):
2060	self.setID['None']=len(self.sets)
2061	self.sets.append(A)
2062	return 'None'
2063
2064	def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
2065	#FIXME Draws over previous triangles - may bloat canvas
2066	Triangles = self.sets[self.setID[setName]]
2067	for triangle in Triangles:
2068	triangle.draw(canvas,1,
2069	scale = SCALE,
2070	colour = colour)
2071
2072	def undrawSet(self,canvas,setName,SCALE,colour='green'):
2073	#FIXME Draws over previous lines - may bloat canvas
2074	Triangles = self.sets[self.setID[setName]]
2075	for triangle in Triangles:
2076	triangle.draw(canvas,1,
2077	scale = SCALE,
2078	colour = colour)
2079
2080
2081	def triangles_to_polySet(self,setName):
2082
2083	Triangles = self.sets[self.setID[setName]]
2084	Triangles_dict = {}
2085	for triangle in Triangles:
2086	Triangles_dict[triangle]=None
2087
2088	userVertices = {}
2089	userSegments = []
2090	index = 0
2091	for triangle in Triangles:
2092	for i in (0,1,2):
2093	if not Triangles_dict.has_key(triangle.neighbors[i]):
2094	a = triangle.vertices[i-1]
2095	b = triangle.vertices[i-2]
2096	userVertices[a]=0
2097	userVertices[b]=0
2098	userSegments.append([a,b])
2099	return userVertices, userSegments
2100
2101	def threshold(self,setName,min=None,max=None,attribute_name = 'elevation'):
2102	"""
2103	threshold using d
2104	"""
2105	triangles = self.sets[self.setID[setName]]
2106	A = []
2107
2108	if attribute_name in self.attributeTitles:
2109	i = self.attributeTitles.index(attribute_name)
2110	if not max == None:
2111	for t in triangles:
2112	if (min<self.av_att(t,i)<max):
2113	A.append(t)
2114	else:
2115	for t in triangles:
2116	if (min<self.av_att(t,i)):
2117	A.append(t)
2118	self.sets[self.setID[setName]] = A
2119
2120	def general_threshold(self,setName,min=None,max=None,attribute_name = 'elevation',function = None):
2121	"""
2122	threshold using d
2123	"""
2124	triangles = self.sets[self.setID[setName]]
2125	A = []
2126
2127	if attribute_name in self.attributeTitles:
2128	i = self.attributeTitles.index(attribute_name)
2129	if not max == None:
2130	for t in triangles:
2131	if (min<function(t,i)<max):
2132	A.append(t)
2133	else:
2134	for t in triangles:
2135	if (min<self.function(t,i)):
2136	A.append(t)
2137	self.sets[self.setID[setName]] = A
2138
2139	def av_att(self,triangle,i):
2140	#evaluates the average attribute of the vertices of a triangle.
2141	V = triangle.getVertices()
2142	a0 = (V[0].attributes[i])
2143	a1 = (V[1].attributes[i])
2144	a2 = (V[2].attributes[i])
2145	return (a0+a1+a2)/3
2146
2147	def Courant_ratio(self,triangle,index):
2148	"""
2149	Not the true Courant ratio, just elevation on area
2150	"""
2151	e = self.av_att(triangle,index)
2152	A = triangle.calcArea()
2153	return e/A
2154
2155	def Gradient(self,triangle,index):
2156	V = triangle.vertices
2157	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]
2158	grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
2159	if ((grad_x2)+(grad_y2))**(0.5)<0:
2160	print ((grad_x2)+(grad_y2))**(0.5)
2161	return ((grad_x2)+(grad_y2))**(0.5)
2162
2163
2164	def append_triangle(self,triangle):
2165	self.meshTriangles.append(triangle)
2166
2167	def replace_triangle(self,triangle,replacement):
2168	i = self.meshTriangles.index(triangle)
2169	self.meshTriangles[i]=replacement
2170	assert replacement in self.meshTriangles
2171
2172
2173
2174
2175	def importUngenerateFile(ofile):
2176	"""
2177	import a file, ofile, with the format
2178	[poly]
2179	poly format:
2180	First line: <# of vertices> <x centroid> <y centroid>
2181	Following lines: <x> <y>
2182	last line: "END"
2183
2184	Note: These are clockwise.
2185	"""
2186	fd = open(ofile,'r')
2187	Dict = readUngenerateFile(fd)
2188	fd.close()
2189	return Dict
2190
2191	def readUngenerateFile(fd):
2192	"""
2193	import a file, ofile, with the format
2194	[poly]
2195	poly format:
2196	First line: <# of polynomial> <x centroid> <y centroid>
2197	Following lines: <x> <y>
2198	last line: "END"
2199	"""
2200	END_DELIMITER = 'END\n'
2201
2202	points = []
2203	segments = []
2204
2205	isEnd = False
2206	line = fd.readline() #not used <# of polynomial> <x> <y>
2207	while not isEnd:
2208	line = fd.readline()
2209	fragments = line.split()
2210	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2211	points.append(vert)
2212	PreviousVertIndex = len(points)-1
2213	firstVertIndex = PreviousVertIndex
2214
2215	line = fd.readline() #Read the next line
2216	while line <> END_DELIMITER:
2217	#print "line >" + line + "<"
2218	fragments = line.split()
2219	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2220	points.append(vert)
2221	thisVertIndex = len(points)-1
2222	segment = [PreviousVertIndex,thisVertIndex]
2223	segments.append(segment)
2224	PreviousVertIndex = thisVertIndex
2225	line = fd.readline() #Read the next line
2226	i =+ 1
2227	# If the last and first segments are the same,
2228	# Remove the last segment and the last vertex
2229	# then add a segment from the second last vert to the 1st vert
2230	thisVertIndex = len(points)-1
2231	firstVert = points[firstVertIndex]
2232	thisVert = points[thisVertIndex]
2233	#print "firstVert",firstVert
2234	#print "thisVert",thisVert
2235	if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
2236	points.pop()
2237	segments.pop()
2238	thisVertIndex = len(points)-1
2239	segments.append([thisVertIndex, firstVertIndex])
2240
2241	line = fd.readline() # read <# of polynomial> <x> <y> OR END
2242	#print "line >>" + line + "<<"
2243	if line == END_DELIMITER:
2244	isEnd = True
2245
2246	#print "points", points
2247	#print "segments", segments
2248	ungenerated_dict = {}
2249	ungenerated_dict['pointlist'] = points
2250	ungenerated_dict['segmentlist'] = segments
2251	return ungenerated_dict
2252
2253	def importMeshFromFile(ofile):
2254	"""returns a mesh object, made from a .xya/.pts or .tsh/.msh file
2255	Often raises SyntaxError, IOError
2256	"""
2257	newmesh = None
2258	if ofile[-4:]== ".xya":
2259	#print "loading " + ofile
2260	try:
2261	dict = load_mesh.loadASCII.load_points_file(ofile, delimiter = ',')
2262	except SyntaxError:
2263	dict = load_mesh.loadASCII.load_points_file(ofile, delimiter = ' ')
2264	#print "dict",dict outline_
2265	dict['segmentlist'] = []
2266	dict['segmenttaglist'] = []
2267	dict['regionlist'] = []
2268	dict['regionattributelist'] = []
2269	dict['regionmaxarealist'] = []
2270	dict['holelist'] = []
2271	newmesh= Mesh()
2272	newmesh.setMesh(dict) #FIXME use IOOutline2Mesh
2273	counter = newmesh.removeDuplicatedUserVertices()
2274	if (counter >0):
2275	print "%i duplicate vertices removed from dataset" % (counter)
2276	elif ofile[-4:]== ".pts":
2277	#print "loading " + ofile
2278	dict = load_mesh.loadASCII.load_points_file(ofile)
2279	#print "dict",dict
2280	dict['points'] = dict['pointlist']
2281	dict['outline_segments'] = []
2282	dict['outline_segment_tags'] = []
2283	dict['regions'] = []
2284	dict['region_tags'] = []
2285	dict['region_max_areas'] = []
2286	dict['holes'] = []
2287	newmesh= Mesh(geo_reference = dict['geo_reference'])
2288	newmesh.IOOutline2Mesh(dict) #FIXME use IOOutline2Mesh
2289	counter = newmesh.removeDuplicatedUserVertices()
2290	if (counter >0):
2291	print "%i duplicate vertices removed from dataset" % (counter)
2292	elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
2293	dict = load_mesh.loadASCII.import_triangulation(ofile)
2294	#print "********"
2295	#print "zq mesh.dict",dict
2296	#print "********"
2297	newmesh= Mesh()
2298	newmesh.IOOutline2Mesh(dict)
2299	newmesh.IOTriangulation2Mesh(dict)
2300	else:
2301	raise RuntimeError
2302
2303	if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
2304	newmesh.geo_reference = dict['geo_reference']
2305	return newmesh
2306
2307	def loadPickle(currentName):
2308	fd = open(currentName)
2309	mesh = pickle.load(fd)
2310	fd.close()
2311	return mesh
2312
2313	def square_outline(side_length = 1,up = "top", left = "left", right = "right",
2314	down = "bottom", regions = False):
2315
2316	a = Vertex (0,0)
2317	b = Vertex (0,side_length)
2318	c = Vertex (side_length,0)
2319	d = Vertex (side_length,side_length)
2320
2321	s2 = Segment(b,d, tag = up)
2322	s3 = Segment(b,a, tag = left)
2323	s4 = Segment(d,c, tag = right)
2324	s5 = Segment(a,c, tag = down)
2325
2326	if regions:
2327	e = Vertex (side_length/2,side_length/2)
2328	s6 = Segment(a,e, tag = down + left)
2329	s7 = Segment(b,e, tag = up + left)
2330	s8 = Segment(c,e, tag = down + right)
2331	s9 = Segment(d,e, tag = up + right)
2332	r1 = Region(side_length/2,3.*side_length/4, tag = up)
2333	r2 = Region(1.*side_length/4,side_length/2, tag = left)
2334	r3 = Region(3.*side_length/4,side_length/2, tag = right)
2335	r4 = Region(side_length/2,1.*side_length/4, tag = down)
2336	mesh = Mesh(userVertices=[a,b,c,d,e],
2337	userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
2338	regions = [r1,r2,r3,r4])
2339	else:
2340	mesh = Mesh(userVertices=[a,b,c,d],
2341	userSegments=[s2,s3,s4,s5])
2342
2343	return mesh
2344
2345
2346
2347	def region_strings2ints(region_list):
2348	"""Given a list of (x_int,y_int,tag_string) lists it returns a list of
2349	(x_int,y_int,tag_int) and a list to convert the tag_int's back to
2350	the tag_strings
2351	"""
2352	# Make sure "" has an index of 0
2353	region_list.reverse()
2354	region_list.append((1.0,2.0,""))
2355	region_list.reverse()
2356	convertint2string = []
2357	for i in xrange(len(region_list)):
2358	convertint2string.append(region_list[i][2])
2359	if len(region_list[i]) == 4: # there's an area value
2360	region_list[i] = (region_list[i][0],
2361	region_list[i][1],i,region_list[i][3])
2362	elif len(region_list[i]) == 3: # no area value
2363	region_list[i] = (region_list[i][0],region_list[i][1],i)
2364	else:
2365	print "The region list has a bad size"
2366	# raise an error ..
2367	raise Error
2368
2369	#remove "" from the region_list
2370	region_list.pop(0)
2371
2372	return [region_list, convertint2string]
2373
2374	def region_ints2strings(region_list,convertint2string):
2375	"""Reverses the transformation of region_strings2ints
2376	"""
2377	if region_list[0] != []:
2378	for i in xrange(len(region_list)):
2379	region_list[i] = [convertint2string[int(region_list[i][0])]]
2380	return region_list
2381
2382	def segment_ints2strings(intlist, convertint2string):
2383	"""Reverses the transformation of segment_strings2ints """
2384	stringlist = []
2385	for x in intlist:
2386	stringlist.append(convertint2string[x])
2387	return stringlist
2388
2389	def segment_strings2ints(stringlist, preset):
2390	"""Given a list of strings return a list of 0 to n ints which represent
2391	the strings and a converting list of the strings, indexed by 0 to n ints.
2392	Also, input an initial converting list of the strings
2393	Note, the converting list starts off with
2394	["internal boundary", "external boundary", "internal boundary"]
2395	example input and output
2396	input ["a","b","a","c"],["c"]
2397	output [[2, 1, 2, 0], ['c', 'b', 'a']]
2398
2399	the first element in the converting list will be
2400	overwritten with "".
2401	?This will become the third item in the converting list?
2402
2403	# note, order the initial converting list is important,
2404	since the index = the int tag
2405	"""
2406	nodups = unique(stringlist)
2407	# note, order is important, the index = the int tag
2408	#preset = ["internal boundary", "external boundary"]
2409	#Remove the preset tags from the list with no duplicates
2410	nodups = [x for x in nodups if x not in preset]
2411
2412	try:
2413	nodups.remove("") # this has to go to zero
2414	except ValueError:
2415	pass
2416
2417	# Add the preset list at the beginning of no duplicates
2418	preset.reverse()
2419	nodups.extend(preset)
2420	nodups.reverse()
2421
2422
2423	convertstring2int = {}
2424	convertint2string = []
2425	index = 0
2426	for x in nodups:
2427	convertstring2int[x] = index
2428	convertint2string.append(x)
2429	index += 1
2430	convertstring2int[""] = 0
2431
2432	intlist = []
2433	for x in stringlist:
2434	intlist.append(convertstring2int[x])
2435	return [intlist, convertint2string]
2436
2437
2438
2439
2440	def unique(s):
2441	"""Return a list of the elements in s, but without duplicates.
2442
2443	For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
2444	unique("abcabc") some permutation of ["a", "b", "c"], and
2445	unique(([1, 2], [2, 3], [1, 2])) some permutation of
2446	[[2, 3], [1, 2]].
2447
2448	For best speed, all sequence elements should be hashable. Then
2449	unique() will usually work in linear time.
2450
2451	If not possible, the sequence elements should enjoy a total
2452	ordering, and if list(s).sort() doesn't raise TypeError it's
2453	assumed that they do enjoy a total ordering. Then unique() will
2454	usually work in O(N*log2(N)) time.
2455
2456	If that's not possible either, the sequence elements must support
2457	equality-testing. Then unique() will usually work in quadratic
2458	time.
2459	"""
2460
2461	n = len(s)
2462	if n == 0:
2463	return []
2464
2465	# Try using a dict first, as that's the fastest and will usually
2466	# work. If it doesn't work, it will usually fail quickly, so it
2467	# usually doesn't cost much to try it. It requires that all the
2468	# sequence elements be hashable, and support equality comparison.
2469	u = {}
2470	try:
2471	for x in s:
2472	u[x] = 1
2473	except TypeError:
2474	del u # move on to the next method
2475	else:
2476	return u.keys()
2477
2478	# We can't hash all the elements. Second fastest is to sort,
2479	# which brings the equal elements together; then duplicates are
2480	# easy to weed out in a single pass.
2481	# NOTE: Python's list.sort() was designed to be efficient in the
2482	# presence of many duplicate elements. This isn't true of all
2483	# sort functions in all languages or libraries, so this approach
2484	# is more effective in Python than it may be elsewhere.
2485	try:
2486	t = list(s)
2487	t.sort()
2488	except TypeError:
2489	del t # move on to the next method
2490	else:
2491	assert n > 0
2492	last = t[0]
2493	lasti = i = 1
2494	while i < n:
2495	if t[i] != last:
2496	t[lasti] = last = t[i]
2497	lasti += 1
2498	i += 1
2499	return t[:lasti]
2500
2501	# Brute force is all that's left.
2502	u = []
2503	for x in s:
2504	if x not in u:
2505	u.append(x)
2506	return u
2507
2508	"""Refines triangles
2509
2510	Implements the #triangular bisection?# algorithm.
2511
2512
2513	"""
2514
2515	def Refine(mesh, triangles):
2516	"""
2517	Given a general_mesh, and a triangle number, split
2518	that triangle in the mesh in half. Then to prevent
2519	vertices and edges from meeting, keep refining
2520	neighbouring triangles until the mesh is clean.
2521	"""
2522	state = BisectionState(mesh)
2523	for triangle in triangles:
2524	if not state.refined_triangles.has_key(triangle):
2525	triangle.rotate_longest_side()
2526	state.start(triangle)
2527	Refine_mesh(mesh, state)
2528
2529	def Refine_mesh(mesh, state):
2530	"""
2531	"""
2532	state.getState(mesh)
2533	refine_triangle(mesh,state)
2534	state.evolve()
2535	if not state.end:
2536	Refine_mesh(mesh,state)
2537
2538	def refine_triangle(mesh,state):
2539	split(mesh,state.current_triangle,state.new_point)
2540	if state.case == 'one':
2541	state.r[3]=state.current_triangle#triangle 2
2542
2543	new_triangle_id = len(mesh.meshTriangles)-1
2544	new_triangle = mesh.meshTriangles[new_triangle_id]
2545
2546	split(mesh,new_triangle,state.old_point)
2547	state.r[2]=new_triangle#triangle 1.2
2548	state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
2549	r = state.r
2550	state.repairCaseOne()
2551
2552	if state.case == 'two':
2553	state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2554
2555	new_triangle = state.current_triangle
2556
2557	split(mesh,new_triangle,state.old_point)
2558
2559	state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
2560	state.r[4]=new_triangle#triangle 2.2
2561	r = state.r
2562
2563	state.repairCaseTwo()
2564
2565	if state.case == 'vertex':
2566	state.r[2]=state.current_triangle#triangle 2
2567	state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2568	r = state.r
2569	state.repairCaseVertex()
2570
2571	if state.case == 'start':
2572	state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2573	state.r[3]=state.current_triangle#triangle 2
2574
2575	if state.next_case == 'boundary':
2576	state.repairCaseBoundary()
2577
2578
2579	def split(mesh, triangle, new_point):
2580	"""
2581	Given a mesh, triangle_id and a new point,
2582	split the corrosponding triangle into two
2583	new triangles and update the mesh.
2584	"""
2585
2586	new_triangle1 = Triangle(new_point,triangle.vertices[0],triangle.vertices[1],attribute = triangle.attribute, neighbors = None)
2587	new_triangle2 = Triangle(new_point,triangle.vertices[2],triangle.vertices[0],attribute = triangle.attribute, neighbors = None)
2588
2589	new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
2590	new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
2591
2592	mesh.meshTriangles.append(new_triangle1)
2593
2594	triangle.vertices = new_triangle2.vertices
2595	triangle.neighbors = new_triangle2.neighbors
2596
2597
2598	class State:
2599
2600	def __init__(self):
2601	pass
2602
2603	class BisectionState(State):
2604
2605
2606	def __init__(self,mesh):
2607	self.len = len(mesh.meshTriangles)
2608	self.refined_triangles = {}
2609	self.mesh = mesh
2610	self.current_triangle = None
2611	self.case = 'start'
2612	self.end = False
2613	self.r = [None,None,None,None,None]
2614
2615	def start(self, triangle):
2616	self.current_triangle = triangle
2617	self.case = 'start'
2618	self.end = False
2619	self.r = [None,None,None,None,None]
2620
2621	def getState(self,mesh):
2622	if not self.case == 'vertex':
2623	self.new_point=self.getNewVertex(mesh, self.current_triangle)
2624	#self.neighbour=self.getNeighbour(mesh, self.current_triangle)
2625	self.neighbour = self.current_triangle.neighbors[0]
2626	if not self.neighbour is None:
2627	self.neighbour.rotate_longest_side()
2628	self.next_case = self.get_next_case(mesh,self.neighbour,self.current_triangle)
2629	if self.case == 'vertex':
2630	self.new_point=self.old_point
2631
2632
2633	def evolve(self):
2634	if self.case == 'vertex':
2635	self.end = True
2636
2637	self.last_case = self.case
2638	self.case = self.next_case
2639
2640	self.old_point = self.new_point
2641	self.current_triangle = self.neighbour
2642
2643	if self.case == 'boundary':
2644	self.end = True
2645	self.refined_triangles[self.r[2]]=1
2646	self.refined_triangles[self.r[3]]=1
2647	if not self.r[4] is None:
2648	self.refined_triangles[self.r[4]]=1
2649	self.r[0]=self.r[2]
2650	self.r[1]=self.r[3]
2651
2652
2653	def getNewVertex(self,mesh,triangle):
2654	coordinate1 = triangle.vertices[1]
2655	coordinate2 = triangle.vertices[2]
2656	a = ([coordinate1.x1.,coordinate1.y1.])
2657	b = ([coordinate2.x1.,coordinate2.y1.])
2658	attributes = []
2659	for i in range(len(coordinate1.attributes)):
2660	att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
2661	attributes.append(att)
2662	new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
2663	newVertex = Vertex(new_coordinate[0],new_coordinate[1], attributes = attributes)
2664	mesh.maxVertexIndex+=1
2665	newVertex.index = mesh.maxVertexIndex
2666	mesh.meshVertices.append(newVertex)
2667	return newVertex
2668
2669	def get_next_case(self, mesh,neighbour,triangle):
2670	"""
2671	Given the locations of two neighbouring triangles,
2672	examine the interior indices of their vertices (i.e.
2673	0,1 or 2) to determine what how the neighbour needs
2674	to be refined.
2675	"""
2676	if (neighbour is None):
2677	next_case = 'boundary'
2678	else:
2679	if triangle.vertices[1].x==neighbour.vertices[2].x:
2680	if triangle.vertices[1].y==neighbour.vertices[2].y:
2681	next_case = 'vertex'
2682	if triangle.vertices[1].x==neighbour.vertices[0].x:
2683	if triangle.vertices[1].y==neighbour.vertices[0].y:
2684	next_case = 'two'
2685	if triangle.vertices[1].x==neighbour.vertices[1].x:
2686	if triangle.vertices[1].y==neighbour.vertices[1].y:
2687	next_case = 'one'
2688	return next_case
2689
2690
2691
2692	def repairCaseVertex(self):
2693
2694	r = self.r
2695
2696
2697	self.link(r[0],r[2])
2698	self.repair(r[0])
2699
2700	self.link(r[1],r[3])
2701	self.repair(r[1])
2702
2703	self.repair(r[2])
2704
2705	self.repair(r[3])
2706
2707
2708	def repairCaseOne(self):
2709	r = self.r
2710
2711
2712	self.link(r[0],r[2])
2713	self.repair(r[0])
2714
2715	self.link(r[1],r[4])
2716	self.repair(r[1])
2717
2718	self.repair(r[4])
2719
2720	def repairCaseTwo(self):
2721	r = self.r
2722
2723	self.link(r[0],r[4])
2724	self.repair(r[0])
2725
2726	self.link(r[1],r[3])
2727	self.repair(r[1])
2728
2729	self.repair(r[4])
2730
2731	def repairCaseBoundary(self):
2732	r = self.r
2733	self.repair(r[2])
2734	self.repair(r[3])
2735
2736
2737
2738	def repair(self,triangle):
2739	"""
2740	Given a triangle that knows its neighbours, this will
2741	force the neighbours to comply.
2742
2743	However, it needs to compare the vertices of triangles
2744	for this implementation
2745
2746	But it doesn't work for invalid neighbour structures
2747	"""
2748	n=triangle.neighbors
2749	for i in (0,1,2):
2750	if not n[i] is None:
2751	for j in (0,1,2):#to find which side of the list is broken
2752	if not (n[i].vertices[j] in triangle.vertices):
2753	#ie if j is the side of n that needs fixing
2754	k = j
2755	n[i].neighbors[k]=triangle
2756
2757
2758
2759	def link(self,triangle1,triangle2):
2760	"""
2761	make triangle1 neighbors point to t
2762	#count = 0riangle2
2763	"""
2764	count = 0
2765	for i in (0,1,2):#to find which side of the list is broken
2766	if not (triangle1.vertices[i] in triangle2.vertices):
2767	j = i
2768	count+=1
2769	assert count == 1
2770	triangle1.neighbors[j]=triangle2
2771
2772
2773
2774	if __name__ == "__main__":
2775	#from mesh import *
2776	# THIS CAN BE DELETED
2777	m = Mesh()
2778	dict = importUngenerateFile("ungen_test.txt")
2779	m.addVertsSegs(dict)
2780	print m

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

Download in other formats: