Context Navigation

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

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

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

Download in other formats: