[2229] | 1 | import sys |
---|
| 2 | from os import sep |
---|
| 3 | sys.path.append('..'+sep+'pyvolution') |
---|
| 4 | |
---|
| 5 | |
---|
| 6 | import least_squares |
---|
| 7 | |
---|
| 8 | alpha = 1 |
---|
| 9 | mesh_file = 'merimbula_10785.tsh' |
---|
| 10 | point_file = 'meri0.xya' |
---|
| 11 | mesh_output_file = 'merimbula_10785_%g.tsh'%alpha |
---|
| 12 | |
---|
| 13 | print mesh_output_file |
---|
| 14 | |
---|
| 15 | from shallow_water import Domain |
---|
| 16 | from pmesh2domain import pmesh_to_domain_instance |
---|
| 17 | from util import file_function, Polygon_function, read_polygon |
---|
| 18 | from Numeric import zeros, Float, maximum, minimum |
---|
| 19 | |
---|
| 20 | #------- |
---|
| 21 | # Domain |
---|
| 22 | filename = mesh_output_file |
---|
| 23 | print 'Creating domain from', filename |
---|
| 24 | domain = pmesh_to_domain_instance(filename, Domain) |
---|
| 25 | print "Number of triangles = ", len(domain) |
---|
| 26 | |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | |
---|
| 30 | |
---|
| 31 | |
---|
| 32 | triangles = domain.triangles |
---|
| 33 | |
---|
| 34 | vertices = domain.coordinates |
---|
| 35 | |
---|
| 36 | tri_vertices = domain.vertex_coordinates |
---|
| 37 | |
---|
| 38 | Bed = domain.quantities['elevation'] |
---|
| 39 | |
---|
| 40 | bed = Bed.vertex_values |
---|
| 41 | |
---|
| 42 | from visual import * |
---|
| 43 | |
---|
| 44 | f = frame() |
---|
| 45 | model = faces(frame=f) |
---|
| 46 | |
---|
| 47 | |
---|
| 48 | max_x = max(vertices[:,0]) |
---|
| 49 | min_x = min(vertices[:,0]) |
---|
| 50 | max_y = max(vertices[:,1]) |
---|
| 51 | min_y = min(vertices[:,1]) |
---|
| 52 | range_x = max_x - min_x |
---|
| 53 | range_y = max_y - min_y |
---|
| 54 | |
---|
| 55 | range_xy = max(range_x, range_y) |
---|
| 56 | |
---|
| 57 | print max_x, min_x, max_y, min_y , range_x, range_y, range_xy |
---|
| 58 | |
---|
| 59 | print 'Calculating triangles' |
---|
| 60 | |
---|
| 61 | pos = zeros( (6*len(domain),3), Float) |
---|
| 62 | colour = zeros( (6*len(domain),3), Float) |
---|
| 63 | normals = zeros( (6*len(domain),3), Float) |
---|
| 64 | |
---|
| 65 | for i in range( len(domain) ): |
---|
| 66 | # v0 = vector(2.0*(vertices[triangles[i,0],0]-min_x)/range_xy-1.0, \ |
---|
| 67 | # 2.0*(vertices[triangles[i,0],1]-min_y)/range_xy-1.0,bed[i,0]/50.0) |
---|
| 68 | # v1 = vector(2.0*(vertices[triangles[i,1],0]-min_x)/range_xy-1.0, \ |
---|
| 69 | # 2.0*(vertices[triangles[i,1],1]-min_y)/range_xy-1.0,bed[i,1]/50.0) |
---|
| 70 | # v2 = vector(2.0*(vertices[triangles[i,2],0]-min_x)/range_xy-1.0, \ |
---|
| 71 | # 2.0*(vertices[triangles[i,2],1]-min_y)/range_xy-1.0,bed[i,2]/50.0) |
---|
| 72 | |
---|
| 73 | |
---|
| 74 | v0 = vector(2.0*(tri_vertices[i,0]-min_x)/range_xy-1.0, \ |
---|
| 75 | 2.0*(tri_vertices[i,1]-min_y)/range_xy-1.0,bed[i,0]/50.0) |
---|
| 76 | v1 = vector(2.0*(tri_vertices[i,2]-min_x)/range_xy-1.0, \ |
---|
| 77 | 2.0*(tri_vertices[i,3]-min_y)/range_xy-1.0,bed[i,1]/50.0) |
---|
| 78 | v2 = vector(2.0*(tri_vertices[i,4]-min_x)/range_xy-1.0, \ |
---|
| 79 | 2.0*(tri_vertices[i,5]-min_y)/range_xy-1.0,bed[i,2]/50.0) |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | try: |
---|
| 83 | normal = norm( cross(v1-v0, v2-v0) ) |
---|
| 84 | except: |
---|
| 85 | normal = vector(0,0,1) |
---|
| 86 | |
---|
| 87 | pos[6*i ,:] = v0 |
---|
| 88 | pos[6*i+1,:] = v1 |
---|
| 89 | pos[6*i+2,:] = v2 |
---|
| 90 | pos[6*i+3,:] = v0 |
---|
| 91 | pos[6*i+4,:] = v2 |
---|
| 92 | pos[6*i+5,:] = v1 |
---|
| 93 | |
---|
| 94 | colour[6*i ,:] = color.blue |
---|
| 95 | colour[6*i+1,:] = color.blue |
---|
| 96 | colour[6*i+2,:] = color.blue |
---|
| 97 | colour[6*i+3,:] = color.blue |
---|
| 98 | colour[6*i+4,:] = color.blue |
---|
| 99 | colour[6*i+5,:] = color.blue |
---|
| 100 | |
---|
| 101 | normals[6*i ,:] = normal |
---|
| 102 | normals[6*i+1,:] = normal |
---|
| 103 | normals[6*i+2,:] = normal |
---|
| 104 | normals[6*i+3,:] = -normal |
---|
| 105 | normals[6*i+4,:] = -normal |
---|
| 106 | normals[6*i+5,:] = -normal |
---|
| 107 | |
---|
| 108 | #FacetedTriangle( v0, v1, v2 ) |
---|
| 109 | #print v0 |
---|
| 110 | #print v1 |
---|
| 111 | #print v2 |
---|
| 112 | |
---|
| 113 | print 'Drawing plot' |
---|
| 114 | model.pos = pos |
---|
| 115 | model.color = colour |
---|
| 116 | model.normal = normals |
---|
| 117 | |
---|
| 118 | #v0 = vector(0.0,0.0,0.0) |
---|
| 119 | #v1 = vector(1.0,0.0,0.1) |
---|
| 120 | #v2 = vector(1.0,1.0,0.1) |
---|
| 121 | |
---|
| 122 | |
---|