"""Example of shallow water wave equation.

Island surrounded by water

"""

######################
# Module imports 
#
from os import sep, path
from mesh_factory import rectangular
from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
     Constant_height
from Numeric import array
from util import Polygon_function, read_polygon
from math import exp


#Create basic mesh
N = 20
points, vertices, boundary = rectangular(N, N, 100, 100)

#Create shallow water domain
domain = Domain(points, vertices, boundary)
domain.store = True
domain.smooth = False
domain.visualise = False
domain.set_name('island')
print 'Output being written to ' + domain.get_datadir() + sep + \
	      domain.get_name() + '.' + domain.format
	

domain.default_order=2


#I tried to introduce this parameter top control the h-limiter,
#but it doesn't remove the 'lapping water'
#NEW (ON): I think it has been fixed (or at least reduced significantly) now
#
# beta_h == 1.0 means that the largest gradients (on h) are allowed
# beta_h == 0.0 means that constant (1st order) gradients are introduced 
# on h. This is equivalent to the constant depth used previously.
domain.beta_h = 0.2


#Initial conditions

def island(x, y):
    z = 0*x
    for i in range(len(x)):
        z[i] = 8*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 )    
    return z	

domain.set_quantity('friction', 0.0)
domain.set_quantity('elevation', island)
domain.set_quantity('stage', 1)
	


######################
# Boundary conditions
Br = Reflective_boundary(domain)

domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
domain.check_integrity()


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

print 'Done'