source: inundation/ga/storm_surge/pmesh/load_mesh/loadASCII.py @ 511

Last change on this file since 511 was 511, checked in by duncan, 20 years ago

remove debug comments

File size: 16.6 KB
Line 
1
2from string import  find, rfind
3
4def mesh_file_to_mesh_dictionary(fileName):
5    """Load a pmesh file.  Returning the mesh dictionary.
6    """
7    try:
8        meshdic = import_trianglulation(fileName)
9    except IOError, e:       
10        msg = 'Could not open file %s ' %fileName
11        raise IOError, msg
12    return meshdic
13
14def import_trianglulation(ofile):
15    """
16    reading triangulation and meshoutline info.
17    import a file, ofile, with the format
18   
19    First line:  <# of vertices> <# of attributes>
20    Following lines:  <vertex #> <x> <y> [attributes]
21    One line:  <# of triangles>
22    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
23    One line:  <# of segments>
24    Following lines:  <segment #> <vertex #>  <vertex #> [boundary marker]
25    Note: This might throw a can't load file error
26    """
27    fd = open(ofile,'r')
28    dict = read_trianglulation(fd)
29    dict_mesh = read_mesh(fd)
30    for element in dict_mesh.keys():
31        dict[element] = dict_mesh[element]
32       
33    fd.close()
34    return dict
35
36
37def read_trianglulation(fd):
38    delimiter = " "
39    ######### loading the point info
40    line = fd.readline()
41    #print line
42    fragments = line.split()
43    #for fragment in fragments:
44    #    print fragment
45    if fragments ==[]:
46        NumOfVertices = 0
47        NumOfVertAttributes = 0
48    else:
49        NumOfVertices = fragments[0]
50        NumOfVertAttributes = fragments[1]
51    points = []
52    pointattributes = []
53    for index in range(int(NumOfVertices)):       
54        #print index
55        fragments = fd.readline().split()
56        #print fragments
57        fragments.pop(0) #pop off the index
58       
59        # pop the x y off so we're left with a list of attributes       
60        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
61        points.append(vert)
62        apointattributes  = []
63        #print fragments
64        for fragment in fragments:
65            apointattributes.append(float(fragment))
66        pointattributes.append(apointattributes)
67       
68    ######### loading the point title info
69    line = fd.readline()
70    #print "point title comments",line
71    vertTitle = []
72    for index in range(int(NumOfVertAttributes)):       
73        #print index
74        fragments = fd.readline().strip() 
75        vertTitle.append(fragments)
76        #print vertTitle
77       
78    ######### loading the triangle info
79    line = fd.readline()
80    #print "triangle comments",line
81    fragments = line.split()
82    #for fragment in fragments:
83    #    print fragment
84    NumOfTriangles = fragments[0]
85    triangles = []
86    triangleattributes = []
87    triangleneighbors = []
88    for index in range(int(NumOfTriangles)):       
89        #print index
90        line = fd.readline()
91        line.strip() # so we can get the region string
92        fragments = line.split()
93        #print "triangle info", fragments
94        fragments.pop(0) #pop off the index
95       
96        tri = [int(fragments[0]),int(fragments[1]),int(fragments[2])]
97        triangles.append(tri)
98        neighbors = [int(fragments[3]),int(fragments[4]),int(fragments[5])]
99        triangleneighbors.append(neighbors)
100        for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
101            line = line[find(line,delimiter):] # remove index
102            line = line.lstrip()
103        stringmarker = line.strip()
104        triangleattributes.append([stringmarker])
105       
106    ######### loading the segment info
107    line = fd.readline()
108    #print "seg comment line",line
109    fragments = line.split()
110    #for fragment in fragments:
111    #    print fragment
112    NumOfSegments = fragments[0]
113    segments = []
114    segmentmarkers = []
115    for index in range(int(NumOfSegments)):       
116        #print index
117        line = fd.readline()
118        line.strip() # to get the segment string
119        fragments = line.split()
120        #print fragments
121        fragments.pop(0) #pop off the index
122        seg = [int(fragments[0]),int(fragments[1])]
123        segments.append(seg)
124        line = line[find(line,delimiter):] # remove index
125        line = line.lstrip()
126        line = line[find(line,delimiter):] # remove x
127        line = line.lstrip()
128        line = line[find(line,delimiter):] # remove y
129        stringmarker = line.strip()
130        segmentmarkers.append(stringmarker)
131    meshDict = {}
132    meshDict['generatedpointlist'] = points
133    meshDict['generatedpointattributelist'] = pointattributes
134    meshDict['generatedtrianglelist'] = triangles
135    meshDict['generatedtriangleattributelist'] = triangleattributes
136    meshDict['generatedtriangleneighborlist'] = triangleneighbors
137    meshDict['generatedsegmentlist'] = segments
138    meshDict['generatedsegmentmarkerlist'] = segmentmarkers
139    meshDict['generatedpointattributetitlelist'] = vertTitle
140    return meshDict
141   
142def import_mesh(ofile):
143    """
144    import a file, ofile, with the format
145   
146    First line:  <# of vertices> <# of attributes>
147    Following lines: <x> <y> [attributes]
148    One line:  <# of segments>
149    Following lines:  <vertex #>  <vertex #> [boundary marker]
150    Note: This might throw a can't load file error
151    """
152    fd = open(ofile,'r')
153    meshDict = read_mesh(fd)
154    fd.close()
155    return meshDict
156
157def read_mesh(fd):
158    """
159    Note, if a file has no mesh info, it can still be read - the meshdic
160    returned will be 'empty'.
161    """
162    delimiter = " " # warning: split() calls are using default whitespace
163   
164    ######### loading the point info
165    line = fd.readline()
166    #print line
167    fragments = line.split()
168    #for fragment in fragments:
169    if fragments ==[]:
170        NumOfVertices = 0
171        NumOfVertAttributes = 0
172    else:
173        NumOfVertices = fragments[0]
174        NumOfVertAttributes = fragments[1]
175    points = []
176    pointattributes = []
177    for index in range(int(NumOfVertices)):       
178        #print index
179        fragments = fd.readline().split() 
180        #print fragments
181        fragments.pop(0) #pop off the index
182        # pop the x y off so we're left with a list of attributes
183        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
184        points.append(vert)
185        apointattributes  = []
186        #print fragments
187        for fragment in fragments:
188            apointattributes.append(float(fragment))
189        pointattributes.append(apointattributes)
190       
191
192    ######### loading the segment info
193    line = fd.readline()
194    #print line
195    fragments = line.split()
196    #for fragment in fragments:
197    #    print fragment
198    if fragments ==[]:
199        NumOfSegments = 0
200    else:
201        NumOfSegments = fragments[0]
202    segments = []
203    segmentmarkers = []
204    for index in range(int(NumOfSegments)): 
205        #print index
206        line = fd.readline()
207        fragments = line.split()
208        #print fragments
209        fragments.pop(0) #pop off the index
210       
211        seg = [int(fragments[0]),int(fragments[1])]
212        segments.append(seg)
213        line = line[find(line,delimiter):] # remove index
214        line = line.lstrip()
215        line = line[find(line,delimiter):] # remove x
216        line = line.lstrip()
217        line = line[find(line,delimiter):] # remove y
218        stringmarker = line.strip()
219        segmentmarkers.append(stringmarker) 
220
221    ######### loading the hole info
222    line = fd.readline()
223    #print line
224    fragments = line.split()
225    #for fragment in fragments:
226    #    print fragment
227    if fragments ==[]:
228        numOfHoles = 0
229    else:
230        numOfHoles = fragments[0]
231    holes = []
232    for index in range(int(numOfHoles)):       
233        #print index
234        fragments = fd.readline().split()
235        #print fragments
236        fragments.pop(0) #pop off the index
237        hole = [float(fragments[0]),float(fragments[1])]
238        holes.append(hole)
239
240   
241    ######### loading the region info
242    line = fd.readline()
243    #print line
244    fragments = line.split()
245    #for fragment in fragments:
246    #    print fragment
247    if fragments ==[]:
248        numOfRegions = 0
249    else:
250        numOfRegions = fragments[0]
251    regions = []
252    regionattributes = []
253    for index in range(int(numOfRegions)):
254        line = fd.readline()
255        fragments = line.split()
256        #print fragments
257        fragments.pop(0) #pop off the index
258        region = [float(fragments[0]),float(fragments[1])]
259        regions.append(region)
260
261        line = line[find(line,delimiter):] # remove index
262        line = line.lstrip()
263        line = line[find(line,delimiter):] # remove x
264        line = line.lstrip()
265        line = line[find(line,delimiter):] # remove y
266        stringmarker = line.strip()
267        regionattributes.append(stringmarker)
268    regionmaxareas = []
269    for index in range(int(numOfRegions)): # Read in the Max area info
270        line = fd.readline()
271        fragments = line.split()
272        #print fragments
273        fragments.pop(0) #pop off the index
274        if len(fragments) == 0: #no max area
275            regionmaxareas.append(None)
276        else:
277            regionmaxareas.append(float(fragments[0]))
278
279       
280    meshDict = {}
281    meshDict['pointlist'] = points
282    meshDict['pointattributelist'] = pointattributes
283    meshDict['segmentlist'] = segments
284    meshDict['segmentmarkerlist'] = segmentmarkers
285    meshDict['holelist'] = holes
286    meshDict['regionlist'] = regions
287    meshDict['regionattributelist'] = regionattributes
288    meshDict['regionmaxarealist'] = regionmaxareas
289
290    return meshDict
291
292def clean_line(line,delimiter):     
293    """Remove whitespace
294    """
295    #print ">%s" %line
296    line = line.strip()
297    #print "stripped>%s" %line
298    numbers = line.split(delimiter)
299    i = len(numbers) - 1
300    while i >= 0:
301        if numbers[i] == '':
302            numbers.pop(i)
303        i += -1
304    #for num in numbers:
305    #    print "num>%s<" %num
306    return numbers
307
308def write_ASCII_trianglulation(fd,
309                               gen_dict):
310    vertices = gen_dict['generatedpointlist']
311    vertices_attributes = gen_dict['generatedpointattributelist']
312    try:
313        vertices_attribute_titles = gen_dict['generatedpointattributetitlelist']
314   
315    except KeyError, e:
316        #FIXME is this the best way?
317        if vertices_attributes == [] or vertices_attributes[0] == []:
318             vertices_attribute_titles = []
319        else:
320            raise KeyError, e
321       
322    triangles = gen_dict['generatedtrianglelist']
323    triangles_attributes = gen_dict['generatedtriangleattributelist']
324    triangle_neighbors = gen_dict['generatedtriangleneighborlist']
325       
326    segments = gen_dict['generatedsegmentlist']
327    segment_markers = gen_dict['generatedsegmentmarkerlist']
328     
329    numVert = str(len(vertices))
330    if (numVert == "0"):
331        numVertAttrib = "0"
332    else:
333        numVertAttrib = str(len(vertices_attributes[0]))
334    fd.write(numVert + " " + numVertAttrib + " # <vertex #> <x> <y> [attributes] ...Triangulation Vertices..." + "\n")
335
336    #<vertex #> <x> <y> [attributes]   
337    index = 0 
338    for vert in vertices:
339        attlist = ""
340        for att in vertices_attributes[index]:
341            attlist = attlist + str(att)+" "
342        attlist.strip()
343        fd.write(str(index) + " "
344                 + str(vert[0]) + " "
345                 + str(vert[1]) + " "
346                 + attlist + "\n")
347        index += 1
348
349    # write comments for title
350    fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
351    for title in vertices_attribute_titles:
352        fd.write(title + "\n")
353       
354    #<# of triangles>
355    n = len(triangles)
356    fd.write(str(n) + " # <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
357       
358    # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
359    for index in range(n):
360        neighbors = ""
361        tri = triangles[index]
362        for neighbor in triangle_neighbors[index]:
363            if neighbor:
364                neighbors += str(neighbor) + " "
365            else:
366                if neighbor == 0:
367                    neighbors +=  "0 "
368                else:
369                    neighbors +=  "-1 "
370        if triangles_attributes[index] == ['']:
371            att = ""
372        else:
373            att = str(triangles_attributes[index])
374        fd.write(str(index) + " "
375                 + str(tri[0]) + " " 
376                 + str(tri[1]) + " " 
377                 + str(tri[2]) + " " 
378                 + neighbors + " "
379                 + att + "\n")
380           
381    #One line:  <# of segments>
382    fd.write(str(len(segments)) + 
383             " # <segment #> <vertex #>  <vertex #> [boundary marker] ...Triangulation Segments..." + "\n")
384       
385    #Following lines:  <segment #> <vertex #>  <vertex #> [boundary marker]
386    for i in range(len(segments)):
387        seg = segments[i]
388        fd.write(str(i) + " "
389                 + str(seg[0]) + " " 
390                 + str(seg[1]) + " " 
391                 + str(segment_markers[i]) + "\n")
392
393
394def export_trianglulation_file(ofile,gen_dict):
395    """
396    write a file, ofile, with the format
397   
398    First line:  <# of vertices> <# of attributes>
399    Following lines:  <vertex #> <x> <y> [attributes]
400    One line:  <# of triangles>
401    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
402    One line:  <# of segments>
403    Following lines:  <segment #> <vertex #>  <vertex #> [boundary marker]
404    """
405    try:
406        fd = open(ofile,'w')
407        write_ASCII_trianglulation(fd,gen_dict)
408        fd.close()
409    except IOError, e:       
410        msg = 'Could not write file %s ' %fileName
411        raise IOError, msg
412
413###
414#  LOADING XYA FILES
415###
416 
417           
418def load_xya_file(ofile,delimiter = ','):
419    """
420    load a file, ofile, with the format
421    x,y, [attributes]
422    """
423    try:
424        fd = open(ofile)
425        xya_dic = read_xya_file(fd, delimiter)
426        fd.close()
427    except IOError, e:       
428        msg = 'Could not open file %s ' %fileName
429        raise IOError, msg
430    return xya_dic
431
432def read_xya_file(fd, delimiter):
433    lines = fd.readlines()
434    points = []
435    pointattributes = []
436    if len(lines) <= 1:
437        raise SyntaxError
438    title = lines.pop(0) # the first (title) line
439    attLength = len(clean_line(lines[0],delimiter))-2
440   
441    #print "initlegth"
442    #print attLength
443    for line in lines:
444        #print "line >%s" %line
445        numbers = clean_line(line,delimiter)
446        #print "numbers >%s<" %numbers
447        if len(numbers) < 2 and numbers != []:
448            raise SyntaxError
449        if numbers != []:
450            try:
451                x = float(numbers[0])
452                y = float(numbers[1])
453                points.append([x,y])
454                numbers.pop(0)
455                numbers.pop(0)
456                attributes = []
457                if attLength != len(numbers):
458                    raise SyntaxError
459                   
460                attLength = len(numbers)
461
462                for num in numbers:
463                    num.strip()
464                    if num != '\n' and num != '':
465                        attributes.append(float(num))
466            except ValueError:
467                raise SyntaxError
468            pointattributes.append(attributes)
469    xya_dict = {}
470    xya_dict['pointlist'] = points
471    xya_dict['pointattributelist'] = pointattributes
472    xya_dict['title'] = title
473    xya_dict['segmentlist'] = []
474    xya_dict['segmentmarkerlist'] = []
475    xya_dict['regionlist'] = []
476    xya_dict['regionattributelist'] = []
477    xya_dict['regionmaxarealist'] = []
478    xya_dict['holelist'] = []
479   
480    return xya_dict
481
482               
483def export_xya_file( file_name, xya_dict, title, delimiter = ','):
484    """
485    export a file, ofile, with the format
486   
487    First line: Title variable
488    Following lines:  <vertex #> <x> <y> [attributes]
489
490    file_name - the name of the new file
491    xya_dict - point and point attribute info in a dictionary
492    title - info to write in the first line
493    """
494    #FIXME, move the test for this from meshharness to loadasciiharness
495    points = xya_dict['pointlist'] 
496    pointattributes = xya_dict['pointattributelist']
497   
498    fd = open(file_name,'w')
499   
500    fd.write(title+"\n")
501    #<vertex #> <x> <y> [attributes]
502    for vert, vertatts in map(None, points, pointattributes):
503        attlist = ""
504        for att in vertatts:
505            attlist = attlist + str(att)+ delimiter
506        attlist = attlist[0:-len(delimiter)] # remove the last delimiter
507        attlist.strip()
508        fd.write( str(vert[0]) + delimiter
509                  + str(vert[1]) + delimiter
510                  + attlist + "\n")
511    fd.close()
512
513if __name__ == "__main__":
514    m = import_mesh("tee.txt")
515    print m
Note: See TracBrowser for help on using the repository browser.