Context Navigation

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

Last change on this file since 669 was 650, checked in by duncan, 20 years ago
internal tags nolonger default to 'internal', they now default to not being a tag.
File size: 63.6 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
	16	import load_mesh.loadASCII
	17
	18	import sys
	19	import math
	20	import triang
	21	import re
	22	import os
	23	import pickle
	24
	25	import types
	26	import exceptions
	27
[650]	28	# 1st and third values must be the same
	29	#initialconversions = ['internal', 'external','internal']
	30	# FIXME: maybe make this a switch that the user can change?
	31	initialconversions = ['', 'external','']
[349]	32
	33	#from os import sep
	34	#sys.path.append('..'+sep+'pmesh')
	35	#print "sys.path",sys.path
	36
	37	class MeshObject:
	38	"""
	39	An abstract superclass for the basic mesh objects, eg vertex, segment,
	40	triangle.
	41	"""
	42	def __init__(self):
	43	pass
	44
	45	class Point(MeshObject):
	46	"""
	47	Define a point in a 2D space.
	48	"""
	49	def __init__(self,X,Y):
	50	__slots__ = ['x','y']
	51	self.x=X
	52	self.y=Y
	53
	54	def DistanceToPoint(self, OtherPoint):
	55	"""
	56	Returns the distance from this point to another
	57	"""
	58	SumOfSquares = ((self.x - OtherPoint.x)2) + ((self.y - OtherPoint.y)2)
	59	return math.sqrt(SumOfSquares)
	60
	61	def IsInsideCircle(self, Center, Radius):
	62	"""
	63	Return 1 if this point is inside the circle,
	64	0 otherwise
	65	"""
	66
	67	if (self.DistanceToPoint(Center)<Radius):
	68	return 1
	69	else:
	70	return 0
	71
	72	def __repr__(self):
	73	return "(%f,%f)" % (self.x,self.y)
	74
	75	def cmp_xy(self, point):
	76	if self.x < point.x:
	77	return -1
	78	elif self.x > point.x:
	79	return 1
	80	else:
	81	if self.y < point.y:
	82	return -1
	83	elif self.y > point.y:
	84	return 1
	85	else:
	86	return 0
	87
	88	def same_x_y(self, point):
	89	if self.x == point.x and self.y == point.y:
	90	return True
	91	else:
	92	return False
	93
	94
	95
	96	class Vertex(Point):
	97	"""
	98	A point on the mesh.
	99	Object attributes based on the Triangle program
	100	"""
	101	def __init__(self,X,Y, attributes = None):
	102	__slots__ = ['x','y','attributes']
	103	self.x=X
	104	self.y=Y
	105	self.attributes=[]
	106
	107	if attributes is None:
	108	self.attributes=[]
	109	else:
	110	self.attributes=attributes
	111
	112
	113	def setAttributes(self,attributes):
	114	"""
	115	attributes is a list.
	116	"""
	117	self.attributes = attributes
	118
	119	VERTEXSQUARESIDELENGTH = 6
	120	def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, yoffset =0, ):
	121	x = scale*(self.x + xoffset)
	122	y = -1scale(self.y + yoffset) # - since for a canvas - is up
	123	#print "draw x:", x
	124	#print "draw y:", y
	125	cornerOffset= self.VERTEXSQUARESIDELENGTH/2
[420]	126
	127	# A hack to see the vert tags
	128	#canvas.create_text(x+ 2*cornerOffset,
	129	# y+ 2*cornerOffset,
	130	# text=tags)
	131
[607]	132	return canvas.create_rectangle(x-cornerOffset,
	133	y-cornerOffset,
	134	x+cornerOffset,
	135	y+cornerOffset,
	136	tags = tags,
	137	outline=colour,
	138	fill = 'white')
[349]	139
	140	def __repr__(self):
	141	return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
	142
	143	class Hole(Point):
	144	"""
	145	A region of the mesh were no triangles are generated.
	146	Defined by a point in the hole enclosed by segments.
	147	"""
	148	HOLECORNERLENGTH = 6
	149	def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ):
	150	x = scale*(self.x + xoffset)
	151	y = -1scale(self.y + yoffset) # - since for a canvas - is up
	152	#print "draw x:", x
	153	#print "draw y:", y
	154	cornerOffset= self.HOLECORNERLENGTH/2
	155	return canvas.create_oval(x-cornerOffset,
	156	y-cornerOffset,
	157	x+cornerOffset,
	158	y+cornerOffset,
	159	tags = tags,
	160	outline=colour,
	161	fill = 'white')
	162
	163	class Region(Point):
	164	"""
	165	A region of the mesh, defined by a point in the region
	166	enclosed by segments. Used to tag areas.
	167	"""
	168	CROSSLENGTH = 6
	169	TAG = 0
	170	MAXAREA = 1
	171
	172	def __init__(self,X,Y, tag = None, maxArea = None):
	173	"""Precondition: tag is a string and maxArea is a real
	174	"""
	175	# This didn't work.
	176	#super(Region,self)._init_(self,X,Y)
	177	self.x=X
	178	self.y=Y
	179	self.attributes =[] # index 0 is the tag string
	180	#optoinal index 1 is the max triangle area
	181	#NOTE the size of this attribute is assumed
	182	# to be 1 or 2 in regionstrings2int
	183	if tag is None:
	184	self.attributes.append("")
	185	else:
	186	self.attributes.append(tag) #this is a string
	187
	188	if maxArea is not None:
	189	self.setMaxArea(maxArea) # maxArea is a number
	190
	191	def getTag(self,):
	192	return self.attributes[self.TAG]
	193
	194	def setTag(self,tag):
	195	self.attributes[self.TAG] = tag
	196
	197	def getMaxArea(self):
	198	""" Returns the Max Area of a Triangle or
	199	None, if the Max Area has not been set.
	200	"""
	201	if self.isMaxArea():
	202	return self.attributes[1]
	203	else:
	204	return None
	205
	206	def setMaxArea(self,MaxArea):
	207	if self.isMaxArea():
	208	self.attributes[self.MAXAREA] = float(MaxArea)
	209	else:
	210	self.attributes.append( float(MaxArea) )
	211
	212	def deleteMaxArea(self):
	213	if self.isMaxArea():
	214	self.attributes.pop(self.MAXAREA)
	215
	216	def isMaxArea(self):
	217	return len(self.attributes)> 1
	218
	219	def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "black"):
	220	"""
	221	Draw a black cross, returning the objectID
	222	"""
	223	x = scale*(self.x + xoffset)
	224	y = -1scale(self.y + yoffset)
	225	cornerOffset= self.CROSSLENGTH/2
	226	return canvas.create_polygon(x,
	227	y-cornerOffset,
	228	x,
	229	y,
	230	x+cornerOffset,
	231	y,
	232	x,
	233	y,
	234	x,
	235	y+cornerOffset,
	236	x,
	237	y,
	238	x-cornerOffset,
	239	y,
	240	x,
	241	y,
	242	tags = tags,
	243	outline = colour,fill = '')
	244
	245	def __repr__(self):
	246	if self.isMaxArea():
	247	area = self.getMaxArea()
	248	return "(%f,%f,%s,%f)" % (self.x,self.y,
	249	self.getTag(), self.getMaxArea())
	250	else:
	251	return "(%f,%f,%s)" % (self.x,self.y,
	252	self.getTag())
	253
	254	class Triangle(MeshObject):
	255	"""
	256	A triangle element, defined by 3 vertices.
	257	Attributes based on the Triangle program.
	258	"""
	259
	260	def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ):
	261	"""
	262	Vertices, the initial arguments, are listed in counterclockwise order.
	263	"""
	264	self.vertices= [vertex1,vertex2, vertex3 ]
	265
	266	if attribute is None:
	267	self.attribute =""
	268	else:
	269	self.attribute = attribute #this is a string
	270
	271	if neighbors is None:
	272	self.neighbors=[]
	273	else:
	274	self.neighbors=neighbors
	275
	276	def getVertices(self):
	277	return self.vertices
	278
	279	def calcArea(self):
	280	ax = self.vertices[0].x
	281	ay = self.vertices[0].y
	282
	283	bx = self.vertices[1].x
	284	by = self.vertices[1].y
	285
	286	cx = self.vertices[2].x
	287	cy = self.vertices[2].y
	288
	289	return abs((bxay-axby)+(cxby-bxcy)+(axcy-cxay))/2
	290
	291	def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
	292	"""
	293	neighbor1 is the triangle opposite vertex1 and so on.
	294	Null represents no neighbor
	295	"""
	296	self.neighbors = [neighbor1, neighbor2, neighbor3]
	297
	298	def setAttribute(self,attribute):
	299	"""
	300	neighbor1 is the triangle opposite vertex1 and so on.
	301	Null represents no neighbor
	302	"""
	303	self.attribute = attribute
	304
	305	def __repr__(self):
	306	return "[%s,%s]" % (self.vertices,self.attribute)
	307
	308
	309	def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"):
	310	"""
	311	Draw a red triagle, returning the objectID
	312	"""
	313	return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
	314	scale-1(self.vertices[1].y + yoffset),
	315	scale*(self.vertices[0].x + xoffset),
	316	scale-1(self.vertices[0].y + yoffset),
	317	scale*(self.vertices[2].x + xoffset),
	318	scale-1(self.vertices[2].y + yoffset),
	319	tags = tags,
	320	outline = colour,fill = '')
	321
	322	class Segment(MeshObject):
	323	"""
	324	Segments are edges whose presence in the triangulation is enforced.
	325
	326	"""
	327	def __init__(self, vertex1, vertex2, marker = None ):
	328	"""
	329	Each segment is specified by listing the vertices of its endpoints
	330	"""
	331
	332	assert(vertex1 != vertex2)
	333	self.vertices = [vertex1,vertex2 ]
	334
	335	if marker is None:
	336	self.marker = self.__class__.default
	337	else:
	338	self.marker = marker #this is a string
	339
	340	def __repr__(self):
	341	return "[%s,%s]" % (self.vertices,self.marker)
	342
	343
	344	def draw(self, canvas, tags,scale=1 , xoffset=0 , yoffset=0,colour='blue' ):
	345	x=[]
	346	y=[]
	347	for end in self.vertices:
	348	#end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
	349	x.append(scale*(end.x + xoffset))
	350	y.append(-1scale(end.y + yoffset)) # - since for a canvas - is up
	351	return canvas.create_line(x[0], y[0], x[1], y[1], tags = tags,fill=colour)
	352	# Class methods
	353	def set_default_tag(cls, default):
	354	cls.default = default
	355
	356	def get_default_tag(cls):
	357	return cls.default
	358
	359	set_default_tag = classmethod(set_default_tag)
	360	get_default_tag = classmethod(get_default_tag)
	361
	362	Segment.set_default_tag("")
