import sys from os import sep sys.path.append('..'+sep+'pyvolution') import least_squares alpha = 1 mesh_file = 'merimbula_10785.tsh' point_file = 'meri0.xya' mesh_output_file = 'merimbula_10785_%g.tsh'%alpha print mesh_output_file from shallow_water import Domain from pmesh2domain import pmesh_to_domain_instance from util import file_function, Polygon_function, read_polygon from Numeric import zeros, Float, maximum, minimum #------- # Domain filename = mesh_output_file print 'Creating domain from', filename domain = pmesh_to_domain_instance(filename, Domain) print "Number of triangles = ", len(domain) triangles = domain.triangles vertices = domain.coordinates tri_vertices = domain.vertex_coordinates Bed = domain.quantities['elevation'] bed = Bed.vertex_values from visual import * f = frame() model = faces(frame=f) max_x = max(vertices[:,0]) min_x = min(vertices[:,0]) max_y = max(vertices[:,1]) min_y = min(vertices[:,1]) range_x = max_x - min_x range_y = max_y - min_y range_xy = max(range_x, range_y) print max_x, min_x, max_y, min_y , range_x, range_y, range_xy print 'Calculating triangles' pos = zeros( (6*len(domain),3), Float) colour = zeros( (6*len(domain),3), Float) normals = zeros( (6*len(domain),3), Float) for i in range( len(domain) ): # v0 = vector(2.0*(vertices[triangles[i,0],0]-min_x)/range_xy-1.0, \ # 2.0*(vertices[triangles[i,0],1]-min_y)/range_xy-1.0,bed[i,0]/50.0) # v1 = vector(2.0*(vertices[triangles[i,1],0]-min_x)/range_xy-1.0, \ # 2.0*(vertices[triangles[i,1],1]-min_y)/range_xy-1.0,bed[i,1]/50.0) # v2 = vector(2.0*(vertices[triangles[i,2],0]-min_x)/range_xy-1.0, \ # 2.0*(vertices[triangles[i,2],1]-min_y)/range_xy-1.0,bed[i,2]/50.0) v0 = vector(2.0*(tri_vertices[i,0]-min_x)/range_xy-1.0, \ 2.0*(tri_vertices[i,1]-min_y)/range_xy-1.0,bed[i,0]/50.0) v1 = vector(2.0*(tri_vertices[i,2]-min_x)/range_xy-1.0, \ 2.0*(tri_vertices[i,3]-min_y)/range_xy-1.0,bed[i,1]/50.0) v2 = vector(2.0*(tri_vertices[i,4]-min_x)/range_xy-1.0, \ 2.0*(tri_vertices[i,5]-min_y)/range_xy-1.0,bed[i,2]/50.0) try: normal = norm( cross(v1-v0, v2-v0) ) except: normal = vector(0,0,1) pos[6*i ,:] = v0 pos[6*i+1,:] = v1 pos[6*i+2,:] = v2 pos[6*i+3,:] = v0 pos[6*i+4,:] = v2 pos[6*i+5,:] = v1 colour[6*i ,:] = color.blue colour[6*i+1,:] = color.blue colour[6*i+2,:] = color.blue colour[6*i+3,:] = color.blue colour[6*i+4,:] = color.blue colour[6*i+5,:] = color.blue normals[6*i ,:] = normal normals[6*i+1,:] = normal normals[6*i+2,:] = normal normals[6*i+3,:] = -normal normals[6*i+4,:] = -normal normals[6*i+5,:] = -normal #FacetedTriangle( v0, v1, v2 ) #print v0 #print v1 #print v2 print 'Drawing plot' model.pos = pos model.color = colour model.normal = normals #v0 = vector(0.0,0.0,0.0) #v1 = vector(1.0,0.0,0.1) #v2 = vector(1.0,1.0,0.1)