"""Simple water flow example using ANUGA

Water flowing down a channel with more complex topography
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
from anuga.shallow_water import Domain
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import Time_boundary
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
from math import tan, sqrt, sin, pi

#------------------------------------------------------------------------------
# Project file
#------------------------------------------------------------------------------
#Set up file structure
anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep
meshes_dir = anuga_dir+'meshes'+sep
meshname = 'circular_mesh.msh'
output_dir = anuga_dir+'output'+sep


##    #------------------------------------------------------------------------------
##    # Copy scripts to time stamped output directory and capture screen
##    # output to file
##    #------------------------------------------------------------------------------
##    print "Processor Name:",get_processor_name()
##
##    #copy script must be before screen_catcher
##    #print kwargs
##
##    print 'output_dir', output_dir
##    if myid == 0:
##        copy_code_files(kwargs['output_dir'],__file__, 
##                 dirname(project.__file__)+sep+ project.__name__+'.py' )
##
##        store_parameters(**kwargs)
##
##    barrier()
##
##    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
##
##    print "Processor Name:",get_processor_name()
##
##    # filenames
###    meshes_dir_name = project.meshes_dir_name+'.msh'
##
##    # creates copy of code in output dir
##    print 'min triangles', project.trigs_min,
##    print 'Note: This is generally about 20% less than the final amount'


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


length = 30.
width = 25.
Cx = 12.96                  # centre of island on the x axis
Cy = 13.8                   # centre of island on the y axis
dx = dy = 1                 # Resolution: Length of subdivisions on both axes
water_depth = 0.32          # Can be 0.32 or 0.42

#boundary
poly_domain = [[0,0],[length,0],[length,width],[0,width]]


# exporting asc grid
xmin = 0
xmax = length   
ymin = 0
ymax = width

#Create interior region
Dome = [[(Cx)-4,(Cy)-4],[(Cx)+4,(Cy)-4,],
        [(Cx)+4,(Cy)+4],[(Cx)-4,(Cy)+4]]

remainder_res=1
Dome_res = .01

interior_dome = [[Dome, Dome_res]]

#Create mesh

print 'start create mesh from regions'
 
create_mesh_from_regions(poly_domain,
                         boundary_tags={'wavemaker': [0], 'right': [1],
                                        'top': [2], 'left': [3]},
                         maximum_triangle_area = remainder_res,
                         filename=meshname,
                         interior_regions = interior_dome,
                         use_cache=False,
                         verbose=False)

# Setup computational domain

domain = Domain(meshes_dir+sep+meshname, use_cache=False, verbose = True)
domain.set_name('circular')                  # Output name
print 'memory usage before del domain',mem_usage()
       
    print domain.statistics()
    print 'triangles',len(domain)


#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
def topography(x,y):
    """Complex topography defined by a function of vectors x and y
    """
    
                                    
    z= 0*x                      # defining z for all values other than the if statements
    r= 3.6                      # radius, provided in document
    angle = 14                  # angle, provided in document
    h= r*tan(angle/57.2957795)  # finding height of cone if not truncated
    

    N = len(x)
    for i in range(N):
       
        #truncated top
        if (x[i]-Cx)**2 + (y[i]-Cy)**2 <1.1**2:
            z[i] += 0.625

        # cone
        if (x[i]-Cx)**2 + (y[i]-Cy)**2 <r**2 and (x[i]-Cx)**2 + (y[i]-Cy)**2 >1.1**2:
           z[i] = -(sqrt(((x[i]-Cx)**2+(y[i]-Cy)**2)/((r/h)**2))-h)
    return z   

domain.set_quantity('elevation', topography, verbose=True)  # Use function for elevation
domain.set_quantity('friction', 0.01)         # Constant friction 
domain.set_quantity('stage',water_depth)   # Dry initial condition


#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
# Create boundary function from timeseries provided in file

##boundary_filename="ts2cnew1_input_20_80sec_new"
##prepare_timeboundary(boundary_filename+'.txt')
##
##function = file_function(boundary_filename+'.tms',
##                         domain, verbose=True)
def wave_form(t):
    return 0.1*sin(2*pi*t/50.)

# Create and assign boundary objects
Bw = Dirichlet_boundary([water_depth, 0, 0])          #wall
Bt = Transmissive_Momentum_Set_Stage_boundary(domain, wave_form) #wavemaker

domain.set_boundary({'left': Bw, 'right': Bw, 'top': Bw, 'wavemaker': Bt})




#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = 0.2, finaltime = 100.0):
    print domain.timestepping_statistics()