Context Navigation

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

Last change on this file since 1035 was 1035, checked in by duncan, 20 years ago
all meshes now have geo-refs. added some assertions
File size: 95.4 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 getUserVertices(self):
628	return self.userVertices
629
630	def getUserSegments(self):
631	allSegments = self.userSegments + self.alphaUserSegments
632	#print "self.userSegments",self.userSegments
633	#print "self.alphaUserSegments",self.alphaUserSegments
634	#print "allSegments",allSegments
635	return allSegments
636
637	def deleteUserSegments(self,seg):
638	if self.userSegments.count(seg) == 0:
639	self.alphaUserSegments.remove(seg)
640	pass
641	else:
642	self.userSegments.remove(seg)
643
644	def clearUserSegments(self):
645	self.userSegments = []
646	self.alphaUserSegments = []
647
648	def getTriangulation(self):
649	return self.meshTriangles
650
651	def getMeshVertices(self):
652	return self.meshVertices
653
654	def getMeshSegments(self):
655	return self.meshSegments
656
657	def getHoles(self):
658	return self.holes
659
660	def getRegions(self):
661	return self.regions
662
663	def isTriangulation(self):
664	if self.meshVertices == []:
665	return False
666	else:
667	return True
668
669	def addUserSegment(self, v1,v2):
670	"""
671	PRECON: A segment between the two vertices is not already present.
672	Check by calling isUserSegmentNew before calling this function.
673
674	"""
675	s=Segment( v1,v2)
676	self.userSegments.append(s)
677	return s
678
679	def clearTriangulation(self):
680
681	#Clear the current generated mesh values
682	self.meshTriangles=[]
683	self.meshSegments=[]
684	self.meshVertices=[]
685
686	def removeDuplicatedUserVertices(self):
687	"""Pre-condition: There are no user segments
688	This function will keep the first duplicate
689	"""
690	assert self.getUserSegments() == []
691	self.userVertices, counter = self.removeDuplicatedVertices(self.userVertices)
692	return counter
693
694	def removeDuplicatedVertices(self, Vertices):
695	"""
696	This function will keep the first duplicate, remove all others
697	Precondition: Each vertex has a dupindex, which is the list
698	index.
699	"""
700	remove = []
701	index = 0
702	for v in Vertices:
703	v.dupindex = index
704	index += 1
705	t = list(Vertices)
706	t.sort(Point.cmp_xy)
707
708	length = len(t)
709	behind = 0
710	ahead = 1
711	counter = 0
712	while ahead < length:
713	b = t[behind]
714	ah = t[ahead]
715	if (b.y == ah.y and b.x == ah.x):
716	remove.append(ah.dupindex)
717	behind += 1
718	ahead += 1
719
720	# remove the duplicate vertices
721	remove.sort()
722	remove.reverse()
723	for i in remove:
724	Vertices.pop(i)
725	pass
726
727	#Remove the attribute that this function added
728	for v in Vertices:
729	del v.dupindex
730	return Vertices,counter
731
732	def thinoutVertices(self, delta):
733	"""Pre-condition: There are no user segments
734	This function will keep the first duplicate
735	"""
736	assert self.getUserSegments() == []
737	#t = self.userVertices
738	#self.userVertices =[]
739	boxedVertices = {}
740	thinnedUserVertices =[]
741	delta = round(delta,1)
742
743	for v in self.userVertices :
744	# tag is the center of the boxes
745	tag = (round(v.x/delta,0)delta,round(v.y/delta,0)delta)
746	#this creates a dict of lists of faces, indexed by tag
747	boxedVertices.setdefault(tag,[]).append(v)
748
749	for [tag,verts] in boxedVertices.items():
750	min = delta
751	bestVert = None
752	tagVert = Vertex(tag[0],tag[1])
753	for v in verts:
754	dist = v.DistanceToPoint(tagVert)
755	if (dist<min):
756	min = dist
757	bestVert = v
758	thinnedUserVertices.append(bestVert)
759	self.userVertices =thinnedUserVertices
760
761
762	def isUserSegmentNew(self, v1,v2):
763	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) ]
764
765	return len(identicalSegs) == 0
766
767	def representedAlphaUserSegment(self, v1,v2):
768	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) ]
769
770	if identicalSegs == []:
771	return None
772	else:
773	# Only return the first one.
774	return identicalSegs[0]
775
776	def representedUserSegment(self, v1,v2):
777	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) ]
778
779	if identicalSegs == []:
780	return None
781	else:
782	# Only return the first one.
783	return identicalSegs[0]
784
785	def deleteSegsOfVertex(self, delVertex):
786	"""
787	Delete this vertex and any segments that connect to it.
788	"""
789	#Find segments that connect to delVertex
790	deletedSegments = []
791	for seg in self.getUserSegments():
792	if (delVertex in seg.vertices):
793	deletedSegments.append(seg)
794	# Delete segments that connect to delVertex
795	for seg in deletedSegments:
796	self.deleteUserSegments(seg)
797	return deletedSegments
798
799
800	def deleteMeshObject(self, MeshObject):
801	"""
802	Returns a list of all objects that were removed
803	"""
804	deletedObs = []
805	if isinstance(MeshObject, Vertex ):
806	deletedObs = self.deleteSegsOfVertex(MeshObject)
807	deletedObs.append(MeshObject)
808	self.userVertices.remove(MeshObject)
809	elif isinstance(MeshObject, Segment):
810	deletedObs.append(MeshObject)
811	self.deleteUserSegments(MeshObject)
812	elif isinstance(MeshObject, Hole):
813	deletedObs.append(MeshObject)
814	self.holes.remove(MeshObject)
815	elif isinstance(MeshObject, Region):
816	deletedObs.append(MeshObject)
817	self.regions.remove(MeshObject)
818	return deletedObs
819
820	def Mesh2triangList(self, userVertices=None,
821	userSegments=None,
822	holes=None,
823	regions=None):
824	"""
825	Convert the Mesh to a dictionary of the lists needed for the triang modul;
826	points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
827	pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
828	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
829	hole list: [(x1,y1),...](Tuples of doubles, one inside each hole region)
830	regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
831
832	Note, this adds an index attribute to the user Vertex objects.
833
834	Used to produce output to triangle
835	"""
836	if userVertices is None:
837	userVertices = self.getUserVertices()
838	if userSegments is None:
839	userSegments = self.getUserSegments()
840	if holes is None:
841	holes = self.getHoles()
842	if regions is None:
843	regions = self.getRegions()
844
845	meshDict = {}
846
847	pointlist=[]
848	pointattributelist=[]
849	index = 0
850	for vertex in userVertices:
851	vertex.index = index
852	pointlist.append((vertex.x,vertex.y))
853	pointattributelist.append((vertex.attributes))
854
855	index += 1
856	meshDict['pointlist'] = pointlist
857	meshDict['pointattributelist'] = pointattributelist
858
859	segmentlist=[]
860	segmenttaglist=[]
861	for seg in userSegments:
862	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
863	segmenttaglist.append(seg.tag)
864	meshDict['segmentlist'] =segmentlist
865	meshDict['segmenttaglist'] =segmenttaglist
866
867	holelist=[]
868	for hole in holes:
869	holelist.append((hole.x,hole.y))
870	meshDict['holelist'] = holelist
871
872	regionlist=[]
873	for region in regions:
874	if (region.getMaxArea() != None):
875	regionlist.append((region.x,region.y,region.getTag(),
876	region.getMaxArea()))
877	else:
878	regionlist.append((region.x,region.y,region.getTag()))
879	meshDict['regionlist'] = regionlist
880	#print "(("
881	#print meshDict
882	#print meshDict['regionlist']
883	#print "(("
884	return meshDict
885
886	def Mesh2MeshList(self):
887	"""
888	Convert the Mesh to a dictionary of lists describing the triangulation variables;
889	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
890	generated point attribute list: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
891	generated point attribute title list:[A1Title, A2Title ...] (list of strings)
892	generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
893	generated segment tag list: [tag,tag,...] list of strings
894
895	generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
896
897	generated triangle attribute list: [s1,s2,...] list of strings
898
899	generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....] tuple of triangles
900
901	Used to produce .tsh file
902	"""
903
904	meshDict = {}
905	pointlist=[]
906	pointattributelist=[]
907
908
909	self.maxVertexIndex=0
910	for vertex in self.meshVertices:
911	vertex.index = self.maxVertexIndex
912	pointlist.append((vertex.x,vertex.y))
913	pointattributelist.append((vertex.attributes))
914	self.maxVertexIndex += 1
915
916	meshDict['generatedpointlist'] = pointlist
917	meshDict['generatedpointattributelist'] = pointattributelist
918	meshDict['generatedpointattributetitlelist'] = self.attributeTitles
919	#segments
920	segmentlist=[]
921	segmenttaglist=[]
922	for seg in self.meshSegments:
923	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
924	segmenttaglist.append(seg.tag)
925	meshDict['generatedsegmentlist'] =segmentlist
926	meshDict['generatedsegmenttaglist'] =segmenttaglist
927
928	# Make sure that the indexation is correct
929	index = 0
930	for tri in self.meshTriangles:
931	tri.index = index
932	index += 1
933
934	trianglelist = []
935	triangleattributelist = []
936	triangleneighborlist = []
937	for tri in self.meshTriangles:
938	trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,tri.vertices[2].index))
939	triangleattributelist.append([tri.attribute])
940	neighborlist = [-1,-1,-1]
941	for neighbor,index in map(None,tri.neighbors,
942	range(len(tri.neighbors))):
943	if neighbor:
944	neighborlist[index] = neighbor.index
945	triangleneighborlist.append(neighborlist)
946
947	meshDict['generatedtrianglelist'] = trianglelist
948	meshDict['generatedtriangleattributelist'] = triangleattributelist
949	meshDict['generatedtriangleneighborlist'] = triangleneighborlist
950
951	#print "mesh.Mesh2MeshList))"
952	#print meshDict
953	#print "mesh.Mesh2MeshList))"
954
955	return meshDict
956
957
958	def Mesh2MeshDic(self):
959	"""
960	Convert the user and generated info of a mesh to a dictionary
961	structure
962	"""
963	dic = self.Mesh2triangList()
964	dic_mesh = self.Mesh2MeshList()
965	for element in dic_mesh.keys():
966	dic[element] = dic_mesh[element]
967	return dic
968
969	def setTriangulation(self, genDict):
970	"""
971	Set the mesh attributes given a dictionary of the lists
972	returned from the triang module
973	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
974	generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
975	generated point attribute title list:[A1Title, A2Title ...] (list of strings)
976	generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
977	generated segment marker list: [S1Tag, S2Tag, ...] (list of ints)
978	triangle list: [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
979	triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
980	triangle attribute list: [(T1att), (T2att), ...] (list of a list of strings)
981	"""
982	#Clear the current generated mesh values
983	self.meshTriangles=[]
984	self.attributeTitles=[]
985	self.meshSegments=[]
986	self.meshVertices=[]
987
988	#print "mesh.setTriangulation@#@#@#"
989	#print genDict
990	#print "@#@#@#"
991
992	self.maxVertexIndex = 0
993	for point in genDict['generatedpointlist']:
994	v=Vertex(point[0], point[1])
995	v.index = self.maxVertexIndex
996	self.maxVertexIndex +=1
997	self.meshVertices.append(v)
998
999	self.attributeTitles = genDict['generatedpointattributetitlelist']
1000
1001	index = 0
1002	for seg,marker in map(None,genDict['generatedsegmentlist'],genDict['generatedsegmentmarkerlist']):
1003	segObject = Segment( self.meshVertices[seg[0]],
1004	self.meshVertices[seg[1]], tag = marker )
1005	segObject.index = index
1006	index +=1
1007	self.meshSegments.append(segObject)
1008
1009	index = 0
1010	for triangle in genDict['generatedtrianglelist']:
1011	tObject =Triangle( self.meshVertices[triangle[0]],
1012	self.meshVertices[triangle[1]],
1013	self.meshVertices[triangle[2]] )
1014	tObject.index = index
1015	index +=1
1016	self.meshTriangles.append(tObject)
1017
1018	index = 0
1019	for att in genDict['generatedtriangleattributelist']:
1020	if att == []:
1021	self.meshTriangles[index].setAttribute("")
1022	else:
1023	# Note, is the first attribute always the region att?
1024	# haven't confirmed this
1025	#Peter - I think so (from manuel)
1026	#...the first such value is assumed to be a regional attribute...
1027	self.meshTriangles[index].setAttribute(att[0])
1028	index += 1
1029
1030	index = 0
1031	for att in genDict['generatedpointattributelist']:
1032	if att == None:
1033	self.meshVertices[index].setAttributes([])
1034	else:
1035	self.meshVertices[index].setAttributes(att)
1036	index += 1
1037
1038	index = 0
1039	for triangle in genDict['generatedtriangleneighborlist']:
1040	# Build a list of triangle object neighbors
1041	ObjectNeighbor = []
1042	for neighbor in triangle:
1043	if ( neighbor != -1):
1044	ObjectNeighbor.append(self.meshTriangles[neighbor])
1045	else:
1046	ObjectNeighbor.append(None)
1047	self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],ObjectNeighbor[1],ObjectNeighbor[2])
1048	index += 1
1049
1050
1051	def setMesh(self, genDict):
1052	"""
1053	Set the user Mesh attributes given a dictionary of the lists
1054	point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1055	point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1056	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1057	segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1058	region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1059	region attribute list: ["","reservoir",""] list of strings
1060	region max area list:[real, None, Real,...] list of None and reals
1061
1062	mesh is an instance of a mesh object
1063	"""
1064	#Clear the current user mesh values
1065	self.clearUserSegments()
1066	self.userVertices=[]
1067	self.Holes=[]
1068	self.Regions=[]
1069
1070	#print "mesh.setMesh@#@#@#"
1071	#print genDict
1072	#print "@#@#@#"
1073
1074	#index = 0
1075	for point in genDict['pointlist']:
1076	v=Vertex(point[0], point[1])
1077	#v.index = index
1078	#index +=1
1079	self.userVertices.append(v)
1080
1081	#index = 0
1082	for seg,tag in map(None,genDict['segmentlist'],genDict['segmenttaglist']):
1083	segObject = Segment( self.userVertices[seg[0]],
1084	self.userVertices[seg[1]], tag = tag )
1085	#segObject.index = index
1086	#index +=1
1087	self.userSegments.append(segObject)
1088
1089	# Remove the loading of attribute info.
1090	# Have attribute info added using least_squares in pyvolution
1091	# index = 0
1092	# for att in genDict['pointattributelist']:
1093	# if att == None:
1094	# self.userVertices[index].setAttributes([])
1095	# else:
1096	# self.userVertices[index].setAttributes(att)
1097	# index += 1
1098
1099	#index = 0
1100	for point in genDict['holelist']:
1101	h=Hole(point[0], point[1])
1102	#h.index = index
1103	#index +=1
1104	self.holes.append(h)
1105
1106	#index = 0
1107	for reg,att,maxArea in map(None,
1108	genDict['regionlist'],
1109	genDict['regionattributelist'],
1110	genDict['regionmaxarealist']):
1111	Object = Region( reg[0],
1112	reg[1],
1113	tag = att,
1114	maxArea = maxArea)
1115	#Object.index = index
1116	#index +=1
1117	self.regions.append(Object)
1118
1119	def addVertsSegs(self, outlineDict):
1120	"""
1121	Add out-line (user Mesh) attributes given a dictionary of the lists
1122	point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1123	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1124	#segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1125
1126	Assume the values are in Eastings and Northings, with no reference
1127	point
1128	"""
1129
1130	localUserVertices = []
1131	#index = 0
1132	for point in outlineDict['pointlist']:
1133
1134	v=Vertex(point[0]-self.geo_reference.xllcorner,
1135	point[1]-self.geo_reference.yllcorner)
1136	#v.index = index
1137	#index +=1
1138	self.userVertices.append(v)
1139	localUserVertices.append(v)
1140
1141	#index = 0
1142	for seg in map(None,outlineDict['segmentlist']):
1143	segObject = Segment( localUserVertices[seg[0]],
1144	localUserVertices[seg[1]] )
1145	#segObject.index = index
1146	#index +=1
1147	self.userSegments.append(segObject)
1148	#DSG!!!
1149
1150	def TestautoSegment(self):
1151	newsegs = []
1152	s1 = Segment(self.userVertices[0],
1153	self.userVertices[1])
1154	s2 = Segment(self.userVertices[0],
1155	self.userVertices[2])
1156	s3 = Segment(self.userVertices[2],
1157	self.userVertices[1])
1158	if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1159	newsegs.append(s1)
1160	if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1161	newsegs.append(s2)
1162	if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1163	newsegs.append(s3)
1164	#DSG!!!
1165	self.userSegments.extend(newsegs)
1166	return newsegs
1167
1168
1169	def savePickle(self, currentName):
1170	fd = open(currentName, 'w')
1171	pickle.dump(self,fd)
1172	fd.close()
1173
1174	def autoSegmentHull(self):
1175	"""
1176	initially work by running an executable
1177	Later compile the c code with a python wrapper.
1178
1179	Precon: There must be 3 or more vertices in the userVertices structure
1180	"""
1181	newsegs = []
1182	inputfile = 'hull_in.txt'
1183	outputfile = inputfile + '-alf'
1184	#write vertices to file
1185	fd = open(inputfile,'w')
1186	for v in self.userVertices:
1187	fd.write(str(v.x))
1188	fd.write(' ')
1189	fd.write(str(v.y))
1190	fd.write('\n')
1191	fd.close()
1192
1193	#run hull executable
1194	#warning need to compile hull for the current operating system
1195	command = 'hull.exe -A -i ' + inputfile
1196	os.system(command)
1197
1198	#read results into this object
1199	fd = open(outputfile)
1200	lines = fd.readlines()
1201	fd.close()
1202	#print "(((*("
1203	#print lines
1204	#print "(((*("
1205	lines.pop(0) #remove the first (title) line
1206	for line in lines:
1207	vertindexs = line.split()
1208	#print 'int(vertindexs[0])', int(vertindexs[0])
1209	#print 'int(vertindexs[1])', int(vertindexs[1])
1210	#print 'self.userVertices[int(vertindexs[0])]' ,self.userVertices[int(vertindexs[0])]
1211	#print 'self.userVertices[int(vertindexs[1])]' ,self.userVertices[int(vertindexs[1])]
1212	v1 = self.userVertices[int(vertindexs[0])]
1213	v2 = self.userVertices[int(vertindexs[1])]
1214
1215	if self.isUserSegmentNew(v1,v2):
1216	newseg = Segment(v1, v2)
1217	newsegs.append(newseg)
1218	#DSG!!!
1219	self.userSegments.extend(newsegs)
1220	return newsegs
1221	def autoSegmentFilter(self,):
1222	pass
1223
1224	def autoSegment(self, alpha = None,
1225	raw_boundary=True,
1226	remove_holes=False,
1227	smooth_indents=False,
1228	expand_pinch=False):
1229	"""
1230	Precon: There must be 3 or more vertices in the userVertices structure
1231	"""
1232
1233	points=[]
1234	for vertex in self.getUserVertices():
1235	points.append((vertex.x,vertex.y))
1236	self.shape = alpha_shape.alpha_shape.Alpha_Shape(points, alpha = alpha)
1237	self.shape.set_boundary_type(raw_boundary=raw_boundary,
1238	remove_holes=remove_holes,
1239	smooth_indents=smooth_indents,
1240	expand_pinch=expand_pinch)
1241	boundary_segs = self.shape.get_boundary()
1242
1243	segs2delete = self.alphaUserSegments
1244
1245	new_segs = []
1246	alpha_segs = []
1247	user_segs = []
1248	for seg in boundary_segs:
1249	v1 = self.userVertices[int(seg[0])]
1250	v2 = self.userVertices[int(seg[1])]
1251	alpha_seg = self.representedAlphaUserSegment(v1, v2)
1252	user_seg = self.representedUserSegment(v1, v2)
1253	#DSG!!!
1254	assert not(not (alpha_seg == None) and not (user_seg == None))
1255	if not alpha_seg == None:
1256	alpha_segs.append(alpha_seg)
1257	elif not user_seg == None:
1258	user_segs.append(user_seg)
1259	else:
1260	unique_seg = Segment(v1, v2)
1261	new_segs.append(unique_seg)
1262
1263	for seg in alpha_segs:
1264	try:
1265	segs2delete.remove(seg)
1266	except:
1267	pass
1268
1269	self.alphaUserSegments = []
1270	self.alphaUserSegments.extend(new_segs)
1271	self.alphaUserSegments.extend(alpha_segs)
1272
1273	optimum_alpha = self.shape.get_alpha()
1274	# need to draw newsegs
1275	return new_segs, segs2delete, optimum_alpha
1276
1277	def joinVertices(self):
1278	"""
1279	Return list of segments connecting all userVertices
1280	in the order they were given
1281
1282	Precon: There must be 3 or more vertices in the userVertices structure
1283	"""
1284
1285	newsegs = []
1286
1287	v1 = self.userVertices[0]
1288	for v2 in self.userVertices[1:]:
1289	if self.isUserSegmentNew(v1,v2):
1290	newseg = Segment(v1, v2)
1291	newsegs.append(newseg)
1292	v1 = v2
1293
1294	#Connect last point to the first
1295	v2 = self.userVertices[0]
1296	if self.isUserSegmentNew(v1,v2):
1297	newseg = Segment(v1, v2)
1298	newsegs.append(newseg)
1299
1300
1301	#Update list of user segments
1302	#DSG!!!
1303	self.userSegments.extend(newsegs)
1304	return newsegs
1305
1306	def normaliseMesh(self,scale, offset, height_scale):
1307	[xmin, ymin, xmax, ymax] = self.boxsize()
1308	[attmin0, attmax0] = self.maxMinVertAtt(0)
1309	#print "[attmin0, attmax0]" ,[attmin0, attmax0]
1310	[attmin1, attmax1] = self.maxMinVertAtt(1)
1311	#print [xmin, ymin, xmax, ymax]
1312	xrange = xmax - xmin
1313	yrange = ymax - ymin
1314	if xrange > yrange:
1315	min,max = xmin, xmax
1316	else:
1317	min,max = ymin, ymax
1318
1319	for obj in self.getUserVertices():
1320	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1321	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1322	if len(obj.attributes) > 0 and attmin0 != attmax0:
1323	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1324	(attmax0-attmin0)*height_scale
1325	if len(obj.attributes) > 1 and attmin1 != attmax1:
1326	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1327	(attmax1-attmin1)*height_scale
1328
1329	for obj in self.getMeshVertices():
1330	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1331	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1332	if len(obj.attributes) > 0 and attmin0 != attmax0:
1333	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1334	(attmax0-attmin0)*height_scale
1335	if len(obj.attributes) > 1 and attmin1 != attmax1:
1336	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1337	(attmax1-attmin1)*height_scale
1338
1339	for obj in self.getHoles():
1340	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1341	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1342	for obj in self.getRegions():
1343	obj.x = (obj.x - xmin)/(max- min)*scale + offset
1344	obj.y = (obj.y - ymin)/(max- min)*scale + offset
1345	[xmin, ymin, xmax, ymax] = self.boxsize()
1346	#print [xmin, ymin, xmax, ymax]
1347
1348	def boxsizeVerts(self):
1349	"""
1350	Returns a list of verts denoting a box or triangle that contains verts on the xmin, ymin, xmax and ymax axis.
1351	Structure: list of verts
1352	"""
1353	# FIXME dsg!!! large is a hack
1354	#You want the kinds package, part of Numeric:
1355	#In [2]: import kinds
1356
1357	#In [3]: kinds.default_float_kind.M
1358	#kinds.default_float_kind.MAX kinds.default_float_kind.MIN
1359	#kinds.default_float_kind.MAX_10_EXP kinds.default_float_kind.MIN_10_EXP
1360	#kinds.default_float_kind.MAX_EXP kinds.default_float_kind.MIN_EXP
1361
1362	#In [3]: kinds.default_float_kind.MIN
1363	#Out[3]: 2.2250738585072014e-308
1364
1365	large = 1e100
1366	xmin= large
1367	xmax=-large
1368	ymin= large
1369	ymax=-large
1370	for vertex in self.userVertices:
1371	if vertex.x < xmin:
1372	xmin = vertex.x
1373	xminVert = vertex
1374	if vertex.x > xmax:
1375	xmax = vertex.x
1376	xmaxVert = vertex
1377
1378	if vertex.y < ymin:
1379	ymin = vertex.y
1380	yminVert = vertex
1381	if vertex.y > ymax:
1382	ymax = vertex.y
1383	ymaxVert = vertex
1384	verts, count = self.removeDuplicatedVertices([xminVert,xmaxVert,yminVert,ymaxVert])
1385
1386	return verts
1387
1388	def boxsize(self):
1389	"""
1390	Returns a list denoting a box that contains the entire structure of vertices
1391	Structure: [xmin, ymin, xmax, ymax]
1392	"""
1393	# FIXME dsg!!! large is a hack
1394	#You want the kinds package, part of Numeric:
1395	#In [2]: import kinds
1396
1397	#In [3]: kinds.default_float_kind.M
1398	#kinds.default_float_kind.MAX kinds.default_float_kind.MIN
1399	#kinds.default_float_kind.MAX_10_EXP kinds.default_float_kind.MIN_10_EXP
1400	#kinds.default_float_kind.MAX_EXP kinds.default_float_kind.MIN_EXP
1401
1402	#In [3]: kinds.default_float_kind.MIN
1403	#Out[3]: 2.2250738585072014e-308
1404
1405	large = 1e100
1406	xmin= large
1407	xmax=-large
1408	ymin= large
1409	ymax=-large
1410	for vertex in self.userVertices:
1411	if vertex.x < xmin:
1412	xmin = vertex.x
1413	if vertex.x > xmax:
1414	xmax = vertex.x
1415
1416	if vertex.y < ymin:
1417	ymin = vertex.y
1418	if vertex.y > ymax:
1419	ymax = vertex.y
1420	return [xmin, ymin, xmax, ymax]
1421
1422	def maxMinVertAtt(self, iatt):
1423	"""
1424	Returns a list denoting a box that contains the entire structure of vertices
1425	Structure: [xmin, ymin, xmax, ymax]
1426	"""
1427	# FIXME dsg!!! large is a hack
1428	#You want the kinds package, part of Numeric:
1429	#In [2]: import kinds
1430
1431	#In [3]: kinds.default_float_kind.M
1432	#kinds.default_float_kind.MAX kinds.default_float_kind.MIN
1433	#kinds.default_float_kind.MAX_10_EXP kinds.default_float_kind.MIN_10_EXP
1434	#kinds.default_float_kind.MAX_EXP kinds.default_float_kind.MIN_EXP
1435
1436	#In [3]: kinds.default_float_kind.MIN
1437	#Out[3]: 2.2250738585072014e-308
1438
1439	large = 1e100
1440	min= large
1441	max=-large
1442	for vertex in self.userVertices:
1443	if len(vertex.attributes) > iatt:
1444	if vertex.attributes[iatt] < min:
1445	min = vertex.attributes[iatt]
1446	if vertex.attributes[iatt] > max:
1447	max = vertex.attributes[iatt]
1448	for vertex in self.meshVertices:
1449	if len(vertex.attributes) > iatt:
1450	if vertex.attributes[iatt] < min:
1451	min = vertex.attributes[iatt]
1452	if vertex.attributes[iatt] > max:
1453	max = vertex.attributes[iatt]
1454	return [min, max]
1455
1456	def scaleoffset(self, WIDTH, HEIGHT):
1457	"""
1458	Returns a list denoting the scale and offset terms that need to be
1459	applied when converting mesh co-ordinates onto grid co-ordinates
1460	Structure: [scale, xoffset, yoffset]
1461	"""
1462	OFFSET = 0.05*min([WIDTH, HEIGHT])
1463	[xmin, ymin, xmax, ymax] = self.boxsize()
1464	SCALE = min([0.9WIDTH, 0.9HEIGHT])/max([xmax-xmin, ymax-ymin])
1465
1466	if SCALE*xmin < OFFSET:
1467	xoffset = abs(SCALE*xmin) + OFFSET
1468	if SCALE*xmax > WIDTH - OFFSET:
1469	xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1470	if SCALE*ymin < OFFSET:
1471	b = abs(SCALE*ymin)+OFFSET
1472	if SCALE*ymax > HEIGHT-OFFSET:
1473	b = -(SCALE*ymax - HEIGHT + OFFSET)
1474	yoffset = HEIGHT - b
1475	return [SCALE, xoffset, yoffset]
1476
1477	def plotMeshTriangle(self,tag = 0,WIDTH = 400,HEIGHT = 400):
1478	"""
1479	Plots all node connections.
1480	tag = 0 (no node numbers), tag = 1 (node numbers)
1481	"""
1482
1483	try:
1484	from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
1485
1486	[SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
1487
1488	root = Tk()
1489	frame = Frame(root)
1490	frame.pack()
1491	button = Button(frame, text="OK", fg="red", command=frame.quit)
1492	button.pack(side=BOTTOM)
1493	canvas = Canvas(frame,bg="white", width=WIDTH, height=HEIGHT)
1494	canvas.pack()
1495	text = Label(frame, width=20, height=10, text='triangle mesh')
1496	text.pack()
1497
1498	#print self.meshTriangles
1499	for triangle in self.meshTriangles:
1500	triangle.draw(canvas,1,
1501	scale = SCALE,
1502	xoffset = xoffset,
1503	yoffset = yoffset )
1504
1505	root.mainloop()
1506
1507	except:
1508	print "Unexpected error:", sys.exc_info()[0]
1509	raise
1510
1511	#print """
1512	#node::plot Failed.
1513	#Most probably, the Tkinter module is not available.
1514	#"""
1515
1516	def plotUserSegments(self,tag = 0,WIDTH = 400,HEIGHT = 400):
1517	"""
1518	Plots all node connections.
1519	tag = 0 (no node numbers), tag = 1 (node numbers)
1520	"""
1521
1522	try:
1523	from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
1524
1525	[SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
1526
1527	root = Tk()
1528	frame = Frame(root)
1529	frame.pack()
1530	button = Button(frame, text="OK", fg="red", command=frame.quit)
1531	button.pack(side=BOTTOM)
1532	canvas = Canvas(frame, bg="white", width=WIDTH, height=HEIGHT)
1533	canvas.pack()
1534	text = Label(frame, width=20, height=10, text='user segments')
1535	text.pack()
1536
1537	for segment in self.getUserSegments():
1538	segment.draw(canvas,SCALE, xoffset, yoffset )
1539
1540	root.mainloop()
1541
1542	except:
1543	print "Unexpected error:", sys.exc_info()[0]
1544	raise
1545
1546	#print """
1547	#node::plot Failed.
1548	#Most probably, the Tkinter module is not available.
1549	#"""
1550	# FIXME let's not use this function,
1551	# use export _mesh_file instead
1552	def export_triangulation_file(self,ofile):
1553	"""
1554	export a file, ofile, with the format
1555
1556	First line: <# of vertices> <# of attributes>
1557	Following lines: <vertex #> <x> <y> [attributes]
1558	One line: <vertex attribute titles>
1559	One line: <# of triangles>
1560	Following lines: <triangle #> <vertex #> <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
1561	One line: <# of segments>
1562	Following lines: <segment #> <vertex #> <vertex #> [boundary tag]
1563	"""
1564	gen_dict = self.Mesh2IODict()
1565	if (ofile[-4:] == ".tsh"):
1566	fd = open(ofile,'w')
1567	load_mesh.loadASCII.write_ASCII_triangulation(fd,gen_dict)
1568	self.writeASCIImesh(fd,
1569	self.userVertices,
1570	self.getUserSegments(),
1571	self.holes,
1572	self.regions)
1573	fd.close()
1574	elif (ofile[-4:] == ".msh"):
1575	#print "mesh gen_dict",gen_dict
1576	load_mesh.loadASCII.write_msh_file(ofile, gen_dict)
1577
1578	# self.writeASCIImesh(fd, does anything use this?
1579
1580	#FIXME this function has a bug..
1581	def exportASCIIsegmentoutlinefile(self,ofile):
1582	"""
1583	export a file, ofile, with no triangulation and only vertices connected to segments.
1584	"""
1585	fd = open(ofile,'w')
1586	meshDict = {}
1587
1588	meshDict['vertices'] = []
1589	meshDict['vertex_attributes'] = []
1590	meshDict['segments'] = []
1591	meshDict['segment_tags'] = []
1592	meshDict['triangles'] = []
1593	meshDict['triangle_tags'] = []
1594	meshDict['triangle_neighbors'] = []
1595
1596	load_mesh.loadASCII.write_ASCII_triangulation(fd,meshDict)
1597	self.writeASCIIsegmentoutline(fd,
1598	self.userVertices,
1599	self.getUserSegments(),
1600	self.holes,
1601	self.regions)
1602	fd.close()
1603
1604	def exportASCIIobj(self,ofile):
1605	"""
1606	export a file, ofile, with the format
1607	lines: v <x> <y> <first attribute>
1608	f <vertex #> <vertex #> <vertex #> (of the triangles)
1609	"""
1610	fd = open(ofile,'w')
1611	self.writeASCIIobj(fd)
1612	fd.close()
1613
1614
1615	def writeASCIIobj(self,fd):
1616	fd.write(" # Triangulation as an obj file\n")
1617	numVert = str(len(self.meshVertices))
1618
1619	index1 = 1
1620	for vert in self.meshVertices:
1621	vert.index1 = index1
1622	index1 += 1
1623
1624	fd.write("v "
1625	+ str(vert.x) + " "
1626	+ str(vert.y) + " "
1627	+ str(vert.attributes[0]) + "\n")
1628
1629	for tri in self.meshTriangles:
1630	fd.write("f "
1631	+ str(tri.vertices[0].index1) + " "
1632	+ str(tri.vertices[1].index1) + " "
1633	+ str(tri.vertices[2].index1) + "\n")
1634
1635	#FIXME I think this has a bug...
1636	def writeASCIIsegmentoutline(self,
1637	fd,
1638	userVertices,
1639	userSegments,
1640	holes,
1641	regions):
1642	"""Write the user mesh info, only with vertices that are connected to segs
1643	"""
1644	verts = []
1645	#dupindex = 0
1646	for seg in self.userSegments:
1647	verts.append(seg.vertices[0])
1648	verts.append(seg.vertices[1])
1649	print 'verts>',verts
1650
1651	verts, count = self.removeDuplicatedVertices(verts)
1652	print 'verts no dups>',verts
1653	self.writeASCIImesh(fd,
1654	verts,
1655	self.getUserSegments(),
1656	self.holes,
1657	self.regions)
1658
1659	# exportASCIImeshfile
1660	def export_mesh_file(self,ofile):
1661	"""
1662	export a file, ofile, with the format
1663	"""
1664
1665	dict = self.Mesh2IODict()
1666	load_mesh.loadASCII.export_mesh_file(ofile,dict)
1667
1668	# not used
1669	def writeASCIImesh(self,
1670	fd,
1671	userVertices,
1672	userSegments,
1673	holes,
1674	regions):
1675
1676	#print "mesh.writeASCIImesh*********"
1677	#print "dict",dict
1678	#print "mesh*********"
1679	load_mesh.loadASCII.write_ASCII_outline(fd, dict)
1680
1681	#FIXME the title is wrong, need more comments
1682	def exportxyafile(self,ofile):
1683	"""
1684	export a file, ofile, with the format
1685
1686	First line: <# of vertices> <# of attributes>
1687	Following lines: <vertex #> <x> <y> [attributes]
1688	"""
1689	#load_mesh.loadASCII
1690	#FIXME, this should be a mesh2io method
1691	if self.meshVertices == []:
1692	Vertices = self.userVertices
1693	else:
1694	Vertices = self.meshVertices
1695
1696	numVert = str(len(Vertices))
1697
1698	if Vertices == []:
1699	raise RuntimeError
1700	numVertAttrib = str(len(Vertices[0].attributes))
1701	title = numVert + " " + numVertAttrib + " # <vertex #> <x> <y> [attributes]"
1702
1703	#Convert the Vertices to pointlist and pointattributelist
1704	xya_dict = {}
1705	pointattributes = []
1706	points = []
1707
1708	for vert in Vertices:
1709	points.append([vert.x,vert.y])
1710	pointattributes.append(vert.attributes)
1711
1712	xya_dict['pointlist'] = points
1713	xya_dict['pointattributelist'] = pointattributes
1714	xya_dict['geo_reference'] = self.geo_reference
1715
1716	load_mesh.loadASCII.export_xya_file(ofile, xya_dict, title, delimiter = " ")
1717
1718
1719	########### IO CONVERTERS ##################
1720	"""
1721	The dict fromat for IO with .tsh files is;
1722	(the triangulation)
1723	vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
1724	vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1725	vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
1726	segments: [[v1,v2],[v3,v4],...] (lists of integers)
1727	segment_tags : [tag,tag,...] list of strings
1728	triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
1729	triangle tags: [s1,s2,...] A list of strings
1730	triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
1731
1732	(the outline)
1733	points: [[x1,y1],[x2,y2],...] (lists of doubles)
1734	point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1735	outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
1736	outline_segment_tags : [tag1,tag2,...] list of strings
1737	holes : [[x1,y1],...](List of doubles, one inside each hole region)
1738	regions : [ [x1,y1],...] (List of 4 doubles)
1739	region_tags : [tag1,tag2,...] (list of strings)
1740	region_max_areas: [ma1,ma2,...] (A list of doubles)
1741	{Convension: A -ve max area means no max area}
1742
1743	"""
1744
1745
1746
1747	def Mesh2IODict(self):
1748	"""
1749	Convert the triangulation and outline info of a mesh to a dictionary
1750	structure
1751	"""
1752	dict = self.Mesh2IOTriangulationDict()
1753	dict_mesh = self.Mesh2IOOutlineDict()
1754	for element in dict_mesh.keys():
1755	dict[element] = dict_mesh[element]
1756
1757	# add the geo reference
1758	dict['geo_reference'] = self.geo_reference
1759	return dict
1760
1761	def Mesh2IOTriangulationDict(self):
1762	"""
1763	Convert the Mesh to a dictionary of lists describing the
1764	triangulation variables;
1765
1766	Used to produce .tsh file
1767	"""
1768
1769	meshDict = {}
1770	vertices=[]
1771	vertex_attributes=[]
1772
1773
1774	self.maxVertexIndex=0
1775	for vertex in self.meshVertices:
1776	vertex.index = self.maxVertexIndex
1777	vertices.append([vertex.x,vertex.y])
1778	vertex_attributes.append(vertex.attributes)
1779	self.maxVertexIndex += 1
1780
1781	meshDict['vertices'] = vertices
1782	meshDict['vertex_attributes'] = vertex_attributes
1783	meshDict['vertex_attribute_titles'] = self.attributeTitles
1784	#segments
1785	segments=[]
1786	segment_tags=[]
1787	for seg in self.meshSegments:
1788	segments.append([seg.vertices[0].index,seg.vertices[1].index])
1789	segment_tags.append(seg.tag)
1790	meshDict['segments'] =segments
1791	meshDict['segment_tags'] =segment_tags
1792
1793	# Make sure that the indexation is correct
1794	index = 0
1795	for tri in self.meshTriangles:
1796	tri.index = index
1797	index += 1
1798
1799	triangles = []
1800	triangle_tags = []
1801	triangle_neighbors = []
1802	for tri in self.meshTriangles:
1803	triangles.append([tri.vertices[0].index,tri.vertices[1].index,tri.vertices[2].index])
1804	triangle_tags.append([tri.attribute])
1805	neighborlist = [-1,-1,-1]
1806	for neighbor,index in map(None,tri.neighbors,
1807	range(len(tri.neighbors))):
1808	if neighbor:
1809	neighborlist[index] = neighbor.index
1810	triangle_neighbors.append(neighborlist)
1811
1812	meshDict['triangles'] = triangles
1813	meshDict['triangle_tags'] = triangle_tags
1814	meshDict['triangle_neighbors'] = triangle_neighbors
1815
1816	#print "mesh.Mesh2IOTriangulationDict))"
1817	#print meshDict
1818	#print "mesh.Mesh2IOTriangulationDict))"
1819
1820	return meshDict
1821
1822
1823	def Mesh2IOOutlineDict(self, userVertices=None,
1824	userSegments=None,
1825	holes=None,
1826	regions=None):
1827	"""
1828	Convert the mesh outline to a dictionary of the lists needed for the
1829	triang module;
1830
1831	Note, this adds an index attribute to the user Vertex objects.
1832
1833	Used to produce .tsh file and output to triangle
1834	"""
1835	if userVertices is None:
1836	userVertices = self.getUserVertices()
1837	if userSegments is None:
1838	userSegments = self.getUserSegments()
1839	if holes is None:
1840	holes = self.getHoles()
1841	if regions is None:
1842	regions = self.getRegions()
1843
1844	meshDict = {}
1845	#print "userVertices",userVertices
1846	#print "userSegments",userSegments
1847	pointlist=[]
1848	pointattributelist=[]
1849	index = 0
1850	for vertex in userVertices:
1851	vertex.index = index
1852	pointlist.append([vertex.x,vertex.y])
1853	pointattributelist.append(vertex.attributes)
1854
1855	index += 1
1856	meshDict['points'] = pointlist
1857	meshDict['point_attributes'] = pointattributelist
1858
1859	segmentlist=[]
1860	segmenttaglist=[]
1861	for seg in userSegments:
1862	segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
1863	segmenttaglist.append(seg.tag)
1864	meshDict['outline_segments'] =segmentlist
1865	meshDict['outline_segment_tags'] =segmenttaglist
1866
1867	holelist=[]
1868	for hole in holes:
1869	holelist.append([hole.x,hole.y])
1870	meshDict['holes'] = holelist
1871
1872	regionlist=[]
1873	regiontaglist = []
1874	regionmaxarealist = []
1875	for region in regions:
1876	regionlist.append([region.x,region.y])
1877	regiontaglist.append(region.getTag())
1878
1879	if (region.getMaxArea() != None):
1880	regionmaxarealist.append(region.getMaxArea())
1881	else:
1882	regionmaxarealist.append( load_mesh.loadASCII.NOMAXAREA)
1883	meshDict['regions'] = regionlist
1884	meshDict['region_tags'] = regiontaglist
1885	meshDict['region_max_areas'] = regionmaxarealist
1886	#print "(("
1887	#print meshDict
1888	#print meshDict['regionlist']
1889	#print "(("
1890	return meshDict
1891
1892
1893	def IOTriangulation2Mesh(self, genDict):
1894	"""
1895	Set the mesh attributes given an tsh IO dictionary
1896	"""
1897	#Clear the current generated mesh values
1898	self.meshTriangles=[]
1899	self.attributeTitles=[]
1900	self.meshSegments=[]
1901	self.meshVertices=[]
1902
1903	#print "mesh.setTriangulation@#@#@#"
1904	#print genDict
1905	#print "@#@#@#"
1906
1907	self.maxVertexIndex = 0
1908	for point in genDict['vertices']:
1909	v=Vertex(point[0], point[1])
1910	v.index = self.maxVertexIndex
1911	self.maxVertexIndex +=1
1912	self.meshVertices.append(v)
1913
1914	self.attributeTitles = genDict['vertex_attribute_titles']
1915
1916	index = 0
1917	for seg,tag in map(None,genDict['segments'],genDict['segment_tags']):
1918	segObject = Segment( self.meshVertices[seg[0]],
1919	self.meshVertices[seg[1]], tag = tag )
1920	segObject.index = index
1921	index +=1
1922	self.meshSegments.append(segObject)
1923
1924	index = 0
1925	for triangle in genDict['triangles']:
1926	tObject =Triangle( self.meshVertices[triangle[0]],
1927	self.meshVertices[triangle[1]],
1928	self.meshVertices[triangle[2]] )
1929	tObject.index = index
1930	index +=1
1931	self.meshTriangles.append(tObject)
1932
1933	index = 0
1934	for att in genDict['triangle_tags']:
1935	if att == []:
1936	self.meshTriangles[index].setAttribute("")
1937	else:
1938	# Note, is the first attribute always the region att?
1939	# haven't confirmed this
1940	#Peter - I think so (from manuel)
1941	#...the first such value is assumed to be a regional attribute...
1942	self.meshTriangles[index].setAttribute(att[0])
1943	index += 1
1944
1945	index = 0
1946	for att in genDict['vertex_attributes']:
1947	if att == None:
1948	self.meshVertices[index].setAttributes([])
1949	else:
1950	self.meshVertices[index].setAttributes(att)
1951	index += 1
1952
1953	index = 0
1954	for triangle in genDict['triangle_neighbors']:
1955	# Build a list of triangle object neighbors
1956	ObjectNeighbor = []
1957	for neighbor in triangle:
1958	if ( neighbor != -1):
1959	ObjectNeighbor.append(self.meshTriangles[neighbor])
1960	else:
1961	ObjectNeighbor.append(None)
1962	self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],ObjectNeighbor[1],ObjectNeighbor[2])
1963	index += 1
1964
1965
1966	def IOOutline2Mesh(self, genDict):
1967	"""
1968	Set the outline (user Mesh attributes) given a IO tsh dictionary
1969
1970	mesh is an instance of a mesh object
1971	"""
1972	#Clear the current user mesh values
1973	self.clearUserSegments()
1974	self.userVertices=[]
1975	self.Holes=[]
1976	self.Regions=[]
1977
1978	#print "mesh.IOOutline2Mesh@#@#@#"
1979	#print "genDict",genDict
1980	#print "@#@#@#"
1981
1982	#index = 0
1983	for point in genDict['points']:
1984	v=Vertex(point[0], point[1])
1985	#v.index = index
1986	#index +=1
1987	self.userVertices.append(v)
1988
1989	#index = 0
1990	for seg,tag in map(None,genDict['outline_segments'],genDict['outline_segment_tags']):
1991	segObject = Segment( self.userVertices[seg[0]],
1992	self.userVertices[seg[1]], tag = tag )
1993	#segObject.index = index
1994	#index +=1
1995	self.userSegments.append(segObject)
1996
1997	# Remove the loading of attribute info.
1998	# Have attribute info added using least_squares in pyvolution
1999	# index = 0
2000	# for att in genDict['point_attributes']:
2001	# if att == None:
2002	# self.userVertices[index].setAttributes([])
2003	# else:
2004	# self.userVertices[index].setAttributes(att)
2005	# index += 1
2006
2007	#index = 0
2008	for point in genDict['holes']:
2009	h=Hole(point[0], point[1])
2010	#h.index = index
2011	#index +=1
2012	self.holes.append(h)
2013
2014	#index = 0
2015	for reg,att,maxArea in map(None,
2016	genDict['regions'],
2017	genDict['region_tags'],
2018	genDict['region_max_areas']):
2019	if maxArea > 0: # maybe I should ref NOMAXAREA? Prob' not though
2020	Object = Region( reg[0],
2021	reg[1],
2022	tag = att,
2023	maxArea = maxArea)
2024	else:
2025	Object = Region( reg[0],
2026	reg[1],
2027	tag = att)
2028
2029	#Object.index = index
2030	#index +=1
2031	self.regions.append(Object)
2032
2033	############################################
2034
2035
2036	def refineSet(self,setName):
2037	Triangles = self.sets[self.setID[setName]]
2038	Refine(self,Triangles)
2039
2040	def selectAllTriangles(self):
2041	A=[]
2042	A.extend(self.meshTriangles)
2043	if not('All' in self.setID.keys()):
2044	self.setID['All']=len(self.sets)
2045	self.sets.append(A)
2046	else:
2047	self.sets[self.setID['All']]=A
2048	return 'All'
2049	# and objectIDs
2050
2051
2052	def clearSelection(self):
2053	A = []
2054	if not('None' in self.setID.keys()):
2055	self.setID['None']=len(self.sets)
2056	self.sets.append(A)
2057	return 'None'
2058
2059	def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
2060	#FIXME Draws over previous triangles - may bloat canvas
2061	Triangles = self.sets[self.setID[setName]]
2062	for triangle in Triangles:
2063	triangle.draw(canvas,1,
2064	scale = SCALE,
2065	colour = colour)
2066
2067	def undrawSet(self,canvas,setName,SCALE,colour='green'):
2068	#FIXME Draws over previous lines - may bloat canvas
2069	Triangles = self.sets[self.setID[setName]]
2070	for triangle in Triangles:
2071	triangle.draw(canvas,1,
2072	scale = SCALE,
2073	colour = colour)
2074
2075
2076	def triangles_to_polySet(self,setName):
2077
2078	Triangles = self.sets[self.setID[setName]]
2079	Triangles_dict = {}
2080	for triangle in Triangles:
2081	Triangles_dict[triangle]=None
2082
2083	userVertices = {}
2084	userSegments = []
2085	index = 0
2086	for triangle in Triangles:
2087	for i in (0,1,2):
2088	if not Triangles_dict.has_key(triangle.neighbors[i]):
2089	a = triangle.vertices[i-1]
2090	b = triangle.vertices[i-2]
2091	userVertices[a]=0
2092	userVertices[b]=0
2093	userSegments.append([a,b])
2094	return userVertices, userSegments
2095
2096	def threshold(self,setName,min=None,max=None,attribute_name = 'elevation'):
2097	"""
2098	threshold using d
2099	"""
2100	triangles = self.sets[self.setID[setName]]
2101	A = []
2102
2103	if attribute_name in self.attributeTitles:
2104	i = self.attributeTitles.index(attribute_name)
2105	if not max == None:
2106	for t in triangles:
2107	if (min<self.av_att(t,i)<max):
2108	A.append(t)
2109	else:
2110	for t in triangles:
2111	if (min<self.av_att(t,i)):
2112	A.append(t)
2113	self.sets[self.setID[setName]] = A
2114
2115	def general_threshold(self,setName,min=None,max=None,attribute_name = 'elevation',function = None):
2116	"""
2117	threshold using d
2118	"""
2119	triangles = self.sets[self.setID[setName]]
2120	A = []
2121
2122	if attribute_name in self.attributeTitles:
2123	i = self.attributeTitles.index(attribute_name)
2124	if not max == None:
2125	for t in triangles:
2126	if (min<function(t,i)<max):
2127	A.append(t)
2128	else:
2129	for t in triangles:
2130	if (min<self.function(t,i)):
2131	A.append(t)
2132	self.sets[self.setID[setName]] = A
2133
2134	def av_att(self,triangle,i):
2135	#evaluates the average attribute of the vertices of a triangle.
2136	V = triangle.getVertices()
2137	a0 = (V[0].attributes[i])
2138	a1 = (V[1].attributes[i])
2139	a2 = (V[2].attributes[i])
2140	return (a0+a1+a2)/3
2141
2142	def Courant_ratio(self,triangle,index):
2143	"""
2144	Not the true Courant ratio, just elevation on area
2145	"""
2146	e = self.av_att(triangle,index)
2147	A = triangle.calcArea()
2148	return e/A
2149
2150	def Gradient(self,triangle,index):
2151	V = triangle.vertices
2152	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]
2153	grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
2154	if ((grad_x2)+(grad_y2))**(0.5)<0:
2155	print ((grad_x2)+(grad_y2))**(0.5)
2156	return ((grad_x2)+(grad_y2))**(0.5)
2157
2158
2159	def append_triangle(self,triangle):
2160	self.meshTriangles.append(triangle)
2161
2162	def replace_triangle(self,triangle,replacement):
2163	i = self.meshTriangles.index(triangle)
2164	self.meshTriangles[i]=replacement
2165	assert replacement in self.meshTriangles
2166
2167
2168
2169
2170	def importUngenerateFile(ofile):
2171	"""
2172	import a file, ofile, with the format
2173	[poly]
2174	poly format:
2175	First line: <# of vertices> <x centroid> <y centroid>
2176	Following lines: <x> <y>
2177	last line: "END"
2178
2179	Note: These are clockwise.
2180	"""
2181	fd = open(ofile,'r')
2182	Dict = readUngenerateFile(fd)
2183	fd.close()
2184	return Dict
2185
2186	def readUngenerateFile(fd):
2187	"""
2188	import a file, ofile, with the format
2189	[poly]
2190	poly format:
2191	First line: <# of polynomial> <x centroid> <y centroid>
2192	Following lines: <x> <y>
2193	last line: "END"
2194	"""
2195	END_DELIMITER = 'END\n'
2196
2197	points = []
2198	segments = []
2199
2200	isEnd = False
2201	line = fd.readline() #not used <# of polynomial> <x> <y>
2202	while not isEnd:
2203	line = fd.readline()
2204	fragments = line.split()
2205	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2206	points.append(vert)
2207	PreviousVertIndex = len(points)-1
2208	firstVertIndex = PreviousVertIndex
2209
2210	line = fd.readline() #Read the next line
2211	while line <> END_DELIMITER:
2212	#print "line >" + line + "<"
2213	fragments = line.split()
2214	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2215	points.append(vert)
2216	thisVertIndex = len(points)-1
2217	segment = [PreviousVertIndex,thisVertIndex]
2218	segments.append(segment)
2219	PreviousVertIndex = thisVertIndex
2220	line = fd.readline() #Read the next line
2221	i =+ 1
2222	# If the last and first segments are the same,
2223	# Remove the last segment and the last vertex
2224	# then add a segment from the second last vert to the 1st vert
2225	thisVertIndex = len(points)-1
2226	firstVert = points[firstVertIndex]
2227	thisVert = points[thisVertIndex]
2228	#print "firstVert",firstVert
2229	#print "thisVert",thisVert
2230	if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
2231	points.pop()
2232	segments.pop()
2233	thisVertIndex = len(points)-1
2234	segments.append([thisVertIndex, firstVertIndex])
2235
2236	line = fd.readline() # read <# of polynomial> <x> <y> OR END
2237	#print "line >>" + line + "<<"
2238	if line == END_DELIMITER:
2239	isEnd = True
2240
2241	#print "points", points
2242	#print "segments", segments
2243	ungenerated_dict = {}
2244	ungenerated_dict['pointlist'] = points
2245	ungenerated_dict['segmentlist'] = segments
2246	return ungenerated_dict
2247
2248	def importMeshFromFile(ofile):
2249	"""returns a mesh object, made from a .xya/.pts or .tsh/.msh file
2250	Often raises SyntaxError, IOError
2251	"""
2252	newmesh = None
2253	if ofile[-4:]== ".xya":
2254	#print "loading " + ofile
2255	try:
2256	dict = load_mesh.loadASCII.load_points_file(ofile, delimiter = ',')
2257	except SyntaxError:
2258	dict = load_mesh.loadASCII.load_points_file(ofile, delimiter = ' ')
2259	#print "dict",dict outline_
2260	dict['segmentlist'] = []
2261	dict['segmenttaglist'] = []
2262	dict['regionlist'] = []
2263	dict['regionattributelist'] = []
2264	dict['regionmaxarealist'] = []
2265	dict['holelist'] = []
2266	newmesh= Mesh()
2267	newmesh.setMesh(dict) #FIXME use IOOutline2Mesh
2268	counter = newmesh.removeDuplicatedUserVertices()
2269	if (counter >0):
2270	print "%i duplicate vertices removed from dataset" % (counter)
2271	elif ofile[-4:]== ".pts":
2272	#print "loading " + ofile
2273	dict = load_mesh.loadASCII.load_points_file(ofile)
2274	#print "dict",dict
2275	dict['points'] = dict['pointlist']
2276	dict['outline_segments'] = []
2277	dict['outline_segment_tags'] = []
2278	dict['regions'] = []
2279	dict['region_tags'] = []
2280	dict['region_max_areas'] = []
2281	dict['holes'] = []
2282	newmesh= Mesh(geo_reference = dict['geo_reference'])
2283	newmesh.IOOutline2Mesh(dict) #FIXME use IOOutline2Mesh
2284	counter = newmesh.removeDuplicatedUserVertices()
2285	if (counter >0):
2286	print "%i duplicate vertices removed from dataset" % (counter)
2287	elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
2288	dict = load_mesh.loadASCII.import_triangulation(ofile)
2289	#print "********"
2290	#print "zq mesh.dict",dict
2291	#print "********"
2292	newmesh= Mesh()
2293	newmesh.IOOutline2Mesh(dict)
2294	newmesh.IOTriangulation2Mesh(dict)
2295	else:
2296	raise RuntimeError
2297
2298	if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
2299	newmesh.geo_reference = dict['geo_reference']
2300	return newmesh
2301
2302	def loadPickle(currentName):
2303	fd = open(currentName)
2304	mesh = pickle.load(fd)
2305	fd.close()
2306	return mesh
2307
2308	def square_outline(side_length = 1,up = "top", left = "left", right = "right",
2309	down = "bottom", regions = False):
2310
2311	a = Vertex (0,0)
2312	b = Vertex (0,side_length)
2313	c = Vertex (side_length,0)
2314	d = Vertex (side_length,side_length)
2315
2316	s2 = Segment(b,d, tag = up)
2317	s3 = Segment(b,a, tag = left)
2318	s4 = Segment(d,c, tag = right)
2319	s5 = Segment(a,c, tag = down)
2320
2321	if regions:
2322	e = Vertex (side_length/2,side_length/2)
2323	s6 = Segment(a,e, tag = down + left)
2324	s7 = Segment(b,e, tag = up + left)
2325	s8 = Segment(c,e, tag = down + right)
2326	s9 = Segment(d,e, tag = up + right)
2327	r1 = Region(side_length/2,3.*side_length/4, tag = up)
2328	r2 = Region(1.*side_length/4,side_length/2, tag = left)
2329	r3 = Region(3.*side_length/4,side_length/2, tag = right)
2330	r4 = Region(side_length/2,1.*side_length/4, tag = down)
2331	mesh = Mesh(userVertices=[a,b,c,d,e],
2332	userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
2333	regions = [r1,r2,r3,r4])
2334	else:
2335	mesh = Mesh(userVertices=[a,b,c,d],
2336	userSegments=[s2,s3,s4,s5])
2337
2338	return mesh
2339
2340
2341
2342	def region_strings2ints(region_list):
2343	"""Given a list of (x_int,y_int,tag_string) lists it returns a list of
2344	(x_int,y_int,tag_int) and a list to convert the tag_int's back to
2345	the tag_strings
2346	"""
2347	# Make sure "" has an index of 0
2348	region_list.reverse()
2349	region_list.append((1.0,2.0,""))
2350	region_list.reverse()
2351	convertint2string = []
2352	for i in xrange(len(region_list)):
2353	convertint2string.append(region_list[i][2])
2354	if len(region_list[i]) == 4: # there's an area value
2355	region_list[i] = (region_list[i][0],
2356	region_list[i][1],i,region_list[i][3])
2357	elif len(region_list[i]) == 3: # no area value
2358	region_list[i] = (region_list[i][0],region_list[i][1],i)
2359	else:
2360	print "The region list has a bad size"
2361	# raise an error ..
2362	raise Error
2363
2364	#remove "" from the region_list
2365	region_list.pop(0)
2366
2367	return [region_list, convertint2string]
2368
2369	def region_ints2strings(region_list,convertint2string):
2370	"""Reverses the transformation of region_strings2ints
2371	"""
2372	if region_list[0] != []:
2373	for i in xrange(len(region_list)):
2374	region_list[i] = [convertint2string[int(region_list[i][0])]]
2375	return region_list
2376
2377	def segment_ints2strings(intlist, convertint2string):
2378	"""Reverses the transformation of segment_strings2ints """
2379	stringlist = []
2380	for x in intlist:
2381	stringlist.append(convertint2string[x])
2382	return stringlist
2383
2384	def segment_strings2ints(stringlist, preset):
2385	"""Given a list of strings return a list of 0 to n ints which represent
2386	the strings and a converting list of the strings, indexed by 0 to n ints.
2387	Also, input an initial converting list of the strings
2388	Note, the converting list starts off with
2389	["internal boundary", "external boundary", "internal boundary"]
2390	example input and output
2391	input ["a","b","a","c"],["c"]
2392	output [[2, 1, 2, 0], ['c', 'b', 'a']]
2393
2394	the first element in the converting list will be
2395	overwritten with "".
2396	?This will become the third item in the converting list?
2397
2398	# note, order the initial converting list is important,
2399	since the index = the int tag
2400	"""
2401	nodups = unique(stringlist)
2402	# note, order is important, the index = the int tag
2403	#preset = ["internal boundary", "external boundary"]
2404	#Remove the preset tags from the list with no duplicates
2405	nodups = [x for x in nodups if x not in preset]
2406
2407	try:
2408	nodups.remove("") # this has to go to zero
2409	except ValueError:
2410	pass
2411
2412	# Add the preset list at the beginning of no duplicates
2413	preset.reverse()
2414	nodups.extend(preset)
2415	nodups.reverse()
2416
2417
2418	convertstring2int = {}
2419	convertint2string = []
2420	index = 0
2421	for x in nodups:
2422	convertstring2int[x] = index
2423	convertint2string.append(x)
2424	index += 1
2425	convertstring2int[""] = 0
2426
2427	intlist = []
2428	for x in stringlist:
2429	intlist.append(convertstring2int[x])
2430	return [intlist, convertint2string]
2431
2432
2433
2434
2435	def unique(s):
2436	"""Return a list of the elements in s, but without duplicates.
2437
2438	For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
2439	unique("abcabc") some permutation of ["a", "b", "c"], and
2440	unique(([1, 2], [2, 3], [1, 2])) some permutation of
2441	[[2, 3], [1, 2]].
2442
2443	For best speed, all sequence elements should be hashable. Then
2444	unique() will usually work in linear time.
2445
2446	If not possible, the sequence elements should enjoy a total
2447	ordering, and if list(s).sort() doesn't raise TypeError it's
2448	assumed that they do enjoy a total ordering. Then unique() will
2449	usually work in O(N*log2(N)) time.
2450
2451	If that's not possible either, the sequence elements must support
2452	equality-testing. Then unique() will usually work in quadratic
2453	time.
2454	"""
2455
2456	n = len(s)
2457	if n == 0:
2458	return []
2459
2460	# Try using a dict first, as that's the fastest and will usually
2461	# work. If it doesn't work, it will usually fail quickly, so it
2462	# usually doesn't cost much to try it. It requires that all the
2463	# sequence elements be hashable, and support equality comparison.
2464	u = {}
2465	try:
2466	for x in s:
2467	u[x] = 1
2468	except TypeError:
2469	del u # move on to the next method
2470	else:
2471	return u.keys()
2472
2473	# We can't hash all the elements. Second fastest is to sort,
2474	# which brings the equal elements together; then duplicates are
2475	# easy to weed out in a single pass.
2476	# NOTE: Python's list.sort() was designed to be efficient in the
2477	# presence of many duplicate elements. This isn't true of all
2478	# sort functions in all languages or libraries, so this approach
2479	# is more effective in Python than it may be elsewhere.
2480	try:
2481	t = list(s)
2482	t.sort()
2483	except TypeError:
2484	del t # move on to the next method
2485	else:
2486	assert n > 0
2487	last = t[0]
2488	lasti = i = 1
2489	while i < n:
2490	if t[i] != last:
2491	t[lasti] = last = t[i]
2492	lasti += 1
2493	i += 1
2494	return t[:lasti]
2495
2496	# Brute force is all that's left.
2497	u = []
2498	for x in s:
2499	if x not in u:
2500	u.append(x)
2501	return u
2502
2503	"""Refines triangles
2504
2505	Implements the #triangular bisection?# algorithm.
2506
2507
2508	"""
2509
2510	def Refine(mesh, triangles):
2511	"""
2512	Given a general_mesh, and a triangle number, split
2513	that triangle in the mesh in half. Then to prevent
2514	vertices and edges from meeting, keep refining
2515	neighbouring triangles until the mesh is clean.
2516	"""
2517	state = BisectionState(mesh)
2518	for triangle in triangles:
2519	if not state.refined_triangles.has_key(triangle):
2520	triangle.rotate_longest_side()
2521	state.start(triangle)
2522	Refine_mesh(mesh, state)
2523
2524	def Refine_mesh(mesh, state):
2525	"""
2526	"""
2527	state.getState(mesh)
2528	refine_triangle(mesh,state)
2529	state.evolve()
2530	if not state.end:
2531	Refine_mesh(mesh,state)
2532
2533	def refine_triangle(mesh,state):
2534	split(mesh,state.current_triangle,state.new_point)
2535	if state.case == 'one':
2536	state.r[3]=state.current_triangle#triangle 2
2537
2538	new_triangle_id = len(mesh.meshTriangles)-1
2539	new_triangle = mesh.meshTriangles[new_triangle_id]
2540
2541	split(mesh,new_triangle,state.old_point)
2542	state.r[2]=new_triangle#triangle 1.2
2543	state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
2544	r = state.r
2545	state.repairCaseOne()
2546
2547	if state.case == 'two':
2548	state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2549
2550	new_triangle = state.current_triangle
2551
2552	split(mesh,new_triangle,state.old_point)
2553
2554	state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
2555	state.r[4]=new_triangle#triangle 2.2
2556	r = state.r
2557
2558	state.repairCaseTwo()
2559
2560	if state.case == 'vertex':
2561	state.r[2]=state.current_triangle#triangle 2
2562	state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2563	r = state.r
2564	state.repairCaseVertex()
2565
2566	if state.case == 'start':
2567	state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2568	state.r[3]=state.current_triangle#triangle 2
2569
2570	if state.next_case == 'boundary':
2571	state.repairCaseBoundary()
2572
2573
2574	def split(mesh, triangle, new_point):
2575	"""
2576	Given a mesh, triangle_id and a new point,
2577	split the corrosponding triangle into two
2578	new triangles and update the mesh.
2579	"""
2580
2581	new_triangle1 = Triangle(new_point,triangle.vertices[0],triangle.vertices[1],attribute = triangle.attribute, neighbors = None)
2582	new_triangle2 = Triangle(new_point,triangle.vertices[2],triangle.vertices[0],attribute = triangle.attribute, neighbors = None)
2583
2584	new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
2585	new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
2586
2587	mesh.meshTriangles.append(new_triangle1)
2588
2589	triangle.vertices = new_triangle2.vertices
2590	triangle.neighbors = new_triangle2.neighbors
2591
2592
2593	class State:
2594
2595	def __init__(self):
2596	pass
2597
2598	class BisectionState(State):
2599
2600
2601	def __init__(self,mesh):
2602	self.len = len(mesh.meshTriangles)
2603	self.refined_triangles = {}
2604	self.mesh = mesh
2605	self.current_triangle = None
2606	self.case = 'start'
2607	self.end = False
2608	self.r = [None,None,None,None,None]
2609
2610	def start(self, triangle):
2611	self.current_triangle = triangle
2612	self.case = 'start'
2613	self.end = False
2614	self.r = [None,None,None,None,None]
2615
2616	def getState(self,mesh):
2617	if not self.case == 'vertex':
2618	self.new_point=self.getNewVertex(mesh, self.current_triangle)
2619	#self.neighbour=self.getNeighbour(mesh, self.current_triangle)
2620	self.neighbour = self.current_triangle.neighbors[0]
2621	if not self.neighbour is None:
2622	self.neighbour.rotate_longest_side()
2623	self.next_case = self.get_next_case(mesh,self.neighbour,self.current_triangle)
2624	if self.case == 'vertex':
2625	self.new_point=self.old_point
2626
2627
2628	def evolve(self):
2629	if self.case == 'vertex':
2630	self.end = True
2631
2632	self.last_case = self.case
2633	self.case = self.next_case
2634
2635	self.old_point = self.new_point
2636	self.current_triangle = self.neighbour
2637
2638	if self.case == 'boundary':
2639	self.end = True
2640	self.refined_triangles[self.r[2]]=1
2641	self.refined_triangles[self.r[3]]=1
2642	if not self.r[4] is None:
2643	self.refined_triangles[self.r[4]]=1
2644	self.r[0]=self.r[2]
2645	self.r[1]=self.r[3]
2646
2647
2648	def getNewVertex(self,mesh,triangle):
2649	coordinate1 = triangle.vertices[1]
2650	coordinate2 = triangle.vertices[2]
2651	a = ([coordinate1.x1.,coordinate1.y1.])
2652	b = ([coordinate2.x1.,coordinate2.y1.])
2653	attributes = []
2654	for i in range(len(coordinate1.attributes)):
2655	att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
2656	attributes.append(att)
2657	new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
2658	newVertex = Vertex(new_coordinate[0],new_coordinate[1], attributes = attributes)
2659	mesh.maxVertexIndex+=1
2660	newVertex.index = mesh.maxVertexIndex
2661	mesh.meshVertices.append(newVertex)
2662	return newVertex
2663
2664	def get_next_case(self, mesh,neighbour,triangle):
2665	"""
2666	Given the locations of two neighbouring triangles,
2667	examine the interior indices of their vertices (i.e.
2668	0,1 or 2) to determine what how the neighbour needs
2669	to be refined.
2670	"""
2671	if (neighbour is None):
2672	next_case = 'boundary'
2673	else:
2674	if triangle.vertices[1].x==neighbour.vertices[2].x:
2675	if triangle.vertices[1].y==neighbour.vertices[2].y:
2676	next_case = 'vertex'
2677	if triangle.vertices[1].x==neighbour.vertices[0].x:
2678	if triangle.vertices[1].y==neighbour.vertices[0].y:
2679	next_case = 'two'
2680	if triangle.vertices[1].x==neighbour.vertices[1].x:
2681	if triangle.vertices[1].y==neighbour.vertices[1].y:
2682	next_case = 'one'
2683	return next_case
2684
2685
2686
2687	def repairCaseVertex(self):
2688
2689	r = self.r
2690
2691
2692	self.link(r[0],r[2])
2693	self.repair(r[0])
2694
2695	self.link(r[1],r[3])
2696	self.repair(r[1])
2697
2698	self.repair(r[2])
2699
2700	self.repair(r[3])
2701
2702
2703	def repairCaseOne(self):
2704	r = self.r
2705
2706
2707	self.link(r[0],r[2])
2708	self.repair(r[0])
2709
2710	self.link(r[1],r[4])
2711	self.repair(r[1])
2712
2713	self.repair(r[4])
2714
2715	def repairCaseTwo(self):
2716	r = self.r
2717
2718	self.link(r[0],r[4])
2719	self.repair(r[0])
2720
2721	self.link(r[1],r[3])
2722	self.repair(r[1])
2723
2724	self.repair(r[4])
2725
2726	def repairCaseBoundary(self):
2727	r = self.r
2728	self.repair(r[2])
2729	self.repair(r[3])
2730
2731
2732
2733	def repair(self,triangle):
2734	"""
2735	Given a triangle that knows its neighbours, this will
2736	force the neighbours to comply.
2737
2738	However, it needs to compare the vertices of triangles
2739	for this implementation
2740
2741	But it doesn't work for invalid neighbour structures
2742	"""
2743	n=triangle.neighbors
2744	for i in (0,1,2):
2745	if not n[i] is None:
2746	for j in (0,1,2):#to find which side of the list is broken
2747	if not (n[i].vertices[j] in triangle.vertices):
2748	#ie if j is the side of n that needs fixing
2749	k = j
2750	n[i].neighbors[k]=triangle
2751
2752
2753
2754	def link(self,triangle1,triangle2):
2755	"""
2756	make triangle1 neighbors point to t
2757	#count = 0riangle2
2758	"""
2759	count = 0
2760	for i in (0,1,2):#to find which side of the list is broken
2761	if not (triangle1.vertices[i] in triangle2.vertices):
2762	j = i
2763	count+=1
2764	assert count == 1
2765	triangle1.neighbors[j]=triangle2
2766
2767
2768
2769	if __name__ == "__main__":
2770	#from mesh import *
2771	# THIS CAN BE DELETED
2772	m = Mesh()
2773	dict = importUngenerateFile("ungen_test.txt")
2774	m.addVertsSegs(dict)
2775	print m

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

Download in other formats: