Context Navigation

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

Last change on this file since 1177 was 1177, checked in by prow, 20 years ago
nearly fixed issues with duplicate segments/vertices in outline draw.
File size: 100.9 KB

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

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

Download in other formats: