Context Navigation

source: inundation/pmesh/mesh.py @ 2378

Last change on this file since 2378 was 2345, checked in by duncan, 19 years ago
fixed error in building command line for generte mesh
File size: 126.8 KB

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

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

Download in other formats: