source: anuga_work/development/analytical solutions/Sydney.py @ 4631

"""Example of the use of the shallow water wave equation
to simulate a bomb blast in an urban area.

   Copyright 2004
   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
   Geoscience Australia
   
"""

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

from shallow_water import Domain, Reflective_boundary, File_boundary,\
     Dirichlet_boundary, Transmissive_boundary
from pmesh2domain import pmesh_to_domain_instance
from anuga.pyvolution.util import Polygon_function,read_polygon


######################
# Domain
filename = 'sydney_ubd.tsh'
print 'Creating domain from', filename
domain = pmesh_to_domain_instance(filename, Domain)
print 'Number of triangles = ', len(domain)

domain.default_order = 2
domain.smooth = True


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

#Reduction operation for get_vertex_values
from anuga.pyvolution.util import mean
domain.reduction = mean 
#domain.reduction = min  #Looks better near steep slopes


######################
#Initial condition
#
print 'Initial condition'

#Set bed-elevation and friction(None)
def x_slope(x,y):
    n = x.shape[0]
    z = 0*x
    return z

domain.set_quantity('elevation', x_slope)

#Set the initial water stage
def stage(x,y):
    z = x_slope(x,y)
    n = x.shape[0]
    h = 0*x
    return h   

domain.set_quantity('stage', stage)
p0 = [[334429.024416, 6251238.11488],[334428.532417, 6251234.42489],
      [334432.468407, 6251233.93289],[334431.730409, 6251238.36088]]
domain.set_quantity('stage',Polygon_function([(p0,10000.0)]))

############
#Boundary
tags = {}
tags['external'] = Transmissive_boundary(domain)
tags[''] = None
domain.set_boundary(tags)
    

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

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