Context Navigation

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

Last change on this file since 1023 was 1023, checked in by duncan, 20 years ago
adding more geo-ref stuff
File size: 95.0 KB

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

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

Download in other formats: