Context Navigation

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

Last change on this file since 1181 was 1181, checked in by duncan, 20 years ago
geo ref stuff
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,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[0])
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
2120	def weed(self,Vertices,Segments):
2121	#weed out existing duplicates
2122	print 'len(self.getUserSegments())'
2123	print len(self.getUserSegments())
2124	print 'len(self.getUserVertices())'
2125	print len(self.getUserVertices())
2126
2127	point_keys = {}
2128	for vertex in Vertices:
2129	point = (vertex.x,vertex.y)
2130	point_keys[point]=vertex
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 triangles_to_polySet(self,setName):
2154	#turn the triangles into a set
2155	Triangles = self.sets[self.setID[setName]]
2156	Triangles_dict = {}
2157	for triangle in Triangles:
2158	Triangles_dict[triangle]=None
2159
2160	userVertices = {}
2161	userSegments = []
2162	point_keys = {}
2163	for vertex in self.getUserVertices():
2164	point = (vertex.x,vertex.y)
2165	point_keys[point]=vertex
2166
2167	line_keys = {}
2168	for segment in self.getUserSegments():
2169	vertex1 = segment.vertices[0]
2170	vertex2 = segment.vertices[1]
2171	point1 = (vertex1.x,vertex1.y)
2172	point2 = (vertex2.x,vertex2.y)
2173	#segment.vertices[0]=point_keys[point1]
2174	#segment.vertices[1]=point_keys[point2]
2175	#vertex1 = segment.vertices[0]
2176	#vertex2 = segment.vertices[1]
2177	#point1 = (vertex1.x,vertex1.y)
2178	#point2 = (vertex2.x,vertex2.y)
2179	line1 = (point1,point2)
2180	line2 = (point2,point1)
2181	if not (line_keys.has_key(line1) \
2182	or line_keys.has_key(line2)):
2183	line_keys[line1]=segment
2184
2185	for triangle in Triangles:
2186	for i in (0,1,2):
2187	#for every triangles neighbour:
2188
2189	if not Triangles_dict.has_key(triangle.neighbors[i]):
2190	#if the neighbour is not in the set:
2191	a = triangle.vertices[i-1]
2192	b = triangle.vertices[i-2]
2193	if not point_keys.has_key((a.x,a.y)):
2194	#if point a does not already exist
2195	#then add it to the points.
2196	userVertices[a]=True
2197	point_keys[(a.x,a.y)]=a
2198	else:
2199	a=point_keys[(a.x,a.y)]
2200	userVertices[a]=point_keys[(a.x,a.y)]
2201	assert userVertices.has_key(a)
2202
2203	if not point_keys.has_key((b.x,b.y)):
2204	#if point b does not already exist
2205	#then add it to the points.
2206	userVertices[b]=True
2207	point_keys[(b.x,b.y)]=b
2208	else:
2209	b=point_keys[(b.x,b.y)]
2210	userVertices[b]=point_keys[(b.x,b.y)]
2211	assert userVertices.has_key(b)
2212
2213	if not (line_keys.has_key(((a.x,a.y),(b.x,b.y)))\
2214	or line_keys.has_key(((b.x,b.y),(a.x,a.y)))):
2215	#if the segment does not already exist then
2216	#add it to the segments
2217	assert ((a.x,a.y)!=(b.x,b.y))
2218	assert a!=b
2219	userSegments.append(Segment(a,b))
2220	line_keys[((a.x,a.y),(b.x,b.y))]=None
2221
2222	userVertices,userSegments = self.weed(userVertices.keys(),userSegments)
2223	self.userVertices.extend(userVertices)
2224	self.userSegments.extend(userSegments)
2225	self.userVertices,self.userSegments = \
2226	self.weed(self.userVertices,self.userSegments)
2227
2228	def threshold(self,setName,min=None,max=None,attribute_name = 'elevation'):
2229	"""
2230	threshold using d
2231	"""
2232	triangles = self.sets[self.setID[setName]]
2233	A = []
2234
2235	if attribute_name in self.attributeTitles:
2236	i = self.attributeTitles.index(attribute_name)
2237	if not max == None:
2238	for t in triangles:
2239	if (min<self.av_att(t,i)<max):
2240	A.append(t)
2241	else:
2242	for t in triangles:
2243	if (min<self.av_att(t,i)):
2244	A.append(t)
2245	self.sets[self.setID[setName]] = A
2246
2247	def general_threshold(self,setName,min=None,max=None,attribute_name = 'elevation',function = None):
2248	"""
2249	threshold using d
2250	"""
2251	triangles = self.sets[self.setID[setName]]
2252	A = []
2253
2254	if attribute_name in self.attributeTitles:
2255	i = self.attributeTitles.index(attribute_name)
2256	if not max == None:
2257	for t in triangles:
2258	if (min<function(t,i)<max):
2259	A.append(t)
2260	else:
2261	for t in triangles:
2262	if (min<self.function(t,i)):
2263	A.append(t)
2264	self.sets[self.setID[setName]] = A
2265
2266	def av_att(self,triangle,i):
2267	#evaluates the average attribute of the vertices of a triangle.
2268	V = triangle.getVertices()
2269	a0 = (V[0].attributes[i])
2270	a1 = (V[1].attributes[i])
2271	a2 = (V[2].attributes[i])
2272	return (a0+a1+a2)/3
2273
2274	def Courant_ratio(self,triangle,index):
2275	"""
2276	Not the true Courant ratio, just elevation on area
2277	"""
2278	e = self.av_att(triangle,index)
2279	A = triangle.calcArea()
2280	return e/A
2281
2282	def Gradient(self,triangle,index):
2283	V = triangle.vertices
2284	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]
2285	grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
2286	if ((grad_x2)+(grad_y2))**(0.5)<0:
2287	print ((grad_x2)+(grad_y2))**(0.5)
2288	return ((grad_x2)+(grad_y2))**(0.5)
2289
2290
2291	def append_triangle(self,triangle):
2292	self.meshTriangles.append(triangle)
2293
2294	def replace_triangle(self,triangle,replacement):
2295	i = self.meshTriangles.index(triangle)
2296	self.meshTriangles[i]=replacement
2297	assert replacement in self.meshTriangles
2298
2299
2300
2301
2302	def importUngenerateFile(ofile):
2303	"""
2304	import a file, ofile, with the format
2305	[poly]
2306	poly format:
2307	First line: <# of vertices> <x centroid> <y centroid>
2308	Following lines: <x> <y>
2309	last line: "END"
2310
2311	Note: These are clockwise.
2312	"""
2313	fd = open(ofile,'r')
2314	Dict = readUngenerateFile(fd)
2315	fd.close()
2316	return Dict
2317
2318	def readUngenerateFile(fd):
2319	"""
2320	import a file, ofile, with the format
2321	[poly]
2322	poly format:
2323	First line: <# of polynomial> <x centroid> <y centroid>
2324	Following lines: <x> <y>
2325	last line: "END"
2326	"""
2327	END_DELIMITER = 'END\n'
2328
2329	points = []
2330	segments = []
2331
2332	isEnd = False
2333	line = fd.readline() #not used <# of polynomial> <x> <y>
2334	while not isEnd:
2335	line = fd.readline()
2336	fragments = line.split()
2337	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2338	points.append(vert)
2339	PreviousVertIndex = len(points)-1
2340	firstVertIndex = PreviousVertIndex
2341
2342	line = fd.readline() #Read the next line
2343	while line <> END_DELIMITER:
2344	#print "line >" + line + "<"
2345	fragments = line.split()
2346	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2347	points.append(vert)
2348	thisVertIndex = len(points)-1
2349	segment = [PreviousVertIndex,thisVertIndex]
2350	segments.append(segment)
2351	PreviousVertIndex = thisVertIndex
2352	line = fd.readline() #Read the next line
2353	i =+ 1
2354	# If the last and first segments are the same,
2355	# Remove the last segment and the last vertex
2356	# then add a segment from the second last vert to the 1st vert
2357	thisVertIndex = len(points)-1
2358	firstVert = points[firstVertIndex]
2359	thisVert = points[thisVertIndex]
2360	#print "firstVert",firstVert
2361	#print "thisVert",thisVert
2362	if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
2363	points.pop()
2364	segments.pop()
2365	thisVertIndex = len(points)-1
2366	segments.append([thisVertIndex, firstVertIndex])
2367
2368	line = fd.readline() # read <# of polynomial> <x> <y> OR END
2369	#print "line >>" + line + "<<"
2370	if line == END_DELIMITER:
2371	isEnd = True
2372
2373	#print "points", points
2374	#print "segments", segments
2375	ungenerated_dict = {}
2376	ungenerated_dict['points'] = points
2377	ungenerated_dict['segments'] = segments
2378	return ungenerated_dict
2379
2380	def importMeshFromFile(ofile):
2381	"""returns a mesh object, made from a .xya/.pts or .tsh/.msh file
2382	Often raises SyntaxError, IOError
2383	"""
2384	newmesh = None
2385	if ofile[-4:]== ".xya":
2386	#print "loading " + ofile
2387	try:
2388	dict = load_mesh.loadASCII.load_points_file(ofile, delimiter = ',')
2389	except SyntaxError:
2390	dict = load_mesh.loadASCII.load_points_file(ofile, delimiter = ' ')
2391	#print "dict",dict outline_
2392	dict['segmentlist'] = []
2393	dict['segmenttaglist'] = []
2394	dict['regionlist'] = []
2395	dict['regionattributelist'] = []
2396	dict['regionmaxarealist'] = []
2397	dict['holelist'] = []
2398	newmesh= Mesh()
2399	newmesh.setMesh(dict) #FIXME use IOOutline2Mesh
2400	counter = newmesh.removeDuplicatedUserVertices()
2401	if (counter >0):
2402	print "%i duplicate vertices removed from dataset" % (counter)
2403	elif ofile[-4:]== ".pts":
2404	#print "loading " + ofile
2405	dict = load_mesh.loadASCII.load_points_file(ofile)
2406	#print "dict",dict
2407	dict['points'] = dict['pointlist']
2408	dict['outline_segments'] = []
2409	dict['outline_segment_tags'] = []
2410	dict['regions'] = []
2411	dict['region_tags'] = []
2412	dict['region_max_areas'] = []
2413	dict['holes'] = []
2414	newmesh= Mesh(geo_reference = dict['geo_reference'])
2415	newmesh.IOOutline2Mesh(dict) #FIXME use IOOutline2Mesh
2416	counter = newmesh.removeDuplicatedUserVertices()
2417	if (counter >0):
2418	print "%i duplicate vertices removed from dataset" % (counter)
2419	elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
2420	dict = load_mesh.loadASCII.import_triangulation(ofile)
2421	#print "********"
2422	#print "zq mesh.dict",dict
2423	#print "********"
2424	newmesh= Mesh()
2425	newmesh.IOOutline2Mesh(dict)
2426	newmesh.IOTriangulation2Mesh(dict)
2427	else:
2428	raise RuntimeError
2429
2430	if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
2431	newmesh.geo_reference = dict['geo_reference']
2432	return newmesh
2433
2434	def loadPickle(currentName):
2435	fd = open(currentName)
2436	mesh = pickle.load(fd)
2437	fd.close()
2438	return mesh
2439
2440	def square_outline(side_length = 1,up = "top", left = "left", right = "right",
2441	down = "bottom", regions = False):
2442
2443	a = Vertex (0,0)
2444	b = Vertex (0,side_length)
2445	c = Vertex (side_length,0)
2446	d = Vertex (side_length,side_length)
2447
2448	s2 = Segment(b,d, tag = up)
2449	s3 = Segment(b,a, tag = left)
2450	s4 = Segment(d,c, tag = right)
2451	s5 = Segment(a,c, tag = down)
2452
2453	if regions:
2454	e = Vertex (side_length/2,side_length/2)
2455	s6 = Segment(a,e, tag = down + left)
2456	s7 = Segment(b,e, tag = up + left)
2457	s8 = Segment(c,e, tag = down + right)
2458	s9 = Segment(d,e, tag = up + right)
2459	r1 = Region(side_length/2,3.*side_length/4, tag = up)
2460	r2 = Region(1.*side_length/4,side_length/2, tag = left)
2461	r3 = Region(3.*side_length/4,side_length/2, tag = right)
2462	r4 = Region(side_length/2,1.*side_length/4, tag = down)
2463	mesh = Mesh(userVertices=[a,b,c,d,e],
2464	userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
2465	regions = [r1,r2,r3,r4])
2466	else:
2467	mesh = Mesh(userVertices=[a,b,c,d],
2468	userSegments=[s2,s3,s4,s5])
2469
2470	return mesh
2471
2472
2473
2474	def region_strings2ints(region_list):
2475	"""Given a list of (x_int,y_int,tag_string) lists it returns a list of
2476	(x_int,y_int,tag_int) and a list to convert the tag_int's back to
2477	the tag_strings
2478	"""
2479	# Make sure "" has an index of 0
2480	region_list.reverse()
2481	region_list.append((1.0,2.0,""))
2482	region_list.reverse()
2483	convertint2string = []
2484	for i in xrange(len(region_list)):
2485	convertint2string.append(region_list[i][2])
2486	if len(region_list[i]) == 4: # there's an area value
2487	region_list[i] = (region_list[i][0],
2488	region_list[i][1],i,region_list[i][3])
2489	elif len(region_list[i]) == 3: # no area value
2490	region_list[i] = (region_list[i][0],region_list[i][1],i)
2491	else:
2492	print "The region list has a bad size"
2493	# raise an error ..
2494	raise Error
2495
2496	#remove "" from the region_list
2497	region_list.pop(0)
2498
2499	return [region_list, convertint2string]
2500
2501	def region_ints2strings(region_list,convertint2string):
2502	"""Reverses the transformation of region_strings2ints
2503	"""
2504	if region_list[0] != []:
2505	for i in xrange(len(region_list)):
2506	region_list[i] = [convertint2string[int(region_list[i][0])]]
2507	return region_list
2508
2509	def segment_ints2strings(intlist, convertint2string):
2510	"""Reverses the transformation of segment_strings2ints """
2511	stringlist = []
2512	for x in intlist:
2513	stringlist.append(convertint2string[x])
2514	return stringlist
2515
2516	def segment_strings2ints(stringlist, preset):
2517	"""Given a list of strings return a list of 0 to n ints which represent
2518	the strings and a converting list of the strings, indexed by 0 to n ints.
2519	Also, input an initial converting list of the strings
2520	Note, the converting list starts off with
2521	["internal boundary", "external boundary", "internal boundary"]
2522	example input and output
2523	input ["a","b","a","c"],["c"]
2524	output [[2, 1, 2, 0], ['c', 'b', 'a']]
2525
2526	the first element in the converting list will be
2527	overwritten with "".
2528	?This will become the third item in the converting list?
2529
2530	# note, order the initial converting list is important,
2531	since the index = the int tag
2532	"""
2533	nodups = unique(stringlist)
2534	# note, order is important, the index = the int tag
2535	#preset = ["internal boundary", "external boundary"]
2536	#Remove the preset tags from the list with no duplicates
2537	nodups = [x for x in nodups if x not in preset]
2538
2539	try:
2540	nodups.remove("") # this has to go to zero
2541	except ValueError:
2542	pass
2543
2544	# Add the preset list at the beginning of no duplicates
2545	preset.reverse()
2546	nodups.extend(preset)
2547	nodups.reverse()
2548
2549
2550	convertstring2int = {}
2551	convertint2string = []
2552	index = 0
2553	for x in nodups:
2554	convertstring2int[x] = index
2555	convertint2string.append(x)
2556	index += 1
2557	convertstring2int[""] = 0
2558
2559	intlist = []
2560	for x in stringlist:
2561	intlist.append(convertstring2int[x])
2562	return [intlist, convertint2string]
2563
2564
2565
2566
2567	def unique(s):
2568	"""Return a list of the elements in s, but without duplicates.
2569
2570	For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
2571	unique("abcabc") some permutation of ["a", "b", "c"], and
2572	unique(([1, 2], [2, 3], [1, 2])) some permutation of
2573	[[2, 3], [1, 2]].
2574
2575	For best speed, all sequence elements should be hashable. Then
2576	unique() will usually work in linear time.
2577
2578	If not possible, the sequence elements should enjoy a total
2579	ordering, and if list(s).sort() doesn't raise TypeError it's
2580	assumed that they do enjoy a total ordering. Then unique() will
2581	usually work in O(N*log2(N)) time.
2582
2583	If that's not possible either, the sequence elements must support
2584	equality-testing. Then unique() will usually work in quadratic
2585	time.
2586	"""
2587
2588	n = len(s)
2589	if n == 0:
2590	return []
2591
2592	# Try using a dict first, as that's the fastest and will usually
2593	# work. If it doesn't work, it will usually fail quickly, so it
2594	# usually doesn't cost much to try it. It requires that all the
2595	# sequence elements be hashable, and support equality comparison.
2596	u = {}
2597	try:
2598	for x in s:
2599	u[x] = 1
2600	except TypeError:
2601	del u # move on to the next method
2602	else:
2603	return u.keys()
2604
2605	# We can't hash all the elements. Second fastest is to sort,
2606	# which brings the equal elements together; then duplicates are
2607	# easy to weed out in a single pass.
2608	# NOTE: Python's list.sort() was designed to be efficient in the
2609	# presence of many duplicate elements. This isn't true of all
2610	# sort functions in all languages or libraries, so this approach
2611	# is more effective in Python than it may be elsewhere.
2612	try:
2613	t = list(s)
2614	t.sort()
2615	except TypeError:
2616	del t # move on to the next method
2617	else:
2618	assert n > 0
2619	last = t[0]
2620	lasti = i = 1
2621	while i < n:
2622	if t[i] != last:
2623	t[lasti] = last = t[i]
2624	lasti += 1
2625	i += 1
2626	return t[:lasti]
2627
2628	# Brute force is all that's left.
2629	u = []
2630	for x in s:
2631	if x not in u:
2632	u.append(x)
2633	return u
2634
2635	"""Refines triangles
2636
2637	Implements the #triangular bisection?# algorithm.
2638
2639
2640	"""
2641
2642	def Refine(mesh, triangles):
2643	"""
2644	Given a general_mesh, and a triangle number, split
2645	that triangle in the mesh in half. Then to prevent
2646	vertices and edges from meeting, keep refining
2647	neighbouring triangles until the mesh is clean.
2648	"""
2649	state = BisectionState(mesh)
2650	for triangle in triangles:
2651	if not state.refined_triangles.has_key(triangle):
2652	triangle.rotate_longest_side()
2653	state.start(triangle)
2654	Refine_mesh(mesh, state)
2655
2656	def Refine_mesh(mesh, state):
2657	"""
2658	"""
2659	state.getState(mesh)
2660	refine_triangle(mesh,state)
2661	state.evolve()
2662	if not state.end:
2663	Refine_mesh(mesh,state)
2664
2665	def refine_triangle(mesh,state):
2666	split(mesh,state.current_triangle,state.new_point)
2667	if state.case == 'one':
2668	state.r[3]=state.current_triangle#triangle 2
2669
2670	new_triangle_id = len(mesh.meshTriangles)-1
2671	new_triangle = mesh.meshTriangles[new_triangle_id]
2672
2673	split(mesh,new_triangle,state.old_point)
2674	state.r[2]=new_triangle#triangle 1.2
2675	state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
2676	r = state.r
2677	state.repairCaseOne()
2678
2679	if state.case == 'two':
2680	state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2681
2682	new_triangle = state.current_triangle
2683
2684	split(mesh,new_triangle,state.old_point)
2685
2686	state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
2687	state.r[4]=new_triangle#triangle 2.2
2688	r = state.r
2689
2690	state.repairCaseTwo()
2691
2692	if state.case == 'vertex':
2693	state.r[2]=state.current_triangle#triangle 2
2694	state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2695	r = state.r
2696	state.repairCaseVertex()
2697
2698	if state.case == 'start':
2699	state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
2700	state.r[3]=state.current_triangle#triangle 2
2701
2702	if state.next_case == 'boundary':
2703	state.repairCaseBoundary()
2704
2705
2706	def split(mesh, triangle, new_point):
2707	"""
2708	Given a mesh, triangle_id and a new point,
2709	split the corrosponding triangle into two
2710	new triangles and update the mesh.
2711	"""
2712
2713	new_triangle1 = Triangle(new_point,triangle.vertices[0],triangle.vertices[1],attribute = triangle.attribute, neighbors = None)
2714	new_triangle2 = Triangle(new_point,triangle.vertices[2],triangle.vertices[0],attribute = triangle.attribute, neighbors = None)
2715
2716	new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
2717	new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
2718
2719	mesh.meshTriangles.append(new_triangle1)
2720
2721	triangle.vertices = new_triangle2.vertices
2722	triangle.neighbors = new_triangle2.neighbors
2723
2724
2725	class State:
2726
2727	def __init__(self):
2728	pass
2729
2730	class BisectionState(State):
2731
2732
2733	def __init__(self,mesh):
2734	self.len = len(mesh.meshTriangles)
2735	self.refined_triangles = {}
2736	self.mesh = mesh
2737	self.current_triangle = None
2738	self.case = 'start'
2739	self.end = False
2740	self.r = [None,None,None,None,None]
2741
2742	def start(self, triangle):
2743	self.current_triangle = triangle
2744	self.case = 'start'
2745	self.end = False
2746	self.r = [None,None,None,None,None]
2747
2748	def getState(self,mesh):
2749	if not self.case == 'vertex':
2750	self.new_point=self.getNewVertex(mesh, self.current_triangle)
2751	#self.neighbour=self.getNeighbour(mesh, self.current_triangle)
2752	self.neighbour = self.current_triangle.neighbors[0]
2753	if not self.neighbour is None:
2754	self.neighbour.rotate_longest_side()
2755	self.next_case = self.get_next_case(mesh,self.neighbour,self.current_triangle)
2756	if self.case == 'vertex':
2757	self.new_point=self.old_point
2758
2759
2760	def evolve(self):
2761	if self.case == 'vertex':
2762	self.end = True
2763
2764	self.last_case = self.case
2765	self.case = self.next_case
2766
2767	self.old_point = self.new_point
2768	self.current_triangle = self.neighbour
2769
2770	if self.case == 'boundary':
2771	self.end = True
2772	self.refined_triangles[self.r[2]]=1
2773	self.refined_triangles[self.r[3]]=1
2774	if not self.r[4] is None:
2775	self.refined_triangles[self.r[4]]=1
2776	self.r[0]=self.r[2]
2777	self.r[1]=self.r[3]
2778
2779
2780	def getNewVertex(self,mesh,triangle):
2781	coordinate1 = triangle.vertices[1]
2782	coordinate2 = triangle.vertices[2]
2783	a = ([coordinate1.x1.,coordinate1.y1.])
2784	b = ([coordinate2.x1.,coordinate2.y1.])
2785	attributes = []
2786	for i in range(len(coordinate1.attributes)):
2787	att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
2788	attributes.append(att)
2789	new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
2790	newVertex = Vertex(new_coordinate[0],new_coordinate[1], attributes = attributes)
2791	mesh.maxVertexIndex+=1
2792	newVertex.index = mesh.maxVertexIndex
2793	mesh.meshVertices.append(newVertex)
2794	return newVertex
2795
2796	def get_next_case(self, mesh,neighbour,triangle):
2797	"""
2798	Given the locations of two neighbouring triangles,
2799	examine the interior indices of their vertices (i.e.
2800	0,1 or 2) to determine what how the neighbour needs
2801	to be refined.
2802	"""
2803	if (neighbour is None):
2804	next_case = 'boundary'
2805	else:
2806	if triangle.vertices[1].x==neighbour.vertices[2].x:
2807	if triangle.vertices[1].y==neighbour.vertices[2].y:
2808	next_case = 'vertex'
2809	if triangle.vertices[1].x==neighbour.vertices[0].x:
2810	if triangle.vertices[1].y==neighbour.vertices[0].y:
2811	next_case = 'two'
2812	if triangle.vertices[1].x==neighbour.vertices[1].x:
2813	if triangle.vertices[1].y==neighbour.vertices[1].y:
2814	next_case = 'one'
2815	return next_case
2816
2817
2818
2819	def repairCaseVertex(self):
2820
2821	r = self.r
2822
2823
2824	self.link(r[0],r[2])
2825	self.repair(r[0])
2826
2827	self.link(r[1],r[3])
2828	self.repair(r[1])
2829
2830	self.repair(r[2])
2831
2832	self.repair(r[3])
2833
2834
2835	def repairCaseOne(self):
2836	r = self.r
2837
2838
2839	self.link(r[0],r[2])
2840	self.repair(r[0])
2841
2842	self.link(r[1],r[4])
2843	self.repair(r[1])
2844
2845	self.repair(r[4])
2846
2847	def repairCaseTwo(self):
2848	r = self.r
2849
2850	self.link(r[0],r[4])
2851	self.repair(r[0])
2852
2853	self.link(r[1],r[3])
2854	self.repair(r[1])
2855
2856	self.repair(r[4])
2857
2858	def repairCaseBoundary(self):
2859	r = self.r
2860	self.repair(r[2])
2861	self.repair(r[3])
2862
2863
2864
2865	def repair(self,triangle):
2866	"""
2867	Given a triangle that knows its neighbours, this will
2868	force the neighbours to comply.
2869
2870	However, it needs to compare the vertices of triangles
2871	for this implementation
2872
2873	But it doesn't work for invalid neighbour structures
2874	"""
2875	n=triangle.neighbors
2876	for i in (0,1,2):
2877	if not n[i] is None:
2878	for j in (0,1,2):#to find which side of the list is broken
2879	if not (n[i].vertices[j] in triangle.vertices):
2880	#ie if j is the side of n that needs fixing
2881	k = j
2882	n[i].neighbors[k]=triangle
2883
2884
2885
2886	def link(self,triangle1,triangle2):
2887	"""
2888	make triangle1 neighbors point to t
2889	#count = 0riangle2
2890	"""
2891	count = 0
2892	for i in (0,1,2):#to find which side of the list is broken
2893	if not (triangle1.vertices[i] in triangle2.vertices):
2894	j = i
2895	count+=1
2896	assert count == 1
2897	triangle1.neighbors[j]=triangle2
2898
2899
2900
2901	if __name__ == "__main__":
2902	#from mesh import *
2903	# THIS CAN BE DELETED
2904	m = Mesh()
2905	dict = importUngenerateFile("ungen_test.txt")
2906	m.addVertsSegs(dict)
2907	print m

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

Download in other formats: