source: inundation/ga/storm_surge/pyvolution/realtime_visualisation.py @ 643

from visual import *

def normal(v1,v2):
    """Safe computation of normalised cross product
    """
    try:
        n = norm(cross(v1, v2))
    except ZeroDivisionError:
        print v1, v2
        raise Exception
        #n = vector(0,0,0)

    return n    


class Triangle:
    def __init__(self,v0,v1,v2,color=(1,1,1), frame=None):
        """Create one flat triangular panel with two sides
        """

        #Outward normal
        
        self.normal=normal(v1-v0, v2-v0)

        self.panel = faces(frame = frame)
        self.backpanel = faces(frame = frame)        
        for v in (v0,v1,v2):
            self.panel.append(pos=v, normal = self.normal, color=color)

        for v in (v2,v1,v0):
            self.backpanel.append(pos=v, normal = -self.normal, color=color)            

            
    def move(self,v):
        """Move panel in direction given by vector v
        """

        for vertex in self.panel.pos:
            vertex += v            

    def update_height(self,d,fix_baseline=False):
        """Change height of triangle by displacing the third vertex by
        d in the direction perpendicular to the baseline (v1-v0) and the
        outward normal vector
        """

        v0 = self.panel.pos[0]
        v1 = self.panel.pos[1]
        v2 = self.panel.pos[2]        
        
        n = normal(v1-v0, self.normal)

        if fix_baseline:
            self.panel.pos[2] -= d*n
        else:
            self.panel.pos[:2] += d*n        

    def update_color(self, c):
        """Change color to c
        """

        self.panel.color = c
        

    def set_vertexheights(self, heights, floor_heights = None):
        from Numeric import zeros, Float

        from config import minimum_allowed_height as hmin
        if floor_heights is None:
            floor_heights = zeros(heights.shape, Float)

        all_vertices_below_threshold = True    
        for k in range(3):
            w = heights[k]
            z = floor_heights[k]

            if w-z >= hmin:
                all_vertices_below_threshold = False
            
            vertex = self.panel.pos[k]
            vertex[2] = w

            vertex = self.backpanel.pos[2-k]
            vertex[2] = w


        #Update color to visualise dry areas 
        if all_vertices_below_threshold: 
            self.panel.color = self.bottom_color
            self.backpanel.color = self.bottom_color                
        else:
            self.panel.color = self.top_color
            self.backpanel.color = self.top_color


        #update normal    
        v0 = self.panel.pos[0]
        v1 = self.panel.pos[1]
        v2 = self.panel.pos[2]


        n = normal(v1-v0, v2-v0)
        self.panel.normal=n
        self.backpanel.normal=-n

        


def create_surface(domain):
    """Link surface of domains to their visual representations
    """

    fr = frame() #Default frame to contain all objects
    s=0

    Q = domain.quantities['level']
    try:
        Z = domain.quantities['elevation']
        bed = True
    except:
        bed = False
        
    N = Q.vertex_values.shape[0]

    domain.visuals = []
    for i in range(N):

        z0 = Q.vertex_values[i, 0]
        z1 = Q.vertex_values[i, 1]
        z2 = Q.vertex_values[i, 2]        

        x0, y0 = domain.get_vertex_coordinate(i,0)
        x1, y1 = domain.get_vertex_coordinate(i,1)
        x2, y2 = domain.get_vertex_coordinate(i,2)        
        
        V0 = vector(x0, y0, z0)
        V1 = vector(x1, y1, z1)
        V2 = vector(x2, y2, z2)


        #Top surface
        c0 = 0.1
        c1 = 0.4
        c2 = 0.99
        if s:
            c2 = 0.7*c2   #To show triangles in slightly different shades

        s = 1-s

        col = (c0,c1,c2)
            
        visual_top = Triangle(V0,V1,V2,color=col,frame=fr)
        visual_top.top_color = col

    
        #Bottom surface
        #v0, v1, v2 = volume.vertices
        #v0 = vector(v0.x, v0.y, v0.field_values[index])
        #v1 = vector(v1.x, v1.y, v1.field_values[index])
        #v2 = vector(v2.x, v2.y, v2.field_values[index])
        if bed:
            z0 = Z.vertex_values[i, 0]
            z1 = Z.vertex_values[i, 1]
            z2 = Z.vertex_values[i, 2]
        else:
            z0 = z1 = z2 = 0.0

        V0 = vector(x0, y0, z0)
        V1 = vector(x1, y1, z1)
        V2 = vector(x2, y2, z2)        

        c0 = 0.3
        c1 = 0.3
        c2 = 0.3
        if s:
            c2 = 0.4*c2   #To show triangles in slightly different shades

        col = (c0,c1,c2)   
        visual_bottom = Triangle(V0, V1, V2, color=col,frame=fr)
        visual_top.bottom_color=col
        
        domain.visuals.append( (visual_top, visual_bottom) )
        
    update(domain)
    #print 'Scale', scene.scale    


def update(domain):    
    """Update vertex heights.
    The argument index refers to which conserved quantity to visualise.    
    If domain.smooth is set True, vertex values will be averaged
    yielding a smoother surface.
    """

    from Numeric import array


    Q = domain.quantities['level']
    N = Q.vertex_values.shape[0]
    
    #print scene.forward
    #FIXME: Use smoother from pyvolution instead
    if domain.smooth:
        #Get all average point values
        vertex_heights = {}

        for k in range(N):
            for i in range(3):
                vertex = domain.triangles[k, i]
                if vertex_heights.has_key(vertex):
                    vertex_heights[vertex].append(
                        Q.vertex_values[k, i])
                else:
                    vertex_heights[vertex] = []
                    vertex_heights[vertex].append(
                        Q.vertex_values[k, i])



    for k in range(N):    
        if domain.smooth:
            #The averages
            x = zeros(3, Float)
            for i in range(3):            
                vertex = domain.triangles[k, i]
                A = array(vertex_heights[vertex])
                x[i] = sum(A)/len(A)
        else:        
            #The true values
            x = [Q.vertex_values[k, 0],
                 Q.vertex_values[k, 1],
                 Q.vertex_values[k, 2]]

            
        #Do it
        floor_heights = array([pos[2] for pos in domain.visuals[k][1].panel.pos])
        
        domain.visuals[k][0].set_vertexheights(x, floor_heights)




scene.width = 1000
scene.height = 800

#Original
scene.center = (0.5,0.5,0)
scene.forward = vector(-0.1, 0.5, -0.5)

#Temporary (for bedslope)
#scene.forward = vector(0.0006, 0.7, -0.03)    


#Temporary for hackett - begin
#scene.autoscale = 0
#scene.scale = (0.002, 0.002, 0.01) #Scale z so that countours stand out more
#scene.center = (300.0,500.0,-10)
#Temporary for hackett - end

    
scene.ambient = 0.4
scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4), (-0.1, 0.1, -0.4),
                (-0.2, 0.2, 0.1)]