Context Navigation

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

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

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

Download in other formats: