"""Script for running a tsunami inundation scenario for Boca do Rio.
Adopted from cairns files

Joaquim Luis, 2007
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------

# Standard modules
import os
import time
import sys

# Related major packages
from anuga.shallow_water import Domain
from anuga.shallow_water import Transmissive_boundary
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import Time_boundary
from anuga.shallow_water import File_boundary
from anuga.shallow_water import Field_boundary

from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
from anuga.shallow_water.data_manager import dem2pts
from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
from anuga.abstract_2d_finite_volumes.util import file_function


#------------------------------------------------------------------------------
# Define scenario as either ...
#------------------------------------------------------------------------------
scenario = 'br'
basename = scenario + 'source'

#------------------------------------------------------------------------------
# Preparation of topographic data
# Convert ASC 2 DEM 2 PTS using source data and store result in source data
#------------------------------------------------------------------------------

# Filenames
dem_name = 'sub_region_bat' 
meshname = 'sub_region.msh'

# Create DEM from asc data
convert_dem_from_ascii2netcdf(dem_name, use_cache=False, verbose=True)

# Create pts file for onshore DEM
dem2pts(dem_name, use_cache=False, verbose=True)

from anuga.utilities.polygon import read_polygon, plot_polygons, \
                                    polygon_area, is_inside_polygon

###############################
# Domain definitions
###############################

# bounding polygon for study area
bounding_polygon = read_polygon('out_rect.csv')

print 'Area of bounding polygon in km?', polygon_area(bounding_polygon)/1000000.0

###############################
# Interior region definitions
###############################

# interior polygons
poly_shallow = read_polygon('shallow.csv')

#------------------------------------------------------------------------------
# Create the triangular mesh based on overall clipping polygon with a tagged
# boundary and interior regions defined in project.py along with
# resolutions (maximal area of per triangle) for each polygon
#------------------------------------------------------------------------------

remainder_res = 2500
shallow_res = 250
interior_regions = [[bounding_polygon, remainder_res], [poly_shallow, shallow_res]]

create_mesh_from_regions(bounding_polygon, 
                         boundary_tags={'west': [0],
                                        'south': [1],
                                        'east': [2],
                                        'north': [3]},
                         maximum_triangle_area=remainder_res,
                         filename=meshname,
                         interior_regions=interior_regions,
                         use_cache=False,
                         verbose=True)

#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------

domain = Domain(meshname, use_cache=False, verbose=True)

print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()
print domain.statistics()

domain.set_name(basename) 
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
domain.set_minimum_storable_height(0.01)

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------

tide = 0.0
domain.set_quantity('stage', tide)
domain.set_quantity('friction', 0.025) 
domain.set_quantity('elevation', 
                    filename=dem_name + '.pts',
                    use_cache=False,
                    verbose=True,
                    alpha=0.1)

#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------

# Create boundary function from timeseries provided in file
#function1 = file_function('mareg.tms', domain, verbose=True)

# Create and assign boundary objects
#Bts1 = Transmissive_Momentum_Set_Stage_boundary(domain, function1) 

print 'Available boundary tags', domain.get_boundary_tags()
#Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])
Bf = Field_boundary('br_swan.sww',
                    domain, time_thinning=1, mean_stage=tide, 
                    use_cache=False, verbose=True)
   
#Boundary condition for sww feed at the east boundary
domain.set_boundary({'west': Br,'south':Bf,'east': Br,'north': Br})
    
#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------

import time
t0 = time.time()

from Numeric import allclose
from anuga.abstract_2d_finite_volumes.quantity import Quantity


    # save every 1sec leading up to wave approaching land
for t in domain.evolve(yieldstep = 1, finaltime = 960): 
        domain.write_time()
        domain.write_boundary_statistics(tags = 'east')
           
print 'That took %.2f seconds' %(time.time()-t0)