[404]	363
[349]	364	class Mesh:
	365	"""
	366	Representation of a 2D triangular mesh.
	367	User attributes describe the mesh region/segments/vertices/attributes
	368
	369	mesh attributes describe the mesh that is produced eg triangles and vertices.
	370
	371	The Mesh holds user information to define the
	372	"""
	373
	374	def __repr__(self):
	375	return """
	376	mesh Triangles: %s
	377	mesh Segments: %s
	378	mesh Vertices: %s
	379	user Segments: %s
	380	user Vertices: %s
	381	holes: %s
	382	regions: %s""" % (self.meshTriangles,
	383	self.meshSegments,
	384	self.meshVertices,
	385	self.userSegments,
	386	self.userVertices,
	387	self.holes,
	388	self.regions)
	389
	390	def __init__(self, userSegments=None, userVertices=None, holes=None, regions=None):
	391	self.meshTriangles=[]
	392	self.meshSegments=[]
	393	self.meshVertices=[]
	394
	395	if userSegments is None:
	396	self.userSegments=[]
	397	else:
	398	self.userSegments=userSegments
	399
	400	if userVertices is None:
	401	self.userVertices=[]
	402	else:
	403	self.userVertices=userVertices
	404
	405	if holes is None:
	406	self.holes=[]
	407	else:
	408	self.holes=holes
	409
	410	if regions is None:
	411	self.regions=[]
	412	else:
	413	self.regions=regions
	414	def __cmp__(self,other):
	415
	416	# A dic for the initial m
	417	dic = self.Mesh2triangList()
	418	dic_mesh = self.Mesh2MeshList()
	419	for element in dic_mesh.keys():
	420	dic[element] = dic_mesh[element]
	421
	422	# A dic for the exported/imported m
	423	dic_other = other.Mesh2triangList()
	424	dic_mesh = other.Mesh2MeshList()
	425	for element in dic_mesh.keys():
	426	dic_other[element] = dic_mesh[element]
	427
	428	#print "dsg************************8"
	429	#print "dic ",dic
	430	#print "*******8"
	431	#print "mesh",dic_other
	432	#print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
	433	#print "dsg************************8"
	434
	435	return (dic.__cmp__(dic_other))
	436
	437	def generateMesh(self, mode = None, maxArea = None, isRegionalMaxAreas = True):
	438	"""
	439	Based on the current user vaules, holes and regions
	440	generate a new mesh
	441	mode is a string that sets conditions on the mesh generations
	442	see triangle_instructions.txt for a definition of the commands
	443	PreCondition: maxArea is a double
	444	"""
	445	if mode == None:
	446	self.mode = ""
	447	else:
	448	self.mode = mode
	449
	450	if not re.match('p',self.mode):
	451	self.mode += 'p' #p - read a planar straight line graph.
	452	#there must be segments to use this switch
	453	# TODO throw an aception if there aren't seg's
	454	# it's more comlex than this. eg holes
	455	if not re.match('z',self.mode):
	456	self.mode += 'z' # z - Number all items starting from zero (rather than one)
	457	if not re.match('n',self.mode):
	458	self.mode += 'n' # n - output a list of neighboring triangles
	459	if not re.match('A',self.mode):
	460	self.mode += 'A' # A - output region attribute list for triangles
	461	if not re.match('V',self.mode):
	462	self.mode += 'V' # V - output info about what Triangle is doing
	463
	464	if maxArea != None:
	465	self.mode += 'a' + str(maxArea)
	466
	467	if isRegionalMaxAreas:
	468	self.mode += 'a'
	469
	470	meshDict = self.Mesh2triangList()
	471	#print "!@!@ This is going to triangle !@!@"
	472	#print meshDict
	473	#print "!@!@ This is going to triangle !@!@"
	474
	475	#print "meshDict['segmentmarkerlist']", meshDict['segmentmarkerlist']
	476	#initialconversions = ['internal boundary', 'external boundary','internal boundary']
	477	[meshDict['segmentmarkerlist'],
	478	segconverter] = segment_strings2ints(meshDict['segmentmarkerlist'],
	479	initialconversions)
	480	#print "regionlist",meshDict['regionlist']
	481	[meshDict['regionlist'],
	482	regionconverter] = region_strings2ints(meshDict['regionlist'])
	483	#print "regionlist",meshDict['regionlist']
	484
	485	#print "meshDict['segmentmarkerlist']", meshDict['segmentmarkerlist'
	486	generatedMesh = triang.genMesh(
	487	meshDict['pointlist'],
	488	meshDict['segmentlist'],
	489	meshDict['holelist'],
	490	meshDict['regionlist'],
	491	meshDict['pointattributelist'],
	492	meshDict['segmentmarkerlist'],
	493	[], # since the trianglelist isn't used
	494	self.mode)
	495	#print "generated",generatedMesh['generatedsegmentmarkerlist']
	496	generatedMesh['generatedsegmentmarkerlist'] = \
	497	segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
	498	segconverter)
	499	#print "processed gen",generatedMesh['generatedsegmentmarkerlist']
	500	generatedMesh['generatedtriangleattributelist'] = \
	501	region_ints2strings(generatedMesh['generatedtriangleattributelist'],
	502	regionconverter)
	503
	504	self.setTriangulation(generatedMesh)
	505
	506	def addUserPoint(self, pointType, x,y):
	507	if pointType == Vertex:
	508	point = self.addUserVertex(x,y)
	509	if pointType == Hole:
	510	point = self.addHole(x,y)
	511	if pointType == Region:
	512	point = self.addRegion(x,y)
	513	return point
	514
	515	def addUserVertex(self, x,y):
	516	v=Vertex(x, y)
	517	self.userVertices.append(v)
	518	return v
	519
	520	def addHole(self, x,y):
	521	h=Hole(x, y)
	522	self.holes.append(h)
	523	return h
	524
	525	def addRegion(self, x,y):
	526	h=Region(x, y)
	527	self.regions.append(h)
	528	return h
	529
	530	def getUserVertices(self):
	531	return self.userVertices
	532
	533	def getUserSegments(self):
	534	return self.userSegments
	535
	536	def getTriangulation(self):
	537	return self.meshTriangles
	538
	539	def getMeshVertices(self):
	540	return self.meshVertices
	541
	542	def getMeshSegments(self):
	543	return self.meshSegments
	544
	545	def getHoles(self):
	546	return self.holes
	547
	548	def getRegions(self):
	549	return self.regions
	550
	551	def isTriangulation(self):
	552	if self.meshVertices == []:
	553	return False
	554	else:
	555	return True
	556
	557	def addUserSegment(self, v1,v2):
	558	"""
	559	PRECON: A segment between the two vertices is not already present.
	560	Check by calling isUserSegmentNew before calling this function.
	561
	562	"""
	563	s=Segment( v1,v2)
	564	self.userSegments.append(s)
	565	return s
	566
	567	def clearTriangulation(self):
	568
	569	#Clear the current generated mesh values
	570	self.meshTriangles=[]
	571	self.meshSegments=[]
	572	self.meshVertices=[]
	573
	574	def removeDuplicatedUserVertices(self):
	575	"""Pre-condition: There are no user segments
	576	This function will keep the first duplicate
	577	"""
	578	assert self.userSegments == []
	579	self.userVertices, counter = self.removeDuplicatedVertices(self.userVertices)
	580	return counter
	581
	582	def removeDuplicatedVertices(self, Vertices):
	583	"""
	584	This function will keep the first duplicate, remove all others
	585	Precondition: Each vertex has a dupindex, which is the list
	586	index.
	587	"""
	588	remove = []
	589	index = 0
	590	for v in Vertices:
	591	v.dupindex = index
	592	index += 1
	593	t = list(Vertices)
	594	t.sort(Point.cmp_xy)
	595
	596	length = len(t)
	597	behind = 0
	598	ahead = 1
	599	counter = 0
	600	while ahead < length:
	601	b = t[behind]
	602	ah = t[ahead]
	603	if (b.y == ah.y and b.x == ah.x):
	604	remove.append(ah.dupindex)
	605	behind += 1
	606	ahead += 1
	607
	608	# remove the duplicate vertices
	609	remove.sort()
	610	remove.reverse()
	611	for i in remove:
	612	Vertices.pop(i)
	613	pass
	614
	615	#Remove the attribute that this function added
	616	for v in Vertices:
	617	del v.dupindex
	618	return Vertices,counter
	619
	620	def thinoutVertices(self, delta):
	621	"""Pre-condition: There are no user segments
	622	This function will keep the first duplicate
	623	"""
	624	assert self.userSegments == []
	625	#t = self.userVertices
	626	#self.userVertices =[]
	627	boxedVertices = {}
	628	thinnedUserVertices =[]
	629	delta = round(delta,1)
	630
	631	for v in self.userVertices :
	632	# marker is the center of the boxes
	633	marker = (round(v.x/delta,0)delta,round(v.y/delta,0)delta)
	634	#this creates a dict of lists of faces, indexed by marker
	635	boxedVertices.setdefault(marker,[]).append(v)
	636
	637	for [marker,verts] in boxedVertices.items():
	638	min = delta
	639	bestVert = None
	640	markerVert = Vertex(marker[0],marker[1])
	641	for v in verts:
	642	dist = v.DistanceToPoint(markerVert)
	643	if (dist<min):
	644	min = dist
	645	bestVert = v
	646	thinnedUserVertices.append(bestVert)
	647	self.userVertices =thinnedUserVertices
	648
	649
	650	def isUserSegmentNew(self, v1,v2):
	651	identicalSegs= [x for x in self.userSegments if (x.vertices[0] == v1 and x.vertices[1] == v2) or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
	652
	653	return len(identicalSegs) == 0
	654
	655	def deleteSegsOfVertex(self, delVertex):
	656	"""
	657	Delete this vertex and any segments that connect to it.
	658	"""
	659	#Find segments that connect to delVertex
	660	deletedSegments = []
	661	for seg in self.userSegments:
	662	if (delVertex in seg.vertices):
	663	deletedSegments.append(seg)
	664	# Delete segments that connect to delVertex
	665	for seg in deletedSegments:
	666	self.userSegments.remove(seg)
	667	# Delete delVertex
	668	# self.userVertices.remove(delVertex)
	669	return deletedSegments
	670
	671
	672	def deleteMeshObject(self, MeshObject):
	673	"""
	674	Returns a list of all objects that were removed
	675	"""
	676	deletedObs = []
	677	if isinstance(MeshObject, Vertex ):
	678	deletedObs = self.deleteSegsOfVertex(MeshObject)
	679	deletedObs.append(MeshObject)
	680	self.userVertices.remove(MeshObject)
	681	elif isinstance(MeshObject, Segment):
	682	deletedObs.append(MeshObject)
	683	self.userSegments.remove(MeshObject)
	684	elif isinstance(MeshObject, Hole):
	685	deletedObs.append(MeshObject)
	686	self.holes.remove(MeshObject)
	687	elif isinstance(MeshObject, Region):
	688	deletedObs.append(MeshObject)
	689	self.regions.remove(MeshObject)
	690	return deletedObs
	691
	692	def Mesh2triangList(self):
	693	"""
	694	Convert the Mesh to a dictionary of the lists needed for the triang modul;
	695	points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	696	pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
	697	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
	698	hole list: [(x1,y1),...](Tuples of doubles, one inside each hole region)
	699	regionlist: [ (x1,y1,index),...] (Tuple of 3 doubles)
	700	Note, this adds an index attribute to the user Vertex objects.
	701	"""
	702
	703	meshDict = {}
	704
	705	pointlist=[]
	706	pointattributelist=[]
[404]	707
[349]	708	index = 0
	709	for vertex in self.userVertices:
	710	vertex.index = index
	711	pointlist.append((vertex.x,vertex.y))
	712	pointattributelist.append((vertex.attributes))
	713
	714	index += 1
	715	meshDict['pointlist'] = pointlist
	716	meshDict['pointattributelist'] = pointattributelist
	717
	718	segmentlist=[]
	719	segmentmarkerlist=[]
	720	for seg in self.userSegments:
	721	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
	722	segmentmarkerlist.append(seg.marker)
	723	meshDict['segmentlist'] =segmentlist
	724	meshDict['segmentmarkerlist'] =segmentmarkerlist
	725
	726	holelist=[]
	727	for hole in self.holes:
	728	holelist.append((hole.x,hole.y))
	729	meshDict['holelist'] = holelist
	730
	731	regionlist=[]
	732	for region in self.regions:
	733	if (region.getMaxArea() != None):
	734	regionlist.append((region.x,region.y,region.getTag(),
	735	region.getMaxArea()))
	736	else:
	737	regionlist.append((region.x,region.y,region.getTag()))
	738	meshDict['regionlist'] = regionlist
	739	#print "(("
	740	#print meshDict
	741	#print meshDict['regionlist']
	742	#print "(("
	743	return meshDict
	744
	745	def Mesh2MeshList(self):
	746	"""
	747	Convert the Mesh to a dictionary of lists describing the triangulation variables;
	748	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	749	generated point attribute list: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
	750	generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
	751	hole list: [(x1,y1),...](Tuples of doubles, one inside each hole region)
	752	regionlist: [ (x1,y1,index),...] (Tuple of 3 doubles)
	753	Note, this adds an index attribute to the user Vertex objects.
	754	"""
	755
	756	meshDict = {}
	757	pointlist=[]
	758	pointattributelist=[]
	759
	760	index = 0
	761	for vertex in self.meshVertices:
	762	vertex.index = index
	763	pointlist.append((vertex.x,vertex.y))
	764	pointattributelist.append((vertex.attributes))
	765	index += 1
	766
	767	meshDict['generatedpointlist'] = pointlist
	768	meshDict['generatedpointattributelist'] = pointattributelist
	769
	770	#segments
	771	segmentlist=[]
	772	segmentmarkerlist=[]
	773	for seg in self.meshSegments:
	774	segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
	775	segmentmarkerlist.append(seg.marker)
	776	meshDict['generatedsegmentlist'] =segmentlist
	777	meshDict['generatedsegmentmarkerlist'] =segmentmarkerlist
	778
	779	# Make sure that the indexation is correct
	780	index = 0
	781	for tri in self.meshTriangles:
	782	tri.index = index
	783	index += 1
	784
	785	trianglelist = []
	786	triangleattributelist = []
	787	triangleneighborlist = []
	788	for tri in self.meshTriangles:
	789	trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,tri.vertices[2].index))
	790	triangleattributelist.append(tri.attribute)
	791	neighborlist = [-1,-1,-1]
	792	for neighbor,index in map(None,tri.neighbors,
	793	range(len(tri.neighbors))):
	794	if neighbor:
	795	neighborlist[index] = neighbor.index
	796	triangleneighborlist.append(neighborlist)
	797
	798	meshDict['generatedtrianglelist'] = trianglelist
	799	meshDict['generatedtriangleattributelist'] = triangleattributelist
	800	meshDict['generatedtriangleneighborlist'] = triangleneighborlist
	801
	802	#print "(("
	803	#print meshDict
	804	#print "(("
