source: anuga_work/development/steve/get_triangle_data.py @ 3929

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 anuga.pyvolution.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)