"""Validation study of Merimbula lake using Pyvolution.
Example of shallow water wave equation applied to
Malpasset dam break simulation.

   Copyright 2004
   Christopher Zoppou, Stephen Roberts
   Australian National University

Specific methods pertaining to the 2D shallow water equation
are imported from shallow_water
for use with the generic finite volume framework

Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
numerical vector named conserved_quantities.

Existence of file 'Hobart_mesh.tsh' is assumed.
"""

###############################
# Setup Path and importing modules
import sys
from os import sep, path
sys.path.append('..'+sep+'pyvolution')

import Numeric
from shallow_water import Domain, Reflective_boundary, File_boundary,\
     Dirichlet_boundary, Transmissive_boundary
from pmesh2domain import pmesh_to_domain_instance
from util import Polygon_function, read_polygon

######################
# Domain
filename = 'hobart_mesh.tsh'
yieldstep = 2
finaltime = 600

print 'Creating domain from', filename
domain = pmesh_to_domain_instance(filename, Domain)
print "Number of triangles = ", len(domain)

domain.default_order = 1
domain.smooth = True
domain.visualise = True

#------------------------------
# Boundary Conditions
tags = {}
tags['external'] = Reflective_boundary(domain)
domain.set_boundary(tags)

#---------------------------------------------------------
#Decide which quantities are to be stored at each timestep
domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']

#-----------------
#Initial condition

p0 = read_polygon('Hobart_recent_source_zone.xya')
domain.set_quantity('stage',Polygon_function([(p0,10.0)]))

stage_vv = domain.quantities['stage'].vertex_values
elevation_vv = domain.quantities['elevation'].vertex_values

new_stage_vv = stage_vv + elevation_vv

domain.set_quantity('stage', new_stage_vv)

#-------------------------------------
# Provide file name for storing output
domain.store = True
domain.format = 'sww'
domain.filename = 'Hobart_first_order'

#----------------------------
# Friction
domain.set_quantity('friction', 0.033)

######################
#Evolution
import time
t0 = time.time()
for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    domain.write_time()

print 'That took %.2f seconds' %(time.time()-t0)