Context Navigation

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

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

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

Download in other formats: