source: anuga_work/development/river2d_integration/River2Dcdg_ANUGAMesh.py @ 6367

""" This script takes a River2D msh file of format:

 node no.,x,y,z,roughness
 no more nodes.
 Vert1,Vert2,Vert3,Tag,TagValue
 no more boundary segments.
 
 
and converts it to a format where it can be used to creat a domain from Points and Vertices.

All you need to do is specify the data file name and output file name


"""
import string

infid = open('Riv2dMesh2.cdg') # input file is a River 2D run file
outfid1 = open('ANUGA_MeshPts.csv', 'w') # output filename
outfid2 = open('ANUGA_MeshVertices.csv', 'w') # output filename
print "\nFILES are open ready for conversion....\n"

lines = infid.readlines()

# Ignore all metadata prior to node information
for i, line in enumerate(lines):
    if line.strip().startswith('Node #'):
        # Found node information
        break


reading_points = True
reading_vertices = False
last_point_index = 0

for L in lines[i+1:]:
    line = L.strip()

    print line
    
    if line == '':
        continue

    if line.startswith('Element Information'):
        continue
        
    if line.startswith('Element #'):
        reading_vertices = True
        reading_points = False
        continue

    if line.startswith('Boundary Element #'):
        print 'Reached boundary info, not processing yet'
        break

    if reading_points is True:   # Read Points
        fields = line.split()

        #print fields, len(fields)

        if len(fields) == 9:
            # No single character flag
            x_index = 1
        elif len(fields) == 10:
            # Single character flag used in in second column
            x_index = 2
        else:
            print fields, len(fields)
            raise Exception, 'Bad file format'
            
        # Check that point indices are consecutive        
        s = fields[0]
        index = int(s)
        assert index-1 == last_point_index
        last_point_index = index

        # Get points
        x=float(fields[x_index])
        y=float(fields[x_index + 1])
        z=float(fields[x_index + 2]) # Elevation
        n=float(fields[x_index + 3]) # Mannings roughness
        outfid1.write('%.3f, %.3f\n' % (x, y))

    elif reading_vertices is True:    # Read Elements
        fields = line.split()

        # Create indices and make them start from 0
        pt1 = int(fields[3]) - 1
        pt2 = int(fields[4]) - 1
        pt3 = int(fields[5]) - 1
        outfid2.write('%d, %d, %d\n' % (pt1, pt2,pt3)) 

    
outfid1.close()
outfid2.close()
print '\nConversion complete.'
raw_input('Press ENTER to exit.')