[369]	805
[349]	806	return meshDict
	807
	808
	809	def Mesh2MeshDic(self):
	810	"""
	811	Convert the user and generated info of a mesh to a dictionary
	812	structure
	813	"""
	814	dic = self.Mesh2triangList()
	815	dic_mesh = self.Mesh2MeshList()
	816	for element in dic_mesh.keys():
	817	dic[element] = dic_mesh[element]
	818	return dic
	819
	820	def setTriangulation(self, genDict):
	821	"""
	822	Set the mesh attributes given a dictionary of the lists
	823	returned from the triang module
	824	generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	825	generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
	826	generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
	827	generated segment marker list: [S1Marker, S2Marker, ...] (list of ints)
	828	triangle list: [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
	829	triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
	830	triangle attribute list: [(T1att), (T2att), ...] (list of a list of strings)
	831	"""
	832
	833	#Clear the current generated mesh values
	834	self.meshTriangles=[]
	835	self.meshSegments=[]
	836	self.meshVertices=[]
	837
	838	#print "@#@#@#"
	839	#print genDict
	840	#print "@#@#@#"
	841
	842	index = 0
	843	for point in genDict['generatedpointlist']:
	844	v=Vertex(point[0], point[1])
	845	v.index = index
	846	index +=1
	847	self.meshVertices.append(v)
	848
	849	index = 0
	850	for seg,marker in map(None,genDict['generatedsegmentlist'],genDict['generatedsegmentmarkerlist']):
	851	segObject = Segment( self.meshVertices[seg[0]],
	852	self.meshVertices[seg[1]], marker = marker )
	853	segObject.index = index
	854	index +=1
	855	self.meshSegments.append(segObject)
	856
	857	index = 0
	858	for triangle in genDict['generatedtrianglelist']:
	859	tObject =Triangle( self.meshVertices[triangle[0]],
	860	self.meshVertices[triangle[1]],
	861	self.meshVertices[triangle[2]] )
	862	tObject.index = index
	863	index +=1
	864	self.meshTriangles.append(tObject)
	865
	866	index = 0
	867	for att in genDict['generatedtriangleattributelist']:
	868	if att == []:
	869	self.meshTriangles[index].setAttribute("")
	870	else:
	871	# Note, is the first attribute always the region att?
	872	# haven't confirmed this
	873	self.meshTriangles[index].setAttribute(att[0])
	874	index += 1
	875
	876	index = 0
	877	for att in genDict['generatedpointattributelist']:
	878	if att == None:
	879	self.meshVertices[index].setAttributes([])
	880	else:
	881	self.meshVertices[index].setAttributes(att)
	882	index += 1
	883
	884	index = 0
	885	for triangle in genDict['generatedtriangleneighborlist']:
	886	# Build a list of triangle object neighbors
	887	ObjectNeighbor = []
	888	for neighbor in triangle:
	889	if ( neighbor != -1):
	890	ObjectNeighbor.append(self.meshTriangles[neighbor])
	891	else:
	892	ObjectNeighbor.append(None)
	893	self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],ObjectNeighbor[1],ObjectNeighbor[2])
	894	index += 1
	895
	896
	897	def setMesh(self, genDict):
	898	"""
	899	Set the user Mesh attributes given a dictionary of the lists
	900	point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	901	point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
	902	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
	903	segment marker list: [S1Marker, S2Marker, ...] (list of ints)
	904	region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	905	region attribute list: ["","reservoir",""] list of strings
	906	region max area list:[real, None, Real,...] list of None and reals
	907
	908	mesh is an instance of a mesh object
	909	"""
	910
	911	#Clear the current user mesh values
	912	self.userSegments=[]
	913	self.userVertices=[]
	914	self.Holes=[]
	915	self.Regions=[]
	916
	917	#print "@#@#@#"
	918	#print genDict
	919	#print "@#@#@#"
	920
[404]	921	#index = 0
[349]	922	for point in genDict['pointlist']:
	923	v=Vertex(point[0], point[1])
[404]	924	#v.index = index
	925	#index +=1
[349]	926	self.userVertices.append(v)
	927
[404]	928	#index = 0
[349]	929	for seg,marker in map(None,genDict['segmentlist'],genDict['segmentmarkerlist']):
	930	segObject = Segment( self.userVertices[seg[0]],
	931	self.userVertices[seg[1]], marker = marker )
[404]	932	#segObject.index = index
	933	#index +=1
[349]	934	self.userSegments.append(segObject)
[364]	935
	936	# Remove the loading of attribute info.
	937	# Have attribute info added using least_squares in pyvolution
	938	# index = 0
	939	# for att in genDict['pointattributelist']:
	940	# if att == None:
	941	# self.userVertices[index].setAttributes([])
	942	# else:
	943	# self.userVertices[index].setAttributes(att)
	944	# index += 1
[349]	945
[404]	946	#index = 0
[349]	947	for point in genDict['holelist']:
	948	h=Hole(point[0], point[1])
[404]	949	#h.index = index
	950	#index +=1
[349]	951	self.holes.append(h)
	952
[404]	953	#index = 0
[349]	954	for reg,att,maxArea in map(None,
	955	genDict['regionlist'],
	956	genDict['regionattributelist'],
	957	genDict['regionmaxarealist']):
	958	Object = Region( reg[0],
	959	reg[1],
	960	tag = att,
	961	maxArea = maxArea)
[404]	962	#Object.index = index
	963	#index +=1
[349]	964	self.regions.append(Object)
[404]	965
	966	def addVertsSegs(self, outlineDict):
	967	"""
	968	Add out-line (user Mesh) attributes given a dictionary of the lists
	969	point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
	970	segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
	971	#segment marker list: [S1Marker, S2Marker, ...] (list of ints)
	972	"""
[349]	973
[404]	974	localUserVertices = []
	975	#index = 0
	976	for point in outlineDict['pointlist']:
	977	v=Vertex(point[0], point[1])
	978	#v.index = index
	979	#index +=1
	980	self.userVertices.append(v)
	981	localUserVertices.append(v)
	982
	983	#index = 0
	984	for seg in map(None,outlineDict['segmentlist']):
	985	segObject = Segment( localUserVertices[seg[0]],
	986	localUserVertices[seg[1]] )
	987	#segObject.index = index
	988	#index +=1
	989	self.userSegments.append(segObject)
[349]	990
	991	def TestautoSegment(self):
	992	newsegs = []
	993	s1 = Segment(self.userVertices[0],
	994	self.userVertices[1])
	995	s2 = Segment(self.userVertices[0],
	996	self.userVertices[2])
	997	s3 = Segment(self.userVertices[2],
	998	self.userVertices[1])
	999	if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
	1000	newsegs.append(s1)
	1001	if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
	1002	newsegs.append(s2)
	1003	if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
	1004	newsegs.append(s3)
	1005	self.userSegments.extend(newsegs)
	1006	return newsegs
	1007
	1008
	1009	def savePickle(self, currentName):
	1010	fd = open(currentName, 'w')
	1011	pickle.dump(self,fd)
	1012	fd.close()
	1013
	1014	def autoSegment(self):
	1015	"""
	1016	initially work by running an executable
	1017	Later compile the c code with a python wrapper.
	1018
	1019	Precon: There must be 3 or more vertices in the userVertices structure
	1020	"""
	1021	newsegs = []
	1022	inputfile = 'hull_in.txt'
	1023	outputfile = inputfile + '-alf'
	1024	#write vertices to file
	1025	fd = open(inputfile,'w')
	1026	for v in self.userVertices:
	1027	fd.write(str(v.x))
	1028	fd.write(' ')
	1029	fd.write(str(v.y))
	1030	fd.write('\n')
	1031	fd.close()
	1032
	1033	# os.system('SET PATH=..\pmesh;%PATH%') # this doesn't work
	1034	#print "os.getcwd()",os.getcwd() # this could be usefull
	1035	#command = "k:\\EQRMsens\\" + s
	1036	#os.chdir(command)
	1037	#command = 'env'
	1038	#os.system(command)
	1039
	1040	#run hull executable
	1041	#warning need to compile hull for the current operating system
	1042	command = 'hull.exe -A -i ' + inputfile
	1043	os.system(command)
	1044
	1045	#read results into this object
	1046	fd = open(outputfile)
	1047	lines = fd.readlines()
	1048	fd.close()
	1049	#print "(((*("
	1050	#print lines
	1051	#print "(((*("
	1052	lines.pop(0) #remove the first (title) line
	1053	for line in lines:
	1054	vertindexs = line.split()
	1055	#print 'int(vertindexs[0])', int(vertindexs[0])
	1056	#print 'int(vertindexs[1])', int(vertindexs[1])
	1057	#print 'self.userVertices[int(vertindexs[0])]' ,self.userVertices[int(vertindexs[0])]
	1058	#print 'self.userVertices[int(vertindexs[1])]' ,self.userVertices[int(vertindexs[1])]
	1059	v1 = self.userVertices[int(vertindexs[0])]
	1060	v2 = self.userVertices[int(vertindexs[1])]
	1061
	1062	if self.isUserSegmentNew(v1,v2):
	1063	newseg = Segment(v1, v2)
	1064	newsegs.append(newseg)
	1065	self.userSegments.extend(newsegs)
	1066	return newsegs
	1067
	1068	def joinVertices(self):
	1069	"""
	1070	Return list of segments connecting all userVertices
	1071	in the order they were given
	1072
	1073	Precon: There must be 3 or more vertices in the userVertices structure
	1074	"""
	1075
	1076	newsegs = []
	1077
	1078	v1 = self.userVertices[0]
	1079	for v2 in self.userVertices[1:]:
	1080	if self.isUserSegmentNew(v1,v2):
	1081	newseg = Segment(v1, v2)
	1082	newsegs.append(newseg)
	1083	v1 = v2
	1084
	1085	#Connect last point to the first
	1086	v2 = self.userVertices[0]
	1087	if self.isUserSegmentNew(v1,v2):
	1088	newseg = Segment(v1, v2)
	1089	newsegs.append(newseg)
	1090
	1091
	1092	#Update list of user segments
	1093	self.userSegments.extend(newsegs)
	1094	return newsegs
	1095
	1096	def normaliseMesh(self,scale, offset, height_scale):
	1097	[xmin, ymin, xmax, ymax] = self.boxsize()
	1098	[attmin0, attmax0] = self.maxMinVertAtt(0)
	1099	#print "[attmin0, attmax0]" ,[attmin0, attmax0]
	1100	[attmin1, attmax1] = self.maxMinVertAtt(1)
	1101	#print [xmin, ymin, xmax, ymax]
	1102	xrange = xmax - xmin
	1103	yrange = ymax - ymin
	1104	if xrange > yrange:
	1105	min,max = xmin, xmax
	1106	else:
	1107	min,max = ymin, ymax
	1108
	1109	for obj in self.getUserVertices():
	1110	obj.x = (obj.x - xmin)/(max- min)*scale + offset
	1111	obj.y = (obj.y - ymin)/(max- min)*scale + offset
	1112	if len(obj.attributes) > 0 and attmin0 != attmax0:
	1113	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
	1114	(attmax0-attmin0)*height_scale
	1115	if len(obj.attributes) > 1 and attmin1 != attmax1:
	1116	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
	1117	(attmax1-attmin1)*height_scale
	1118
	1119	for obj in self.getMeshVertices():
	1120	obj.x = (obj.x - xmin)/(max- min)*scale + offset
	1121	obj.y = (obj.y - ymin)/(max- min)*scale + offset
	1122	if len(obj.attributes) > 0 and attmin0 != attmax0:
	1123	obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
	1124	(attmax0-attmin0)*height_scale
	1125	if len(obj.attributes) > 1 and attmin1 != attmax1:
	1126	obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
	1127	(attmax1-attmin1)*height_scale
	1128
	1129	for obj in self.getHoles():
	1130	obj.x = (obj.x - xmin)/(max- min)*scale + offset
	1131	obj.y = (obj.y - ymin)/(max- min)*scale + offset
	1132	for obj in self.getRegions():
	1133	obj.x = (obj.x - xmin)/(max- min)*scale + offset
	1134	obj.y = (obj.y - ymin)/(max- min)*scale + offset
	1135	[xmin, ymin, xmax, ymax] = self.boxsize()
	1136	#print [xmin, ymin, xmax, ymax]
	1137
	1138	def boxsize(self):
	1139	"""
	1140	Returns a list denoting a box that contains the entire structure of vertices
	1141	Structure: [xmin, ymin, xmax, ymax]
	1142	"""
	1143	# FIXME dsg!!! large is a hack
	1144	#You want the kinds package, part of Numeric:
	1145	#In [2]: import kinds
	1146
	1147	#In [3]: kinds.default_float_kind.M
	1148	#kinds.default_float_kind.MAX kinds.default_float_kind.MIN
	1149	#kinds.default_float_kind.MAX_10_EXP kinds.default_float_kind.MIN_10_EXP
	1150	#kinds.default_float_kind.MAX_EXP kinds.default_float_kind.MIN_EXP
	1151
	1152	#In [3]: kinds.default_float_kind.MIN
	1153	#Out[3]: 2.2250738585072014e-308
	1154
	1155	large = 1e100
	1156	xmin= large
	1157	xmax=-large
	1158	ymin= large
	1159	ymax=-large
	1160	for vertex in self.userVertices:
	1161	if vertex.x < xmin:
	1162	xmin = vertex.x
	1163	if vertex.x > xmax:
	1164	xmax = vertex.x
	1165
	1166	if vertex.y < ymin:
	1167	ymin = vertex.y
	1168	if vertex.y > ymax:
	1169	ymax = vertex.y
	1170	return [xmin, ymin, xmax, ymax]
	1171
	1172	def maxMinVertAtt(self, iatt):
	1173	"""
	1174	Returns a list denoting a box that contains the entire structure of vertices
	1175	Structure: [xmin, ymin, xmax, ymax]
	1176	"""
	1177	# FIXME dsg!!! large is a hack
	1178	#You want the kinds package, part of Numeric:
	1179	#In [2]: import kinds
	1180
	1181	#In [3]: kinds.default_float_kind.M
	1182	#kinds.default_float_kind.MAX kinds.default_float_kind.MIN
	1183	#kinds.default_float_kind.MAX_10_EXP kinds.default_float_kind.MIN_10_EXP
	1184	#kinds.default_float_kind.MAX_EXP kinds.default_float_kind.MIN_EXP
	1185
	1186	#In [3]: kinds.default_float_kind.MIN
	1187	#Out[3]: 2.2250738585072014e-308
	1188
	1189	large = 1e100
	1190	min= large
	1191	max=-large
	1192	for vertex in self.userVertices:
	1193	if len(vertex.attributes) > iatt:
	1194	if vertex.attributes[iatt] < min:
	1195	min = vertex.attributes[iatt]
	1196	if vertex.attributes[iatt] > max:
	1197	max = vertex.attributes[iatt]
	1198	for vertex in self.meshVertices:
	1199	if len(vertex.attributes) > iatt:
	1200	if vertex.attributes[iatt] < min:
	1201	min = vertex.attributes[iatt]
	1202	if vertex.attributes[iatt] > max:
	1203	max = vertex.attributes[iatt]
	1204	return [min, max]
	1205
	1206	def scaleoffset(self, WIDTH, HEIGHT):
	1207	"""
	1208	Returns a list denoting the scale and offset terms that need to be
	1209	applied when converting mesh co-ordinates onto grid co-ordinates
	1210	Structure: [scale, xoffset, yoffset]
	1211	"""
	1212	OFFSET = 0.05*min([WIDTH, HEIGHT])
	1213	[xmin, ymin, xmax, ymax] = self.boxsize()
	1214	SCALE = min([0.9WIDTH, 0.9HEIGHT])/max([xmax-xmin, ymax-ymin])
	1215
	1216	if SCALE*xmin < OFFSET:
	1217	xoffset = abs(SCALE*xmin) + OFFSET
	1218	if SCALE*xmax > WIDTH - OFFSET:
	1219	xoffset= -(SCALE*xmax - WIDTH + OFFSET)
	1220	if SCALE*ymin < OFFSET:
	1221	b = abs(SCALE*ymin)+OFFSET
	1222	if SCALE*ymax > HEIGHT-OFFSET:
	1223	b = -(SCALE*ymax - HEIGHT + OFFSET)
	1224	yoffset = HEIGHT - b
	1225	return [SCALE, xoffset, yoffset]
	1226
	1227	def plotMeshTriangle(self,tag = 0,WIDTH = 400,HEIGHT = 400):
	1228	"""
	1229	Plots all node connections.
	1230	tag = 0 (no node numbers), tag = 1 (node numbers)
	1231	"""
	1232
	1233	try:
	1234	from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
	1235
	1236	[SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
	1237
	1238	root = Tk()
	1239	frame = Frame(root)
	1240	frame.pack()
	1241	button = Button(frame, text="OK", fg="red", command=frame.quit)
	1242	button.pack(side=BOTTOM)
	1243	canvas = Canvas(frame,bg="white", width=WIDTH, height=HEIGHT)
	1244	canvas.pack()
	1245	text = Label(frame, width=20, height=10, text='triangle mesh')
	1246	text.pack()
	1247
	1248	#print self.meshTriangles
	1249	for triangle in self.meshTriangles:
	1250	triangle.draw(canvas,1,
	1251	scale = SCALE,
	1252	xoffset = xoffset,
	1253	yoffset = yoffset )
	1254
	1255	root.mainloop()
	1256
	1257	except:
	1258	print "Unexpected error:", sys.exc_info()[0]
	1259	raise
	1260
	1261	#print """
	1262	#node::plot Failed.
	1263	#Most probably, the Tkinter module is not available.
	1264	#"""
	1265
	1266	def plotUserSegments(self,tag = 0,WIDTH = 400,HEIGHT = 400):
	1267	"""
	1268	Plots all node connections.
	1269	tag = 0 (no node numbers), tag = 1 (node numbers)
	1270	"""
	1271
	1272	try:
	1273	from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
	1274
	1275	[SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
	1276
	1277	root = Tk()
	1278	frame = Frame(root)
	1279	frame.pack()
	1280	button = Button(frame, text="OK", fg="red", command=frame.quit)
	1281	button.pack(side=BOTTOM)
	1282	canvas = Canvas(frame, bg="white", width=WIDTH, height=HEIGHT)
	1283	canvas.pack()
	1284	text = Label(frame, width=20, height=10, text='user segments')
	1285	text.pack()
	1286
	1287	for segment in self.userSegments:
	1288	segment.draw(canvas,SCALE, xoffset, yoffset )
	1289
	1290	root.mainloop()
	1291
	1292	except:
	1293	print "Unexpected error:", sys.exc_info()[0]
	1294	raise
	1295
	1296	#print """
	1297	#node::plot Failed.
	1298	#Most probably, the Tkinter module is not available.
	1299	#"""
	1300
	1301	def exportASCIItrianglulationfile(self,ofile):
	1302	"""
	1303	export a file, ofile, with the format
	1304
	1305	First line: <# of vertices> <# of attributes>
	1306	Following lines: <vertex #> <x> <y> [attributes]
	1307	One line: <# of triangles>
	1308	Following lines: <triangle #> <vertex #> <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
	1309	One line: <# of segments>
	1310	Following lines: <segment #> <vertex #> <vertex #> [boundary marker]
	1311	"""
	1312	fd = open(ofile,'w')
	1313	gen_dict = self.Mesh2MeshList()
	1314	load_mesh.loadASCII.write_ASCII_trianglulation(fd,gen_dict)
	1315	self.writeASCIImesh(fd,
	1316	self.userVertices,
	1317	self.userSegments,
	1318	self.holes,
	1319	self.regions)
	1320	fd.close()
	1321
	1322	def exportASCIIsegmentoutlinefile(self,ofile):
	1323	"""
	1324	export a file, ofile, with no triangulation and only vertices connected to segments.
	1325	"""
	1326	fd = open(ofile,'w')
[364]	1327	meshDict = {}
	1328
	1329	meshDict['generatedpointlist'] = []
	1330	meshDict['generatedpointattributelist'] = []
	1331	meshDict['generatedsegmentlist'] = []
	1332	meshDict['generatedsegmentmarkerlist'] = []
	1333
	1334	meshDict['generatedtrianglelist'] = []
	1335	meshDict['generatedtriangleattributelist'] = []
	1336	meshDict['generatedtriangleneighborlist'] = []
	1337
	1338	load_mesh.loadASCII.write_ASCII_trianglulation(fd,meshDict)
[349]	1339	self.writeASCIIsegmentoutline(fd,
	1340	self.userVertices,
	1341	self.userSegments,
	1342	self.holes,
	1343	self.regions)
	1344	fd.close()
	1345
	1346	def exportASCIIobj(self,ofile):
	1347	"""
	1348	export a file, ofile, with the format
	1349	lines: v <x> <y> <first attribute>
	1350	f <vertex #> <vertex #> <vertex #> (of the triangles)
	1351	"""
	1352	fd = open(ofile,'w')
	1353	self.writeASCIIobj(fd)
	1354	fd.close()
	1355
	1356
	1357	def writeASCIIobj(self,fd):
	1358	fd.write(" # Triangulation as an obj file\n")
	1359	numVert = str(len(self.meshVertices))
	1360
	1361	index1 = 1
	1362	for vert in self.meshVertices:
	1363	vert.index1 = index1
	1364	index1 += 1
	1365
	1366	fd.write("v "
	1367	+ str(vert.x) + " "
	1368	+ str(vert.y) + " "
	1369	+ str(vert.attributes[0]) + "\n")
	1370
	1371	for tri in self.meshTriangles:
	1372	fd.write("f "
	1373	+ str(tri.vertices[0].index1) + " "
	1374	+ str(tri.vertices[1].index1) + " "
	1375	+ str(tri.vertices[2].index1) + "\n")
	1376
	1377	def writeASCIIsegmentoutline(self,
	1378	fd,
	1379	userVertices,
	1380	userSegments,
	1381	holes,
	1382	regions):
	1383	"""Write the user mesh info, only with vertices that are connected to segs
	1384	"""
	1385	verts = []
	1386	#dupindex = 0
	1387	for seg in self.userSegments:
	1388	verts.append(seg.vertices[0])
	1389	verts.append(seg.vertices[1])
	1390	#print 'verts>',verts
	1391
	1392	verts, count = self.removeDuplicatedVertices(verts)
	1393	#print 'verts no dups>',verts
	1394	self.writeASCIImesh(fd,
	1395	verts,
	1396	self.userSegments,
	1397	self.holes,
	1398	self.regions)
[364]	1399
[349]	1400	def exportASCIImeshfile(self,ofile):
	1401	"""
	1402	export a file, ofile, with the format
	1403
	1404	First line: <# of user vertices> <# of attributes>
	1405	Following lines: <x> <y> [attributes]
	1406	One line: <# of segments>
	1407	Following lines: <segment #> <vertex #> <vertex #> [boundary marker]
	1408	"""
	1409
	1410	fd = open(ofile,'w')
	1411	self.writeASCIImesh(fd,
	1412	self.userVertices,
	1413	self.userSegments,
	1414	self.holes,
	1415	self.regions)
	1416	fd.close()
	1417
	1418	def writeASCIImesh(self,
	1419	fd,
	1420	userVertices,
	1421	userSegments,
	1422	holes,
	1423	regions):
	1424	numVert = str(len(userVertices))
	1425	if (numVert == "0"):
	1426	numVertAttrib = "0"
	1427	else:
	1428	numVertAttrib = str(len(userVertices[0].attributes))
	1429	fd.write(numVert + " " + numVertAttrib + " # <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
	1430
	1431	# <x> <y> [attributes]
	1432	index = 0
	1433	for vert in userVertices:
	1434	vert.index = index
	1435	index += 1
	1436	attlist = ""
	1437	for att in vert.attributes:
	1438	attlist = attlist + str(att)+" "
	1439	attlist.strip()
	1440	fd.write(str(vert.index) + " "
	1441	+str(vert.x) + " "
	1442	+ str(vert.y) + " "
	1443	+ attlist + "\n")
	1444
	1445	#One line: <# of segments>
	1446	fd.write(str(len(userSegments)) +
	1447	" # <segment #> <vertex #> <vertex #> [boundary marker] ...Mesh Segments..." + "\n")
	1448
	1449	#Following lines: <vertex #> <vertex #> [boundary marker]
	1450	index = 0
	1451	for seg in userSegments:
	1452	fd.write(str(index) + " "
	1453	+ str(seg.vertices[0].index) + " "
	1454	+ str(seg.vertices[1].index) + " "
	1455	+ str(seg.marker) + "\n")
	1456	index += 1
	1457
	1458	#One line: <# of holes>
	1459	fd.write(str(len(holes)) +
	1460	" # <Hole #> <x> <y> ...Mesh Holes..." + "\n")
	1461	# <x> <y>
	1462	index = 0
	1463	for h in holes:
	1464	fd.write(str(index) + " "
	1465	+ str(h.x) + " "
	1466	+ str(h.y) + "\n")
	1467	index += 1
	1468
	1469	#One line: <# of regions>
	1470	fd.write(str(len(regions)) +
	1471	" # <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")
	1472	# <index> <x> <y> <tag>
	1473	index = 0
	1474	for r in regions:
	1475	fd.write(str(index) + " "
	1476	+ str(r.x) + " "
	1477	+ str(r.y)+ " "
	1478	+ str(r.getTag()) + "\n")
	1479	index += 1
	1480	index = 0
	1481	# <index> [<MaxArea>\|""]
	1482	for r in regions:
	1483	area = r.getMaxArea()
	1484	if area == None:
	1485	area = ""
	1486	else:
	1487	area = str(area)
	1488
	1489	fd.write(str(index) + " " + area + "\n")
	1490	index += 1
	1491
	1492	def exportxyafile(self,ofile):
	1493	"""
	1494	export a file, ofile, with the format
	1495
	1496	First line: <# of vertices> <# of attributes>
	1497	Following lines: <vertex #> <x> <y> [attributes]
	1498	"""
[355]	1499	#load_mesh.loadASCII
	1500
[349]	1501	if self.meshVertices == []:
	1502	Vertices = self.userVertices
	1503	else:
	1504	Vertices = self.meshVertices
	1505
	1506	numVert = str(len(Vertices))
	1507
	1508	if Vertices == []:
	1509	raise RuntimeError
	1510	numVertAttrib = str(len(Vertices[0].attributes))
[355]	1511	title = numVert + " " + numVertAttrib + " # <vertex #> <x> <y> [attributes]"
	1512
	1513	#Convert the Vertices to pointlist and pointattributelist
	1514	xya_dict = {}
	1515	pointattributes = []
	1516	points = []
	1517
	1518	for vert in Vertices:
	1519	points.append([vert.x,vert.y])
	1520	pointattributes.append(vert.attributes)
[349]	1521
[355]	1522	xya_dict['pointlist'] = points
	1523	xya_dict['pointattributelist'] = pointattributes
[349]	1524
[364]	1525	load_mesh.loadASCII.export_xya_file(ofile, xya_dict, title, delimiter = " ")
[404]	1526
	1527	def importUngenerateFile(ofile):
	1528	"""
	1529	import a file, ofile, with the format
	1530	[poly]
	1531	poly format:
	1532	First line: <# of vertices> <x centroid> <y centroid>
	1533	Following lines: <x> <y>
	1534	last line: "END"
[420]	1535
	1536	Note: These are clockwise.
[404]	1537	"""
	1538	fd = open(ofile,'r')
	1539	Dict = readUngenerateFile(fd)
	1540	fd.close()
	1541	return Dict
[349]	1542
[404]	1543	def readUngenerateFile(fd):
	1544	"""
	1545	import a file, ofile, with the format
	1546	[poly]
	1547	poly format:
	1548	First line: <# of polynomial> <x centroid> <y centroid>
	1549	Following lines: <x> <y>
	1550	last line: "END"
	1551	"""
	1552	END_DELIMITER = 'END\n'
	1553
	1554	points = []
	1555	segments = []
	1556
	1557	isEnd = False
	1558	line = fd.readline() #not used <# of polynomial> <x> <y>
	1559	while not isEnd:
	1560	line = fd.readline()
	1561	fragments = line.split()
	1562	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
	1563	points.append(vert)
	1564	PreviousVertIndex = len(points)-1
	1565	firstVertIndex = PreviousVertIndex
	1566
	1567	line = fd.readline() #Read the next line
	1568	while line <> END_DELIMITER:
	1569	#print "line >" + line + "<"
	1570	fragments = line.split()
	1571	vert = [float(fragments.pop(0)),float(fragments.pop(0))]
	1572	points.append(vert)
[431]	1573	thisVertIndex = len(points)-1
	1574	segment = [PreviousVertIndex,thisVertIndex]
[404]	1575	segments.append(segment)
[431]	1576	PreviousVertIndex = thisVertIndex
[404]	1577	line = fd.readline() #Read the next line
	1578	i =+ 1
[431]	1579	# If the last and first segments are the same,
[404]	1580	# Remove the last segment and the last vertex
	1581	# then add a segment from the second last vert to the 1st vert
[431]	1582	thisVertIndex = len(points)-1
	1583	firstVert = points[firstVertIndex]
	1584	thisVert = points[thisVertIndex]
	1585	#print "firstVert",firstVert
	1586	#print "thisVert",thisVert
	1587	if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
	1588	points.pop()
	1589	segments.pop()
	1590	thisVertIndex = len(points)-1
	1591	segments.append([thisVertIndex, firstVertIndex])
[404]	1592
	1593	line = fd.readline() # read <# of polynomial> <x> <y> OR END
	1594	#print "line >>" + line + "<<"
	1595	if line == END_DELIMITER:
	1596	isEnd = True
	1597
	1598	#print "points", points
	1599	#print "segments", segments
	1600	ungenerated_dict = {}
	1601	ungenerated_dict['pointlist'] = points
	1602	ungenerated_dict['segmentlist'] = segments
	1603	return ungenerated_dict
[355]	1604
[349]	1605	def importMeshFromFile(ofile):
	1606	"""returns a mesh object, made from a .xya or .txh file
	1607	Often raises SyntaxError, IOError
	1608	"""
	1609	newmesh = None
	1610	if ofile[-4:]== ".xya":
	1611	print "loading " + ofile
	1612	try:
	1613	dict = load_mesh.loadASCII.load_xya_file(ofile, delimiter = ',')
	1614	except SyntaxError:
	1615	dict = load_mesh.loadASCII.load_xya_file(ofile, delimiter = ' ')
	1616	newmesh= Mesh()
	1617	newmesh.setMesh(dict)
	1618	counter = newmesh.removeDuplicatedUserVertices()
	1619	if (counter >0):
	1620	print "%i duplicate vertices removed from dataset" % (counter)
	1621	elif ofile[-4:]== ".tsh":
	1622	dict = load_mesh.loadASCII.import_trianglulation(ofile)
	1623	newmesh= Mesh()
	1624	newmesh.setMesh(dict)
	1625	newmesh.setTriangulation(dict)
	1626	else:
	1627	raise RuntimeError
	1628	return newmesh
	1629
	1630	def loadPickle(currentName):
	1631	fd = open(currentName)
	1632	mesh = pickle.load(fd)
	1633	fd.close()
	1634	return mesh
	1635
	1636	def square_outline(side_length = 1,up = "top", left = "left", right = "right",
	1637	down = "bottom"):
	1638
	1639	a = Vertex (0,0)
	1640	b = Vertex (0,side_length)
	1641	c = Vertex (side_length,0)
	1642	d = Vertex (side_length,side_length)
	1643
	1644	s2 = Segment(b,d, marker = up)
	1645	s3 = Segment(b,a, marker = left)
	1646	s4 = Segment(d,c, marker = right)
	1647	s5 = Segment(a,c, marker = down)
	1648
	1649	return Mesh(userVertices=[a,b,c,d],
	1650	userSegments=[s2,s3,s4,s5])
	1651
	1652	def region_strings2ints(region_list):
	1653	"""Given a list of (x_int,y_int,marker_string) lists it returns a list of
	1654	(x_int,y_int,marker_int) and a list to convert the marker_int's back to
	1655	the marker_strings
	1656	"""
	1657	# Make sure "" has an index of 0
	1658	region_list.reverse()
	1659	region_list.append((1.0,2.0,""))
	1660	region_list.reverse()
	1661	convertint2string = []
	1662	for i in xrange(len(region_list)):
	1663	convertint2string.append(region_list[i][2])
	1664	if len(region_list[i]) == 4: # there's an area value
	1665	region_list[i] = (region_list[i][0],
	1666	region_list[i][1],i,region_list[i][3])
	1667	elif len(region_list[i]) == 3: # no area value
	1668	region_list[i] = (region_list[i][0],region_list[i][1],i)
	1669	else:
	1670	print "The region list has a bad size"
	1671	# raise an error ..
	1672	raise Error
	1673
	1674	#remove "" from the region_list
	1675	region_list.pop(0)
	1676
	1677	return [region_list, convertint2string]
	1678
	1679	def region_ints2strings(region_list,convertint2string):
	1680	"""Reverses the transformation of region_strings2ints
	1681	"""
	1682	if region_list[0] != []:
	1683	for i in xrange(len(region_list)):
	1684	region_list[i] = [convertint2string[int(region_list[i][0])]]
	1685	return region_list
	1686
	1687	def segment_ints2strings(intlist, convertint2string):
	1688	"""Reverses the transformation of segment_strings2ints """
	1689	stringlist = []
	1690	for x in intlist:
	1691	stringlist.append(convertint2string[x])
	1692	return stringlist
	1693
	1694	def segment_strings2ints(stringlist, preset):
	1695	"""Given a list of strings return a list of 0 to n ints which represent
	1696	the strings and a converting list of the strings, indexed by 0 to n ints.
	1697	Also, input an initial converting list of the strings
	1698	Note, the converting list starts off with
[650]	1699	["internal boundary", "external boundary", "internal boundary"]
[349]	1700	example input and output
	1701	input ["a","b","a","c"],["c"]
	1702	output [[2, 1, 2, 0], ['c', 'b', 'a']]
	1703
[650]	1704	the first element in the converting list will be
	1705	overwritten with "".
	1706	?This will become the third item in the converting list?
[349]	1707
	1708	# note, order the initial converting list is important,
	1709	since the index = the int marker
	1710	"""
	1711	nodups = unique(stringlist)
	1712	# note, order is important, the index = the int marker
	1713	#preset = ["internal boundary", "external boundary"]
	1714	#Remove the preset tags from the list with no duplicates
	1715	nodups = [x for x in nodups if x not in preset]
	1716
	1717	try:
	1718	nodups.remove("") # this has to go to zero
	1719	except ValueError:
	1720	pass
	1721
	1722	# Add the preset list at the beginning of no duplicates
	1723	preset.reverse()
	1724	nodups.extend(preset)
	1725	nodups.reverse()
	1726
	1727
	1728	convertstring2int = {}
	1729	convertint2string = []
	1730	index = 0
	1731	for x in nodups:
	1732	convertstring2int[x] = index
	1733	convertint2string.append(x)
	1734	index += 1
	1735	convertstring2int[""] = 0
	1736
	1737	intlist = []
	1738	for x in stringlist:
	1739	intlist.append(convertstring2int[x])
	1740	return [intlist, convertint2string]
	1741
	1742
	1743	def unique(s):
	1744	"""Return a list of the elements in s, but without duplicates.
	1745
	1746	For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
	1747	unique("abcabc") some permutation of ["a", "b", "c"], and
	1748	unique(([1, 2], [2, 3], [1, 2])) some permutation of
	1749	[[2, 3], [1, 2]].
	1750
	1751	For best speed, all sequence elements should be hashable. Then
	1752	unique() will usually work in linear time.
	1753
	1754	If not possible, the sequence elements should enjoy a total
	1755	ordering, and if list(s).sort() doesn't raise TypeError it's
	1756	assumed that they do enjoy a total ordering. Then unique() will
	1757	usually work in O(N*log2(N)) time.
	1758
	1759	If that's not possible either, the sequence elements must support
	1760	equality-testing. Then unique() will usually work in quadratic
	1761	time.
	1762	"""
	1763
	1764	n = len(s)
	1765	if n == 0:
	1766	return []
	1767
	1768	# Try using a dict first, as that's the fastest and will usually
	1769	# work. If it doesn't work, it will usually fail quickly, so it
	1770	# usually doesn't cost much to try it. It requires that all the
	1771	# sequence elements be hashable, and support equality comparison.
	1772	u = {}
	1773	try:
	1774	for x in s:
	1775	u[x] = 1
	1776	except TypeError:
	1777	del u # move on to the next method
	1778	else:
	1779	return u.keys()
	1780
	1781	# We can't hash all the elements. Second fastest is to sort,
	1782	# which brings the equal elements together; then duplicates are
	1783	# easy to weed out in a single pass.
	1784	# NOTE: Python's list.sort() was designed to be efficient in the
	1785	# presence of many duplicate elements. This isn't true of all
	1786	# sort functions in all languages or libraries, so this approach
	1787	# is more effective in Python than it may be elsewhere.
	1788	try:
	1789	t = list(s)
	1790	t.sort()
	1791	except TypeError:
	1792	del t # move on to the next method
	1793	else:
	1794	assert n > 0
	1795	last = t[0]
	1796	lasti = i = 1
	1797	while i < n:
	1798	if t[i] != last:
	1799	t[lasti] = last = t[i]
	1800	lasti += 1
	1801	i += 1
	1802	return t[:lasti]
	1803
	1804	# Brute force is all that's left.
	1805	u = []
	1806	for x in s:
	1807	if x not in u:
	1808	u.append(x)
	1809	return u
	1810
	1811
	1812
	1813	if __name__ == "__main__":
[404]	1814	#from mesh import *
	1815	# THIS CAN BE DELETED
	1816	m = Mesh()
	1817	dict = importUngenerateFile("ungen_test.txt")
	1818	m.addVertsSegs(dict)
	1819	print m

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

Download in other formats: