Context Navigation

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

Last change on this file since 1221 was 1221, checked in by prow, 20 years ago
fixed. point_on_line needed an allclose.
File size: 119.4 KB

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

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

Download in other formats: